Software for Calculating Trap Induced Broadening in Potential Hydrogen Lattice Clocks

Introduction

There is growing interest in optically trapping atomic (anti-)hydrogen in a quest for improved spectroscopic precision [1, 2, 3, 4]. A full picture of optical trapping requires knowledge of the atomic polarisability, which determines optical confinement [5]; and atom-photon scattering rates, which limit trap lifetime and coherence times [5]. The usual approach to calculating these quantities in atomic hydrogen is analytic, expressing them in terms of hyper-geometric functions [6, 7, 8]. Despite the advantages inherent in an analytic approach, these calculations are extremely involved and inaccessible to those unfamiliar with the necessary techniques. In addition, a given solution is specific to a particular state or pair of states (in the case of atomic scattering), and considering a different state or pair requires re-calculation.

There are a number of available AMO physics packages which include routines for calculating atom polarisability of atom-photon scattering rates, such as ARC 3.0 [9] or atomphys [10] for Python. However, these routines are often designed to be general and calculations proceed via summation over a finite number of bound states, neglecting contributions from unbound states. Whilst these contributions are small in alkali atoms, they are known to be significant in hydrogen and cannot be ignored [1, 2]. One can express these sums in closed form via the implicit summation of Dalgarno and Lewis [11]. They can then be computed exactly for hydrogen, in a basis of Laguerre functions known as Sturmian functions, which span both the discrete and continuous parts of the energy spectrum. This approach to these calculations is well established, beginning with the pioneering calculations of Zernik [12] and later work by Gontier and Trahin [13]. The STRFLO program [14] applies similar methods to the calculation of photo-ionisation rates and interactions of hydrogen with bi-chromatic fields (we make use of this software in [1] to calculate multi-photon ionisation rates). Whilst these functions can be extended to calculations of scattering rates [15, 16], such code is unpublished. Recently, these techniques have been applied to efficient calculation of scattering rates in frequency regimes beyond the ionisation threshold [17]. In addition, R-matrix codes have previously been used to calculate polarisability from photo-ionisation rates [18]. Despite this, the authors know of no generally released software which allows for the calculation of polarisability and atom-photon scattering rates in atomic hydrogen which correctly includes the continuum contributions.

In this paper, we present software for the calculation of the atomic polarisability and atom-photon scattering rates of S (zero orbital angular momentum) state atomic hydrogen in far off-resonant optical fields. We apply an implicit summation method over a finite basis of radial Sturmian functions to naturally include all contributions from both bound and unbound states, while angular terms are calculated from analytic expressions [1, 2]. The software also allows for the calculation of magic wavelengths (where the atomic polarisability of two states are equal) for pairs of S-states. This software was developed for the calculations of trap potentials and scattering rates in [1] and [2]. As such, it is limited to calculations involving initial S states and assumes a linearly polarised optical field. Although developed with low lying states in mind, the software is applicable to any hydrogen (or anti-hydrogen) S state. It can be reasonably extended to other hydrogen-like systems or to processes involving other orbital angular momentum states.

Implementation and architecture

This software calculates matrix elements,

between a given pair of initial and final non-relativistic hydrogen states: a and b. States and operators are expressed as vectors and matrices across a finite basis of Sturmian functions with user-defined dimension n_max and free parameter k. For example, b is the vector describing the final state b in the Sturmian basis, whilst 𝒟^± is the matrix of the electric dipole operator (split into l-“raising” and “lowering” parts denoted +⁣∕⁣–). The vector $\Psi$ arises from the implicit summation method, it is a solution to the matrix equation,

with an initial state energy E_a and a user-defined value for the external field frequency ω. Here ℋ is the atomic Hamiltonian and 𝒯 is the overlap function for the Sturmian basis. The software assumes that the initial state a is an S state (zero orbital angular momentum, l = 0).

The module “main_calculation_functions.py” makes extensive use of the arrays of NumPy [19] and linear algebra routines of SciPy [20] to calculate these matrix elements for a given pair of states and external field wavelength. Matrices are constructed from Python dictionaries of analytic matrix elements as required, whilst state vectors are found by solving the Schrödinger equation as a generalised eigenvalue problem,

The correct vector is selected from the spectrum according to its eigenvalue — the binding energy E_a. The lower sign in equation (1) relates to processes where both b and a are S states, the upper signs are only relevant for S to D state scattering.

These matrix elements are then used to calculate the polarisability

and two photon scattering rates,

Here ω_s is the scattered photon angular momentum, α_FS is the fine structure constant, and A_ba is an analytical angular term, equal to $\frac{8\pi }{27}$ if b is an S state and $\frac{16\pi }{27}$ if it is a D state. The derivations of these expressions and further discussion of the numerical method can be found in [1] and [2]. Magic wavelengths are then calculated using the Newton-Raphson method in SciPy to find zeros in the differential polarisability between two states. Relevant constants are defined in this module and take their CODATA values [21]. Please note that all calculations use Hartree atomic units (see [22] for details) but this is usually converted to S.I. for outputs. The module also includes expressions for the reduced mass of hydrogen, deuterium, and tritium. The user can perform calculations for each isotope by selecting the appropriate reduced mass.

