# Complete solution to the two-particle scattering problem on a line graph

Luna Lima e Silva  Instituto de Física, Universidade Federal Fluminense,
Niterói, RJ, 24210-340, Brazil  Daniel J. Brod  Instituto de Física,
Universidade Federal Fluminense, Niterói, RJ, 24210-340, Brazil

###### Abstract

Quantum walks have been used to develop quantum algorithms since their
inception, and can be seen as an alternative to the usual circuit model;
combining single-particle quantum walks on sparse graphs with two-particle
scattering on a line lattice is sufficient to perform universal quantum
computation. In this work we solve the problem of two-particle scattering on
the line lattice for a family of interactions without translation invariance,
recovering the Bose-Hubbard interaction as the limiting case. Due to its
generality, our systematic approach lays the groundwork to solve the more
general problem of multi-particle scattering on general graphs, which in turn
can enable design of different or simpler quantum gates and gadgets. As a
consequence of this work, we show that a CPHASE gate can be achieved with high
fidelity when the interaction acts only on a small portion of the line graph.

##  I Introduction

There are two well-known models of quantum walks, discrete- Ambainis  _et al._
(2001); Nayak and Vishwanath (2000) and continuous- Childs  _et al._ (2002)
time quantum walks. The former is closely related to classical random walks—it
is composed of a quantum particle and a quantum coin, and discrete time
unitary dynamics mimics a coin toss that determines the direction of a step
taken by the particle. Continuous-time quantum walks, on the other hand,
consist of a quantum particle that evolves by the Schrodinger equation, with a
Hamiltonian on a lattice or graph that is analogous to the Laplacian in free
space. This model of quantum walk is the focus of our work.

Quantum walks were first studied as a tool for building quantum algorithms.
They were shown to exponentially outperform their classical counterparts Farhi
and Gutmann (1998); Childs  _et al._ (2002, 2003), and an important milestone
was the discovery that they can be used to perform universal quantum
computation, first by Childs Childs (2009), later improved in collaboration
with Gosset and Webb Childs  _et al._ (2013).

In Childs  _et al._ (2013), the quantum gates and other gadgets operate mostly
in the regime where one particle is propagating far from any others,
essentially performing a single-particle quantum walk on a complicated graph.
Interactions are required only for two-qubit controlled-phase gates, where two
particles scatter off of each other on the simplest graph, a very long line,
via a translation-invariant Hamiltonian. This (approximate) translation
invariance implies conservation of the momentum which, together with energy
conservation, severely restricts the outcome of the scattering, as the momenta
of the particles are individually preserved.

In this work, we solve the more general problem of two-particle quantum walks
on the infinite line with an interaction that does not display translation
symmetry. More concretely, we consider two particles that propagate on an
infinite line graph, but only interact in a finite set of contiguous vertices
via the Bose-Hubbard Hamiltonian Valiente and Petrosyan (2008). Inspired by
scattering theory on continuous space, we employ the Lippmann-Schwinger
formalism Sakurai (1994) to find the S-matrix for this scenario. Crucially, in
contrast to previous works Childs and Gosset (2012); Varbanov and Brun (2009);
Childs  _et al._ (2013) that only compute the S-matrix for fixed particle
momenta, our approach explicitly treats the S-matrix as an operator on the
full Hilbert space, allowing us to describe momentum exchanges between the
particles.

After solving the finite-region scattering problem, we also investigate how
the solution converges to the previously-known results for the infinite line
Childs  _et al._ (2013). We analytically prove this convergence in the
asymptotic limit, and numerically investigate the rate of convergence.
Surprisingly, our results suggest that a few dozen interaction vertices might
already be sufficient to obtain a high-fidelity controlled-phase gate.

By not relying on translation symmetry, our results provide new tools to
analyse the scattering of particles on complex graphs, adding to the toolbox
of Childs  _et al._ (2013). In particular, we lay the groundwork to solve the
more general problem of multi-particle scattering on general graphs, which
might lead to better and more efficient quantum gates and gadgets for quantum
computation. We also note that multi-particle scattering theory has broad
applications, e.g. in condensed matter and particle physics, and its discrete-
space counterpart has been much less explored, which suggests applications for
our results beyond the field of quantum computing.

This work is organized as follows. We review some basic notions about
continuous time quantum walks in Sec. II, and we solve the scattering problem
in Sec. III.2. In Sec. III.3 we show analytically that in the limit of
interaction on all sites we recover known results, and in Sec. IV we show
numerical results indicating the corresponding rate of convergence.

##  II Single-particle scattering

Let us begin by briefly motivating the general setup for continuous-time
quantum walks Childs (2009); Childs  _et al._ (2013); Childs and Gosset
(2012). Single-particle quantum walks on a one dimensional lattice can be seen
as the discrete analogue of free-particle evolution on continuous space. The
latter is described by the Laplacian (up to some constants):

$H^{\textit{continuous
space}}=-\dfrac{\hbar^{2}}{2m}\nabla^{2}=-\dfrac{\hbar^{2}}{2m}\dfrac{\partial^{2}}{\partial
x^{2}}.$ (1)

To discretize the one dimensional Laplacian, we can look to the natural
definition of the derivative as a limit

$\dfrac{\partial^{2}f}{\partial x^{2}}=\lim_{\epsilon\rightarrow
0}\dfrac{[f(x+\epsilon)-f(x)]-[f(x)-f(x-\epsilon)]}{\epsilon^{2}}=\lim_{\epsilon\rightarrow
0}\dfrac{f(x+\epsilon)+f(x-\epsilon)-2f(x)}{\epsilon^{2}}.$ (2)

If we consider a particle on a one dimensional lattice, with spacing
$\epsilon>0$ and amplitudes $\braket{j}{\psi}=\psi(j)$ defined at every
lattice point $j\epsilon$, with $j\in\mathbb{Z}$, a natural choice for the
Hamiltonian is

$\bra{j}H^{\textit{lattice}}\ket{\psi}=-\dfrac{\hbar^{2}}{2m\epsilon^{2}}(\psi(j+1)+\psi(j-1)-2\psi(j)),$
(3)

or

$H^{\textit{lattice}}=-\dfrac{\hbar^{2}}{2m\epsilon^{2}}\sum_{j\in\mathbb{Z}}\ket{j}\bra{j+1}+\ket{j+1}\bra{j}-2\ket{j}\bra{j}.$
(4)

We can adopt some simplifications: this Hamiltonian has a term proportional to
the identity which adds an unimportant phase factor to the evolved state, so
we drop it. We choose units such that $\hbar^{2}/2m\epsilon^{2}=1$. Also,
following a convention of previous literature we choose the negative of this
Hamiltonian. With these changes, the Hamiltonian now reads

$H^{\textit{lattice}}=\sum_{j\in\mathbb{Z}}\ket{j}\bra{j+1}+\ket{j+1}\bra{j}.$
(5)

Notice that this is exactly the adjacency matrix for the infinite line graph,
i.e., a matrix whose $ij$ entry is $1$ if $\\{i,j\\}$ is an edge of the graph
and $0$ otherwise. The same procedure works for higher dimensional lattices,
which suggests that the adjacency matrix is the natural choice of Hamiltonian
for quantum walks on graphs.

