# 2.8. Rydberg Atom Ion Interaction¶

Here we look at the interacton between a single ion and a Rydberg atom. More details regarding the theory can be found in A. Duspayev et al., Phys. Rev. Research 3, 023114 (2021) or M. Deiß et al., Atoms 9, 2 (2021). The inernuclear axis is aligned with the z-direction and the interaction is expanded in multipolar terms. With this choice of coordinate system, the interaction potential is rotationally symmetric around the quantization axis and only states of same $$m_J$$ couple with each other, reducing the number of basis states.

We will reproduce here Fig. 2 (a) and (b) from Phys. Rev. Research 3, 023114 (2021), showing the influence of higher multipole order terms

## 2.8.1. Import the Library¶

# We call an IPython magic function to make the output of plotting commands displayed inline.
%matplotlib inline

# Arrays
import numpy as np

# Plotting
import matplotlib.pyplot as plt

# Operating system interfaces
import os

# pairinteraction :-)
from pairinteraction import pireal as pi

# Create cache for matrix elements
if not os.path.exists("./cache"):
os.makedirs("./cache")
cache = pi.MatrixElementCache("./cache")


## 2.8.2. Application: Potential Curves¶

The SystemOne class defines the Rydberg atom and the effect of the electric field of the ion can be included into the calculations. The charge of the ion in units of the elementary charge can be passed by SystemOne.setIonCharge(charge) and the distance between the ion and the Rydberg core with SystemOne.setRydIonDistance(distance) in units of micrometer. The orientatin is fixed such that the internuclear axis points along the z-axis.

As an example we show how to calculate the energy of the Rydberg state for different internuclear distances and including a different number of multipole orders. Choosing a maximum multipole order of one is equivalent to calculating a StarkMap in a homogenious electric field.

# Define the Rydberg state for which the interaction with the ion should be calculated
state = pi.StateOne("Rb", 45, 1, 1.5, 0.5)

# Setup system, considering states only with similar energy and principle quantum number.
# (if a high accuracy is required, the energy and principle quantum number ranges must be increased)
LowN = 42
HighN = 48
LowEnergy = -100
HighEnergy = 100
system = pi.SystemOne(state.getSpecies(), cache)
system.restrictEnergy(state.getEnergy() + LowEnergy, state.getEnergy() + HighEnergy)
system.restrictN(LowN, HighN)

# Since the ion and Rydberg atom are placed on the z-axis, the magnetic momentum is conserved
system.restrictM(0.5, 0.5)

# Define the charge of the ion in units of the elementary charge
system.setIonCharge(1)

# Define the maximum considered order of the multipole expansion.
# 1: monopole-dipole, 2: monopole-quadrupole, ...
system.setRydIonOrder(1)

# Loop over different distances to the ion
array_distances = np.linspace(1.5, 1.1, 100)  # um
array_eigenvalues = []
array_overlaps = []

for distance in array_distances:
# Set the ion Rydberg-atom distance in units of um
system.setRydIonDistance(distance)
# Diagonalize the system
system.diagonalize()
# Store the eigenenergies
array_eigenvalues.append(system.getHamiltonian().diagonal())
# Store the overlap of the eigenstates with the defined state
array_overlaps.append(system.getOverlap(state))

array_eigenvalues = np.ravel(array_eigenvalues)
array_overlaps = np.ravel(array_overlaps)
array_distances = np.repeat(array_distances, system.getNumBasisvectors())
array_eigenvalues = (
array_eigenvalues / 29.9792458
)  # Convert GHz into inverse cm wavenumber

# Plot the interaction potential, the color code visualizes overlap of eigenstates with defined state
plt.scatter(array_distances, array_eigenvalues, 8, array_overlaps)
plt.xlabel("Distance (um)")
plt.ylabel("Energy (1/cm)")
plt.colorbar(label=f"Overlap with {state}")
plt.title("Up to 1st order in the multipole expansion")
plt.grid()
plt.ylim(-61.5, -61.2)
plt.show()

# Calculate the same but including orders up to 6
system.setRydIonOrder(6)

# Loop over different distances to the ion
array_distances = np.linspace(1.5, 1.1, 100)  # um
array_eigenvalues = []
array_overlaps = []

for distance in array_distances:
# Set the ion Rydberg-atom distance in units of um
system.setRydIonDistance(distance)
# Diagonalize the system
system.diagonalize()
# Store the eigenenergies
array_eigenvalues.append(system.getHamiltonian().diagonal())
# Store the overlap of the eigenstates with the defined state
array_overlaps.append(system.getOverlap(state))

array_eigenvalues = np.ravel(array_eigenvalues)
array_overlaps = np.ravel(array_overlaps)
array_distances = np.repeat(array_distances, system.getNumBasisvectors())
array_eigenvalues = (
array_eigenvalues / 29.9792458
)  # Convert GHz into inverse cm wavenumber

# Plot the interaction potential, the color code visualizes overlap of eigenstates with defined state
plt.scatter(array_distances, array_eigenvalues, 8, array_overlaps)
plt.xlabel("Distance (um)")
plt.ylabel("Energy (1/cm)")
plt.colorbar(label=f"Overlap with {state}")
plt.title("Up to 6th order in the multipole expansion")
plt.grid()
plt.ylim(-61.5, -61.2);