ARC 3.0: An expanded Python toolbox for atomic physics calculations

ARC 3.0 is a modular, object-oriented Python library combining data and algorithms to enable the calculation of a range of properties of alkali and divalent atoms. Building on the initial version of the ARC library [N. \v{S}ibali\'c et al, Comput. Phys. Commun. 220, 319 (2017)], which focused on Rydberg states of alkali atoms, this major upgrade introduces support for divalent atoms. It also adds new methods for working with atom-surface interactions, for modelling ultracold atoms in optical lattices and for calculating valence electron wave functions and dynamic polarisabilities. Such calculations have applications in a variety of fields, e.g., in the quantum simulation of many-body physics, in atom-based sensing of DC and AC fields (including in microwave and THz metrology) and in the development of quantum gate protocols. ARC 3.0 comes with an extensive documentation including numerous examples. Its modular structure facilitates its application to a wide range of problems in atom-based quantum technologies.


Introduction
Neutral atoms are ideal building blocks for both fundamental research and technological applications exploiting quantum mechanics. The stability of their physical properties makes them an ideal choice for a sensing medium, virtually eliminating the need for recalibration [1,2]. By promoting atoms to highly-excited, long-lived Rydberg states, the sensitivity to applied electric fields can be enhanced by many orders of magnitude compared to the ground state [3], which has applications to electric field metrology over a wide range of frequencies from DC to THz fields [2,4,5,6,7].
The reproducible properties of a given atomic species allows for the creation of hundreds of identical atomic qubits [8,9]. Together with the dramatic increase in the interatomic interaction strength associated with Rydberg excitation, this enables applications in quantum simulation [8,10,11,12] and quantum optics [13,14,15,16,17].
Having precise and accurate values of the relevant atomic parameters (e.g., transition dipole moments, Stark shifts, etc.) is crucial for developing new experiments, comparing to theory, or simply for interpreting the data. As the number of potentially relevant parameters vastly exceeds what can be tabulated, there is a need for an atom-calculator computing these parameters on demand. ARC 1.0 [18] aimed to provide such a research tool for alkali metals by combining the best available algorithms and experimental data into a modular, object-oriented Python library. The open source nature of the package and the popularity of this programming language facilitated its rapid adoption by the community and the inclusion of newly developed methods.
However, a growing number of research groups are now working on Rydberg states of atoms with two valence electrons [19], such as the alkaline earth elements (e.g., Ca, Sr) and some of the lanthanides (e.g., Yb, Ho [20]). While several groups have performed calculations of Stark maps [21,22,23] and interaction potentials [24,25,26], open-source codes like those provided by ARC 1.0 [18] and Pairinteraction [27] for the alkalis are not available for these species. This motivated the major upgrade of the ARC library we are presenting here. Specifically, we have extended ARC to divalent atoms, modifying the algorithms as necessary and including supporting data for 88 Sr, 40 Ca and 174 Yb. Methods for calculating long-range van der Waals interactions have been extended to degenerate perturbative calculations and pair-state calculations to arbitrary interspecies calculations. This new version of ARC also adds functionalities for atom-surface interactions, optical trapping and visualisation of results. The miniaturisation of atom-based sensors [28], and the strong long-wavelength transitions between Rydberg states [3], make atom-surface interaction effects prominent [29]. We have therefore included non-retarded van der Waals atom-surface potential calculations in ARC 3.0. Data are currently included only for sapphire and perfectly reflective surfaces, but other surfaces can easily be incorporated.
The importance of optical trapping stimulated the addition of two modules. OpticalLattice1D is dedicated to optical lattices and allows easy calculation of Bloch bands, Bloch states and Wannier states. DynamicPolarizability calculates the wavelength-dependent atomic polarizability relevant for optical trapping, and enables searches for "magic-wavelengths" where two different atomic states have the same polarizability.
Finally, to allow easy visualisation of different atomic states, a module Wavefunction provides sectional views of the atomic wave function for arbitrary atomic states. These are important both for pedagogical and research purposes, especially since the size of Rydberg electron orbitals can be large enough to encompass other atoms, and may even approach the typical lengthscale over which external trapping potentials vary [30,31].
How to start using the library is explained in Section 2. The rest of the paper gives an overview of the newly implemented calculations, with comments on restrictions and implementation details.

Installation and getting started
ARC 3.0 is available from the online repository Python Package Index (PyPI) and can be installed from the command line simply by invoking pip install ARC-Alkali-Rydberg-Calculator This installs the package correct for the user's environment, which can be based on Python 3 for Windows, Mac OS or Linux. All the methods discussed in the following can be used after importing them from the arc library: The examples of code given in the paper assume that this line was included at the beginning of the program to import the library.
The detailed documentation of the ARC library is available online [32] and is updated at each upgrade of the code. New users are recommended to consult the accompanying Jupyter iPython interactive notebooks listed on the "Getting started with ARC" page of the online documentation [32], which provide examples of calculations and describe the relevant physics. Finally, we note that a selection of the functionality of the ARC library is available at https://atomcalc.jqc.org.uk as a part of the web-app Atom Calculator. This web page will also generate code that can be used as example-on-demand to help users start their own calculations. Bug reporting, questions and further code development are tracked on the project's GitHub page [33].

Overview of the new functions
A high-level view of the ARC 3.0 library is shown in Fig. 1. This upgrade generalises methods originally included in ARC 1.0, which were developed for AlkaliAtom with electron spin S = 1/2, to the singlet (S = 0) and triplet (S = 1) states of divalent atoms. This approach is underpinned by the broad validity of the single active electron approximation for highly-excited states of divalent atoms, as discussed in Appendix A. Thus most of the divalent_atom_functions methods are directly inherited from alkali_atom_functions, upgraded as necessary to support divalent atoms. The spin state is now specified through an optional keyword parameter s, which defaults to 0.5 for alkali atoms and must be set to 0 or 1 for divalent atoms. To support these calculations, a semi-classical method for calculating dipole and quadrupole matrix elements has been implemented that does not require the use of a model potential.
The range of divalent atom data currently included in ARC 3.0 is discussed in Appendix Appendix B. New species/series may easily be added by generating the appropriate data files. Inter-species calculations for both PairStateInteractions and StarkMapResonances are possible.
In addition, PairStateInteractions.getC6perturbatively now supports degenerate perturbation theory for the calculation of C 6 coefficients, The possibility of taking into account a weak magnetic field directed along the quantization axis has also been introduced (only paramagnetic terms linear in field strength which shift energy levels but do not mix states are included in the calculation  Finally, we have also created a group of advanced calculations which are likely to be too specialized for a core toolbox. The corresponding modules are intended for expert use and can be imported with # from arc.advanced.<mn> import * # where <mn> is module name, for example: from arc.advanced.population_lifetime import getPopulationLifetime The modules in this advanced library will be built on top of the ARC core library and will provide solutions for specialised research questions. The first of these modules is population_lifetime, which has been contributed by the authors of Ref. [34]. It gives access to a getPopulationLifetime function calculating population lifetimes taking into account the redistribution of population within a Rydberg manifold driven by black-body radiation, including repopulation processes (such population lifetimes thus differ from state lifetimes).
In the following we outline the calculations newly implemented in ARC 3.0 and provide examples of code for each of them.
As noted above, the calculations are performed within the single active electron approximation. The radial matrix element Table 1: Energy level data for divalent atoms currently included in diva-lent_atom_data, with references. for a dipole transition between states |n, L, J and |n , L , J of hydrogen or an alkali atom are calculated as in ARC 1.0, i.e., by numerical integration between suitably chosen bounds r i and r o : with S = S = 1/2 and the wave functions R nLJ (r) obtained as solutions of the Schrödinger equation for a model potential.
Since model potential methods are problematic for divalent atoms (Appendix A) the dipole radial matrix elements for these species are obtained in a semi-classical approximation [50]. This approach is not based on numerical wave functions. Instead, the dipole radial matrix elements take on the form with S = S = 0 or 1. In this equation n * is the reduced principal quantum number (n * ≡ n − δ(n, L, S , J) with δ(n, L, S , J) being the quantum defect for the |n, L, S , J state) and c , n * c , γ, ∆ and ∆n * are defined as Moreover, where the J s (−∆n * )'s are Anger functions: The radial quadrupole matrix elements are calculated as in the case of hydrogen and alkali atoms, and by using the corresponding semi-classical formulae for divalent atoms. The latter differ between different values of |∆L|. For ∆L = ±2, whereas for ∆L = 0, The expansion coefficients Q p (∆n * ) are the same in both cases: Q 0 (∆n * ) = − 6 5(∆n * ) 2 g 1 (∆n * ), sin π∆n * π(∆n * ) 2 , We have compared the model potential and semiclassical methods for calculating dipole matrix elements using rubidium as a test case. The two methods give results in close agreement when |n − n | ≈ 0, even for n, n as low as 10 (below which multi-electron effects can be expected to become important). While the agreement deteriorates when |n − n | increases, the semi-classical results do not differ by more than 5% from the model potential results in the range 0.65 n n 1.5 n. Given that outside this range the dipole matrix elements are less than 1% of their values at |n − n | ≈ 0 (except for very small values of n), the semi-classical approach is normally appropriate for any value of n − n for which the dipole matrix element is large enough to be significant in calculations of Stark maps or dispersion coefficients. A multi-channel quantum defect approach would be more appropriate in the regions where perturbers mix states of different symmetries [51]; however, such calculations are beyond the scope of the current version of ARC.
All dipole and quadrupole matrix elements for calcium, strontium and ytterbium currently used or produced by ARC 3.0 are obtained as described above, as no literature values of these quantities are currently available for these atoms. Matrix elements obtained in the future from more accurate calculations or from measurements can be added to a literature file, as described in Section 4.3. The library will use the values found in this literature file, should there be any, rather than recalculate them.

3.2.
Pair-state calculation of atom-atom C 6 interactions in degenerate perturbation theory In the single-active electron approximation, the interaction between two atoms at inter-atomic distance R arises from the interaction between the valence electrons and the interaction of each of these electrons with the screened nucleus of the other atom. Denoting the coordinates of the valence electrons relative to the respective nuclei by r 1 and r 2 , this interaction is given by the following multipolar expansion (c.f. Eq. (19) in Ref. [18] for details): propagating in z-direction (c) J = 1 Figure 2: Two examples of cases where degenerate perturbation theory has to be used for C 6 calculations. (a): Two atoms with inter-atomic axis parallel to the quantisation axis may exchange virtual photons that drive two σ + transitions [corresponding to Y 1,+1 in Eq. (7)] in one atom and two σ − transitions in the other atom, thus coupling the initial |r 1 and |r 2 states to energetically permitted |r 3 and |r 4 states. The projection of the total angular momentum on the quantisation axis is conserved in this process (i.e., m J 1 + m J 2 = m J 3 + m J 4 ). (b): Alternatively, when the inter-atomic axis is perpendicular to the quantisation axis, virtual photons driving a σ − transition in one atom may be viewed as having linear polarization perpendicular to the quantisation axis and can thus drive both σ + and σ − transitions in the other atom. The projection of the total angular momentum on the quantisation axis is then not conserved and m J 3 + m J 4 may differ from m J 1 + m J 2 . Using degenerate perturbation theory, the energy eigenstates can be obtained with their corresponding C 6 coefficients, as shown in (c) for the case of two 88 Sr atoms initially in a pair state with n 1 = n 2 = 40, L 1 = L 2 , J 1 = J 2 and S 1 = S 2 (the red circles correspond to the stretched states m J 1 = J 1 = m J 2 = J 2 , assuming that the inter-atomic axis is parallel to the quantization axis). with Here the n m are binomial coefficients and the Y k,p (r) are spherical harmonics, and it is assumed that the quantization axis of both atoms is directed along the internuclear axis [as in Fig. 2(a)]. Terms with k 1 + k 2 = 2, 3, 4, . . . correspond respectively to dipole-dipole, dipole-quadrupole, quadrupole-quadrupole interactions and so on.
In the large R limit, the interaction between Rydberg atoms is typically dominated by dipole-dipole interactions. Such an interaction couples an initial pair-state |r 1 r 2 to other states |r r whose energy differs by an energy defect ∆ r ,r = E r + E r − E r 1 − E r 2 . To second order in this interaction, the interaction energy of two atoms in the initial pair-state |r 1 r 2 is given by a van der Waals interaction potential of the form −C 6 /R 6 , where C 6 is defined by the following equation if the pair-state energy E r 1 + E r 2 is non-degenerate: Here the sum goes over all the pair states |r r that are coupled by electric dipole transitions to the original pair-state |r 1 r 2 , and V dd (R) is the dipole-dipole part of V(R): However, there are situations in which the pair-state energy E r 1 + E r 2 is degenerate and the dipole-dipole interaction couples the initial pair-state |r 1 r 2 to some other pair-states |r 3 r 4 of the same energy -for instance, in the absence of external fields, to pair-states differing in magnetic quantum numbers only. Two examples of such off-diagonal couplings mixing pair-states are given in Fig. 2. In such situations, a perturbative calculation of the C 6 coefficients requires the diagonalization of the matrix C describing the second-order coupling between the energydegenerate pair-states [24]. To this end, we work in the basis {| r i r j } of degenerate pair-states, i.e., where (n 1 , L 1 , S 1 , J 1 , m J 1 ) and (n 2 , L 2 , S 2 , J 2 , m J 2 ) are the quantum numbers associated with the states | r 1 of atom 1 and | r 2 of atom 2, respectively. The element of C corresponding to the second-order coupling between the q-th and the p-th basis states is Diagonalising this matrix yields a C 6 coefficient for each of the energy curves the degenerate pair-state energy splits into due to the dipole-dipole interaction, these coefficients being the eigenvalues C (i) 6 of C [see, e.g., Fig. 2(c)]. These eigenvalues are independent of the orientation of the inter-atomic axis relative to the quantization axis, in the absence of external fields, although the composition of the corresponding energy eigenstates in terms of the basis states defined above depends on the choice of the quantization axis [24].
Degenerate perturbation theory can be used within the getC6perturbatively function by setting the flag degener-atePerturbation=True. The function will then return the C (i) 6 eigenvalues and the corresponding eigenvectors of the relevant C matrix. For example, the following fragment of code produces the results displayed in the third column of Fig. 2 Fig. 2

(c).
Here theta and phi are, respectively, the polar angle (Θ) and azimuthal angle (Φ) defining the orientation of the inter-atomic axis in a referential whose z-axis is parallel to the axis of quantisation of the angular momenta. The dependence on Θ and Φ of the elements of C is taken into account by rotating the atomic basis states using Wigner D-matrices, as in ARC 1.0 [18]. I.e., V dd (R) defined above is replaced by the angle-dependent V dd (R, Θ, Φ), with, in a simplified notation, where the D(J m J , Θ, Φ) represent the relevant rotation matrices. Although the elements of the matrix C and the composition of its eigenstates in terms of the basis states defined above depend on Θ and Φ, this is not the case for its eigenvalues C (i) 6 (hence for the C 6 coefficients resulting from this calculation) [24]. Invoking getC6perturbatively without the flag degener-atePerturbation=True or with degeneratePerturbation=False will only return the individual element of the matrix C corresponding to the values of the quantum numbers specified in the call, rather than the eigenvalues and eigenvectors of C. These individual elements normally depend on Θ and Φ and can be taken to be effective C 6 coefficients for particular combinations of magnetic quantum numbers.
We note that the interaction energies can also be obtained non-perturbatively by full diagonalisation of the Hamiltonian using the function diagonalise. An example of results obtained Sr 60 3 S 1 m j = 1, Sr 60 3 S 1 m j = 1 Figure 3: Example interaction strength calculations using the full diagonalisation method for a pair of divalent atoms in identical states ( 88 Sr 5s60s 3 S 1 m J = 1). The plot shows how the energy levels of the atom pair vary as functions of the inter-atomic distance, R. The energies are measured relative to the energy of the initial pair-state in the R → ∞ limit, and only pair-states coupled to this initial pair-state are represented. The fraction of the initial pair state present in each eigenstate is colour-coded as per the colour bar.
in this way is given by Fig. 3, which is produced by running the following code: We also note that the implementation of degenerate perturbation theory made in ARC does not take into account energy degeneracies between states differing in L or S . It is therefore not appropriate for atomic hydrogen or for high angular momentum states. We recommend using the full diagonalisation method for such cases.

Inter-species interaction calculations
PairStateInteractions supports inter-species calculations. Users can initialize such calculations using the keyword argument atom2 to explicitly state the species of the second atom.
Note that to specify the spin state of the second atom, the keyword argument s2 should also be set. Setting s2 can also be used for calculations where the atoms are from the same atomic species but have different spin.
For example, pair-state calculations between rubidium atoms in the | 5s60s 2 S 1/2 m J = 1/2 state and ytterbium atoms in the | 6s54s 1 S 0 m J = 0 state [ Fig. 4] are initialized as follows: 3.4. Single-atom properties for divalent atoms ARC 3.0 extends most of the single atom methods available in ARC 1.0 to divalent atoms. For example, Stark maps can be obtained for divalent atoms by setting the additional key argument s to define the spin state. Note that in the single electron approximation, states of different total spin are not coupled. Example results from such a calculation are shown in Fig. 5. Static electric fields are often used to adjust pair-state energies. StarkMapResonances allows the user to search for electric field strengths where two pair-states have same energies (Förster resonances). Lastly LevelPlot allows the plotting and interactive exploration of energy level diagrams. These diagrams may be opened as interactive stand-alone plots (from a command line Python call or in a Jupyter notebook with %matplotlib qt); then will then display the transition wavelength and transition frequency for pairs of states selected interactively by clicking on energy levels.

Electronic wave functions
Visualisations of atomic wave functions are a useful pedagogical tool enabling visual interpretation of effects such as  dipole moments. In addition, the detailed shape of Rydberg electron wave functions plays an important role in a number of effects. For example, the wave functions can become so spread out that they encompass other atoms. The modulations of the electron probability density may then induce a significant variation in the potential energy of the encompassed atoms, which may lead to the formation of Rydberg molecules [53]. At the same time, the optical potentials used for atom trapping may vary substantially over the length-scale of the Rydberg electron wave functions, giving rise to energy shifts [30] and affecting the trap lifetime [31].
Wavefunction enables the calculation of atomic wave functions for arbitrary superposition states. Quick 2D and 3D visualisations are possible, with a choice of units (atomic or SI). For example, the following code can be used to obtain and plot the probability density function for the 10f 2 F 7/2 m J = 7/2 state of rubidium [Figs. 6(a and b)]: atom = Rubidium() n = 10; l=3; j=3.5; mj=3.5; stateBasis = [[n, l, j, mj]] stateCoef = [1] # pure 10 F_7/2 mj=7/2 state wf = Wavefunction(atom, stateBasis, stateCoef) wf.plot2D(plane="x-z", units="atomic"); plt.show() wf.plot2D(plane="x-y", units="atomic"); plt.show() Wavefunction can be integrated with other ARC functions. For example, one can find the electronic wave function for an atom perturbed by an electric field by getting the state from the corresponding StarkMap and using Wavefunction as per the following code to plot the result. An example is shown in Figs. 6(c-d). wf.plot2D(plane="x-z", units="nm", pointsPerAxis =400, axisLength=2800) plt.show() The density and scale of the mesh on which the wave functions are calculated can be adjusted with optional keyword parameters. The probability density functions can be provided in Cartesian as well as in spherical coordinates (respectively through getRtimesPsi and getRtimesPsiSpherical), and it is also possible to obtain arrays of wave functions for all different possible spin states (using the getPsi method).
Note that Wavefunction is currently only supported for species in the alkali_atom_data class, as calculations based on model potentials are not currently supported for divalent atoms.

Atom-surface interactions: van der Waals potentials
In the vicinity of a surface, an atomic dipole interacts with its image in the surface, leading to shifts of the atomic energy levels. For small atom-surface distances z < λ/(2π), where λ is the wavelength associated with the strongest transition, the interaction potential V(r) for an atom in the state a is of the non-retarded, van der Waals form [54] Here n(ω ab ) is the frequency dependent refractive index associated with the surface. The summation covers all the states b dipole coupled to state a. The corresponding transition frequencies and dipole matrix elements in the x, y and z directions are respectively ω ab , µ ab x , µ ab y and µ ab z . The z-axis is taken to be perpendicular to the surface. Note that different states will have different C 3 values, which leads to a modification of the atomic transition frequencies near the surface.
To specify the surface material, ARC provides an abstract class materials.OpticalMaterial, with a method getN returning the refractive index n for a specified wavelength. A subclass Sapphire is provided as an example. The AtomSurfaceVdW uses information on the optical properties of the surface and the atomic transition frequencies to calculate C 3 . For example, the following code will return the energy shift with error for the 6s 2 S 1/2 state of caesium in the proximity of a sapphire surface (the result is 1.259 (2)  Such calculations are possible for all the atomic species supported by ARC, to the extent that the required dipole matrix elements are available.

Optical lattices: Bloch bands, Bloch states, Wannier states
Atoms in optical lattices are important in many areas of science and technology including atomic clocks and gravimeters, quantum gas microscopes and quantum simulators. A laser standing wave gives rise, through the AC Stark shift (Sec. 3.8), to a spatially periodic potential for the atoms. As is well-known from solid state physics, such a periodic potential can also be considered in reciprocal space. In momentum space (k−basis), a potential with spatial period λ/2 couples free particle states that are separated by an integer multiple of δk = 4π/λ = 2k ≡ k l where λ is the wavelength of the optical field, and k l = 2k is the lattice momentum. Therefore, from the Bloch theorem, the eigenfunctions of the Hamiltonian (called Bloch states or Bloch wave functions in this context) can be parametrized by a quasimomentum q (|q| < k), and these eigenfunctions are of the form exp(iqr) times a periodic function of r whose period is the same as that of the lattice. The Hamiltonian can thus be diagonalised on a discrete basis of free-particle states {e iqr , e i(q+k l )r , e i(q+2k l )r , . . .} separated by an integer multiple of the lattice momentum in momentum space. Plotting the resulting energy levels for different values of the quasimomentum gives Bloch bands [ Fig. 7(a)].
Both Bloch states and Bloch bands can be calculated using OpticalLattice1D. For example, the following code produces a Bloch band diagram for rubidium atoms trapped in an optical lattice formed by a 1064 nm standing wave with a maximal depth of 5 E R , using a basis which includes states up to q + 35k l in momentum (E R = h 2 /(2mλ 2 ) is the recoil energy):  The Bloch wave functions are delocalised across the lattice sites [ Fig. 7(b) bottom panel]. For many calculations and a more intuitive mapping to atomic physics experiments, it is convenient to switch to a basis of localised functions, namely Wannier functions. For each Bloch band, the Wannier functions w i,R (x) are defined (up to a normalisation factor) as a complete orthogonal set of functions localised at lattice points defined by a lattice vector R: where the sum goes over all values of the quasimomentum q and b i,q are the Bloch wave functions for a given Bloch band index i and quasimomentum q. Values of the Wannier function can be obtained after diagonalisation of the interaction potential Hamiltonian for which we defined Bloch band index to be saved by setting saveBandIndex keyword argument in diagonalise method. We can then call getWannierFunction to obtain values of the Wannier function in a given Bloch band [ Fig. 7 Note that saveBandIndex selects the band index, and latti-ceIndex sets the localisation of the function at the site with the corresponding index. The Wannier functions returned by the program should be normalised on the relevant lattice by the user.

Dynamic polarisabilities and magic wavelengths
The dynamic (AC) polarisability α(ω) of an atom exposed to an oscillating electric field of angular frequency ω can be expressed as the sum of a contribution from the polarisability of the valence electron(s) α v and the core polarisability α c . The valence polarisability often dominates. For an electron in state |a with total angular momentum J and projection m J the valence polarisability for linearly polarised light can be written as [55] α (ω); (14) that is as the sum of a scalar polarisability and a tensor polarisability Here b || er || a are reduced dipole matrix elements and the summation runs over all the bound states |b , with total angular momentum j b , dipole-coupled to the state |a of interest. Finally, there is also a term α (cont.) 0 contributed by the continuum of unbound states (this contribution will be discussed in more detail below).
For example, the scalar and tensor polarisabilities of the caesium 100p 2 P 1/2 state in the AC field given by a 1100 nm wavelength optical trapping laser can be obtained using the following code (in this calculation, the sum over intermediate states b is limited to all the bound states with n ≤ 115, including the ground state and the low lying excited states): atom = Caesium() n = 100 calc = DynamicPolarizability(atom, n, 1, 1.5) calc.defineBasis(atom.groundStateN, n+15) alpha0, alpha1, alpha2, alphaC, alphaP, closestState = calc.getPolarizability(1100e-9, units="SI") closestState saves the state whose transition frequency is closest to that of the driving field. In addition to alpha0 and alpha2 the method getPolarizability also returns the core polarizability, α c (alphaC), and the ponderomotive polarizability, α p (alphaP). The core polarizability is approximated by its static value (saved in atom.alphaC) which is appropriate when the driving laser is far from resonance with core transitions. The ponderomotive polarizability can be linked to the contribution from the continuum of unbound states α (cont.) 0 (which is currently not calculated by ARC). Close to a bound-state resonance, the contributions from the bound intermediate states dominate α v (ω). However, for Rydberg states one is often in the opposite limit far from boundstate resonances. Away from strong resonances, the total scalar polarisability α 0 +α (cont.) 0 then approaches the free-electron ponderomotive polarisability α P associated with the time-averaged motion of a free electron in an oscillating electric field [α P = e 2 /(2hm e ω 2 ), where ω is frequency of driving field and e and m e are electron charge and mass respectively]. This result is obtained by applying the limit where the Keppler frequency of the electron orbit around the core is much smaller than ω as is typical for Rydberg states [56], and applying a Born-Oppenheimer approximation to separate the fast quiver motion of the loosely bound electron driven by the electric field component and the slower relative motion of electron around the ionic core.
The code also returns alpha1 which is a vector polarisability relevant in driving with non-linearly polarised light, see, e.g., Ref. [57]. This feature is not exploited in the current implementation of DynamicPolarizability, which focuses on the simplest case of driving under linearly polarized light; however it provides a path for future support of other polarizations. Other future extenstions could include adding bound-state resonances to the core polarizability and an extended treatment of the continuum contribution near the ionization threshold.
DynamicPolarizability calculations can involve dipole matrix elements between low-lying states or between low-lying states and highly excited states. For low-lying states dipole matrix elements calculated both with the semiclassical approximation and the model potential method become less accurate. That is why ARC uses literature values for these states, where available. Users should check the availability of literature values in the corresponding files (Sec. 4.3). If values for low-lying transitions are available, the expected accuracy can be of the order of ∼ 1%, otherwise the accuracy is limited to ∼ 10% for typical input parameters. Currently, there are significant compiled literature sources for dipole matrix elements in alkali atoms, but the data for divalent atoms is more scarce.

Handling of the spin quantum number
To allow the use of both alkali atoms and divalent atoms, most of the methods have a new keyword argument s specifying the spin quantum number The default value of this keyword argument, s=0.5, is appropriate for alkali atoms. As such, s does not have to be specified for these atoms, making the API of ARC 3.0 completely backwards compatible with ARC 1.0. For divalent atoms, users should specify whether they are working with singlet or triplet states by setting the corresponding spin quantum number to s=0 or s=1 respectively. To state that the second atom has a different spin in PairStateInteractions, the optional keyword argument s2 should be set explicitly.

Use of fitted quantum defects for divalent atoms
The available experimental energy level data is much sparser and less precise for divalent atoms than for the alkalis. Inaccuracies arising from experimental errors may be reduced by using energies derived from a fit of the data to the Rydberg-Ritz formula over a range of principal quantum numbers (see Appendix Appendix B), rather than using the experimental energies directly. However, in many cases the calculations will require energies outside the range of values of n used in the fitting to the Rydberg-Ritz formula. Using this formula is normally justified for values of n above this range, but might be invalid for values of n below this range (e.g., because the Rydberg-Ritz formula does not apply to the ground state and low excited states, or because the energy levels do not vary smoothly and monotonically with n due to perturbers). Using the experimental energies is thus often preferable for small values of n.
For clarity, we have tabulated the range of values of n over which the quantum defects have been fitted using a dictionary variable defectFittingRange indexed by terms. For each term the dictionary returns the smallest and largest principal quantum numbers defining the range of principal quantum numbers used in the calculation of the Rydberg-Ritz coefficients provided in ARC. For example, the following code will do this for the 1  For principal quantum numbers above the fitted range, the energies are calculated by extrapolation of the Rydberg-Ritz formula to outside this range. No extrapolation is done for principal quantum numbers below the fitted range. Instead, ARC tries to use tabulated energies if they exist -either energies provided by the user or the energies provided by ARC, which are literature values as described in the online documentation [32]. In the case of S, P or D series, the calculation is aborted and a value error raised, explaining the problem, if the program requires a missing tabulated energy. For L > 2, the calculation is not aborted but a value of zero is assumed for the quantum defect. These tabulated energies must be stored in the local data directory (Sec. 4.3) under a file name defined for each atom in the levelDataFromNIST variable. For example for strontium they are stored in sr_level_data.csv. The energyLevelsExtrapolated variable of the used atom will be set to True if the calculation uses energies obtained by extrapolation of the Rydberg-Ritz formula to outside the range of values of n used in the fitting.

Local data directory and updates to literature values
Dipole matrix elements from the literature can be added to ARC by setting the literatureDMEfilename variable in the DivalentAtom class to the appropriate file name and ensuring the format of the file matches that specified in the documentation.
When ARC 3.0 is first used, a local hidden directory named .arc-data is created in the user's home folder. The results of calculations of, e.g, dipole matrix elements are saved in humanreadable files in this folder, forming a look-up table that can be used to speed up future calculations. There is an option to add dipole matrix elements from the literature (or obtained using other calculation methods) by modifying the corresponding files. For example, to add a new literature value for a 88 Sr dipole matrix element, the user should modify the file stron-tium_literature_dme.csv as per the table header. This can be conveniently done by loading the appropriate .csv file in any spreadsheet program, as long as the file format is preserved in terms of separators and header comments. To change other atomic properties, like quantum defects, users should make their own subclasses which inherit classes defined in the ARC. An example is given in the following code, which could be used to update the quantum defects for calcium: class MyCalcium40(Calcium40): quantumDefect = ... # write updated quantum defects # list in order specified in documentation dipoleMatrixElementFile = "my_ca_dme.npy" Once defined, MyCalcium40 can be used instead of Calcium40 in all the ARC calculations. To ensure the program does not use old cached values, new names of caching files should be provided when redefining atomic species in this way. Also, note that if a new version of arc-data is released (changing the version number in version.txt in the data directory), the data entered manually in the .arc-data folder may be overwritten automatically if no new file names had been specified. We encourage users to submit new experimental data and parameters for use by the community via a pull request on the ARC GitHub page [33].

Looking under the hood with debugOutput=True
As for ARC 1.0, setting the keyword argument debugOut-put=True results in verbose output, which may be useful for checking basis states and other intermediate results. Additionally, many methods have a progressOutput=True option which can be used for tracking the progress of the calculations.

Outlook
In summary, we have presented ARC 3.0, a major new release of the ARC Python library that extends the library to divalent atoms and adds a number of new methods of general interest in Rydberg physics and beyond. We believe this common code base and consistent interface for many different atomic properties can speed up the development of new applications and lower knowledge barriers, e.g, in quantum technologies based on neutral atoms. The library also offers rich possibilities of advanced educational projects for students in atomic and quantum physics.
Future improvements of ARC could include the addition of the calculation of wave functions for divalent atoms and additional methods for the accurate calculation of low-lying wave functions [58], support for multi-channel quantum defect theory calculations [51], and the calculation of scattering properties. New experimental data can be straightforwardly added to the existing base. ARC is a community-oriented open source package, and the authors welcome contributions of new core data or algorithms to the main library as well as contributions of more specialized codes to arc.advanced.

Acknowledgements
We thank Tsz-Chun Tsui and Trey Porto for tabulating the ytterbium data. We also thank Hei Yin Andrea Kam for the use of her code to check the calculations of strontium dipole matrix elements, Paul Huillery for help with Stark map calculations, and Ifan Hughes for suggesting the bootstrap method. N. Š. is supported by the H2020 Marie Skłodowska-Curie Actions (CO-QUDDE, H2020-MSCA-IF-2017 grant agreement No.786702). The project was supported by the EPSRC Platform grant EP/R002061/1 and Standard grant EP/R035482/1. We also acknowledge funding from the project EMPIR-USOQS (EMPIR projects are cofunded by the European Union's Horizon2020 research and innovation programme and the EMPIR Participating States).

Appendix A. The single active electron approximation
Atomic structure calculations for Rydberg states are simplified by the rapid scaling of the size of the wave function of the Rydberg electron with principal quantum number n. In atoms with a single valence electron, the calculation of energy levels and wave functions is achieved by solving the Schrödinger equation with a modified "model" potential that accounts for the screening effects of the closed-shell core. Each Rydberg series labelled by the quantum numbers L, S and J is characterized by a quantum defect δn that quantifies the deviation of each energy level from its hydrogenoic equivalent. This method works well for alkali atoms and is used by alkali_atom_functions.
For divalent atoms, the situation is complicated by two effects. The first is that the ionic core is no longer a closed shell, since it contains the remaining valence electron. The core is therefore more strongly polarized by the Rydberg electron, leading to an additional n dependence in the quantum defects. The second effect is the existence of compact states where both electrons are excited (e.g., a 5p 2 configuration), known as perturbers. These perturber states exhibit strong interelectronic correlation effects and lead to perturbations of nearby Rydberg levels. Even away from the energy of a perturber, Rydberg states may acquire a small admixture of doubly excited states [59].
Nevertheless, under many circumstances the properties of Rydberg states in divalent atoms can be described under a single active electron approximation [24]. Under this approximation, Russell-Saunders L − S coupling is assumed to hold, and L, S and J are regarded as good quantum numbers. The effects of core polarizability and perturbers are partially included via the energy dependence of the quantum defects. Calculations can then be made in a similar way to those performed in alkali atoms, with the appropriate generalization of the angular momentum algebra to integer spin. Previous work has shown that this treatment gives good agreement with experiment for quantities that depend on the long-range part of the wave function, such as the DC Stark effect [22,60].
The interaction between two divalent Rydberg atoms is also dominated by the long-range character of the wave function. Therefore a single active electron approach may be used here also [24]. A study of the effect of perturbers on the interactions [61] demonstrated that the single-electron treatment is valid to a high degree of accuracy (<2%) except for Rydberg states in the immediate vicinity of perturbers.
Note that other observables that depend on the short-range properties of the wave function, such as the coupling to lowlying states (eg radiative lifetimes), are significantly modified by the presence of even small amounts of perturbing states, and are not well treated in the single active electron approximation [51]. Here other methods such as multi-channel quantum defect theory (MQDT) [51] that explicitly include the effects of perturbers must be used.

Appendix B. Atomic Data for Divalent Atoms
As outlined in section 3.1, calculations involving divalent atoms rely on fitted values for the quantum defects. In the case of Ca, no new data was available, and so values compiled in [24,52] are used. For Sr and Yb, we provide a new analysis of the available data for Sr and Yb that takes into account recent improvements in the spectroscopic data.
Experimental energy levels were fitted to the modified Rydberg-Ritz formula using δ 0 , δ 2 and δ 4 as free parameters. The ionization energy I a was constrained to the values obtained from the analysis of the best available spectroscopic data. Ry a is the atom-specific Rydberg constant: Ry a = R ∞ m a /(m a + m e ), where R ∞ is the Rydberg constant, and m a and m e are the mass of the considered isotope and the electron respectively. The fine-structure splitting of the ionization threshold was neglected. A least-squares fitting method was used, implemented via the curve_fit function from the scipy Python package.
The results, along with references to the experimental dataset are provided in Table B.1 The choice of the range of n used for each series was a compromise between maximising the number of data points, and reducing the effect of series perturbations not described by equation B.1 and experimental systematic uncertainties. Particular care should be taken in extrapolating to values of n below the stated range, where the fits are often poor.
The uncertainties on the Rydberg-Ritz parameters are 68% confidence limits obtained using a "bootstrap" method based on resampling with replacement [62], such that each confidence interval includes 68% of the results falling above the quoted value of the corresponding parameter as well as 68% of the results falling below it. For each series the fitted range was resampled 150 times. The asymmetry of the confidence limits reflects the asymmetric dependence of the value of the energy on the quantum defects encapsulated in Eq. (B.1), as well as the limitations of the experimental data. Note that correlations in the uncertainties are expected to be strong and are not explicitly considered here; users seeking to set rigorous error bounds on derived quantities should take this into account.