Finding the corresponding eigenstates requires solving the Schrödinger
equation:

$H^{\textit{lattice}}\ket{\psi}=E\ket{\psi}.$ (6)

In the position basis, we have the following difference equation:

$\displaystyle\psi(j+1)+\psi(j-1)$ $\displaystyle=E\psi(j)$
$\displaystyle\forall j\in\mathbb{Z}.$ (7)

Considering the ansatz $\psi(j)=z^{j}$, we have a solution if $z+z^{-1}=E$. If
$|z|\neq 1$ the state is not normalizable, so we may consider $z=e^{-ip}$ with
$p\in[-\pi,\pi)$. From now on, the eigenstate with momentum $p$ is denoted by
$\ket{p}$. The energy of this state is $E=z+z^{-1}=2\cos{p}$ and its
amplitudes are given by

$\braket{j}{p}=\frac{1}{\sqrt{2\pi}}e^{-ijp},$ (8)

with normalization $\braket{k}{p}=\delta(k-p)$. Due to our choice of phase
convention, the velocity of a wavepacket with momentum centered on $p$ is
$\frac{dE}{d(-p)}=2\sin{p}$, which is positive for $p\in(0,\pi)$ and negative
for $p\in(-\pi,0)$.

##  III Two-particle scattering

###  III.1 Two-particle free states

The usual Hamiltonian describing two noninteracting particles evolving with
respect to Hamiltonians, say, $H_{A}$ and $H_{B}$, is the Kronecker sum
$H^{(2)}=H_{A}\otimes\mathds{1}+\mathds{1}\otimes H_{B}$. In this case, given
eigenstates for both particles $\ket{\psi}_{A}$ and $\ket{\phi}_{B}$ with
energies $E_{\psi}$ and $E_{\phi}$ respectively, the tensor product
$\ket{\psi}_{A}\otimes\ket{\phi}_{B}$ is a solution of the two-particle
Schrodinger equation, with energy $E=E_{\psi}+E_{\phi}$. In our case, for
distinguishable particles on a lattice, we write
$\ket{p_{1}}\ket{p_{2}}=\ket{p_{1}p_{2}^{\textit{dist}}}$ with amplitudes

$\displaystyle\braket{j_{1}j_{2}}{p_{1}p_{2}^{\textit{dist}}}=\frac{1}{2\pi}e^{-i(p_{1}j_{1}+p_{2}j_{2})}.$
(9)

For bosons and fermions we have, respectively111For the special case where
$p_{1}=p_{2}$, the correct states are
$\ket{p_{1}p_{2}^{\textit{boson}}}=\ket{p_{1}p_{2}^{\textit{dist}}}$ and
$\ket{p_{1}p_{2}^{\textit{fermion}}}=0$.

$\displaystyle\braket{j_{1}j_{2}}{p_{1}p_{2}^{\textit{boson}}}=\frac{1}{2\pi}\dfrac{e^{-i(p_{1}j_{1}+p_{2}j_{2})}+e^{-i(p_{2}j_{1}+p_{1}j_{2})}}{\sqrt{2}},$
(10)
$\displaystyle\braket{j_{1}j_{2}}{p_{1}p_{2}^{\textit{fermion}}}=\frac{1}{2\pi}\dfrac{e^{-i(p_{1}j_{1}+p_{2}j_{2})}-e^{-i(p_{2}j_{1}+p_{1}j_{2})}}{\sqrt{2}}.$

To avoid double counting, we always consider that $p_{1}<p_{2}$ for the
symmetric states. From now on, the notation $\ket{p_{1}p_{2}}$ will be used to
represent any of the states above when the distinction is unimportant.

For the next section, it is convenient to define a constant $b$ such that

$\displaystyle\braket{jj}{p_{1}p_{2}}=\frac{b}{2\pi}e^{-i(p_{1}+p_{2})j},$
(11)

where $b=1$ for distinguishable particles, $b=\sqrt{2}$ for bosons and $b=0$
for fermions.

###  III.2 Two-particle interacting quantum walk

We now consider the scattering of two particles propagating on an infinite
line with a Bose-Hubbard type interaction on $2L+1$ contiguous vertices,
namely:

$H=H^{(2)}+V:=\sum_{w=1}^{2}\sum_{j\in\mathbb{Z}}\left(\ket{j}\bra{j+1}_{w}+\ket{j+1}\bra{j}_{w}\right)+U\sum_{j=-L}^{L}\ket{jj}\bra{jj},$

where $\ket{a}\bra{b}_{1}\equiv\ket{a}\bra{b}\otimes\mathds{1}$,
$\ket{a}\bra{b}_{2}\equiv\mathds{1}\otimes\ket{a}\bra{b}$, so that $H^{(2)}$
is the Hamiltonian for two free particles on the line and $V$ is the
interaction, which is proportional to some constant $U$. In other words, the
particles interact with each other when they are on the same site, but only
for the central $2L+1$ sites. The Lippmann-Schwinger equation relates the
scattering states in the presence of an interaction to the free-particle
scattering states Sakurai (1994):

$\ket{{p_{1}p_{2}}^{\pm}}=\ket{p_{1}p_{2}}+(E-H^{(2)}\pm
i\epsilon)^{-1}V\ket{{p_{1}p_{2}}^{\pm}},$

where take the limit $\epsilon\rightarrow 0$ at the end of the computation.
The choice of positive (negative) sign corresponds to eigenstates of the
interacting Hamiltonian which in the distant past (future) asymptotically
converge to the free-particle states $\ket{p_{1}p_{2}}$.

We are interested in computing the S-matrix, defined by Weinberg (1995):

$\displaystyle
S_{k_{1}k_{2};p_{1}p_{2}}=\braket{k_{1}k_{2}^{-}}{p_{1}p_{2}^{+}},$ (12)

in other words, $S_{k_{1}k_{2};p_{1}p_{2}}$ is the probability amplitude of
preparing the state $\ket{p_{1}p_{2}}$ in the distant past and measuring the
state $\ket{k_{1}k_{2}}$ in the distant future. After some manipulations it is
possible to deduce a simpler expression for the S-matrix:

$\displaystyle
S_{k_{1}k_{2};p_{1}p_{2}}=\delta(k_{1}-p_{1})\delta(k_{2}-p_{2})-2\pi
i\delta(E_{k_{1}k_{2}}-E_{p_{1}p_{2}})\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}},$
(13)

where $E_{q_{1}q_{2}}=2\cos{q_{1}}+2\cos{q_{2}}$. To compute this expression,
due to the nature of the interaction, we only require the amplitudes when
particles are on the same site, namely, terms of the type
$\braket{jj}{p_{1}p_{2}}$. We can write in general, from the definition of $V$

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}=\frac{bU}{2\pi}\sum_{l=-L}^{L}e^{i(k_{1}+k_{2})l}\braket{ll}{{p_{1}p_{2}}^{+}}.$
(14)

