# Building 1D Rydberg Crystals

The following notebook shows a study of many-body dynamics on a 1D system. It is based on 1707.04344. The authors of that paper studied the preparation of symmetry-breaking states in antiferromagnetic Ising chains, by tuning the interaction and driving parameters accross the phase diagram. In this notebook, we reproduce some results of this paper. Since this is a particular experiment not based on certified devices, we will use the MockDevice class to allow for a wide range of configuration settings.

[1]:

import numpy as np
import matplotlib.pyplot as plt
import qutip

from pulser import Pulse, Sequence, Register
from pulser_simulation import QutipEmulator
from pulser.waveforms import CompositeWaveform, RampWaveform, ConstantWaveform
from pulser.devices import MockDevice


## 1. Rydberg Blockade at Resonant Driving

We first consider clusters of $$1, 2$$ and $$3$$ atoms under resonant ($$\delta = 0$$) driving. If all the atoms are placed within each other’s blockade volume, only one excitation per group will be possible at a time. The Rabi frequency will be enhanced by $$\sqrt{N}$$

[2]:

def occupation(reg, j):
r = qutip.basis(2, 0)
N = len(reg.qubits)
prod = [qutip.qeye(2) for _ in range(N)]
prod[j] = r * r.dag()
return qutip.tensor(prod)


Given a value of the maximum Rabi Frequency applied to the atoms, we can calculate the corresponding blockade radius using the rydberg_blockade_radius() method from MockDevice. We use this to arrange clusters of atoms which will experience this blockade effect:

[3]:

Omega_max = 2 * 2 * np.pi
groups = 3

atom_coords = [
((R_blockade / N) * x + 40 * group, 0)
for group in range(groups)
for x in range(1, N + 1)
]
reg = Register.from_coordinates(atom_coords, prefix="q")
reg.draw(
)

resonant_pulse = Pulse.ConstantPulse(1500, Omega_max, 0.0, 0.0)

seq = Sequence(reg, MockDevice)
seq.declare_channel("ising", "rydberg_global")

simul = QutipEmulator.from_sequence(seq, sampling_rate=0.2)

obs = [
sum(occupation(reg, j) for j in range(i, i + N))
for i in range(0, groups * N, N)
]

res = simul.run(progress_bar=True, method="bdf")
return res.expect(obs)

Blockade Radius is: 8.692279598222772µm.


Next we run blockade_cluster(N), which runs the simulation, for clusters of sizes $$N \in \{1,2,3\}$$:

[4]:

data = [blockade_cluster(N) for N in [1, 2, 3]]

10.0%. Run time:   0.01s. Est. time left: 00:00:00:00
20.0%. Run time:   0.02s. Est. time left: 00:00:00:00
30.0%. Run time:   0.02s. Est. time left: 00:00:00:00
40.0%. Run time:   0.03s. Est. time left: 00:00:00:00
50.0%. Run time:   0.04s. Est. time left: 00:00:00:00
60.0%. Run time:   0.05s. Est. time left: 00:00:00:00
70.0%. Run time:   0.07s. Est. time left: 00:00:00:00
80.0%. Run time:   0.07s. Est. time left: 00:00:00:00
90.0%. Run time:   0.08s. Est. time left: 00:00:00:00
Total run time:   0.08s

10.0%. Run time:   0.02s. Est. time left: 00:00:00:00
20.0%. Run time:   0.05s. Est. time left: 00:00:00:00
30.0%. Run time:   0.07s. Est. time left: 00:00:00:00
40.0%. Run time:   0.13s. Est. time left: 00:00:00:00
50.0%. Run time:   0.17s. Est. time left: 00:00:00:00
60.0%. Run time:   0.22s. Est. time left: 00:00:00:00
70.0%. Run time:   0.27s. Est. time left: 00:00:00:00
80.0%. Run time:   0.34s. Est. time left: 00:00:00:00
90.0%. Run time:   0.39s. Est. time left: 00:00:00:00
Total run time:   0.45s

10.0%. Run time:   1.02s. Est. time left: 00:00:00:09
20.0%. Run time:   1.92s. Est. time left: 00:00:00:07
30.0%. Run time:   2.74s. Est. time left: 00:00:00:06
40.0%. Run time:   3.70s. Est. time left: 00:00:00:05
50.0%. Run time:   4.74s. Est. time left: 00:00:00:04
60.0%. Run time:   5.74s. Est. time left: 00:00:00:03
70.0%. Run time:   7.16s. Est. time left: 00:00:00:03
80.0%. Run time:   8.58s. Est. time left: 00:00:00:02
90.0%. Run time:  10.25s. Est. time left: 00:00:00:01
Total run time:  12.29s


We now plot the probability that a Rydberg state withing the cluster is occupied (by summing the expectation values of the $$|r\rangle\langle r|_i$$ operators for each cluster) as it evolves in time, revealing the Rabi frequency of each configuration:

