# 2.7. Rydberg Pair Potentials Near Surfaces¶

This tutorial is based on results that were published by J. Block and S. Scheel, “van der Waals interaction potential between Rydberg atoms near surfaces” Phys. Rev. A 96, 062509 (2017). We will reproduce the pair potentials shown in Figure 4. The final result is that for states around the $$|70p_{3/2};70p_{3/2}\rangle$$-asymptote of Rubidium the strength of the pair interaction is reduced when placing the atoms in front of a perfect mirror (perfectly conducting plate) compared to the vacuum interaction.

As described in the introduction, we start our code with some preparations and load the necessary modules.

%matplotlib inline

# Arrays
import numpy as np

# Plotting
import matplotlib.pyplot as plt
from itertools import product

# Operating system interfaces
import os, sys

# Parallel computing
from multiprocessing import Pool

# 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")


The plate lies in the $$xy$$-plane with the surface at $$z = 0$$. The atoms lie in the $$xz$$-plane with $$z>0$$.

We can set the angle between the interatomic axis and the z-axis theta and the center of mass distance from the surface distance_surface. distance_atom defines the interatomic distances for which the pair potential is plotted. The units of the respective quantities are given as comments.

Be careful: theta = np.pi/2 corresponds to horizontal alignment of the two atoms with respect to the surface. For different angles, large interatomic distances distance_atom might lead to one of the atoms being placed inside the plate. Make sure that distance_surface is larger than distance_atom*np.cos(theta)/2.

theta = np.pi/2 # rad
distance_atom = np.linspace(6, 1.5, 50) # µm
distance_surface = 1 # µm


Next we define the state that we are interested in using pairinteraction’s StateOne class . As shown in Figures 4 and 5 of Phys. Rev. A 96, 062509 (2017) we expect large changes for the $$C_6$$ coefficient of the $$|69s_{1/2},m_j=1/2;72s_{1/2},m_j=1/2\rangle$$ pair state, so this provides a good example.

We set up the one-atom system using restrictions of energy, main quantum number n and angular momentum l. This is done by means of the restrict... functions in SystemOne.

state_one1 = pi.StateOne("Rb", 69, 0, 0.5, 0.5)
state_one2 = pi.StateOne("Rb", 72, 0, 0.5, 0.5)

# Set up one-atom system
system_one = pi.SystemOne(state_one1.getSpecies(), cache)
system_one.restrictEnergy(min(state_one1.getEnergy(),state_one2.getEnergy()) - 50, \
max(state_one1.getEnergy(),state_one2.getEnergy()) + 50)
system_one.restrictN(min(state_one1.getN(),state_one2.getN()) - 2, \
max(state_one1.getN(),state_one2.getN()) + 2)
system_one.restrictL(min(state_one1.getL(),state_one2.getL()) - 2, \
max(state_one1.getL(),state_one2.getL()) + 2)


The pair state state_two is created from the one atom states state_one1 and state_one2 using the StateTwo class.

From the previously set up system_one we define system_two using SystemTwo class. This class also contains methods set.. to set angle, distance, surface distance and to enableGreenTensor in order implement a surface.

# Set up pair state
state_two = pi.StateTwo(state_one1, state_one2)

# Set up two-atom system
system_two = pi.SystemTwo(system_one, system_one, cache)
system_two.restrictEnergy(state_two.getEnergy() - 5, state_two.getEnergy() + 5)

system_two.setAngle(theta)
system_two.enableGreenTensor(True)
system_two.setDistance(distance_atom[0])
system_two.setSurfaceDistance(distance_surface)
system_two.buildInteraction()


Next, we diagonalize the system for the given interatomic distances in distance_atom and compare the free space system to a system at distance_surface away from the perfect mirror. The energy is calculated with respect to a Rubidium $$|70p_{3/2},m_j=3/2;70p_{3/2},m_j=3/2\rangle$$ two atom state, defined in energyzero.

# Diagonalize the two-atom system for different surface and interatomic distances
def getDiagonalizedSystems(distances):
system_two.setSurfaceDistance(distances[0])
system_two.setDistance(distances[1])
system_two.diagonalize(1e-3)
return system_two.getHamiltonian().diagonal()

if sys.platform != "win32":
with Pool() as pool:
energies = pool.map(getDiagonalizedSystems, product([1e12, distance_surface], distance_atom))
else:
energies = list(map(getDiagonalizedSystems, product([1e12, distance_surface], distance_atom)))

energyzero = pi.StateTwo(["Rb", "Rb"], [70, 70], [1, 1], [1.5, 1.5], [1.5, 1.5]).getEnergy()
y = np.array(energies).reshape(2, -1)-energyzero
x = np.repeat(distance_atom, system_two.getNumBasisvectors())

# Plot pair potentials
fig = plt.figure()
ax.set_xlabel(r"Distance ($\mu m$)")