Molecular binding in interacting quantum walks

The easiest way to see the molecular binding effect is by direct simulation of the dynamics described in the previous section, starting from the standard initial state (internal singlet at the origin). The resulting joint probability distribution of the positions (x₁,x₂) after t = 50 steps is shown in figure 1, which allows a direct comparison between the non-interacting case g = 0 and the maximal interaction phase g = π. Both pictures are symmetric with respect to the main diagonal x₁ = x₂. In the non-interacting case, the dominant features are two peaks at x₁ = −x₂ with $\vert x_i\vert \approx t/\sqrt 2$. These peaks are well known from the Hadamard walk of a single particle, and to facilitate the comparison, a dashed square with sides at $x_1,x_2=\pm t/\sqrt 2$ is shown in the left panel. There is a correlation effect resulting from the antisymmetric initial state [10], which is expressed by the absence of peaks at $x_1=x_2=\pm t/\sqrt 2$.

Zoom In Zoom Out Reset image size

Turning to the interacting case we observe the following main features:

• (a)  
  The peaks near $x_1=-x_2\approx t/\sqrt 2$ are still visible, but smaller. Indeed, they will continue to move as in the non-interacting case because particle pairs with a long distance do not see the interaction anyway. It will be shown later in section 3.6 that the total weight of these 'free' peaks is 1/3.
• (b)  
  The rest of the probability is concentrated near the diagonal x₁ = x₂. Along the diagonal the spreading is proportional to the time, which is called ballistic behavior in the context of quantum walks. The center-of-mass motion of the molecules thus exhibits the same behavior as other walking particles, including two leading peaks. These are slower than the single-atom peaks. We will later show them to be at x₁ = x₂ = ± t/3 for large number of steps t.
• (c)  
  The width around the diagonal is constant in time. More precisely, apart from the free peaks that move away, the probability to find |x₁ − x₂| > r decreases exponentially in r, with constants independent of t. This is the dynamical statement of molecular binding.