[5]:

fig, ax = plt.subplots(1, 3, figsize=(18, 3))
for N, expectation in enumerate(data):
ax[N].set_xlabel(r"Time ($µs$)", fontsize=10)
ax[N].set_ylabel(r"Probability of $|r\rangle$", fontsize=10)
ax[N].set_title(f"Atoms per cluster N = {N+1}", fontsize=12)
avg = sum(expectation) / groups
ax[N].plot(np.arange(len(avg)) / 1000, avg)


Only one excitation will be shared between the atoms on each cluster. Notice how the Rabi frequency increases by a factor of $$\sqrt{N}$$

## 2. Ordered Crystalline phases

The pulse sequence that we will prepare is based on the following parameters:

[6]:

# Parameters in rad/µs and ns
delta_0 = -6 * 2 * np.pi
delta_f = 10 * 2 * np.pi
Omega_max = 2 * 2 * np.pi
t_rise = 500
t_stop = 4500


[7]:

R_blockade = MockDevice.rydberg_blockade_radius(Omega_max)
a = 7.0

reg = Register.rectangle(1, 11, spacing=a, prefix="q")

Blockade Radius is: 8.692279598222772µm.


Create the pulses using Pulser objects:

[8]:

hold = ConstantWaveform(t_rise, delta_0)
excite = RampWaveform(t_stop - t_rise, delta_0, delta_f)
sweep = Pulse.ConstantAmplitude(
Omega_max, CompositeWaveform(hold, excite), 0.0
)

[9]:

seq = Sequence(reg, MockDevice)
seq.declare_channel("ising", "rydberg_global")

seq.draw()


The pulse sequence we just created corresponds a path in the Phase space of the ground state, which we represent schematically with the following function:

[10]:

def phase_diagram(seq):
ratio = []
for x in seq._schedule["ising"]:
if isinstance(x.type, Pulse):
ratio += list(x.type.detuning.samples / Omega_max)

fig, ax = plt.subplots()
ax.grid(True, which="both")

ax.set_ylabel(r"Interaction Range $R_b/a$", fontsize=14)
ax.set_xlabel(r"Detuning $\delta/\Omega$", fontsize=14)
ax.set_xlim(-4, 6)
ax.set_ylim(0, 3.2)
ax.axhline(y=0, color="k")
ax.axvline(x=0, color="k")

y = np.arange(0.0, 5, 0.01)
x = 2 * (0.6 + 8 * (y - 1.2) ** 2)
ax.fill_between(x, y, alpha=0.4)

y = np.arange(0.0, 5, 0.01)
x = 2 * (0.8 + 50 * (y - 2.45) ** 2)
ax.fill_between(x, y, alpha=0.4)

y = np.arange(0.0, 5, 0.01)
x = 2 * (1.0 + 170 * (y - 3.06) ** 2)
ax.fill_between(x, y, alpha=0.4)

ax.plot(np.array(ratio), np.full(len(ratio), R_blockade / a), "red", lw=2)
plt.show()

[11]:

phase_diagram(seq)


### 2.1 Simulation

We run our simulation, for a list of observables corresponding to $$|r\rangle \langle r|_j$$ for each atom in the register:

[12]:

simul = QutipEmulator.from_sequence(seq, sampling_rate=0.1)

occup_list = [occupation(reg, j) for j in range(len(reg.qubits))]

res = simul.run(progress_bar=True)
occupations = res.expect(occup_list)

10.0%. Run time:   0.76s. Est. time left: 00:00:00:06
20.0%. Run time:   1.61s. Est. time left: 00:00:00:06
30.0%. Run time:   2.55s. Est. time left: 00:00:00:05
40.0%. Run time:   3.27s. Est. time left: 00:00:00:04
50.0%. Run time:   3.86s. Est. time left: 00:00:00:03
60.0%. Run time:   4.40s. Est. time left: 00:00:00:02
70.0%. Run time:   4.99s. Est. time left: 00:00:00:02
80.0%. Run time:   5.60s. Est. time left: 00:00:00:01
90.0%. Run time:   6.35s. Est. time left: 00:00:00:00
Total run time:   7.21s


The following function plots the evolution of the expectation values with respect to time:

[13]:

def plot_evolution(results):
plt.figure(figsize=(10, 5))
plt.xlabel("Time (µs)", fontsize=14)
plt.ylabel("Rydberg Occupation Probabilty", fontsize=14)
for expv in results:
plt.plot(np.arange(len(expv)) / 1000, expv)

plt.show()

[14]:

plot_evolution(occupations)


We finally plot the probability of occupation of the Rydberg level with respect to the values of detuning, for each atom in the array:

[15]:

def heat_detuning(data, start, end):
N = len(reg.qubits)
time_window = []
x = []
detunings = simul.samples_obj.to_nested_dict()["Global"]["ground-rydberg"][
"det"
][[int(1000 * t) for t in simul.evaluation_times[:-1]]]