The terms $\braket{ll}{{p_{1}p_{2}}^{+}}$ can be found by taking the inner
product of the Lippmann-Schwinger equation with $\bra{ll}$ for $-L\leq l\leq
L$:

$\displaystyle\braket{ll}{{p_{1}p_{2}}^{+}}=\frac{b}{2\pi}e^{-i(p_{1}+p_{2})l}+\frac{U}{4\pi^{2}}\sum_{m=-L}^{L}\braket{mm}{{p_{1}p_{2}}^{+}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\frac{e^{-i(k_{1}+k_{2})(l-m)}}{E-2\cos{k_{1}}-2\cos{k_{2}}+i\epsilon}dk_{1}dk_{2},$
(15)

where we used the resolution of the identity
$\mathds{1}=\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\ket{k_{1}k_{2}^{\textit{dist}}}\bra{k_{1}k_{2}^{\textit{dist}}}dk_{1}dk_{2}$.
Let us define the family of integrals (to be computed in Appendix A), and the
new variables $k_{\pm}=\frac{k_{1}\pm k_{2}}{2}$

$\displaystyle J(E,n)$ $\displaystyle=\dfrac{1}{2\pi}\lim_{\epsilon\rightarrow
0}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\frac{e^{-i(k_{1}+k_{2})n}}{E-2\cos{k_{1}}-2\cos{k_{2}}+i\epsilon}dk_{1}dk_{2}$
(16) $\displaystyle=\dfrac{1}{2\pi}\lim_{\epsilon\rightarrow
0}\int_{-\pi}^{\pi}e^{-2ik_{+}n}\int_{-\pi}^{\pi}\frac{1}{E-4\cos{k_{+}}\cos{k_{-}}+i\epsilon}dk_{-}dk_{+}.$

To find the S-matrix, we have to solve the system of $2L+1$ linear equations
(15). Simplifying notation, we can write

$\displaystyle\sum_{m=-L}^{L}(2\pi\delta_{lm}-UJ(E,l-m))x_{m}=b(c_{p_{1}p_{2}})_{l}\,,$
(17)

where $x_{m}=\braket{mm}{{p_{1}p_{2}}^{+}}$ and
$(c_{p_{1}p_{2}})_{l}:=e^{-i(p_{1}+p_{2})l}$. Note that $-L\leq m,l\leq L$.
The coefficients $t(n):=2\pi\delta_{n0}-UJ(E,n)$ form an $N\times N$ Toeplitz
matrix

$\displaystyle T_{N}=\begin{bmatrix}t(0)&t(-1)&t(-2)&\cdots&t(-(N-1))\\\
t(1)&t(0)&t(-1)&\cdots&t(-(N-2))\\\ t(2)&t(1)&t(0)&\cdots&t(-(N-3))\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
t(N-1)&t(N-2)&t(N-3)&\cdots&t(0)\end{bmatrix},$ (18)

such that Eq. (17) can be written as the matrix equation
$T_{N}x=bc_{p_{1}p_{2}}$. Looking back at Eq. (14), the factor related to the
interaction is

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}=\frac{bU}{2\pi}c_{k_{1}k_{2}}^{\dagger}x=\frac{b^{2}U}{2\pi}c_{k_{1}k_{2}}^{\dagger}T_{N}^{-1}c_{p_{1}p_{2}},$
(19)

and finally the S-matrix for finite $N$ is

$\displaystyle
S^{N}_{k_{1}k_{2};p_{1}p_{2}}=\delta(k_{1}-p_{1})\delta(k_{2}-p_{2})-ib^{2}Uc_{k_{1}k_{2}}^{\dagger}T_{N}^{-1}c_{p_{1}p_{2}}\delta(E_{k_{1}k_{2}}-E_{p_{1}p_{2}}).$
(20)

This expression shows that, in general, the particles can exchange momentum as
a result of scattering, for any finite region of interaction. We show in the
next section that, in the limit $L\rightarrow\infty$, that is no longer the
case and the momentum of each particle is individually conserved222For
distinguishable particles, they actually switch their momenta, but they still
cannot trade an arbitrary amount of momentum..

###  III.3 Asymptotic behaviour for $L\rightarrow\infty$

So far, our results are exact for the Bose-Hubbard interaction restricted to a
finite region. In this Section we explore the limit $L\rightarrow\infty$.

It is known Zhu and Wakin (2016) that Toeplitz matrices can be asymptotically
approximated by circulant matrices. One possible circulant approximation is
the matrix

$\displaystyle C_{N}=\begin{bmatrix}c(0)&c(1)&c(2)&\cdots&c(N-1)\\\
c(N-1)&c(0)&c(1)&\cdots&c(N-2)\\\ c(N-2)&c(N-1)&c(0)&\cdots&c(N-2)\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
c(1)&c(2)&c(3)&\cdots&c(0)\end{bmatrix},$ (21)

with coefficients given by

$\displaystyle c(n)=\begin{cases}t(0)\,&n=0,\\\
t(n)+t(n-N)\,&n=1,\cdots,N-1.\end{cases}$ (22)

Every circulant matrix is diagonalized by the unitary Fourier matrix
$F_{lm}=\frac{1}{\sqrt{N}}w^{-lm}$, where $w=e^{\frac{2\pi i}{N}}$ Gray
(2006). In our case, it is convenient to define both indices running from $-L$
to $L$. The eigenvalues of $C_{N}$ are

$(\lambda_{N})_{m}=\sum_{n=0}^{N-1}c(n)w^{nm},$

for $-L\leq m\leq L$, such that we can write
$C_{N}=F^{\dagger}diag(\lambda_{N})F$.

Now we can approximate Eq. (19) with a much simpler expression:

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}=\frac{b^{2}U}{2\pi}c_{k_{1}k_{2}}^{\dagger}T_{N}^{-1}c_{p_{1}p_{2}}\approx\frac{b^{2}U}{2\pi}c_{k_{1}k_{2}}^{\dagger}F^{\dagger}diag(\lambda_{N})^{-1}Fc_{p_{1}p_{2}}\text{,
as }N\rightarrow\infty.$ (23)

Let us obtain the asymptotic behavior of the $\lambda_{l}$. Noting that
$J(E,-n)=J(E,n)$, we have:

$\displaystyle(\lambda_{N})_{l}=2\pi-U\sum_{n=-2L}^{2L}J(E,n)w^{nl}.$ (24)

Using the definition in Eq. (16), switching the sum and the integral and then
taking the limit of large $N$, we can write

$\displaystyle(\lambda_{N})_{l}$
$\displaystyle=2\pi-\frac{U}{2}\int_{-\pi}^{\pi}\frac{1}{E-4\cos{\frac{\pi
l}{N}}\cos{k_{-}}+i\epsilon}dk_{-}-\frac{U}{2}\int_{-\pi}^{\pi}\frac{1}{E+4\cos{\frac{\pi
l}{N}}\cos{k_{-}}+i\epsilon}dk_{-},$ (25)

