This page was generated from the Jupyter notebook pair_potentials.ipynb. Open in Google Colab.

Pair Potentials

We calculate a pair potential plot for Rubidium \(|60S,m=1/2;60S,m=1/2\rangle\). For doing so, we start with calculating a system pi.SystemAtom for the first atom and a system for the second atom. Then, we combine these systems to a pair basis and construct the system of the pair of atoms pi.SystemPair. Our approach provides a lot of flexibility:

• We can use different atomic species for the first and the second atom.

• We can set arbitrary fields using pi.SystemAtom.set_electric_field and pi.SystemAtom.set_magnetic_field.

• We can set the distance between atoms arranged along the \(z\) axis with pi.SystemPair.set_distance or set an arbitrary distance vector with pi.SystemPair.set_distance_vector.

We start by importing libraries and initializing PairInteraction’s database.

%pip install -q matplotlib numpy pairinteraction

import matplotlib.pyplot as plt
import numpy as np
import pairinteraction as pi
from pairinteraction.visualization.colormaps import alphamagma

if pi.Database.get_global_database() is None:
    pi.Database.initialize_global_database(download_missing=True)

Afterwards, we construct a pi.BasisAtom and pi.SystemAtom. We use it to construct a pi.BasisPair and pi.SystemPair. The latter contains the Hamiltonian that describes the pair of atoms. By diagonalizing it, we obtain the eigenenergies that make up the pair potentials.

ket = pi.KetAtom("Rb", n=60, l=0, m=0.5)
basis = pi.BasisAtom("Rb", n=(56, 64), l=(0, 3))
print(f"Number of single-atom basis states: {basis.number_of_states}")
system = pi.SystemAtom(basis)

pair_energy = 2 * ket.get_energy(unit="GHz")
min_energy, max_energy = pair_energy - 10, pair_energy + 10

pair_basis = pi.BasisPair(
    [system, system], energy=(min_energy, max_energy), energy_unit="GHz", m=(1, 1)
)
print(f"Number of two-atom basis states: {pair_basis.number_of_states}")

distances = np.linspace(1, 4, 100)
pair_systems = [
    pi.SystemPair(pair_basis).set_distance(d, unit="micrometer") for d in distances
]

# Diagonalize the systems in parallel
pi.diagonalize(pair_systems, diagonalizer="eigen", sort_by_energy=True)

eigenenergies = [system.get_eigenenergies(unit="GHz") - pair_energy for system in pair_systems]
overlaps = [system.get_eigenbasis().get_overlaps([ket, ket]) for system in pair_systems]

Number of single-atom basis states: 288
Number of two-atom basis states: 547

We plot the pair potentials using matplotlib.

fig, ax = plt.subplots()

ax.set_xlabel(r"Distance [$\mu$m]")
ax.set_ylabel("Energy [GHz]")

try:
    ax.plot(distances, np.array(eigenenergies), c="0.9", lw=0.25, zorder=-10)
except ValueError:  # inhomogeneous shape -> no simple line plot possible
    for x, es in zip(distances, eigenenergies):
        ax.plot([x] * len(es), es, c="0.9", ls="None", marker=".", zorder=-10)

x_repeated = np.hstack([val * np.ones_like(es) for val, es in zip(distances, eigenenergies)])
energies_flattened = np.hstack(eigenenergies)
overlaps_flattened = np.hstack(overlaps)
sorter = np.argsort(overlaps_flattened)

scat = ax.scatter(
    x_repeated[sorter],
    energies_flattened[sorter],
    c=overlaps_flattened[sorter],
    s=15,
    vmin=0,
    vmax=1,
    cmap=alphamagma,
)

fig.colorbar(scat, ax=ax, label="Overlap with state of interest")

plt.show()