Computational AMO and Quantum Chemistry

Computational AMO connects a physical model to quantities such as an atomic energy, a molecular vibration, or a photon-counting probability. The examples in this chapter show how to choose a representation and recognize whether a discrepancy comes from solving the equations, truncating their state space, or approximating the physics. Begin with a system whose answer you can calculate by hand, then change one assumption at a time.

Helpful background. Numerical mathematics explains matrix solvers and numerical integration. Computational quantum mechanics develops grids, basis sets, propagation, and convergence tests in detail. The examples below apply those ideas to AMO systems.

Enter this chapter

Choose the calculation closest to your question:

Lesson	What you will learn
Computational atomic structure	Separate radial, electron-correlation, and relativistic approximations.
Radial Schrödinger solvers	Turn origin and large-distance boundary conditions into a reliable bound-state calculation.
Electronic-structure methods	Choose a basis and a correlation method for the property you need.

For a first calculation, try variational helium or molecular orbitals, then move to the helium Hartree–Fock or vibrational-spectra lab. For dynamics, follow the driven two-level system, optical Bloch equations, and cavity QED labs. The laser-cooling lab adds motion and recoil; AMO calculation benchmarks brings these examples together.

What does a converged calculation mean?

Suppose a radial grid produces a matrix $H_N$. An eigensolver returns $\widetilde E$ and a normalized vector $\widetilde c$. The residual

$$r=H_N\widetilde c-\widetilde E\widetilde c$$

tests how accurately that matrix equation was solved. For Hermitian $H_N$, some matrix eigenvalue lies within $\lVert r\rVert$ of $\widetilde E$. A small residual alone does not identify which eigenvalue it is, guarantee an accurate eigenvector near a degeneracy, or show that the grid represents the continuum problem accurately. Eigenvector accuracy also depends on spectral separation; see Trefethen and Bau in the references.

Three distinct comparisons are therefore useful:

Comparison	Question answered
Tighten the solver tolerance with the same matrix.	Have I solved the finite problem accurately?
Refine the grid or basis and enlarge its spatial extent.	Is the chosen observable stable as the representation improves?
Add finite-mass, relativistic, correlation, or environmental effects.	Does the physical model include effects relevant to this observable?

For hydrogen, the infinite-nuclear-mass Coulomb energy $E_n=-Z^2E_{\mathrm h}/(2n^2)$ is exact within its model. Agreement with that formula does not test a fine-structure or isotope-shift calculation. Likewise, full configuration interaction solves the electron problem within a specified finite orbital space; it does not remove that space’s incompleteness.

Check the quantity you intend to use. Energies can converge well before polarizabilities, weak transition amplitudes, or wave-function tails. An absolute error is more informative than a relative error near a selection-rule zero. If two large energies are subtracted to obtain a small binding energy, study convergence of the difference itself.

For a self-adjoint Hamiltonian bounded below, the Rayleigh–Ritz principle gives $E_0\leq\langle\psi_T|H|\psi_T\rangle/\langle\psi_T|\psi_T\rangle$, with the expectation interpreted as a quadratic form when necessary. The nonzero trial state must lie in the form domain. This upper bound concerns the same Hamiltonian; inaccurate quadrature or a nonvariational approximation can spoil it. A small SCF energy change also need not mean self-consistency: check the orbital or density residual and stability of the solution.

Atomic Benchmarks

Atomic calculations provide a progression from analytic one-electron tests to correlated and relativistic many-electron structure. Johnson’s textbook develops this progression, including radial orbitals and transition amplitudes.

Hydrogenic radial equation

For a central potential, write

$$\psi_{nlm}(\mathbf r) = \frac{u_{nl}(r)}{r} Y_{lm}(\theta,\phi).$$

In atomic units, the radial Coulomb equation is

$$\left[ -\frac{1}{2} \frac{d^2}{dr^2} + \frac{l(l+1)}{2r^2} - \frac{Z}{r} \right] u_{nl}(r) = E_n u_{nl}(r),$$

with

$$u_{nl}(0)=0, \qquad u_{nl}(r)\to0 \quad \text{as } r\to\infty.$$

A radial solver should test:

• the $-Z^2/(2n^2)$ energy sequence;
• normalization $\int_0^\infty |u_{nl}(r)|^2\,dr=1$;
• node count $n-l-1$;
• near-origin behavior $u_{nl}\propto r^{l+1}$;
• exponential tail behavior;
• degeneracy in $l$ for the pure Coulomb Hamiltonian;
• the Coulomb virial identity $2\langle T\rangle+\langle V\rangle=0$;
• convergence with grid spacing, outer radius, and origin treatment.

