Jupyter notebook main.ipynb

³⁷³ views

Kernel: SageMath 6.10

Quantum walks in curved spacetime

This Jupyter/SageMath worksheet implements some computations of the article Quantum walks in curved spacetime, from Pablo Arrighi, Stefano Facchini and Marcelo Forets. It is available as an arXiv preprint arXiv:1505.07023 . It is accepted for publication in Quantum Information Processing.

This document consists of a commented version of the raw SageMath worksheet, that can be downloaded here. First we tell ipython shell to interpret sage commands:

In [ ]:

%load_ext sage

In [ ]:

%display latex

1. Setup

General setup for the simulation, you don't need to modify this cell.

In [ ]:

import numpy as np
from scipy import linalg

var('t, x_d, t_d, k, q')

class FastMatrix:
    def __init__(self, matrix, **kwds):
        self._matrix = []

        self._width, self._height = matrix.dimensions()

        self._m = [[0 for j in xrange(self._width)] for i in xrange(self._height)]

        self._matrix = [[fast_callable(matrix[i,j], **kwds) for j in xrange(self._width)] for i in xrange(self._height)]

    def __call__(self, *args):
        for i in xrange(self._height):
            for j in xrange(self._width):
                self._m[i][j] = self._matrix[i][j](*args)
        return self._m
        

# Some useful matrices

SigmaX = matrix([[0,1],[1,0]])
SigmaZ = matrix([[1,0],[0,-1]])
X = SigmaX.tensor_product(matrix.identity(2))
Z = SigmaZ.tensor_product(matrix.identity(2))
H = matrix([[1,0,1,0],[0,1,0,1],[-1,0,1,0],[0,-1,0,1]])/sqrt(2)

Parameters of the simulation, please modify.

In [ ]:

# Boundaries of the simulation, expressed in continuous coordinates
XMINcont = 0
TMINcont = 0
XMAXcont = 10.0
TMAXcont = 10.0

# Number of points in the simulation lattice, for the X coordinate.
# Choose multiple of 4.
XMAX = 8000

In [ ]:

# The lattice spacing is calculated automatically
epsilon = (XMAXcont - XMINcont) / XMAX

# The number of points for the T coordinate is obtained accordingly
TMAX = (TMAXcont - TMINcont) / epsilon

epsilon;TMAX

General solution parameters. Please modify.

In [ ]:

# The spacetime metric is expressed in terms of the dyad fields, e^\mu_\nu.
# Because of symmetry, in order to specify the components of the dyads
# only three variables are needed: e00, e11, e10. They must be expressions
# of the variables x and t.

# In the following, we use the Scwartzschild metric with mass parameter
# given by the MASS variable.

MASS = 0.5
e00 = (1-2*MASS/x)^(-1/2)
e11 = (1-2*MASS/x)^(1/2)
e10 = 0

# The mass of the Dirac particle.
mass = 50.0

General solution. You don't need to modify this cell.

In [ ]:

# Eigenvalues of B1
d_1 = -(e11/e00+e10)
d_2 = -(-e11/e00+e10)

# Definition of C
C = -SigmaX*mass/e00

E0 = matrix([
[sqrt(1+d_1), 0, -sqrt(1+d_1)*(1-d_1)-i*sqrt(1-d_1)*d_1, 0],
[0, sqrt(1+d_2), 0, -sqrt(1+d_2)*(1-d_2)-i*sqrt(1-d_2)*d_2],
[sqrt(1-d_1), 0, sqrt(1-d_1)*(1+d_1)+i*sqrt(1+d_1)*d_1, 0],
[0, sqrt(1-d_2), 0, sqrt(1-d_2)*(1+d_2)+i*sqrt(1+d_2)*d_2]
])/sqrt(2)

B = E0.H * Z * E0
B2 = B[2:4,0:2]
U = matrix.identity(2) + 2*B2.H
W0 = E0 * matrix.block_diagonal(matrix.identity(2), U) * E0.H * X

N = diff(E0,t).H * E0
M = E0.H * Z * diff(E0,x)

N1 = N[0:2,0:2]   # identically 0 for the above E0
N2 = N[2:4,0:2]
M1 = M[0:2,0:2]   # M1.H == M1 identically with above E0
M2 = M[2:4,0:2]

T2 = -2*N2 - 2*U*M2

T1 = 2 * I * C + M1.H - M1 - 2 * N1

T = matrix.block([ [T1, -T2.H * U], [T2, 0] ])

W1 = E0 * T * E0.H * X

Boundaries of the simulation. Please modify.

In [ ]:

# These functions specify, for each discrete time step, the minimal and
# the maximal value of the discrete X coordinate that we want to simulate.
# They are useful to avoid singularities, or regions where the metric
# is not defined, etc...
Xmin(t_d) = 2*MASS / epsilon
Xmax(t_d) = XMAX