where $\frac{1}{T}\sum_{n\in\mathbb{Z}}e^{2\pi
i\frac{n}{T}}=\delta_{T}(x)=\sum_{n\in\mathbb{Z}}\delta(x-nT)$ is the Dirac
comb, the periodic version of the Dirac delta, and we used the identity
$2\delta_{T}(2x)=\delta_{T}(x)+\delta_{T}(x+\frac{T}{2})$, proved in Appendix
B. The computation of these integrals is left to Appendix A. The result, as
$\epsilon\rightarrow 0$, is

$\displaystyle(\lambda_{N})_{l}=2\pi-\dfrac{2\pi
U}{\sqrt{E^{2}-16\cos^{2}{\frac{\pi l}{N}}}},$ (26)

where the square root is defined to have positive imaginary part if the
radicand is negative.

Going back to Eq. (23), we have

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}=\frac{b^{2}U}{2\pi
N}\sum_{l,m,n=-L}^{L}w^{l(n-m)}e^{-i(k_{1}+k_{2})m}e^{i(p_{1}+p_{2})n}\left(2\pi-\frac{2\pi
U}{\sqrt{E^{2}-16\cos^{2}{\frac{\pi l}{N}}}}\right)^{-1}.$ (27)

Note that the variables $m$ and $n$ are uncoupled. Summing these variables in
the limit of large $L$, we have

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}=\frac{2\pi
b^{2}U}{N}\sum_{l=-L}^{L}\delta_{2\pi}\left(k_{1}+k_{2}+\frac{2\pi
l}{N}\right)\delta_{2\pi}\left(p_{1}+p_{2}+\frac{2\pi
l}{N}\right)\left(2\pi-\frac{2\pi U}{\sqrt{E^{2}-16\cos^{2}{\frac{\pi
l}{N}}}}\right)^{-1}.$ (28)

Changing the sum to an integral with variable $\ell=\frac{2\pi}{N}l$,
$d\ell=\frac{2\pi}{N}$ we have

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}$
$\displaystyle=b^{2}U\delta_{2\pi}(k_{1}+k_{2}-p_{1}-p_{2})\int_{-\pi}^{\pi}\delta_{2\pi}(p_{1}+p_{2}+\ell)\left(2\pi-\frac{2\pi
U}{\sqrt{E^{2}-16\cos^{2}{\frac{\ell}{2}}}}\right)^{-1}$ (29)
$\displaystyle=\dfrac{b^{2}U}{2\pi-\frac{2\pi
U}{\sqrt{E^{2}-16\cos^{2}{\frac{p_{1}+p_{2}}{2}}}}}\delta_{2\pi}(k_{1}+k_{2}-p_{1}-p_{2}).$

Since $E^{2}=16\cos^{2}{p_{+}}\cos^{2}{p_{-}}\leq 16\cos^{2}{p_{+}}$, where
$p_{\pm}=\frac{p_{1}\pm p_{2}}{2}$, and recalling that this square root was
defined to have positive imaginary part, we can write

$\displaystyle\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}=\dfrac{b^{2}U}{2\pi+\frac{\pi
iU}{|\sin{p_{1}}-\sin{p_{2}}|}}\delta_{2\pi}(k_{1}+k_{2}-p_{1}-p_{2}),$ (30)

where we have used the identity
$\sin{p_{1}}-\sin{p_{2}}=2\cos{p_{+}}\sin{p_{-}}$.

To compute the S-matrix of Eq. (13), we have to use the following identity,
proven in Appendix B:

$\displaystyle\delta(E_{k_{1}k_{2}}-E_{p_{1}p_{2}})$
$\displaystyle\delta_{2\pi}(k_{1}+k_{2}-p_{1}-p_{2})=\frac{1}{2|\sin
k_{1}-\sin
k_{2}|}(\delta_{2\pi}(k_{1}-p_{1})\delta_{2\pi}(k_{2}-p_{2})+\delta_{2\pi}(k_{1}-p_{2})\delta_{2\pi}(k_{2}-p_{1})).$
(31)

Finally, plugging this into Eq. (13), the S-matrix is

$\displaystyle S_{k_{1}k_{2};p_{1}p_{2}}$
$\displaystyle=\delta(k_{1}-p_{1})\delta(k_{2}-p_{2})-2\pi
i\delta(E_{k_{1}k_{2}}-E_{p_{1}p_{2}})\bra{k_{1}k_{2}}V\ket{{p_{1}p_{2}}^{+}}$
(32)
$\displaystyle=\delta(k_{1}-p_{1})\delta(k_{2}-p_{2})\left(1-\dfrac{2ib^{2}U}{4|\sin{p_{1}}-\sin{p_{2}}|+2iU}\right)-$
$\displaystyle-\delta(k_{1}-p_{2})\delta(k_{2}-p_{1})\left(\dfrac{2ib^{2}U}{4|\sin{p_{1}}-\sin{p_{2}}|+2iU}\right).$

This expression shows us that when we recover translation invariance, the
particles cannot exchange any amount of momentum; they either conserve both
momenta, or they exchange their momenta.

For distinguishable particles, we have $b=1$ and the reflection and
transmission coefficients given by, respectively:

$\displaystyle R(p_{1},p_{2})$
$\displaystyle=\dfrac{-iU}{2|\sin{p_{1}}-\sin{p_{2}}|+iU},$ (33)
$\displaystyle T(p_{1},p_{2})$
$\displaystyle=\dfrac{2|\sin{p_{1}}-\sin{p_{2}}|}{2|\sin{p_{1}}-\sin{p_{2}}|+iU}.$

For bosons we have $b=\sqrt{2}$, and since we considered ordered momenta for
the bosonic states, namely, $p_{1}<p_{2}$ and $k_{1}<k_{2}$, we are left with
the following S-matrix

$\displaystyle
S_{k_{1}k_{2};p_{1}p_{2}}=\delta(k_{1}-p_{1})\delta(k_{2}-p_{2})\dfrac{2|\sin{p_{1}}-\sin{p_{2}}|-iU}{2|\sin{p_{1}}-\sin{p_{2}}|+iU}.$
(34)

This expression coincides with the two-particle result on Valiente and
Petrosyan (2008) and in Appendix B of Childs  _et al._ (2013), without the
need to restrict the particles to be counterpropagating [i.e.
$p_{1}\in(-\pi,0)$ and $p_{2}\in(0,\pi)$], as in the latter. Convergence of
operators in infinite dimensional Hilbert spaces is not straightforward, so
the fact that we can recover the correct S-matrices by this limiting process
is not entirely obvious and is quite useful, as the limiting problem might be
much more difficult to solve.

In the next section we analyse two natural questions about this result; how to
measure the convergence and what is the rate of convergence.

##  IV Numerical analysis

In this section we numerically analyse the behaviour of the S-matrix for a
finite number of interaction sites. Up until Eq. (20), every computation was
exact for finite $L$, so that is our starting point.