Hydrogen reveals several common bugs at once: a missing radial Jacobian, an incorrect centrifugal factor, an origin singularity, a wrong sign in the potential, or normalization with Cartesian rather than radial weights.

Helium as a correlation benchmark

For an infinitely massive nucleus in atomic units,

$$H = \sum_{i=1}^{2} \left( -\frac{1}{2}\nabla_i^2 - \frac{2}{r_i} \right) + \frac{1}{r_{12}}.$$

Use a normalized product of hydrogenic $1s$ spatial orbitals with variational effective charge $\zeta$, multiplied by the normalized spin singlet. Its energy is

$$E(\zeta)/E_{\mathrm h} = \zeta^2 - \frac{27}{8}\zeta.$$

Stationarity gives

$$\zeta_\star = \frac{27}{16}, \qquad E(\zeta_\star) = -\frac{729}{256} E_{\mathrm h} = -2.84765625 E_{\mathrm h}.$$

This small calculation is useful because every integral is analytic and the variational inequality is known. Its error has two sources: restricting each orbital to a single exponential, and restricting the two-electron spatial state to a product. More flexible Hartree–Fock orbitals improve the energy before explicit electron correlation is introduced.

The explicitly correlated calculation of Sims and Hagstrom (2002) gives a nonrelativistic infinite-mass ground-state reference near

$$E_0 \approx -2.9037243770341196 E_{\mathrm h}.$$

The value is a benchmark for that Hamiltonian, not a prediction including finite nuclear mass, relativity, or QED. Helium Atom explains the variational derivation and the meaning of correlation energy.

Central-field and many-electron checks

For self-consistent atomic structure, useful checks include:

• one-electron ions, where electron–electron terms must vanish;
• two-electron closed shells, where restricted formulas are transparent;
• invariance under rotations within a fully occupied orbital subspace;
• orthonormality of orbitals;
• self-consistency residuals;
• total-energy stationarity;
• term, parity, and total-angular-momentum labels;
• agreement between length and velocity forms of transition amplitudes as a diagnostic, with due care about model consistency;
• convergence in radial grid, orbital set, configuration space, and partial waves.

Length–velocity agreement is not a theorem that the answer is physically exact. Two forms can agree because they share an approximation, and they can disagree when basis completeness or Hamiltonian consistency is insufficient.

Comparison with measured and evaluated spectra

The NIST Atomic Spectra Database provides critically compiled atomic levels, wavelengths, and transition data. A responsible comparison records:

• database version or access date;
• species, ionization stage, isotope when relevant, and level labels;
• observed versus Ritz value;
• vacuum versus standard-air wavelength convention;
• uncertainty and accuracy code;
• whether the entry is experimental, theoretical, or derived;
• the original source cited by the database.

Evaluated data validate a physical prediction only when the calculation and entry refer to the same quantity. Agreement of a nonrelativistic energy with a measured line can be accidental if omitted corrections cancel.

Molecular Benchmarks

Molecular computation adds moving nuclei, geometry, nonorthogonal basis functions, multiple electronic structures, and observables obtained from derivatives of an energy surface.

Helgaker, Jørgensen, and Olsen develop the electronic-structure methods used here. The first difficulty is often the basis metric, even before electron correlation is considered.

The generalized eigenproblem

For basis functions $|\chi_\mu\rangle$ that are not orthonormal,

$$|\phi_i\rangle = \sum_\mu C_{\mu i} |\chi_\mu\rangle.$$

Projection gives

$$FC = SC\varepsilon,$$

where

$$F_{\mu\nu} = \langle\chi_\mu|\widehat F|\chi_\nu\rangle, \qquad S_{\mu\nu} = \langle\chi_\mu|\chi_\nu\rangle.$$

The metric normalization is

$$C^\dagger S C = I.$$

A solver must not silently treat $S$ as the identity. Small eigenvalues of $S$ signal near-linear dependence and can amplify roundoff. Removing directions requires a documented threshold and a test of observable sensitivity.

Hydrogen molecular ion

The one-electron ion $\mathrm{H}_2^+$ is the molecular analogue of the hydrogen benchmark. For fixed internuclear separation $R$,

$$H_{\mathrm e}(R) = -\frac{1}{2}\nabla^2 - \frac{1}{r_A} - \frac{1}{r_B} + \frac{1}{R}.$$