# The discrete starting time of the simulation.
Tmin = 0

In [ ]:

EH_d = E0(t=epsilon*t_d+TMINcont, x=epsilon*x_d+XMINcont)*H.n()
Wtilde_d = epsilon * W0(t=epsilon*t_d+TMINcont, x=epsilon*x_d+XMINcont).H * W1(t=epsilon*t_d+TMINcont, x=epsilon*x_d+XMINcont)
W0_d = W0(t=epsilon*t_d+TMINcont, x=epsilon*x_d+XMINcont)

fastE = FastMatrix(EH_d, vars=[t_d,x_d], domain=CDF)
fastWtilde = FastMatrix(Wtilde_d, vars=[t_d,x_d], domain=CDF)
fastW0 = FastMatrix(W0_d, vars=[t_d,x_d], domain=CDF)

2. Initial condition

In the following cells we introduce the initial condition, given as Gaussian wavepacket centered around $k_0$ and $x_{\text{mean}}$ , for momentum and space, respectively. The spread in momentum is given by sigma.

In [ ]:

psi = np.zeros((TMAX, XMAX, 2), dtype=complex)
chi = np.zeros((TMAX, XMAX, 2), dtype=complex)
prob = np.zeros((TMAX, XMAX))

In [ ]:

k0 = 100

x_mean = XMINcont + (XMAXcont - XMINcont) / 2.0
sigma_min = (2*pi*sqrt(1/mass^2 + 1/k0^2)).n()

sigma = 4*sigma_min

In [ ]:

KMIN = -4/sigma+k0
KMAX = 4/sigma+k0


u_pos = vector((sqrt(1+k/sqrt(k^2+mass^2)), sqrt(1-k/sqrt(k^2+mass^2))))/sqrt(2.0)
u_neg = vector((-sqrt(1-k/sqrt(k^2+mass^2)), sqrt(1+k/sqrt(k^2+mass^2))))/sqrt(2.0)


phi_pos = e^(-sigma^2*(k-k0)^2/2)
phi_neg = e^(-sigma^2*(k-k0)^2/2)

f0re = (u_pos[0] * phi_pos + u_neg[0] * phi_neg) * cos(k*(x-x_mean))
f0im = (u_pos[0] * phi_pos + u_neg[0] * phi_neg) * sin(k*(x-x_mean))
f1re = (u_pos[1] * phi_pos + u_neg[1] * phi_neg) * cos(k*(x-x_mean))
f1im = (u_pos[1] * phi_pos + u_neg[1] * phi_neg) * sin(k*(x-x_mean))

# Functions to perform the Fourier transform numerically.
def F0re(x):
    return numerical_integral(f0re(x=x), KMIN, KMAX)[0]
    
def F0im(x):
    return numerical_integral(f0im(x=x), KMIN, KMAX)[0]

def F1re(x):
    return numerical_integral(f1re(x=x), KMIN, KMAX)[0]

def F1im(x):
    return numerical_integral(f1im(x=x), KMIN, KMAX)[0]

In [ ]:

plot(phi_pos, KMIN, KMAX)+plot(-phi_neg, KMIN, KMAX, color='red')

In [ ]:

start_time = walltime()

init_cond0 = vector([(F0re(epsilon*j+XMINcont)+I*F0im(epsilon*j+XMINcont)) for j in xrange(0,XMAX)])
init_cond1 = vector([(F1re(epsilon*j+XMINcont)+I*F1im(epsilon*j+XMINcont)) for j in xrange(0,XMAX)])

init_cond_norm = sqrt(init_cond0.norm()^2 + init_cond1.norm()^2)

init_cond0 = init_cond0 / init_cond_norm
init_cond1 = init_cond1 / init_cond_norm

walltime(start_time)

In [ ]:

# Here we initialize the first row of the simulation array with the initial condition.
chi[Tmin] = np.array([[init_cond0[j], init_cond1[j]]  for j in xrange(0,XMAX)])

In [ ]:

phi_pos_norm = sqrt(numerical_integral(phi_pos*phi_pos.conjugate(), KMIN, KMAX)[0])
if phi_pos_norm > 0:
    phi_pos_n = phi_pos/phi_pos_norm

    # Verify normalization
    print 'phi_pos_n: %f' % numerical_integral(phi_pos_n*phi_pos_n.conjugate(), KMIN, KMAX)[0]

    # Compute k_mean
    k_mean_pos = numerical_integral(k*(phi_pos_n.conjugate()*phi_pos_n)/sqrt(k^2+mass^2), KMIN, KMAX)[0]
    print 'k_mean_pos: %f' % k_mean_pos

