Emulation of sequences: AFM state sequence¶

This example is taken from Pulser documentation.

To illustrate the simulation of sequences, let us study a simple one-dimensional system with periodic boundary conditions (a ring of atoms):

In [11]:

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

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

In [12]:

# Setup
L = 10

Omega_max = 2.3 * 2 * np.pi
U = Omega_max / 2.3

delta_0 = -3 * U
delta_f = 1 * U

t_rise = 2000
t_fall = 2000
t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 5000

# Define a ring of atoms distanced by a blockade radius distance:
R_interatomic = pulser.MockDevice.rydberg_blockade_radius(U)
coords = (
    R_interatomic
    / (2 * np.sin(np.pi / L))
    * np.array(
        [
            (np.cos(theta * 2 * np.pi / L), np.sin(theta * 2 * np.pi / L))
            for theta in range(L)
        ]
    )
)

# ring, periodic register
reg = pulser.Register.from_coordinates(coords, prefix="q")
# or try open boundaries
#reg = pulser.Register.rectangle(1,L,spacing=R_interatomic/1.2, prefix="q")

reg.draw(blockade_radius=R_interatomic, draw_half_radius=True, draw_graph=True)

# Setup L = 10 Omega_max = 2.3 * 2 * np.pi U = Omega_max / 2.3 delta_0 = -3 * U delta_f = 1 * U t_rise = 2000 t_fall = 2000 t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 5000 # Define a ring of atoms distanced by a blockade radius distance: R_interatomic = pulser.MockDevice.rydberg_blockade_radius(U) coords = ( R_interatomic / (2 * np.sin(np.pi / L)) * np.array( [ (np.cos(theta * 2 * np.pi / L), np.sin(theta * 2 * np.pi / L)) for theta in range(L) ] ) ) # ring, periodic register reg = pulser.Register.from_coordinates(coords, prefix="q") # or try open boundaries #reg = pulser.Register.rectangle(1,L,spacing=R_interatomic/1.2, prefix="q") reg.draw(blockade_radius=R_interatomic, draw_half_radius=True, draw_graph=True)

We use the drawing capabilities of the Register class to highlight the area half the blockade radius away from each atom, which makes it so that overlapping circles correspond to interacting atoms. This is further fleshed out by the graph edges drawn using the draw_graph option.

In this register, we shall act with the following pulser.Sequence, which is designed to reach a state with antiferromagnetic order:

In [13]:

rise = pulser.Pulse.ConstantDetuning(
    pulser.RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0
)
sweep = pulser.Pulse.ConstantAmplitude(
    Omega_max, pulser.RampWaveform(t_sweep, delta_0, delta_f), 0.0
)
fall = pulser.Pulse.ConstantDetuning(
    pulser.RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0
)

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

seq.add(rise, "global")
seq.add(sweep, "global")
seq.add(fall, "global")

seq.draw()

rise = pulser.Pulse.ConstantDetuning( pulser.RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0 ) sweep = pulser.Pulse.ConstantAmplitude( Omega_max, pulser.RampWaveform(t_sweep, delta_0, delta_f), 0.0 ) fall = pulser.Pulse.ConstantDetuning( pulser.RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0 ) seq = pulser.Sequence(reg, pulser.MockDevice) seq.declare_channel("global", "rydberg_global") seq.add(rise, "global") seq.add(sweep, "global") seq.add(fall, "global") seq.draw()

1. Running an emulation¶

First, we initialize the state vector backend:

In [14]:

from emu_sv import SVBackend
bknd = SVBackend()

from emu_sv import SVBackend bknd = SVBackend()

Then, we need to initialize the desired observables, their evaluation times, and the basis. In this example we chose the final state and the magnetization profile.

In [15]:

from emu_sv import StateResult, SVConfig, QubitDensity

dt = 10  # timestep in ns
seq_duration = seq.get_duration()

eval_times = list(range(0, seq_duration + dt, dt))
#magnetization  is given by -1+2*qubit_density
density = QubitDensity(evaluation_times=eval_times, basis=("r","g"), nqubits=L)
state = StateResult(evaluation_times=[seq_duration])

config = SVConfig(dt=dt, observables = [density, state], log_level=2000)