It tests:

• two-center integrals;
• overlap and generalized diagonalization;
• gerade and ungerade symmetry;
• bonding and antibonding combinations;
• the dissociation and united-atom limits;
• the shape and minimum of a potential-energy curve;
• interpolation and differentiation of that curve.

At large $R$, the electronic state should approach a hydrogen atom plus a proton, with the chosen energy zero stated. At $R\to0$, the electronic part approaches the united-atom problem while the nuclear repulsion diverges. These limits help identify an incorrect energy zero or nuclear-repulsion term.

Basis and method are independent axes

A molecular result depends on both the electronic-structure approximation $M$ and finite basis $X$:

$$E_{M,X}.$$

A basis sequence at fixed method estimates

$$\delta_{\mathrm{basis}} \sim E_{M,X} - E_{M,\infty},$$

while a method sequence at a large basis probes a different error:

$$\delta_{\mathrm{method}} \sim E_{M,\infty} - E_{\mathrm{exact},\infty}.$$

Increasing basis size does not turn Hartree–Fock into a correlated method. Increasing correlation rank in a poor basis does not reach the complete-basis result. Report both axes.

Size consistency and dissociation

For noninteracting fragments at infinite separation, a size-consistent method should satisfy

$$E(A+B) \longrightarrow E(A)+E(B).$$

Truncated configuration interaction can violate this property even when its energy is variational. Restricted single-reference methods can also give a qualitatively wrong dissociation curve when one determinant ceases to describe the separated fragments.

The lesson is broader than one method: a near-equilibrium energy benchmark does not validate bond breaking, excited states, charge transfer, or multireference regions.

Basis-set superposition

In a dimer calculation, each fragment can borrow basis functions centered on the other fragment. The raw interaction energy

$$\Delta E_{\mathrm{int}} = E_{AB}^{AB} - E_A^A - E_B^B$$

then compares calculations in unequal one-particle spaces. A counterpoise comparison instead evaluates monomers in the dimer basis:

$$\Delta E_{\mathrm{int}}^{\mathrm{CP}} = E_{AB}^{AB} - E_A^{AB} - E_B^{AB}.$$

The superscript specifies the basis and the subscript the electrons and nuclei actually present. Hold the fragment geometries fixed at their dimer values for this interaction energy; relaxation energies are separate. Boys and Bernardi (1970) introduced this counterpoise comparison. It helps estimate basis-set superposition, while a basis sequence still tests the remaining error.

Vibrational benchmarks

Near a stationary geometry, expand the potential:

$$V(\mathbf R_0+\Delta\mathbf R) \approx V(\mathbf R_0) + \frac{1}{2} \Delta\mathbf R^{\mathsf T} K \Delta\mathbf R.$$

The mass-weighted Hessian is

$$D_{A\alpha,B\beta} = \frac{ K_{A\alpha,B\beta} }{ \sqrt{M_A M_B} }.$$

Its positive nonzero eigenvalues are $\omega_k^2$ in the harmonic approximation. A molecular-vibration implementation should check:

• three translational and three rotational zero modes for an isolated, unconstrained nonlinear molecule at a stationary geometry, or two rotational zero modes for a linear molecule;
• invariance under rigid translation and rotation;
• Hessian symmetry;
• finite-difference step convergence when derivatives are numerical;
• unit conversion from angular frequency to wavenumber;
• analytic oscillator cases;
• sensitivity to geometry and electronic-structure method.

Imaginary frequencies can identify a saddle point, but they can also result from an unconverged geometry, insufficient integration grid, noisy finite differences, or inconsistent constraints.

Quantum-Optical Simulations

Quantum optics supplies small models with strong analytic structure and large models where truncation and open-system dynamics become decisive.

Gardiner and Zoller explain the open-system models; Dalibard, Castin, and Mølmer give the quantum-trajectory construction. These equations can be implemented directly or with a toolbox such as QuTiP (Lambert et al., 2026).

Closed two-level dynamics

For a time-independent finite Hermitian Hamiltonian,

$$U(t) = \exp \left( -\frac{iHt}{\hbar} \right)$$

is unitary. A propagation benchmark should test:

$$U^\dagger U = I, \qquad \lVert\psi(t)\rVert = 1, \qquad \langle H\rangle_t = \langle H\rangle_0.$$

The last identity applies only when $H$ is time independent. For a driven Hamiltonian, energy exchange with the drive is expected, while norm conservation remains.