The software also contains a number of secondary modules. Four of these are used to produce the plots available in [1]:

• “1S_and_2S_polarisability.py” — plot showing the polarisability of the 1S and 2S states for wavelengths from 395 nm to 700 nm.

• “Variations_with_depth.py” — plot comparing the total inelastic scattering rate to the 2-photon ionisation rate at 514.6 nm for lattice depths varying from 1 to 100 times the recoil energy.

• “Compare_scattering_rates.py” — plot showing both the elastic and total inelastic scattering rates out of the 2S state from 395 nm to 1000 nm. Also compares to 2 and 3 photon ionisation rates.

• “Continuum_corrections.py” — plot showing the cumulative contribution of low lying, bound P states to the inelastic scattering rates at 4 reported magic wavelengths.

These modules utilise the representations defined in the core module, but often have their own calculation routines which are optimised for multiple iterative calculations across a range of wavelengths and final states. These are not essential parts of the software and are included in the repository: 1. for the sake of transparency, 2. as examples of how the software can be used. Very similar implementations were used to create the plots in [2].

The final module “Convergence_and_correctness_tests.py” contains a series of functions used to test the numerical consistency and convergence of calculations (see the next section). A complete list of additional modules and the included functions can be found in the “README” document included in the GitHub repository identified below.

Choice of basis parameters

For any given calculation, the user must define two basis parameters: n_max and k. The first defines the size of the Sturmian basis used in the calculations. Larger basis sizes result in improved computational accuracy at the expense of significantly increased computational load. Additionally, calculations involving higher lying states (states with larger values of n) require larger basis sets to properly converge; this is indicated in Figure 1. The second is a free parameter which appears in the definition of the Sturmian basis (see [1] and [2]). Whilst any positive definite value may be chosen, it is advantageous to choose k to be small, since larger values require larger basis sets to properly converge.

Figure 1

For the calculations in [1], the parameters n_max = 300 and k = 0.3 were chosen. This is sufficient for the calculations of low n states presented in [1], but for calculations involving higher lying states, a larger basis set or smaller k is necessary (see Figure 1).

Accuracy of calculations

Calculations in this software are performed in the framework of non-relativistic quantum mechanics. Leading order relativistic and field configuration corrections are O(10^–4) [1, 21], so results are not physically meaningful beyond the fourth significant figure. The computational accuracy — the agreement between computation and analytic theory — is well beyond this level (e.g. Figure 1). So, results can be considered as exact to the physical accuracy of the non-relativistic theory. This has been verified by comparison with existing analytic literature. Indeed there is good agreement with the results of Adhikari et. al. [23], Klarsfeld [24], and Heno et. al. [25]. Functions to assess the computational accuracy are available in the additional module “Convergence_and_correctness_tests.py”.

Quality control

There is no general error handling for the input of forbidden values. Some cases are handled (such as trying to calculate the scattering rate on an energetically or dipole forbidden transition), but otherwise, it is the responsibility of the user to ensure that any attempted values are physically meaningful and conform to the requirements of the non-relativistic theory (such as only integer values of n and l).

It is the user’s responsibility to ensure that they are using a sufficient basis set for their purposes. The best assessment of this is to check that calculated values do not vary significantly with small changes to the basis parameters. Figure 2 is an example of such an assessment. It consists of heat maps showing the fractional variation of certain calculated quantities as n_max varies about 300 and k varies about 0.3. In this case, it is clear that the choice of basis is appropriate for desired calculations. The functions used to produce this figure are made available in the additional module “Convergence_and_correctness_tests.py”.

Figure 2

Operating system

This software was developed on Windows 10 and has since been tested on Windows 11.

Programming language

This software was developed in Python 3.7.3 but has also been tested in Python 3.13.7. Any version of Python 3 after 3.7 will be compatible if compatible versions of NumPy and SciPy are available.

Additional system requirements

The only requirement is the ability to run and edit .py files. The entire software is 92kB, whilst the central module “main_calculation_functions.py” is 16kB.

Dependencies

The primary module “main_calculation_functions.py” requires “numpy 1.19.5” and “scipy 1.6.0” only. All other modules require “matplotlib 3.1.0”, whilst “Convergence_and_correctness_tests.py” also requires “seaborn 0.11.1”. both it and “Compare_scattering_rates.py” require “pandas 1.3.0”.

The versions listed above were used in the development of the software, compatibility has since been tested with: “numpy 2.3.3”, “scipy 1.16.2”, “matplotlib 3.10.6”, “seaborn 0.13.2”, and “pandas 2.3.3”.

List of contributors

We thank Dr. Robert Potvliege of the Department of Physics at Durham University, who introduced us to the calculation method and assisted with the underlying mathematics.

Software location

Archive

• Name: Hydrogen-S-state-2-photon

• Persistent identifier: https://doi.org/10.5281/zenodo.10034256

• Licence: CC0-1.0

• Publisher: Joseph P. Scott

• Version published: v2.0

• Date published: 25/10/25

Code repository

GitHub

• Name: Hydrogen-S-state-2-photon

• Persistent identifier: https://github.com/JosephPScott/Hydrogen-S-state-2-photon/releases/tag/v2.0

• Licence: CC0-1.0

• Date published: 25/10/25

Language

English