# 2.6. Dispersion Coefficients Near Surfaces¶

In this tutorial we reproduce the results depicted in Figure 5 from J. Block and S. Scheel “van der Waals interaction potential between Rydberg atoms near surfaces” Phys. Rev. A 96, 062509 (2017). We calculate the van der Waals \(C_6\)-coefficient between two Rubidium Rydberg atoms that are equidistantly placed in front of a perfect mirror (i.e. in horizontal alignment in front of a perfectly conducting plate). One finds that the relevant length scale is interatomic distance devided by distance from surface and that for decreasing surface distance the \(C_6\) coefficient is significantly reduced.

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

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_atoms = 10 # µm
distance_surface = np.linspace(
distance_atoms * np.abs(np.cos(theta)) / 2, 2 * distance_atoms, 30
) # µ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 changes of about 50% for the \(C_6\) coefficient of the
\(|69s_{1/2},m_j=1/2;72s_{1/2},m_j=1/2\rangle\) pair state of
Rubidium, 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()) - 30,
max(state_one1.getEnergy(), state_one2.getEnergy()) + 30,
)
system_one.restrictN(
min(state_one1.getN(), state_one2.getN()) - 3,
max(state_one1.getN(), state_one2.getN()) + 3,
)
system_one.restrictL(
min(state_one1.getL(), state_one2.getL()) - 1,
max(state_one1.getL(), state_one2.getL()) + 1,
)
```

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() - 3, state_two.getEnergy() + 3)
system_two.setAngle(theta)
system_two.setDistance(distance_atoms)
system_two.setSurfaceDistance(distance_surface[0])
system_two.enableGreenTensor(True)
system_two.buildInteraction()
```

We calculate the \(C_6\) coefficients. The `energyshift`

is given
by the difference between the interaction energy at given
`surface_distance`

and the unperturbed energy of the two atom state
`state_two.getEnergy()`

. The \(C_6\) coefficient is then given by
the product of `energyshift`

and `distance_atoms**6`

.

`idx`

is the index of the two atom state. The command
`getOverlap(state_two, 0, -theta, 0)`

rotates the quantisation axis of
`state_two`

by `theta`

around the y-axis. The rotation is given by
the Euler angles `(0, -theta, 0)`

in zyz convention. The negative sign
of theta is needed because the Euler angles used by pairinteraction
represent a rotation of the coordinate system. Thus, the quantisation
axis has to be rotated by the inverse angle.

```
# Calculate C6 coefficients
C6 = []
for d in distance_surface:
system_two.setSurfaceDistance(d)
system_two.diagonalize()
idx = np.argmax(system_two.getOverlap(state_two, 0, -theta, 0))
energyshift = system_two.getHamiltonian().diagonal()[idx] - state_two.getEnergy()
C6.append(energyshift * distance_atoms**6)
```

```
# Plot results
plt.plot(distance_surface / distance_atoms, np.abs(C6))
plt.xlim(min(distance_surface / distance_atoms), max(distance_surface / distance_atoms))
plt.xlabel("distance to surface / interatomic distance")
plt.ylabel(r"|C$_6$| (GHz $\mu m^6$)")
plt.show()
```