For an atom initially in $|g\rangle$, a constant real Rabi frequency $\Omega$, constant detuning $\Delta=\omega_L-\omega_0$, and no decay, the rotating-wave model gives

$P_e(t)=\frac{\Omega^2}{\Omega^2+\Delta^2} \sin^2\!\left(\frac{t}{2}\sqrt{\Omega^2+\Delta^2}\right).$

This tests frequency, amplitude, and basis ordering. Useful cases include:

• zero coupling, where populations remain fixed;
• resonance, where a $\pi$ pulse inverts the ideal system;
• large detuning, where excitation is suppressed;
• a drive phase change, which rotates the Bloch-sphere axis;
• comparison of rotating-wave and full driven Hamiltonians.

Optical Bloch steady state

Let $\Gamma_1>0$ be decay from the excited state to the ground state and $\Gamma_2=\Gamma_1/2+\gamma_\phi$ be coherence decay, with nonnegative pure dephasing rate $\gamma_\phi$. Take $\sigma_z=|e\rangle\langle e|-|g\rangle\langle g|$ and a real, constant drive. With Rabi frequency convention

$$H_{\mathrm{rot}} = -\frac{\hbar\Delta}{2}\sigma_z + \frac{\hbar\Omega}{2}\sigma_x,$$

the steady excited-state population can be written

$$\rho_{ee}^{(\mathrm{ss})} = \frac{s}{2(1+s)},$$

where

$$s = \frac{ \Omega^2\Gamma_2 }{ \Gamma_1 \left( \Delta^2+\Gamma_2^2 \right) }.$$

For pure radiative broadening, $\Gamma_2=\Gamma_1/2$. This population tests factors of two in the drive, decay conventions, and the saturation limit

$$\rho_{ee}^{(\mathrm{ss})} \longrightarrow \frac{1}{2}.$$

The population is even in $\Delta$, so it cannot detect a reversed detuning sign. Compare the complex optical coherence or a phase-sensitive transient to check that sign. Optical Bloch equations derives both population and coherence.

Master-equation structure

A Lindblad-form model is

$$\dot\rho = -\frac{i}{\hbar} [H,\rho] + \sum_k \left( L_k\rho L_k^\dagger - \frac{1}{2} \left\{ L_k^\dagger L_k, \rho \right\} \right).$$

The represented evolution should preserve trace and Hermiticity and, for a valid implementation of the stated generator, positivity. Numerical checks include:

$$\left| \operatorname{Tr}\rho-1 \right|, \qquad \lVert\rho-\rho^\dagger\rVert, \qquad \min_j\lambda_j(\rho).$$

A tiny negative eigenvalue can be roundoff; a persistent negative population can indicate an integrator, generator, or approximation problem. Clipping negative eigenvalues hides the symptom and changes the dynamics.

Quantum trajectories

An unraveling estimates the density operator by an ensemble:

$$\rho(t) \approx \frac{1}{M} \sum_{j=1}^{M} |\psi_j(t)\rangle \langle\psi_j(t)|.$$

The ensemble mean should converge to the master-equation result for the same generator. This cross-method comparison is stronger than either solver alone. Report confidence intervals and compare multiple seeds; individual trajectories are not expected to match deterministic density-matrix curves.

Jaynes–Cummings blocks

Choose the phase of the states so that $g\geq0$, and define the cavity detuning $\Delta_c=\omega_a-\omega_c$. For

$$H = \hbar\omega_c a^\dagger a + \frac{\hbar\omega_a}{2}\sigma_z + \hbar g \left( a^\dagger\sigma_- + a\sigma_+ \right),$$

the excitation number

$$N_{\mathrm{exc}} = a^\dagger a + \sigma_+\sigma_-$$

is conserved. In the manifold $\{|e,n\rangle,|g,n+1\rangle\}$,

$$E_{n,\pm} = \hbar\omega_c \left( n+\frac{1}{2} \right) \pm \frac{\hbar}{2} \sqrt{ \Delta_c^2 + 4g^2(n+1) }.$$

At resonance, the splitting is

$$E_{n,+}-E_{n,-} = 2\hbar g\sqrt{n+1}.$$

This gives an exact finite-block benchmark for matrix construction. A full-space simulation still needs a photon cutoff. Convergence should be checked through both observables and the occupation near the highest retained number state.

The finite photon ladder also changes an operator identity. On $|0\rangle,\ldots,|n_{\max}\rangle$, with the raising operator set to zero at the top,