It is not trivial to measure how close the finite-$L$ S-matrix from Eq. (20)
is to the limiting case described in Eq. (34). For example, let us consider
narrow wavepackets in momentum space; in the limiting case, they scatter off
of each other and acquire a phase factor depending on both momenta; on the
other hand, for any finite region of interaction, these wavepackets do not
acquire a phase after scattering, as described later. Thus, there always
exists states for which both S-matrices give consistently different results,
therefore the finite-case S-matrix does not converge to the asymptotic
S-matrix in operator norm, even though we showed they converge element-wise.

To avoid this issue, we take a more pragmatic approach and measure the
distance between the S-matrices only when acting on a subspace of states that
would arise in quantum walk based computation.

There are three families of states that might propagate in a long linear graph
in the universal construction of Childs  _et al._ (2013): there might be no
particles, a single particle or two particles present. Within the dual-rail
encoding used in Childs  _et al._ (2013), for the purposes of implementing a
CPHASE gate, the presence (or absence) of a particle in the linear graph is
equivalent to the logical state $\ket{1}$ (or $\ket{0}$).

The mapping now is as follows: the logical two-qubit state $\ket{00}$ is
represented by the absence of any particles. The logical states $\ket{01}$ and
$\ket{10}$ are represented by one particle moving to the left or to the right,
respectively—this choice is arbitrary. So the logical states $\ket{01}$ and
$\ket{10}$ are both represented by a particle on some state of the form

$\displaystyle\ket{k^{\sigma}}:=\dfrac{1}{\sqrt{2\sigma}}\int_{k-\sigma}^{k+\sigma}\ket{k^{\prime}}dk^{\prime},$
(35)

for momentum $k$ and width $2\sigma$. Since we want to deal with bosons, for
the case of both particles present, we chose the logical state $\ket{11}$ as
the (normalized) state

$\displaystyle\dfrac{\ket{k_{1}^{\sigma}}\ket{k_{2}^{\sigma}}+\ket{k_{2}^{\sigma}}\ket{k_{1}^{\sigma}}}{\sqrt{2+2|\braket{k_{1}^{\sigma}}{k_{2}^{\sigma}}|^{2}}}.$
(36)

For concreteness, we chose one particle with momentum centered around
$k_{1}=\frac{\pi}{4}$ and the other one centered around
$k_{2}=-\frac{\pi}{2}$, following the choice made in Childs  _et al._ (2013).
We also chose the coupling constant to be $U=2+\sqrt{2}$, so that, for the
eigenstate $\ket{p_{1}p_{2}^{\textit{boson}}}$, which is related to the limit
$\sigma\rightarrow 0$, the phase in equation (34) is $-i$. Thus, we expect to
recover this behaviour for large $L$ and small $\sigma$.

The average gate fidelity between the finite-case S-matrix and a CPHASE gate
with phase $\phi$ is computed as in Brod and Combes (2016):

$\displaystyle F(\phi):=\int
d\psi\bra{\psi}C_{\phi}^{\dagger}S_{N}\ket{\psi}\bra{\psi}S_{N}^{\dagger}C_{\phi}\ket{\psi},$
(37)

where $S_{N}$ is the S-matrix given in Eq. (13) for $N=2L+1$ interaction
sites, the integration is carried over the two-qubit state space with the Haar
measure and $C_{\phi}$ is the ideal two qubit gate

$\displaystyle\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&e^{i\phi}\end{bmatrix}.$ (38)

Similarly, we can compute the fidelity between the finite-case S-matrix and
the limit S-matrix directly, but as we can see in Figure 1, it gives us the
same qualitative results, when we fix $\phi=-\frac{\pi}{2}$. Therefore, we are
going to focus our analysis on the fidelity described in Eq. (37), as it gives
us more flexibility by allowing us to choose different values for the
parameter $\phi$.

Figure 1: Solid line: Fidelity between $S_{N}$ and the CPHASE gate with the
ideal phase $\phi=-\frac{\pi}{2}$ as a function of the wavepacket half width
$\sigma$. Dashed line: Fidelity between $S_{N}$ and the limit S-matrix as a
function of $\sigma$.

Using the fact that the S-matrix acts as the identity if there are less than
two particles present we can write the average gate fidelity as

$\displaystyle
F(\phi)=\dfrac{1}{10}(6+3Re(e^{-i\phi}\mathcal{F})+|\mathcal{F}|^{2}),$ (39)

where $\mathcal{F}=\bra{11}S_{N}\ket{11}$.

It is possible to show that, in the limit $\sigma\rightarrow 0$ for any finite
$L$, the S-matrix acts as the identity on the wavepackets considered. This can
be interpreted as follows: since the interaction is proportional to the
amplitude of both particles being at the same site, which vanishes for
$\sigma\rightarrow 0$, any finite number of vertices of interaction is not
enough to produce an effect asymptotically. As a consequence, the fidelity for
the ideal phase $\phi=-\frac{\pi}{2}$ converges to $7/10$ as
$\sigma\rightarrow 0$, seen in Fig. 1. Also, as we can observe in both graphs
of Fig. 2, for small $\sigma$, the best phase that approximates the result of
the interaction is $\phi=0$, with fidelity close to one.

Figure 2: Left: Maximized fidelity over all phases as a function of $\sigma$.
Right: Corresponding phase $\phi$ for which the fidelity is maximized.

Nevertheless, we observe that the fidelity for the phase $\phi=-\frac{\pi}{2}$
for small $\sigma$ increases for larger values of $L$, which tells us that in
the limit $L\rightarrow\infty$, narrow momentum wavepackets scatter with the
correct phase.

Since the S-matrix acts as the identity for small $\sigma$, in other words, as
the $C_{\phi}$ gate with $\phi=0$, we can ask ourselves what is the target
value of $\phi$ for which the fidelity is maximum, for any given value of $L$
and $\sigma$. In Fig. 2 we can see that, in fact, except for very small
$\sigma$, the best phase approaches $\phi=-\dfrac{\pi}{2}$ for large $L$.

From Fig. 1, we can see that the maximal fidelity is approaching $1$ as a
function of $L$. In Fig. 3 we plot the infidelity, $1-F(-\frac{\pi}{2})$
(minimized over choice of $\sigma$) as a function of $L$ to extract a rate of
convergence. The data can be approximated by the power law
$1-F=0.24L^{-0.76}$. Though we have few data points to extrapolate the
results, we can infer that the fidelity is indeed converging to unity
relatively fast.

Figure 3: Loglog plot of infidelity as a funcion of $L$.

Two important features stand out in the graphs in Fig. 1 and Fig. 2. First,
that all of them display a plateau of behaviour above a certain wavepacket
width, which suggests that a gate implemented this way would be robust to
variations of the wavepacket preparation (contrast this with the results of
Brod and Combes (2016) for two-photon scattering on a line of interaction
sites, where the maximum only occurs for a single well-calibrated value of
$\sigma$). Second, very high fidelities (over 95%) are already achieved for
relatively small number of interaction sites (under 30). This suggests that
these two-qubit gates might be feasible with a modest amount of resources in a
real-world implementation of this computational model. Furthermore, since the
fidelity achieved is close to unity, it means that the scattering does not
change appreciably the shape of the wavepackets, which is surprising given
that the particles have many degrees of freedom that can become entangled.

##  V Conclusion