from emu_sv import StateResult, SVConfig, QubitDensity dt = 10 # timestep in ns seq_duration = seq.get_duration() eval_times = list(range(0, seq_duration + dt, dt)) #magnetization is given by -1+2*qubit_density density = QubitDensity(evaluation_times=eval_times, basis=("r","g"), nqubits=L) state = StateResult(evaluation_times=[seq_duration]) config = SVConfig(dt=dt, observables = [density, state], log_level=2000)

To run the simulation we simply use the method run() of the backend, passing the sequence we wish to emulate and the config file. It returns a Results object, which will allow the study or post-processing of the states for each time step in our simulation.

In [16]:

results = bknd.run(seq, config)

results = bknd.run(seq, config)

2. Using the Results¶

The Results object returned at the end of an emulation run, is a dictionary that contains the observables defined above, at each time step. We can access them using their name and time. The Results object that we created contains the final quantum state for example:

In [17]:

results["state"][seq_duration].vector[:10]

results["state"][seq_duration].vector[:10]

Out[17]:

tensor([-2.2128e-05+0.0002j,  8.0325e-04-0.0002j,  8.0325e-04-0.0002j,
        -7.4170e-04+0.0004j,  8.0325e-04-0.0002j, -1.6267e-03-0.0034j,
        -7.4170e-04+0.0004j,  7.1285e-04-0.0003j,  8.0325e-04-0.0002j,
        -2.7328e-03-0.0029j], dtype=torch.complex128)

In [18]:

#getting the observables names
results.get_result_names()

#getting the observables names results.get_result_names()

Out[18]:

['qubit_density', 'state']

We can sample the final state directly, using its sample() method from the Results object. For example, we sample it and discard the less frequent bitstrings:

In [19]:

final_state = results["state"][seq_duration]
counts = final_state.sample(1000)

large_counts = {k: v for k, v in counts.items() if v > 5}

plt.figure(figsize=(15, 4))
plt.xticks(rotation=90, fontsize=14)
plt.title("Most frequent observations")
plt.bar(large_counts.keys(), large_counts.values())

final_state = results["state"][seq_duration] counts = final_state.sample(1000) large_counts = {k: v for k, v in counts.items() if v > 5} plt.figure(figsize=(15, 4)) plt.xticks(rotation=90, fontsize=14) plt.title("Most frequent observations") plt.bar(large_counts.keys(), large_counts.values())

Out[19]:

<BarContainer object of 67 artists>

Notice how the most frequent bitstrings correspond to the antiferromagnetic order states.

We can also get the expectation values of operators for the states in the evolution. Since we asked for the qubit density (the Rydberg state population) in the config, we can plot its time evolution on each atom.

In [20]:

times = np.array(eval_times)
mag = np.stack([-1.0+2.0*np.array(results["qubit_density"][t]) for t in eval_times[1:]], axis=1)

fig, axs = plt.subplots(1,2, figsize=(8,4))

for i in range(L):
    axs[0].plot(eval_times[1:], mag[i,:])
axs[0].set_xlabel("time [ns]")
axs[0].set_ylabel(r"$\langle Z_i\rangle$")

pc = axs[1].pcolor(eval_times[1:], range(L), mag)
axs[1].set_xlabel("time [ns]")
axs[1].set_ylabel("qubit")
fig.colorbar(pc, ax=axs[1], label = r"$\langle Z_i\rangle$")

times = np.array(eval_times) mag = np.stack([-1.0+2.0*np.array(results["qubit_density"][t]) for t in eval_times[1:]], axis=1) fig, axs = plt.subplots(1,2, figsize=(8,4)) for i in range(L): axs[0].plot(eval_times[1:], mag[i,:]) axs[0].set_xlabel("time [ns]") axs[0].set_ylabel(r"$\langle Z_i\rangle$") pc = axs[1].pcolor(eval_times[1:], range(L), mag) axs[1].set_xlabel("time [ns]") axs[1].set_ylabel("qubit") fig.colorbar(pc, ax=axs[1], label = r"$\langle Z_i\rangle$")

Out[20]:

<matplotlib.colorbar.Colorbar at 0x71bf4ac1a830>

Notice how the local qubit density on each atom goes in the same way from 0 (which corresponds to the ground state) to 1/2. This is expected since as we saw above, the state after the evolution has antiferromagnetic-order, so at each site, there is a compensation of magnetization. The parity (even) and the boundary conditions (periodic) allow for two lowest-energy states, whose superposition is similar to that of the perfectly antiferromagnetic state: $\left( |rgrg\dots\rangle+|grgr\dots\rangle \right)/\sqrt{2}$