$$[a,a^\dagger]=I-(n_{\max}+1)|n_{\max}\rangle\langle n_{\max}|.$$

Demanding an identity matrix here would incorrectly reject the properly truncated operators. Instead, keep occupied states well away from the boundary and test the observables as the cutoff increases.

AMO Platform Simulations

Platform calculations combine internal quantum dynamics with motion, electromagnetic fields, geometry, feedback, and noise. Their outputs are often effective parameters or measured signals rather than bare eigenvalues.

Low-saturation Doppler force

Take $k>0$, angular detuning $\Delta=\omega_L-\omega_0$, population linewidth $\Gamma$, and equal on-resonance saturation parameter $s_0\ll1$ for each beam. For two counterpropagating beams, the elementary two-level model in Metcalf and van der Straten gives

$$F(v) = \frac{\hbar k\Gamma s_0}{2} \Bigg[ \frac{1}{ 1+4(\Delta-kv)^2/\Gamma^2 } - \frac{1}{ 1+4(\Delta+kv)^2/\Gamma^2 } \Bigg].$$

For red detuning $\Delta<0$, the small-velocity slope should oppose motion:

$$F(v) \approx -\alpha v, \qquad \alpha>0.$$

This benchmark catches beam-direction, Doppler-shift, detuning, and recoil sign errors. It does not validate a sub-Doppler, multilevel, polarization gradient, or high-saturation model.

Separate internal and motional approximations

A semiclassical simulation may evolve

$$m\ddot{\mathbf r} = \mathbf F \left( \mathbf r, \mathbf v, \rho, t \right) + \boldsymbol\xi(t)$$

while a master equation evolves the internal state $\rho$. The assumptions behind classical motion, local steady state, recoil diffusion, and white noise are independent. State them separately.

For a stochastic recoil model, checks include:

• zero mean recoil in an isotropic emission model;
• the expected momentum-diffusion growth;
• convergence in trajectory count and time step;
• detailed treatment of absorption and emission directions;
• consistency between friction and diffusion in the regime where an equilibrium temperature is inferred.

Traps and fields

Near a stable trap center,

$$V(\mathbf r) \approx V_0 + \frac{m}{2} \sum_i \omega_i^2 x_i^2.$$

A field or trajectory solver should test:

• Maxwell constraints appropriate to the field representation;
• symmetry points and known zero-field or zero-force locations;
• Hessian eigenvalues and trap frequencies;
• energy conservation for static conservative motion;
• secular-frequency limits for driven traps;
• convergence in mesh, multipole order, and integration step;
• sensitivity to electrode, beam, and alignment parameters.

The trap depth, harmonic frequency, scattering rate, and coherence time are different observables and need not converge at the same rate.

Pulse envelopes and phase

For an exactly resonant isolated two-level system driven along a fixed Bloch sphere axis, a real signed envelope defines the rotation angle $\Theta=\int\Omega(t)\,dt$. If the atom starts in its ground state, the final excitation probability is $\sin^2(\Theta/2)$.

A changing optical phase rotates the instantaneous drive axis. Hamiltonians at different times then generally do not commute: the integral of a complex $\Omega(t)$ is not a rotation angle. Retain the full envelope and phase and use time-ordered propagation. Detuning, Stark shifts, leakage, and decoherence can also make pulses with the same scalar area produce different results.

Exercises

1. Classify four exactness claims

Classify each statement and give the missing qualification:

1. “The hydrogen ground-state energy is exact.”
2. “Dense diagonalization gives the exact molecular energy.”
3. “The helium value is exact to fifteen digits.”
4. “The Rabi curve is exact.”

Solution

1. The value $-Z^2E_{\mathrm h}/2$ is analytically exact for the specified nonrelativistic Coulomb Hamiltonian, nuclear-mass convention, and energy zero. It is not the complete physical energy of a real hydrogenic atom.
2. Dense diagonalization can solve a specified finite matrix to floating-point accuracy. The molecular model, one-particle basis, active space, and correlation approximation remain finite or approximate.
3. A high-precision helium value is a numerically converged reference for a stated Hamiltonian. The digits require an uncertainty and convergence argument; finite-mass, relativistic, and QED terms are separate.
4. The Rabi curve is exact for a declared effective two-level Hamiltonian and pulse convention. Level truncation, classical driving, rotating-wave treatment, and absence of decoherence are model assumptions.

The common repair is to say exact for what mathematical object and under which physical assumptions.