We have fully solved the scattering problem for two particles propagating on
an infinite line interacting via the Bose-Hubbard Hamiltonian on a finite
region of contiguous vertices. We analyzed the scattering in the asymptotic
limit of infinite interaction vertices both analytically and numerically. The
S-matrix we obtained reduces to previously-known results Valiente and
Petrosyan (2008); Childs  _et al._ (2013) in the asymptotic limit (though the
full S-matrix is not described explicitly in those papers, it can be
inferred). By numerically analyzing the scattering in the finite region
regime, we observed that a two-qubit gate with high fidelity can be already be
achieved for very few interaction sites, and with wavepackets that may be
easier to prepare than approximations to momentum eigenstates (e.g. with the
aid of momentum filters Childs (2009)).

Our work showcases how the Lippmann-Schwinger formalism can be used to solve
the more general problem of multiparticle scattering on graphs, paving the way
for a variety of future research directions. The more obvious next steps are
to develop a full framework to treat multiparticle scattering on graphs with
semi-infinite paths attached, such as those considered in Childs  _et al._
(2013). This might be used, for example, to develop new graphs that implement
multi-qubit gates directly, or to simplify gadgets such as the momentum
switches Childs  _et al._ (2014). A primitive in this direction would be to
find nontrivial graphs with perfect transmission from input to output paths
for two or more interacting particles.

Another goal is to filter out which ingredients make multiparticle scattering
universal for quantum computing, developing new tools to study the
computational complexity of quantum walks. For example, it is interesting to
consider whether it is possible to perform universal computation where all
particles have the same momenta (dispensing with the use of momentum
switches), and for which values of momentum this holds. Conversely, one might
ask whether there are values of momenta for which multiparticle scattering is
classically simulable or otherwise nonuniversal Aaronson and Arkhipov (2011).

Finally, we envision applications of these results beyond the field quantum
computation. Scattering theory is a successful formalism used in other fields
such as quantum optics, quantum field theory, condensed matter, and so on, and
its discrete-space analogue is not so well-understood. In that sense, our
result is analogous to that of Brod and Combes (2016); Brod  _et al._ (2016),
and it is an interesting question to try to extend this analogy further.

###### Acknowledgements.

This work was supported by Brazilian funding agencies CNPq and FAPERJ.

## References

  * Ambainis  _et al._ (2001) A. Ambainis, E. Bach, A. Nayak, A. Vishwanath, and J. Watrous, in  _Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing_, STOC ’01 (ACM, New York, 2001) pp. 37–49.
  * Nayak and Vishwanath (2000) A. Nayak and A. Vishwanath, arXiv:quant-ph/0010117 (2000).
  * Childs  _et al._ (2002) A. Childs, E. Farhi, and S. Gutmann, Quantum Information Processing 1, 35 (2002).
  * Farhi and Gutmann (1998) E. Farhi and S. Gutmann, Phys. Rev. A 58, 915 (1998).
  * Childs  _et al._ (2003) A. M. Childs, R. Cleve, E. Deotto, E. Farhi, S. Gutmann, and D. A. Spielman, in  _Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing_, STOC ’03 (ACM, New York, 2003) pp. 59–68.
  * Childs (2009) A. M. Childs, Phys. Rev. Lett. 102, 180501 (2009).
  * Childs  _et al._ (2013) A. M. Childs, D. Gosset, and Z. Webb, Science 339, 791 (2013).
  * Valiente and Petrosyan (2008) M. Valiente and D. Petrosyan, J. Phys. B: At. Mol. Opt. Phys. 41, 161002 (2008).
  * Sakurai (1994) J. J. Sakurai,  _Modern quantum mechanics_ (Addison-Wesley, Reading, MA, 1994).
  * Childs and Gosset (2012) A. M. Childs and D. Gosset, Journal of Mathematical Physics 53, 102207 (2012).
  * Varbanov and Brun (2009) M. Varbanov and T. A. Brun, Phys. Rev. A 80, 052330 (2009).
  * Weinberg (1995) S. Weinberg,  _The Quantum Theory of Fields, Volume I Foundations_ (Cambridge University Press, 1995).
  * Zhu and Wakin (2016) Z. Zhu and M. B. Wakin, arXiv:1608.04820 [cs.IT] (2016).
  * Gray (2006) R. M. Gray,  _Toeplitz and Circulant Matrices: A review_ (Foundations and Trends in Communications and Information Theory, Vol 2, Issue 3, pp 155-239, 2006).
  * Brod and Combes (2016) D. J. Brod and J. Combes, Phys. Rev. Lett. 117, 080502 (2016).
  * Childs  _et al._ (2014) A. Childs, D. Gosset, D. Nagaj, M. Raha, and Z. Webb, Quantum Information and Computation 15 (2014), 10.26421/QIC15.7-8-5.
  * Aaronson and Arkhipov (2011) S. Aaronson and A. Arkhipov, in  _Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing_, STOC ’11 (Association for Computing Machinery, New York, NY, USA, 2011) pp. 333–342.
  * Brod  _et al._ (2016) D. J. Brod, J. Combes, and J. Gea-Banacloche, Phys. Rev. A 94, 023833 (2016).
  * Byrd and Friedman (1971) P. F. Byrd and M. D. Friedman,  _Handbook of Elliptic Integrals for Engineers and Scientists_ (Springer Berlin, Heidelberg, 1971).

##  Appendix A Computing $J(E,n)$

Let us consider the family of integrals

$\displaystyle J(E,n)=\dfrac{1}{2\pi}\lim_{\epsilon\rightarrow
0}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\frac{e^{-i(k_{1}+k_{2})n}}{E-2\cos{k_{1}}-2\cos{k_{2}}+i\epsilon}dk_{1}dk_{2}.$
(40)

It is easy to see that $J(E,n)=J(E,-n)=-\overline{J(-E,n)}$, so we only need
to compute the integral for $n\geq 0$ and $E\geq 0$. Let us assume $E\geq 0$
in the following steps. Changing to the variables $k_{\pm}=\frac{k_{1}\pm
k_{2}}{2}$, we have

$\displaystyle J(E,n)=\dfrac{1}{2\pi}\lim_{\epsilon\rightarrow
0}\int_{-\pi}^{\pi}e^{-2ik_{+}n}\int_{-\pi}^{\pi}\frac{1}{E-4\cos{k_{+}}\cos{k_{-}}+i\epsilon}dk_{-}dk_{+}.$
(41)

Introducing a complex variable $z=e^{ik_{-}}$, the inner integral is written
as

$\displaystyle
I=i\oint\frac{1}{2\cos(k_{+})z^{2}-(E+i\epsilon)z+2\cos{k_{+}}}dz,$ (42)

where the contour is the unit circle in the positive orientation. The poles of
the integrand are the roots of the denominator,
$z_{\pm}=\frac{E+i\epsilon}{4\cos{k_{+}}}\pm\sqrt{\left(\frac{E+i\epsilon}{4\cos{k_{+}}}\right)^{2}-1}$,
where the square root is chosen to have positive real part. Let us assume for
the moment that $\cos{k_{+}}>0$.