Another signature of molecular binding can be identified in the eigenvalue spectrum of the walk operator W. Since, by construction, it commutes with the joint translations (x₁,x₂)↦(x₁ + a,x₂ + a) we will jointly diagonalize the walk operator W with the unitary translation operator T. Both operators are unitary and hence have eigenvalues on the unit circle. For graphical representation we unwrap these circles and hence represent a pair of eigenvalues $({\mathrm { exp}}({\rm i} p),{\rm \exp }({\rm i}\omega )$ of T and W by a point (p,ω) in the square [ − π,π] × [ − π,π]. This introduction of quasi-momentum p in the Brillouin zone [ − π,π] is very familiar from solid state theory. Here, as generally in discrete time dynamics and Floquet theory, we do the same for the quasi-energy ω.

In order to obtain the joint eigenvalues numerically we close the system to a ring of length L. The resulting Hilbert space will have dimension 4L², so the diagonalization can be carried out numerically for fairly large ring sizes. A typical result (again with interaction phase g = π, and for a small ring size L = 32) is shown in the right panel of figure 2.

Zoom In Zoom Out Reset image size

The densely filled regions, which we call bands, can be understood completely from the non-interacting case. To this end consider the left panel in figure 2. The free walk operator commutes not only with joint translations but also with separate translations of the particles. Therefore, we can jointly diagonalize the two separate translation operators and the walk operator, resulting in a cube with axes total quasi-momentum p = p₁ + p₂, relative quasi-momentum k = (p₁ − p₂)/2 and quasi-energy ω, respectively. In the figure, k is taken as depth, roughly perpendicular to the paper plane. Now, from the known spectrum of the single-particle walk W₁ we can immediately determine the joint spectrum of W₂ explicitly. This gives the blue surface in the left panel. The points on that surface are the discrete eigenvalue triples (p,k,ω) one obtains by closing the system to a ring of the same size as in the right panel. Consider now the image obtained by a parallel projection of the left panel figure to the paper plane: this gives the bands also observed in the interacting case. To be precise, for finite L the continuum eigenvalues do not fall exactly on the interacting ones inside the bands: there is a slight deformation due to the interaction. However, as we will show in section 5, the continuous spectrum obtained in the limit L → ∞ is independent of the interaction phase g.

The new features due to the interaction are the lines of eigenvalues appearing in the gaps between the bands (see the left panel). These are the molecular bound states that will be determined analytically in the following section.

We will only sketch the main ideas, because the derivation will be given with more detail and in a much more general setting in section 5.

For two-particle systems with continuous translational symmetry, and an interaction depending only on the distance of the particles, the obvious first step for identifying molecular states is to introduce center-of-mass and relative coordinates. The same is true in the lattice case, but with an important difference: in continuous space the decoupling is complete, and the problem to be solved for the relative motion is not only independent of the center-of-mass position but, by Galilean boost invariance also independent of the total momentum. In discrete space, this symmetry is meaningless and the relative motion in general depends on the total quasi-momentum p = p₁ + p₂. In the momentum representation, this means that the diagonalization of W can be done separately for each value of the conserved quantity p, and one is left with a problem in a space of wave functions depending on the relative momentum k = (p₁ − p₂)/2, with p treated as a fixed parameter. The ranges of p and k can be taken to be [ − π,π] by a suitable reshuffling of Brillouin zones (see section 5.1).

The free walk thus leads to the operator

$$\begin{equation} W_p(k)=W_1(p/2+k)\otimes W_1(p/2-k) \end{equation} \tag{ 4 }$$

acting on the tensor product of the two internal state spaces. The collision space is given by (x₁ − x₂) = 0, which in the momentum representation corresponds to the constant functions of k. We denote by N the projection onto this subspace, so $N\Psi =(2\pi )^{-1}\int \!\!\mathrm {d}k\ \Psi (k)$. Then

$$\begin{equation} (G\Psi)(k)=\Psi(k)+({\rm e}^{{\rm i} g}-1)N\Psi, \end{equation} \tag{ 5 }$$

and Ψ is an eigenvector of W = W₂G with eigenvalue exp(iω) if

$$\begin{equation} W_{p}(k)\Bigl(\Psi(k)+({\rm e}^{{\rm i} g}-1)N\Psi\Bigr) = {\rm e}^{{\rm i} \omega}\Psi(k). \end{equation} \tag{ 6 }$$

One can show by a transfer matrix argument that in the case at hand there are no eigenvalues embedded in the continuous spectrum. Therefore, from now on exp(iω) is understood to be in the spectral gap of W_p, i.e. W_p(k) − exp(iω) is invertible for all k. According to our choice of the antisymmetric subspace in section 2, NΨ must be proportional to the singlet ψ₋, and we can choose NΨ = ψ₋. We can thus solve (6) for Ψ:

$$\begin{equation} \Psi(k)=(W_{p}(k)-{\rm e}^{{\rm i}\omega})^{-1}W_{p}(k)(1-{\rm e}^{{\rm i}g})\psi_-, \end{equation} \tag{ 7 }$$

What determines the eigenvalues ω is the consistency condition NΨ = ψ₋. This means that we have to evaluate the integral over k of the right-hand side of (7). Since W_p depends on k as a polynomial in ${\mathrm { \exp }}({\rm i} k)$, this integral is a complex integral over the unit circle, which can be evaluated via the residue theorem. The resulting equation can be solved, and gives

$$\begin{equation} {\rm e}^{{\rm i}\omega} =\frac{{\rm e}^{{\rm i} g}}{2\rm\,e^{{\rm i} g}-1}\left(\cos\!\,p\pm \mathrm{i}\sqrt{\sin^2p+4(1-\cos\!\,g)}\right) \end{equation} \tag{ 8 }$$

$$\begin{equation} \mbox{ provided}\quad \sin\omega\cdot\sin(g-\omega)>0. \end{equation} \tag{ 9 }$$

Condition (9) results from picking only poles inside the unit circle for residue evaluation. The two branches of relation (8) are plotted in figure 3, and directly confirm the numerical result in figure 2. Some of the eigenvalue curves end at the edge of the continuous spectrum, because (9) has to be taken into account. Ignoring that condition leads to unphysical 'virtual' molecule states, shown in the right panel of figure 3.

Zoom In Zoom Out Reset image size

Of course, for admissible triples (g,p,ω), (7) gives an explicit expression for the bound state wave function. These functions or rather their Fourier transforms from the variable k to the variable (x₁ − x₂) are shown in figure 4 for two values of g. Clearly, they are decaying exponentially on the scale of a few lattice sites. See section 5.4 for a explicit construction of the eigenfunctions.

Zoom In Zoom Out Reset image size

Let us now turn to the dynamical characteristics of the molecules. Choosing, for each p, normalized eigenfunctions χ^±_p(k) in the upper and in the lower gap wherever they exist, we can write the overall wave function of a molecule as Ψ(p,k) = ϕ₊(p)χ⁺_p(k) + ϕ₋(p)χ⁻_p(k). Here, it is understood that we put ϕ_±(p) = 0, whenever the corresponding branch is ruled out by (9). The time evolution of such wave functions will be completely determined by the branches ω_±(p) of the quasi-energy. After t time steps we get

$$\begin{equation} (W^t\Psi)(p,k)=\mathchoice{{\rm e}^{{\rm i}t\omega_+(p)}}{\exp({\rm i}t\omega_+(p))}{{\rm e}^{{\rm i}t\omega_+(p)}}{\cdot}\,\chi^+_p(k)\ \phi_+(p)+\mathchoice{{\rm e}^{{\rm i}t\omega_-(p)}}{\exp({\rm i}t\omega_-(p))}{{\rm e}^{{\rm i}t\omega_-(p)}}{\cdot}\,\chi^-_p(k)\ \phi_-(p). \end{equation} \tag{ 10 }$$

The factors χ just play the role of a 'spinor basis' for the space of internal states of the molecules. With this abstraction we have a discrete time evolution W_eff on a lattice, with two internal degrees of freedom, although some internal states may be forbidden for some values of p. This is the effective time evolution of the molecules.

However, it is far from clear whether W_eff is again a quantum walk, i.e. whether it allows for a description in position space that involves only strictly finite jumps. The position representation of an effective walk is obtained from the momentum representation by Fourier transform which can, however, be preceded by a p-dependent basis change. Up to such a change the walk is obviously given completely in terms of the dispersion relations ω_±(p). Hence, the question boils down to deciding whether these functions are the dispersion relations of some quantum walk. Since the full interacting walk has strictly finite jumps for pairs of particles, it would seem natural that the subspace of molecule states inherits this property. However, the projection to the molecule subspace (given by the eigenfunctions χ_±) is not local, since the eigenfunctions are only exponentially localized.

It turns out that for the family of interacting Hadamard walks studied here W_eff is indeed a quantum walk. This family even allows for a simple decomposition into a single shift and coin operation: $W_{\mathrm { eff}}=S(\mathbbm{1}\otimes C)$ with the shift (1) and a 2 × 2 coin matrix C. Indeed, we exactly obtain the dispersion relations (8) with the coin

$$\begin{equation} C=\frac{1}{2{\rm e}^{{\rm i} g}-1}\left(\begin{array}{@{}cc@{}}{\rm e}^{{\rm i}g} & \sqrt{2} ({\rm e}^{{\rm i}g} -1) \\ \sqrt{2} ({\rm e}^{{\rm i}g} -1) {\rm e}^{{\rm i}g} & {\rm e}^{{\rm i}g}\end{array}\right)\!.\end{equation}\noindent \tag{ 11 }$$

The walk uses both branches in (8), though only 'virtually' if one or both of them are forbidden by the constraint.

The asymptotic position distribution of W_eff can be computed directly from the dispersion relation ω_±(p). In fact, for any walk the probability distribution of Q(t)/t, i.e. the position in 'ballistic scaling', converges to the distribution of the group velocity operator V [23]. V has the same p-dependent eigenvectors as W, but eigenvalues ω_±'(p) = dω_±(p)/dp instead of $\mathchoice {{\mathrm { e}}^{{\rm i}\,\omega _\pm (p)}}{\exp ({\rm i}\,\omega _\pm (p))}{{\rm e}^{{\rm i}\,\omega _\pm (p)}}{\cdot }$. Thus, we can compute the asymptotic position distribution directly for every initial state, and this depends only on the joint distribution of p and the branch index ±. For an initial state at the origin the quasi-momentum distribution is flat. Therefore, when ω_± has an inflection point at some p_*, the group velocity ω_±'(p_*) has a maximum, and consequently the probability density has a square root singularity ('caustic'). This explains all the peaks seen in the simulated position density diagrams of this paper. For the free Hadamard walk, the peaks are at $\omega _\pm '(\pi /2)=1/\sqrt 2$. For the effective walks with coin (11) the inflection point is always at p_* = ± π/2 (see figure 3). The left panel of figure 5 shows the asymptotic velocity distributions as a function of the interaction parameter g. For g < π/3 and g > 2π − π/3, the constraint (9) forbids momenta around p_*. Hence, the quasi-momenta producing the caustic peaks in the effective walk starting from the origin fall into the virtual branch of the dispersion relation, and have probability density zero. Hence, we expect not to see any peaks in the walk of molecules. This is confirmed by the right panel in figure 5, and the corresponding asymptotic distribution (blue line) in the right panel.

Zoom In Zoom Out Reset image size

The single-particle character of the molecular state can also be seen from its reaction to an additional static force. Such a force can be modeled by an additional position-dependent phase factor that the particles acquire in each time step. This results in additional shifts in the total momentum variable, and is seen as an oscillation of the molecule's position distribution. Such oscillations occurring at a frequency proportional to the gradient of the external potential are well known as Bloch oscillations in condensed matter physics [24]. As can be seen in figure 6 the oscillation of the bound state exhibits a doubled frequency, i.e. the molecule experiences twice the force as compared to the single particles. Bloch oscillations lead to rapid dissociation of the molecule, if they can shift momentum to a region without bound state, e.g. g = π/4 in figure 3. Even for the most stable case g = π, however, the static force will eventually dissociate the molecule, which is foreshadowed by the slow increase of the variance in figure 6.

Zoom In Zoom Out Reset image size

Molecules have to be initially prepared. In order to make them detectable in an actual experiment, their preparation efficiency has to be sufficiently high. Since the dynamics of the walk W itself has no dissipation, an initial state with two widely separated atoms would have negligible overlap with any bound state, and this remains true for all times. However, by initially using a dynamics that is different than W it is possible to prepare two atoms in the same site (see section 4). This state then has a non-zero overlap with molecule states. The resulting probability of forming a molecule, i.e. the modulus squared of the overlap, will be studied in this section showing its dependence on the interaction parameter g.

Starting from two atoms in the same lattice site, the formation of molecule states is apparent from figure 1: by just letting the interacting walk run for a few steps the unbound pairs will move away, and the remainder is in the molecular state. The critical parameter for this second stage of the preparation is the modulus squared of the overlap of the initial state with the molecular subspace P_ovl(g).

As before, the computation of P_ovl(g) is done separately for every quasi-momentum p and every branch ω_± of the dispersion relation. The result P_ovl(p,g) is plotted in figure 7. Since we consider the situation of a localized initial state, P_ovl(g) can be determined from P_ovl(p,g) by integration over p. This integrated overlap (see right panel in figure 7) has a maximum at g = π; so this case, also shown in figure 1, has the highest preparation efficiency P_ovl(π) = 2/3.

Zoom In Zoom Out Reset image size

Here, we will give a short account of the introduction of relative coordinates in a two-particle quantum walk. We will look at a slightly more general setting allowing the walkers to move on an s-dimensional lattice and carrying a d-dimensional internal state.

In analogy to section 2, we can then pass via Fourier transform from the tensor product of the two particles position and internal states in the Hilbert space $(\ell ^2({\mathbb Z}^s)\otimes {{\mathbb C}}^d)^{\otimes 2}$ to the quasi-momentum representation in $L^2([-\pi ,\pi ]^2,\mathrm {d}p_1 \mathrm {d}p_2)\otimes {{\mathbb C}}^{2d}$. This corresponds to a change of variables from the lattice positions x₁, x₂ of the two particles to their quasi-momenta p₁, p₂.

The generator of joint particle translations is the total quasi-momentum p = p₁ + p₂, which hence commutes with the walk. The remaining variable is the relative quasi-momentum k = (p₁ − p₂)/2. The factor 1/2 in the definition of k allows us to define the first Brillouin zone again as [ − π,π]². Figure 9 shows how this is done by identifying points outside [ − π,π]² with points inside this interval in the 1D lattice case (or separately for every component of the quasi-momentum vector in higher dimensions).

Zoom In Zoom Out Reset image size

The core of the computation in section 3.3 was done in relative coordinates, and effectively reduced to studying a quantum walk perturbed at the origin. This problem is also of independent interest, and we will here treat such systems in generality. We allow a lattice of arbitrary dimension s, i.e. ${\mathbb Z}^s$. The Hilbert space of the walk is then ${\mathcal H}= {\ell ^2({{\mathbb Z}}^s)\otimes {{\mathcal K}}}$ for some finite-dimensional internal space ${\mathcal K}$. The 'free' evolution will be an arbitrary unitary operator W₀ commuting with translations. By Fourier transform we pass to the momentum representation, so ${\mathcal H}$ is seen as the space of square integrable functions Ψ(k) on the Brillouin zone [ − π,π]^s with values in ${\mathcal K}$. In this representation W₀ thus acts by multiplication with k-dependent unitary matrix W(k).

We usually take it as part of the definition of a quantum walk that W₀ has a strictly finite propagation region. In this case each entry of the corresponding matrix W(k) is given by a trigonometric polynomial, i.e. a finite linear combination of functions $\mathchoice {{\mathrm { e}}^{{\rm i}n\cdot k}}{\exp ({\rm i}n\cdot k)}{{\rm e}^{{\rm i}n\cdot k}}{\cdot }$ with n a vector with integer components. We will not use much of this structure, but only that the eigenvalues $\mathchoice {{\mathrm { e}}^{{\rm i}\omega _\alpha (k)}}{\exp ({\rm i}\omega _\alpha (k))}{{\rm e}^{{\rm i}\omega _\alpha (k)}}{\cdot }$ of W(k) can be taken as a set of differentiable functions in the neighborhood of almost every k.

The subspace of ${\ell ^2({{\mathbb Z}}^s)\otimes {{\mathcal K}}}$ of functions, which are non-zero only at the origin, becomes in momentum space the space of k-independent functions. The corresponding projection N is

$$\begin{equation} (N\Psi)(k)=N\Psi=\int\frac{\mathrm{d}^sk'}{(2\pi)^2}\,\Psi(k'). \end{equation} \tag{ 12 }$$

Here, we have indicated that we denote by NΨ both the constant function in ${\mathcal H}$ and its value in ${\mathcal K}$. By modifying the unitary operator W₀ exactly on this subspace by a unitary operator Γ acting on ${\mathcal K}$ we get the perturbed walk as

$$\begin{equation} W = W_0\left((\mathbbm{1}-N)\otimes\mathbbm{1}_{\mathcal K} + N\otimes\Gamma\right). \end{equation} \tag{ 13 }$$

Note that as long as Γ is unitary, the second factor on the right-hand side of (13) and therefore also W have the same property.

Some qualitative features are immediate. Since $W-W_0=W_0N\otimes (\Gamma -\mathbbm{1})$ is a finite rank operator, the operators W and W₀ have the same absolutely continuous spectrum. This result is well known from scattering theory [45] and can be rigorously proved, e.g. see theorem IV.5.35 in [46]. The spectrum of W₀ is easy to determine: it consists of all numbers $\mathchoice {{\mathrm { e}}^{{\rm i}\omega _\alpha (k)}}{\exp ({\rm i}\omega _\alpha (k))}{{\rm e}^{{\rm i}\omega _\alpha (k)}}{\cdot }$ where k ranges over the Brillouin zone and α over the band indices. Eigenvalues can only occur when one of these functions is constant near some k, all other values belong to the absolutely continuous spectrum, and hence also to the continuous spectrum of W. We call the set of points $z={\mathrm { \exp }}({\rm i}\omega )$ not in the spectrum the gap set of W₀. These points are characterized by the property that the operator $W_0-z\mathbbm{1}$ has a bounded inverse. (We leave out the identity operator $\mathbbm{1}$ for notational convenience in the following).

The most prominent new spectral feature of W and one of the main interests in this paper are eigenvalues in the gap set of W₀. In general, allowing larger collisional Hilbert space there can also be embedded eigenvalues in the absolutely continuous spectrum, one example being the Bose case of an interacting two particle quantum walk, but we will limit the discussion here to eigenvalues in the gap set. Suppose that Ψ is an eigenvector of W with eigenvalue z in the gap set. Then, the eigenvalue equation in momentum space reads

$$\begin{equation} W(k)(\mathbbm{1}+N\otimes(\Gamma-\mathbbm{1}))\Psi(k) = z\Psi(k). \end{equation} \tag{ 14 }$$

Since N projects Ψ onto its value $\psi =N\Psi \in {\mathcal K}$ at the origin, we can take this as an equation for Ψ(k) in terms of the finite-dimensional vector $\psi \in {\mathcal K}$. Since z is in the gap set of W₀, and hence $W(k)-z\mathbbm{1}$ is invertible for all k, we can solve this equation and get

$$\begin{equation} \Psi(k) = \bigl(W(k)-z\bigl)^{-1}W(k)(\mathbbm{1}-\Gamma)\psi. \end{equation} \tag{ 15 }$$

This determines the eigenfunction Ψ, but we have to satisfy the consistency requirement that its projection NΨ is indeed the vector ψ. Hence, z is an eigenvalue if and only if the non-zero vector ψ satisfies the consistency condition

$$\begin{equation} \psi = \int\frac{\mathrm{d}^sk}{(2\pi)^2}\, \bigl(W(k)-z\bigl)^{-1}W(k)(\mathbbm{1}-\Gamma)\psi = R(z) (\mathbbm{1}-\Gamma)\psi, \end{equation} \tag{ 16 }$$

where we defined the operator R(z) by the integral. This self-consistency requirement is well known for the Lippmann–Schwinger equation, of which (16) is the homogeneous version [47]. We can further rewrite (16) in the following form:

$$\begin{equation} \Gamma\psi = (\mathbbm{1}-R(z)^{-1})\psi . \end{equation} \tag{ 17 }$$

As will be shown in the following subsection for z in the gap set the operator on the right-hand side of this equation is always unitary. Hence, every z in the gap set is an eigenvalue of a point perturbed walk W for suitable Γ. Indeed, we can choose $\Gamma =(\mathbbm{1}-R(z)^{-1})$ and find an eigenvalue that is k-fold degenerate with $k=\dim {\mathcal K}$.

The question is, how many solutions of (17) can we get for possibly different z. An easy general answer is not to be expected, since in the example of the antisymmetric Hadamard walk we can get 0, 1 or 2 bound states, although at most one per gap. 'At most one per gap' is indeed a general result for rank 1 perturbations of unitary operators [48]. This can be generalized to rank n perturbations where, in our case, $n=\dim N{\mathcal H}=\dim {\mathcal K}$. Via Cayley transform [49] we can import an elementary argument [50, proposition 2.1] from the perturbation theory of self-adjoint operators, concluding that there can be at most n eigenvalues per gap. From the previous argument we know that this upper bound is tight in each gap separately, but random examples often have significantly fewer bound states.

In this section, we prove the claim that the operator $(\mathbbm{1}-R(z)^{-1})$ appearing in (16) is unitary. We will do this by proving the following, more general claim:

Lemma 5.1. Consider a unitary operator W on a Hilbert space ${\mathcal H}$ and a point z on the unit circle which is not in the spectrum of W. Now let N be an arbitrary projection and consider $R(z)=N(W-z\mathbbm{1})^{-1}WN$ as an operator on $N{\mathcal H}$. Then R(z) invertible and $U=\mathbbm{1}-R(z)^{-1}$ is unitary.

The crucial step for the proof is to show the following relation:

$$\begin{equation} R(z)+R(z)^\dagger=\mathbbm{1}. \end{equation} \tag{ 18 }$$

We can show this for the special case $N=\mathbbm{1}$, and then multiply the above equation with N from both sides. Since W is unitary it has a spectral decomposition $W=\int \!\Omega \ E_W(\mathrm {d}\Omega )$, where E_W is the spectral measure on the unit circle. We can then express the operator (18) directly in the functional calculus (we use † for Hermitian adjoint and an overbar for complex conjugate):

$$\begin{eqnarray*} R(z)+R(z)^{\dagger} &=& \int\left(\frac\Omega{\Omega-z}+ \frac{\overline\Omega}{\overline\Omega-\overline z}\right)E_W(\mathrm{d}\Omega)\\ &=& \int\frac{\Omega({\overline\Omega}-\overline z)+ {\overline\Omega}(\Omega-z)} {\vert\Omega-z\vert^2}\,E_W(\mathrm{d}\Omega) \\ &=& \int\frac{2 -2 \Re({\overline\Omega}z)}{2 -2 \Re({\overline\Omega}z)}\,E_W(\mathrm{d}\Omega) \\ &=& \mathbbm{1}. \end{eqnarray*}$$

This determines the Hermitian part of R(z), and we can write $R(z)=\frac 12(\mathbbm{1}+\mathrm {i}K)$ for some Hermitian operator K on $N{\mathcal H}$.

For the next we use the functional calculus of $K=\int k\,E_K(\mathrm {d}k)$, with the spectral measure of K, which is now supported in the real axis. We note that due to the projection N, which need not commute with W, this measure cannot be easily obtained from the spectral measure E_W . Then

$$\begin{eqnarray*} U &=& \mathbbm{1}- R(z)^{-1}\\ &=&\int\left(1-\frac1{\frac12(1+\mathrm{i}k)}\right)\,E_K(\mathrm{d}k)\nonumber\\ &=&\int\frac{-1+\mathrm{i}k}{1+\mathrm{i}k}\,E_K(\mathrm{d}k). \end{eqnarray*}$$

Now, since | − 1 + ik| = |1 + ik| the integrand has absolute value 1, and so the integral represents a unitary operator.

The goal of this section is to derive an explicit formula for the bound states of a quantum walk with a point defect at the origin. We then apply this procedure to our standard example, the interacting Hadamard walk. This allows us to derive an explicit form of the p-dependent transformation needed to construct the virtual walk (11) of the molecule as explained at the end of section 3.2.

The components Ψ_x of the unnormalized eigenvector Ψ corresponding to relative distance x can be determined from (7), via inverse Fourier transformation. Substituting v for ${\mathrm { \exp }}({\rm i} k)$ this gives for the interacting Hadamard walk

$$\begin{equation} \Psi_x =\frac1{2\pi {\rm i}}\int\!\!\frac{\mathrm{d}v}{v^{x+1}}\left(\mathbbm{1} +{\rm e}^{{\rm i} \omega}(W_{p}(v)-{\rm e}^{{\rm i}\omega})^{-1}\right)({\rm e}^{{\rm i} g} -1) \psi_{-}. \end{equation} \tag{ 19 }$$

We have to substitute the allowed form (8) of the eigenvalue e^iω. Using residual calculus and the Fermi symmetry of the problem we can compute the normalized eigenvector $\tilde \Psi _x$ as

$$\begin{equation} \tilde\Psi_x =P_{\rm ovl}\cdot\left\{\begin{array}{@{}ll@{}} (1-{\rm e}^{{\rm i} g} )v_1^{|x|}R_{v_1}\psi_- , & x<0,\\ \psi_- , & x=0,\\ -(1-{\rm e}^{{\rm i} g} )v_1^{|x|}{\mathbb F} R_{v_1}\psi_- , & x>0. \end{array}\right. \end{equation} \tag{ 20 }$$

Here, P_ovl is the overlap of ψ₋ with the bound state (see section 3.6) and v₁ is the ω-dependent singularity of the integral kernel in (19) inside the unit circle, which has the explicit form

$$\begin{equation} v_1=-\cos \omega +\cot \frac{g}{2} \sin \omega,\quad \mbox{provided} \sin(\omega-g)\cdot\sin(\omega)<0, \end{equation} \tag{ 21 }$$

R_v1 is given by

$$\begin{eqnarray*} &&R_{v_1}=\mathop{{\mathrm{Res}}}\limits_{v \rightarrow v_1}\left( v^{-1}z(W(v)-z)^{-1}\right) \end{eqnarray*}$$

and ${\mathbb F}$ is the flip operator ( ${\mathbb F}\phi \otimes \psi =\psi \otimes \phi$). Note that since the modulus of v₁ is smaller than 1 the state $\tilde \Psi$ decays exponentially in |x|. This implies that the distance of the two particles in the molecule state is exponentially concentrated around zero.

In the case that the singularities of (W(v) − z)⁻¹ are at v₁ = 0, the integral kernel in (19) is analytic for |x| > 1, hence $\tilde \Psi _x=0$ for these cases. This means that the eigenstate $\tilde \Psi$ is strictly localized on a finite set of lattice sites.

In our example of the interacting Hadamard walk this happens exactly at the points where the modulus of the quasi-momentum p equals the quasi-energy ω (p = ± ω(p)). This property leads to an interesting feature of the considered quantum walk: by engineering the initial state of the quantum walk on a few lattice sites it is possible to generate states that are close to a momentum eigenstate at the points p = ± ω(p). Here, the spread of the initial state in the relative coordinate is bounded by the spread of the eigenstate $\tilde \Psi _x$ and in the center-of-mass direction by the desired accuracy of the momentum preparation.

With the help of the eigenstates $\tilde \Psi _\pm$ corresponding to the two branches ω_±, we can construct the 1D quantum walk mimicking the time evolution of the molecule explicitly. This is done by identifying the states $\tilde \Psi _\pm$ with the p-dependent eigenvectors of some 1D quantum walk with the same dispersion relation ω_± (see also the discussion in section 3.4).