2. Reduce a generalized eigenproblem

Assume $F=F^\dagger$ and $S=S^\dagger>0$. Show how to transform

$$FC = SC\varepsilon$$

into an ordinary Hermitian eigenproblem, and state the numerical warning when $S$ has a very small eigenvalue.

Solution

Because $S$ is positive definite, it has a Hermitian square root. Write

$$S = U s U^\dagger, \qquad S^{-1/2} = U s^{-1/2}U^\dagger.$$

Set

$$C = S^{-1/2}Y.$$

Multiplying on the left by $S^{-1/2}$ gives

$$\left( S^{-1/2}FS^{-1/2} \right) Y = Y\varepsilon.$$

The transformed matrix is Hermitian. If $Y^\dagger Y=I$, then

$$C^\dagger S C = I.$$

A very small eigenvalue of $S$ makes $S^{-1/2}$ large and amplifies roundoff and integral errors. Removing that direction changes the represented space, so the threshold and observable sensitivity must be reported.

3. Diagnose a photon cutoff

A cavity simulation with cutoff $n_{\max}=10$ reports $P(n=10)=10^{-9}$ at the final time. Is cutoff convergence established? Design a stronger test.

Solution

No. Final-time boundary population can miss transient occupation, coherences with the boundary, or an observable that is sensitive to a small high-number tail.

A stronger test records the maximum boundary and near-boundary occupation over the full trajectory, then repeats the calculation for a sequence such as $n_{\max}=8,10,12,14$. Compare every claimed observable, trace, positivity, and conservation law where applicable. The cutoff is controlled only when the observables and diagnostics are stable within a predeclared tolerance.

4. Verify the Doppler-force sign

For the low-saturation force on this page, expand to first order in $v$ and show that red detuning produces damping.

Solution

Define

$$f(x) = \frac{1}{ 1+4x^2/\Gamma^2 }.$$

Then

$$F(v) = \frac{\hbar k\Gamma s_0}{2} \left[ f(\Delta-kv) - f(\Delta+kv) \right].$$

For small $v$,

$$f(\Delta\mp kv) \approx f(\Delta) \mp kv f'(\Delta).$$

Hence

$$F(v) \approx -\hbar k^2\Gamma s_0 f'(\Delta)v.$$

Since

$$f'(\Delta) = -\frac{ 8\Delta/\Gamma^2 }{ \left( 1+4\Delta^2/\Gamma^2 \right)^2 },$$

$f'(\Delta)>0$ for $\Delta<0$. Therefore

$$F(v) \approx -\alpha v, \qquad \alpha = \hbar k^2\Gamma s_0 f'(\Delta) >0.$$

This verifies the sign in the low-saturation two-level model, not the validity of that model for multilevel or sub-Doppler cooling.

References

• S. F. Boys and F. Bernardi, “The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors,” Molecular Physics 19, 553–566 (1970), doi:10.1080/00268977000101561.
• J. Dalibard, Y. Castin, and K. Mølmer, “Wave-function approach to dissipative processes in quantum optics,” Physical Review Letters 68, 580–583 (1992), doi:10.1103/PhysRevLett.68.580.
• C. W. Gardiner and P. Zoller, Quantum Noise, 3rd ed., Springer, 2004, publisher record.
• T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory, Wiley, 2000, doi:10.1002/9781119019572.
• W. R. Johnson, Atomic Structure Theory: Lectures on Atomic Physics, Springer, 2007, doi:10.1007/978-3-540-68013-0.
• A. Kramida, Yu. Ralchenko, J. Reader, and the NIST ASD Team, NIST Atomic Spectra Database, National Institute of Standards and Technology, version 5.12 (2024), accessed 2026-10-14.
• N. Lambert et al., “QuTiP 5: The Quantum Toolbox in Python,” Physics Reports 1153, 1–62 (2026), doi:10.1016/j.physrep.2025.10.001.
• H. J. Metcalf and P. van der Straten, Laser Cooling and Trapping, Springer, 1999, doi:10.1007/978-1-4612-1470-0.
• J. S. Sims and S. A. Hagstrom, “High precision Hy-CI variational calculations for the ground state of neutral helium and helium-like ions,” International Journal of Quantum Chemistry 90, 1600–1609 (2002), NIST publication record.
• L. N. Trefethen and D. Bau III, Numerical Linear Algebra, twenty-fifth anniversary ed., SIAM, 2022, doi:10.1137/1.9781611977165.