If $4\cos{k_{+}}<E$, then the roots are close to real and as
$\epsilon\rightarrow 0$ we have $z_{-}<1<z_{+}$ as we can see in figure 4. So
for sufficiently small $\epsilon$, the only root inside the unit circle is
$z_{-}$.

If $4\cos{k_{+}}>E$, we can write the roots as
$z_{\pm}=\frac{E+i\epsilon}{4\cos{k_{+}}}\pm
i\sqrt{1-\left(\frac{E+i\epsilon}{4\cos{k_{+}}}\right)^{2}}$, where the square
root is chosen to have positive real part, similar to before. then the roots
are not on the real line, and approach the unit circle as $\epsilon\rightarrow
0$, as pictured in figure 4. But for $\epsilon\neq 0$ and $|a|<1$, we have

$\displaystyle\left|a+i\epsilon\pm
i\sqrt{1-(a+i\epsilon)^{2}}\right|^{2}=1\pm\frac{2\epsilon}{\sqrt{1-a^{2}}}+O(\epsilon^{2}),$
(43)

so, in both cases, the only root inside the contour is $z_{-}$. For the case
$\cos{k_{+}}<0$, the roots $z_{\pm}$ switch roles, specifically, only the root
$z_{+}$ is inside the contour.

Figure 4: Left: Real roots $z_{\pm}$ for the case $0<4\cos{k_{+}}<E$ as
$\epsilon\rightarrow 0$. Right: Complex roots $z_{\pm}$ for the case
$4\cos{k_{+}}>E$.

By residues, we have

$\displaystyle I$ $\displaystyle=2\pi
i\begin{cases}\mathop{\mathrm{Res}}_{z=z_{-}}\dfrac{i}{2\cos{k_{+}}(z-z_{-})(z-z_{+})}\,&\text{if}\,\cos{k_{+}}>0\\\
\mathop{\mathrm{Res}}_{z=z_{+}}\dfrac{i}{2\cos{k_{+}}(z-z_{-})(z-z_{+})}\,&\text{if}\,\cos{k_{+}}<0\end{cases}$
(44)
$\displaystyle=\dfrac{2\pi}{\sqrt{\left(E+i\epsilon\right)^{2}-16\cos^{2}{k_{+}}}},$

where the square root is chosen as before; in the limit $\epsilon\rightarrow
0$, the square root is $\sqrt{E^{2}-16\cos^{2}{k_{+}}}$ for $4|\cos{k_{+}}|<E$
and $i\sqrt{16\cos^{2}{k_{+}}-E^{2}}$ for $4|\cos{k_{+}}|>E$. Now we can take
the $\epsilon\rightarrow 0$ limit and we have

$J(E,n)=\int_{0}^{\frac{\pi}{2}}\frac{\cos{2np}}{\sqrt{\left(\frac{E}{4}\right)^{2}-\cos^{2}{p}}}dp.$

By decomposing $\cos{2np}=T_{2n}(\cos{p})$, where $T_{m}$ is the $m$-th
Chebyshev polynomial of the first kind, we only need to compute the integrals

$\int_{0}^{\frac{\pi}{2}}\frac{\cos^{2n}{p}}{\sqrt{\left(\frac{E}{4}\right)^{2}-\cos^{2}{p}}}dp=\frac{1}{i\sqrt{1-\left(\frac{E}{4}\right)^{2}}}\left(\int_{0}^{\frac{\pi}{2}}\frac{\cos^{2n}{p}}{\sqrt{1-k^{2}\sin^{2}{p}}}dp\right)^{*}=:\frac{1}{i\sqrt{1-\left(\frac{E}{4}\right)^{2}}}(C_{n})^{*},$

where $k=\sqrt{\frac{16}{16-E^{2}}}$ is called the elliptic modulus, and the
square root of a negative number is taken to have positive imaginary part at
all steps. The integrals $C_{n}$ are related to the complete elliptic
integrals of the first and the second kind, $K(k)$ and $E(k)$ respectively.
They are defined as

$\displaystyle K(k)$
$\displaystyle=\int_{0}^{\frac{\pi}{2}}\frac{dp}{\sqrt{1-k^{2}\sin^{2}{p}}},$
(45) $\displaystyle E(k)$
$\displaystyle=\int_{0}^{\frac{\pi}{2}}\sqrt{1-k^{2}\sin^{2}{p}}\,dp.$

The sequence of integrals $C_{n}$ can be computed with the following
recurrence relation (see equation 312.05 of Byrd and Friedman (1971)):

$\displaystyle C_{0}$ $\displaystyle=K(k),$ (46) $\displaystyle C_{1}$
$\displaystyle=\frac{1}{k^{2}}\left(E(k)-(1-k^{2})K(k)\right),$ $\displaystyle
C_{n}$
$\displaystyle=\frac{(2n-2)(2k^{2}-1)C_{n-1}-(2n-3)(1-k^{2})C_{n-2}}{(2n-1)k^{2}}\qquad\forall
n\geq 2.$

Finally, the integral $J(E,n)$ for $E\geq 0$ can be computed as

$J(E,n)=\frac{4^{n}n}{i\sqrt{1-\left(\frac{E}{4}\right)^{2}}}\sum_{m=0}^{n}\frac{(-1)^{m}(2n-m-1)!}{4^{m}m!(2n-2m)!}(C_{n-m})^{*}.$

##  Appendix B Some Dirac comb identities

The Dirac comb is just a periodic extension of the Dirac delta, and one
possible representation is given by its Fourier series:

$\displaystyle\delta_{T}(x)=\frac{1}{T}\sum_{n\in\mathbb{Z}}e^{i2\pi
n\frac{x}{T}},$ (47)

where $T$ is the period of the Dirac comb.

One intuitive identity we used (in (25)) is the scaling of the Dirac comb by a
factor of two,
$2\delta_{T}(2x)=\delta_{\frac{T}{2}}(x)=\delta_{T}(x)+\delta_{T}(x+\frac{T}{2})$,
similar to the related scaling identity for the Dirac delta, but also changing
the frequency. This can be easily proved using the Fourier series
representation:

$\displaystyle\delta_{T}(x)+\delta_{T}(x+\frac{T}{2})$
$\displaystyle=\frac{1}{T}\sum_{n\in\mathbb{Z}}e^{i2\pi
n\frac{x}{T}}+\frac{1}{T}\sum_{n\in\mathbb{Z}}e^{i2\pi
n\frac{x+\frac{T}{2}}{T}}$ (48)
$\displaystyle=\frac{1}{T}\sum_{n\in\mathbb{Z}}e^{i2\pi
n\frac{x}{T}}+\frac{1}{T}\sum_{n\in\mathbb{Z}}(-1)^{n}e^{i2\pi n\frac{x}{T}}$
$\displaystyle=\frac{2}{T}\sum_{n\,even}e^{i2\pi n\frac{x}{T}}$
$\displaystyle=\frac{2}{T}\sum_{n\in\mathbb{Z}}e^{i4\pi
n\frac{x}{T}}=2\delta_{T}(2x).$