for t, d in enumerate(detunings):
if start <= d <= end:
time_window.append(t)
x.append(d / (2 * np.pi))

y = np.arange(1, N + 1)

X, Y = np.meshgrid(x, y)
Z = np.array(data)[:, time_window]

plt.figure(figsize=(14, 3))
plt.xlabel("Detuning/2π (MHz)", fontsize=14)
plt.ylabel("Atom in array", fontsize=14)
plt.yticks(range(1, N + 1), [f"q{i}" for i in range(N)], va="center")

plt.show()

[16]:

heat_detuning(occupations, delta_0, delta_f)


### 2.2 Rydberg Crystals: $$Z_3$$ Order

To arrive at a different phase, we reduce the interatomic distance $$a$$, thus increasing the interaction range between the atoms. This will lead to a $$Z_3$$ ordered phase:

[17]:

a = 3.5
reg = Register.rectangle(1, 10, spacing=a, prefix="q")

delta_0 = -4 * 2 * np.pi
delta_f = 10 * 2 * np.pi
Omega_max = 2.0 * 2 * np.pi  # btw 1.8-2 * 2pi MHz
t_rise = 600
t_stop = 2500

#
hold = ConstantWaveform(t_rise, delta_0)
excite = RampWaveform(t_stop - t_rise, delta_0, delta_f)
sweep = Pulse.ConstantAmplitude(
Omega_max, CompositeWaveform(hold, excite), 0.0
)
#
seq = Sequence(reg, MockDevice)
seq.declare_channel("ising", "rydberg_global")

phase_diagram(seq)

simul = QutipEmulator.from_sequence(seq, sampling_rate=0.1)

occup_list = [occupation(reg, j) for j in range(len(reg.qubits))]

res = simul.run(progress_bar=True, method="bdf")

occupations = res.expect(occup_list)

plot_evolution(occupations)
heat_detuning(occupations, delta_0, delta_f)

10.0%. Run time:   2.20s. Est. time left: 00:00:00:19
20.0%. Run time:   4.31s. Est. time left: 00:00:00:17
30.0%. Run time:   6.38s. Est. time left: 00:00:00:14
40.0%. Run time:   8.46s. Est. time left: 00:00:00:12
50.0%. Run time:  10.56s. Est. time left: 00:00:00:10
60.0%. Run time:  12.80s. Est. time left: 00:00:00:08
70.0%. Run time:  14.89s. Est. time left: 00:00:00:06
80.0%. Run time:  16.99s. Est. time left: 00:00:00:04
90.0%. Run time:  19.08s. Est. time left: 00:00:00:02
Total run time:  21.14s


### 2.3 Rydberg Crystals: $$Z_4$$ Order

Decreasing even more the interatomic distance leads to a $$Z_4$$ order. The magnitude of the Rydberg interaction with respect to that of the applied pulses means our solver has to control terms with a wider range, thus leading to longer simulation times:

[18]:

a = 2.8
reg = Register.rectangle(1, 9, spacing=a, prefix="q")

# Parameters in rad/µs and ns

delta_0 = -4 * 2 * np.pi
delta_f = 10 * 2 * np.pi
Omega_max = 2.0 * 2 * np.pi  # btw 1.8-2 2pi*MHz
t_rise = 600
t_stop = 2500

#
hold = ConstantWaveform(t_rise, delta_0)
excite = RampWaveform(t_stop - t_rise, delta_0, delta_f)
sweep = Pulse.ConstantAmplitude(
Omega_max, CompositeWaveform(hold, excite), 0.0
)

#
seq = Sequence(reg, MockDevice)
seq.declare_channel("ising", "rydberg_global")

phase_diagram(seq)

simul = QutipEmulator.from_sequence(seq, sampling_rate=0.4)

occup_list = [occupation(reg, j) for j in range(len(reg.qubits))]

#
res = simul.run(progress_bar=True, method="bdf")
occupations = res.expect(occup_list)

plot_evolution(occupations)
heat_detuning(occupations, delta_0, delta_f)

10.0%. Run time:   4.70s. Est. time left: 00:00:00:42
20.0%. Run time:  10.75s. Est. time left: 00:00:00:42
30.0%. Run time:  16.96s. Est. time left: 00:00:00:39
40.0%. Run time:  24.73s. Est. time left: 00:00:00:37
50.0%. Run time:  32.96s. Est. time left: 00:00:00:32
60.0%. Run time:  39.45s. Est. time left: 00:00:00:26
70.0%. Run time:  44.43s. Est. time left: 00:00:00:19
80.0%. Run time:  49.02s. Est. time left: 00:00:00:12
90.0%. Run time:  53.25s. Est. time left: 00:00:00:05
Total run time:  58.02s