phi_neg_norm = sqrt(numerical_integral(phi_neg*phi_neg.conjugate(), KMIN, KMAX)[0])
if phi_neg_norm > 0:
    phi_neg_n = phi_neg/phi_neg_norm

    # Verify normalization
    print 'phi_neg_n: %f' % numerical_integral(phi_neg_n*phi_neg_n.conjugate(), KMIN, KMAX)[0]

    # Compute k_mean
    k_mean_neg = numerical_integral(k*(phi_neg_n.conjugate()*phi_neg_n)/sqrt(k^2+mass^2), KMIN, KMAX)[0]
    print 'k_mean_neg: %f' % k_mean_neg

In [ ]:

k_mean = k_mean_pos

3. Simulation

This is main loop of the simulation. In general you don't need to modify it.

In [ ]:

start_time = walltime()

aggressive_optimization = True
opt_threshold = 1e-6

# Apply first layer of enconding
for n in xrange(0,XMAX,4):
    E_num = np.array(fastE(Tmin,n))
    v = np.dot(E_num, np.concatenate((chi[Tmin,n], chi[Tmin,n+2])))
    psi[Tmin,n] = v[:2]
    psi[Tmin,n+2] = v[2:]
    
    E_num = np.array(fastE(Tmin,n+1))
    v = np.dot(E_num, np.concatenate((chi[Tmin,n+1], chi[Tmin,n+3])))
    psi[Tmin,n+1] = v[:2]
    psi[Tmin,n+3] = v[2:]


# Main loop
for m in xrange (Tmin,TMAX-2,2):

    for n in xrange(XMAX):
        prob[m,n] = np.dot(chi[m,n].conjugate(), chi[m,n])

    print '%d/%d: %.3fs, norm: %.3f' % (m, TMAX, walltime(start_time), sum(prob[m]))

    offset = 2*((m/2+1) % 2)
    n_min = 4*(ceil(Xmin(t_d=m).n()/4)+1)
    n_max = min(4*floor(Xmax(t_d=m).n()/4), XMAX-2*offset)
    for n in xrange(n_min, n_max,4):
        #print("N:%s" % n)
        v_in = np.concatenate((psi[m,n+offset], psi[m,n+offset+2]))
        if aggressive_optimization and (np.dot(v_in.conj(), v_in) < opt_threshold):
            v = v_in
            w = v_in
        else:
            Wt_num = np.array(fastWtilde(m,n), dtype=complex)
            W_num = np.dot(np.array(fastW0(m,n), dtype=complex), linalg.expm(Wt_num))
            v = np.dot(W_num, v_in)
            E_num = np.array(fastE(m,n), dtype=complex)
            w = np.dot(E_num.conj().T, v)
        psi[m+2,n+offset] = v[:2]
        psi[m+2,n+offset+2] = v[2:]
        chi[m+2,n+offset] = w[:2]
        chi[m+2,n+offset+2] = w[2:]

        v_in = np.concatenate((psi[m,n+offset+1], psi[m,n+offset+3]))
        if aggressive_optimization and (np.dot(v_in.conj(), v_in) < opt_threshold):
            v = v_in
            w = v_in
        else:
            Wt_num = np.array(fastWtilde(m,n+1), dtype=complex)
            W_num = np.dot(np.array(fastW0(m,n+1), dtype=complex), linalg.expm(Wt_num))
            v = np.dot(W_num, v_in)
            E_num = np.array(fastE(m,n+1), dtype=complex)
            w = np.dot(E_num.conj().T, v)
        psi[m+2,n+offset+1] = v[:2]
        psi[m+2,n+offset+2+1] = v[2:]
        chi[m+2,n+offset+1] = w[:2]
        chi[m+2,n+offset+2+1] = w[2:]
        
walltime(start_time)

4. Plots

We generate some plots corresponding to the numerical simulation computed in the previous step.

In [ ]:

# We reverse the order of variables to have time on the vertical axis and space on horizontal axis
def prob_cont(x,t):
    t_disc = 2*floor(((t-TMINcont)/epsilon)/2)
    x_disc = floor((x-XMINcont)/epsilon)
    return (prob[t_disc][x_disc])

Plot boundary parameters. They can be modified to "zoom" in a particular region of the simulation.

In [ ]:

PlotXMIN = XMINcont
PlotXMAX = XMAXcont - epsilon
PlotTMIN = TMINcont
PlotTMAX = TMAXcont - epsilon

In [ ]:

start_time = walltime()

plot_points = round((PlotXMAX-PlotXMIN)/epsilon)
WalkPlot = contour_plot(prob_cont, (x,PlotXMIN,PlotXMAX), (t,PlotTMIN,PlotTMAX), cmap='jet', contours=50, plot_points=plot_points)

walltime(start_time)

In [ ]:

# Show the figure
WalkPlot

Quantum walks in curved spacetime

1. Setup

2. Initial condition

3. Simulation

4. Plots

Product

Resources

Company