Another identity was used (in (32)) :

$\displaystyle\delta(E_{k_{1}k_{2}}-E_{p_{1}p_{2}})\delta_{2\pi}(k_{1}+k_{2}-p_{1}-p_{2})=\frac{1}{2|\sin
k_{1}-\sin
k_{2}|}(\delta_{2\pi}(k_{1}-p_{1})\delta_{2\pi}(k_{2}-p_{2})+\delta_{2\pi}(k_{1}-p_{2})\delta_{2\pi}(k_{2}-p_{1})),$

where $E_{q_{1}q_{2}}=2\cos{q_{1}}+2\cos{q_{2}}=4\cos{q_{+}}\cos{q_{-}}$ with
$q_{\pm}=\frac{q_{1}\pm q_{2}}{2}$.

Starting with the left side, we can simplify it

$\displaystyle\delta(E_{k_{1}k_{2}}-E_{p_{1}p_{2}})$
$\displaystyle\delta_{2\pi}(k_{1}+k_{2}-p_{1}-p_{2})=\delta(4\cos{k_{+}}\cos{k_{-}}-4\cos{p_{+}}\cos{p_{-}})\delta_{2\pi}(2k_{+}-2p_{+})$
(49)
$\displaystyle=\delta(4\cos{k_{+}}\cos{k_{-}}-4\cos{p_{+}}\cos{p_{-}})\frac{1}{2}(\delta_{2\pi}(k_{+}-p_{+})+\delta_{2\pi}(k_{+}-p_{+}+\pi))$
$\displaystyle=\frac{1}{8|\cos{k_{+}}|}(\delta(\cos{k_{-}}-\cos{p_{-}})\delta_{2\pi}(k_{+}-p_{+})+\delta(\cos{k_{-}}+\cos{p_{-}})\delta_{2\pi}(k_{+}-p_{+}+\pi)),$

where we used the identity (47) on the second line.

The first delta of cosines can be further simplified to a sum of Dirac combs:

$\displaystyle\delta(\cos{k_{-}}-\cos{p_{-}})=\frac{1}{|\sin{k_{-}}|}(\delta_{2\pi}(k_{-}-p_{-})+\delta_{2\pi}(k_{-}+p_{-})),$
(50)

and similarly for the second product. So we have the full expression

$\displaystyle\frac{1}{8|\cos{k_{+}}\sin{k_{-}}|}($
$\displaystyle\delta_{2\pi}(k_{-}-p_{-})\delta_{2\pi}(k_{+}-p_{+})+\delta_{2\pi}(k_{-}+p_{-})\delta_{2\pi}(k_{+}-p_{+})+$
(51)
$\displaystyle\delta_{2\pi}(k_{-}-p_{-}+\pi)\delta_{2\pi}(k_{+}-p_{+}+\pi)+\delta_{2\pi}(k_{-}+p_{-}+\pi)\delta_{2\pi}(k_{+}-p_{+}+\pi)),$

and we want to go back to the original variables, $k_{1},k_{2},p_{1}$ and
$p_{2}$.

We compute in detail the first product of the above expression.

######  Lemma 1.

$\frac{1}{2}\delta_{2\pi}(k_{-}-p_{-})\delta_{2\pi}(k_{+}-p_{+})=\delta_{4\pi}(k_{1}-p_{1})\delta_{4\pi}(k_{2}-p_{2})+\delta_{4\pi}(k_{1}-p_{1}+2\pi)\delta_{4\pi}(k_{2}-p_{2}+2\pi)$.

###### Proof.

Let us expand the factors on the right hand side as Fourier series:

$\displaystyle\frac{1}{(4\pi)^{2}}\sum_{n,m\in\mathbb{Z}}e^{in\frac{k_{1}-p_{1}}{2}}$
$\displaystyle
e^{im\frac{k_{2}-p_{2}}{2}}(1+(-1)^{n+m})=\frac{1}{(4\pi)^{2}}\sum_{n,l\in\mathbb{Z}}e^{in\frac{k_{1}-p_{1}}{2}}e^{i(l-n)\frac{k_{2}-p_{2}}{2}}(1+(-1)^{l})$
(52)
$\displaystyle=\frac{2}{(4\pi)^{2}}\sum_{\begin{subarray}{c}n\in\mathbb{Z}\\\
l\,even\end{subarray}}e^{in(k_{-}-p_{-})}e^{il\frac{k_{+}-k_{-}-p_{+}+p_{-}}{2}}$
$\displaystyle=\frac{2}{(4\pi)^{2}}\sum_{n,l\in\mathbb{Z}}e^{in(k_{-}-p_{-})}e^{il(k_{+}-k_{-}-p_{+}+p_{-})}$
$\displaystyle=\frac{1}{2}\frac{1}{(2\pi)^{2}}\sum_{q,l\in\mathbb{Z}}e^{iq(k_{-}-p_{-})}e^{il(k_{+}-p_{+})}=\frac{1}{2}\delta_{2\pi}(k_{-}-p_{-})\delta_{2\pi}(k_{+}-p_{+}),$

where we set $l=m+n$ in the first line and $q=n-l$ in the last line. ∎

Similarly, by substituting $k_{1}\rightarrow k_{1}+2\pi$ in Lemma 1, we get an
identity for the third delta product of Eq. (51):

$\frac{1}{2}\delta_{2\pi}(k_{-}-p_{-}+\pi)\delta_{2\pi}(k_{+}-p_{+}+\pi)=\delta_{4\pi}(k_{1}-p_{1}+2\pi)\delta_{4\pi}(k_{2}-p_{2})+\delta_{4\pi}(k_{1}-p_{1})\delta_{4\pi}(k_{2}-p_{2}+2\pi).$

So, summing the first and third delta products in Eq. (51), we get a $2\pi$
periodic expression instead of a $4\pi$ periodic one:

$\frac{1}{2}\delta_{2\pi}(k_{-}-p_{-}+\pi)\delta_{2\pi}(k_{+}-p_{+}+\pi)+\frac{1}{2}\delta_{2\pi}(k_{-}-p_{-})\delta_{2\pi}(k_{+}-p_{+})=\delta_{2\pi}(k_{1}-p_{1})\delta_{2\pi}(k_{2}-p_{2}).$

Similarly, the same steps can be done for the second and fourth delta products
in Eq. (51), by exchanging $k_{1}$ and $k_{2}$. Finally, we can simplify Eq.
(51) to

$\displaystyle\frac{1}{8|\cos{k_{+}}\sin{k_{-}}|}(2\delta_{2\pi}(k_{1}-p_{1})\delta_{2\pi}(k_{2}-p_{2})+2\delta_{2\pi}(k_{1}-p_{2})\delta_{2\pi}(k_{2}-p_{1}))$
(53) $\displaystyle=\frac{1}{2|\sin k_{1}-\sin
k_{2}|}(\delta_{2\pi}(k_{1}-p_{1})\delta_{2\pi}(k_{2}-p_{2})+\delta_{2\pi}(k_{1}-p_{2})\delta_{2\pi}(k_{2}-p_{1})).$