---
abstract: 'We investigate numerically and theoretically the effect of spatial disorder on two-dimensional split-step discrete-time quantum walks with two internal “coin” states. Spatial disorder can lead to Anderson localization, inhibiting the spread of quantum walks, putting them at a disadvantage against their diffusively spreading classical counterparts. We find that spatial disorder of the most general type, i.e., position-dependent Haar random coin operators, does not lead to Anderson localization, but to a diffusive spread instead. This is a delocalization, which happens because disorder places the quantum walk to a critical point between different anomalous Floquet-Anderson insulating topological phases. We base this explanation on the relationship of this general quantum walk to a simpler case more studied in the literature, and for which disorder-induced delocalization of a topological origin has been observed. We review topological delocalization for the simpler quantum walk, using time-evolution of the wavefunctions and level spacing statistics. We apply scattering theory to two-dimensional quantum walks, and thus calculate the topological invariants of disordered quantum walks, substantiating the topological interpretation of the delocalization, and finding signatures of the delocalization in the finite-size scaling of transmission. Our results showcase how theoretical ideas and numerical tools from solid-state physics can help us understand spatially random quantum walks.'
author:
- 'János K. Asbóth'
- Arindam Mallick
bibliography:
- 'walkbib.bib'
title: 'Topological delocalization in the completely disordered two-dimensional quantum walk'
---

Introduction
============

Discrete-Time Quantum Walks (or quantum walks for short) are the quantum generalizations of random walks[@kempe2003quantum; @genske2013electric]. They are promising components of quantum algorithms because they spread faster than their classical counterparts[@ambainis2019quadratic]. These quantum walks can be described as periodic sequences of unitary “coin toss” and “shift” operations applied to a particle on a lattice. If the particle is initially on one of the lattice sites, during the quantum walk, its wavefunction spreads out over the lattice to infinity in a ballistic way, i.e., with root-mean-square distance from the origin proportional to time elapsed. This is faster than the diffusive spread characteristic of classical random walks. Spatial disorder in the quantum walk parameters can be detrimental to the fast spread of the walk. Intuition from solid-state physics suggest that the combination of disorder and coherence generically leads to Anderson localization [@nagaoka1985theory]: all eigenstates assume an exponentially localized envelope. As a consequence, a particle starting from a single site cannot spread off to infinity, its wavefunction remains within a bounded region, to exponential accuracy. Anderson localization can be tested numerically from the spectrum of the Hamiltonian, alternatively from the conductance, or from the spreading of the wavefunction[@sepehrinia2011numerical]. Anderson localization indeed happens in quantum walks [@joye2012dynamical; @vakulchyk2017anderson], and has already been observed in experiment[@schreiber2011decoherence; @pankov2019anderson]. We note this is very different than the dynamical effect of fluctuating disorder: that, by inducing a loss of coherence, would lead to diffusive spreading, i.e., a loss of the “quantum advantage”[@dur2002quantum; @schreiber2011decoherence]. This is also different from the so-called trapping effect in quantum walks[@kollar2015strongly; @machida2015localization; @kollar2020complete], which is a form of localization without disorder, related to the presence of flat bands in the quasienergy spectrum.

Another interesting feature of quantum walks is that they can have topological phases, of the kind known from the physics of topological insulators and superconductors[@hasan2010colloquium; @kitagawa2010exploring; @tarasinski2014scattering; @cedzich2018topological]. On the one hand, this allows quantum walks to be used as simulators for these solid state physics systems. On the other hand, quantum walks also have topological phases beyond those of topological insulators, which are prototypical of periodically driven, Floquet systems[@rudner2013anomalous; @asboth2015edge]. Thus, they can be used as versatile toy models for periodically driven systems in the nonperturbative limit of strong driving.

Recently, a striking disorder-induced *delocalization* phenomenon was observed for the simplest two-dimensional quantum walk[@edge2015localization], related to its topological phases. This quantum walk is the split-step walk on a square lattice, with two internal (coin) states, and a real-valued coin operator. Disorder, added to the walk via an onsite complex phase factor (mimicking onsite potential disorder), leads to an Anomalous Floquet-Anderson Insulator (AFAI) state[@titum2016anomalous], with Anderson localization of all bulk eigenstates yet with topologically protected edge states present in the spectrum. If the coin operator is finetuned to a critical value at the transition between two topological insulating phases, disorder does not lead to Anderson localization, but rather to a diffusive spread [@edge2015localization]. This route to delocalization also explains the lack of Anderson localization observed elsewhere for some two-dimensional walks[@zeng2017discrete]. Intriguingly, this delocalization can also be achieved if instead of finetuning the coin operator parameters, they are chosen to be maximally disordered[@edge2015localization]. This happens because coin disorder can tune the bulk topological invariant, and hence also put the system to a critical point – an example is the case where the coin rotation angle is uniformly random. We note that this localization-delocalization transition has a completely different physical origin than that induced by correlated disorder[@mendes2019localization].

In this article we ask whether we see Anderson localization or topological delocalization in the *completely disordered* general two-state quantum walk on the square lattice, where the coin is picked from all $U(2)$ operators in a Haar uniform random way. We employ a numerical tool used in previous work[@edge2015localization], namely, time evolution of the wavefunction, which is more efficient to compute for quantum walks than for static Hamiltonians. We also use other numerical tools, which have been applied to quantum walks less often: (1) analysis of the level spacing distribution to help identify Anderson localization and (2) calculation of the transmission matrix both to calculate the topological invariants with disorder and to find signatures in the finite-size scaling of Anderson localization and diffusion. We make use of the fact that when there is a complete phase disorder – as in the case when coin operators are Haar uniformly random – all properties of the spectrum should be quasienergy independent, and thus disorder averaging can be replaced/supplemented by quasienergy averaging, which results in a substantial numerical advantage. We find that the completely disordered (Haar random) two-dimensional quantum walk is topologically delocalized, for a similar reason as the more restricted quantum walk.

Our article is organized as follows. In Sect. \[sec:model\] we introduce the split-step two-dimensional quantum walk with two internal states, and discuss the topological properties in the clean, translational invariant setting. In Sect. \[sec:topo\_scattering\] we introduce the scattering matrix of the quantum walk, and use it to calculate the topological invariants in the case with disorder. In Sect. \[sec:phase\_disorder\] we compute time evolution, the level spacing statistics, and the finite-size scaling of the transmission, to show that disorder in the phase parameter leads to Anderson localization, except if the system is in a critical state. In Sect. \[sec:haar\_disorder\] we show using the numerical tools mentioned above that turning the disorder to the “most random” case, i.e., coin operator from Haar uniform random distribution, leads to delocalization. We conclude in Sect. \[sec:conclusion\], and show some preliminary investigation of the case of binary disorder in the coin angle parameter $\theta$ in Appendix \[sec:binary\_disorder\].

Two-dimensional quantum walk with two internal states {#sec:model}
=====================================================

The quantum walk we consider in this paper is the time evolution of a particle with two internal (spin) states on a square lattice. Its wavefunction reads, $$\begin{aligned}
{|{\Psi(t)}\rangle} &= \sum_{x} \sum_{y} \sum_{s=\pm 1} 
\Psi_{x,y,s}(t) {|{x,y,s}\rangle},\end{aligned}$$ where $x,y\in\mathbb{Z}$ are lattice coordinates and we are interested in the state at discrete times $t\in \mathbb{N}$. The dynamics is given by a periodic sequence of internal rotations (coin operator) and spin-dependent displacements (shift operator) on the lattice.

The internal rotations are general $U(2)$ operations acting on the internal degree of freedom, in a position-dependent way. We split off a position-dependent phase operator, $$\begin{aligned}
  \hat{F} &= \sum_{x,y}\sum_{s=\pm 1} 
  e^{i \phi(x,y)}{|{x, y, s}\rangle}{\langle{x,y, s}|},
  \label{eq:def_F}\end{aligned}$$ and $\hat{R}_j$ denotes the remaining SU(2) rotation operator, $$\begin{aligned}
  \hat{R}_j &= \sum_{x,y} {|{x, y}\rangle}{\langle{x,y}|} \otimes
  \hat{R}\left(\alpha_j(x,y), \beta_j(x,y),\theta_j(x,y)\right),\end{aligned}$$ where $j$ is used to differentiate between different rotations, and $$\begin{aligned}
  \label{eq:rotation_matrix}
  R(\alpha, \beta,\theta) &=
  \begin{pmatrix} e^{-i(\alpha+\beta)} \cos\theta & -e^{i(\alpha-\beta)} \sin\theta \\
    e^{-i(\alpha-\beta)} \sin\theta & e^{i(\alpha+\beta)} \cos\theta 
  \end{pmatrix}.\end{aligned}$$ Note that this can be written succinctly in Euler angle representation as $\hat{R} = e^{-i\beta \hat{\sigma}_z} e^{-i\theta \hat{\sigma}_y}
e^{-i\alpha \hat{\sigma}_z}$, where $\alpha$, $\theta$, and $\beta$, respectively, are the first, second and third Euler angles. Setting $\alpha=\beta=0$ reduces this coin operator to that used in several quantum walk papers[@kitagawa2010exploring; @kitagawa2012topological; @edge2015localization].

The shifts are spin-dependent translations on the lattice,

\[eq:shifts\_def\] $$\begin{aligned}
  \hat{S}_{x} &= \sum_{x,y}\sum_{s=\pm 1} 
  {|{x+s, y}\rangle}{\langle{x,y}|} \otimes{|{s}\rangle}{\langle{s}|};\\
  \hat{S}_{y} &= \sum_{x,y}\sum_{s=\pm 1} 
  {|{x, y+s}\rangle}{\langle{x,y}|} \otimes{|{s}\rangle}{\langle{s}|}.
  \label{eq:Sy_def}\end{aligned}$$

The wavefunction of the walk after $t\in\mathbb{N}$ timesteps reads $$\begin{aligned}
  {|{\Psi(t)}\rangle} &= \hat{U}^t {|{\Psi(0)}\rangle},\end{aligned}$$ where the timestep operator, $\hat{U}$, represents the effect of one period of the quantum walk, $$\begin{aligned}
  \label{eq:def_timestep}
  \hat{U} &= \hat{F} \hat{S}_{y} \hat{R}_2
  \hat{S}_{x} \hat{R}_1.\end{aligned}$$ As initial state we took ${|{\Psi(0}\rangle} = {|{x=0,y=0,s=+1}\rangle}$ – due to the disorder, the internal state of the initial condition plays no role in the time evolution, so any choice of initial state will give qualitatively the same results The parameters of the timestep operator are the position-dependent angle variables: a phase $\phi(x,y)$, and the rotation parameters $\theta_{1,2}(x,y)$, $\alpha_{1,2}(x,y)$, $\beta_{1,2}(x,y)$. In the formulas that follow we will often suppress the explicit position and spin dependence for better readability.

Sublattice symmetry of the quantum walk {#subsec:sublattice}
---------------------------------------

The two-dimensional quantum walk we consider here has a sublattice structure. There are four sublattices on the square lattice of $x,y\in\mathbb{Z}$, according to whether $x$ and $y$ are even ($e$) or odd ($o$), with corresponding projectors defined as, e.g., $$\begin{aligned}
  \label{eq:sublattice_projector_def}
  \hat{\Pi}_{e,o} = \sum_{x \text{ even}} \sum_{y \text{ odd}}
    {|{x,y}\rangle}{\langle{x,y}|},\end{aligned}$$ with $\hat{\Pi}_{o,e}$, $\hat{\Pi}_{o,o}$, and $\hat{\Pi}_{e,e}$ defined similarly. Since the quantum walk we consider has one shift along both $x$ and $y$ in every timestep, $\hat{U}$ switches sublattices, according to $(e,e) \leftrightarrow (o,o)$ and $(e,o)
\leftrightarrow (o,e)$. Therefore, a quantum walk started from one of the 4 sublattices never interferes with a walk started from any other sublattice – assuming that boundary conditions also respect the sublattice structure, e.g., are periodic with both $L_x$ and $L_y$ even. (We will later also use absorbing boundary conditions, which always respect the sublattice structure independent of system size.)

The sublattice structure has important consequences for the spectrum of $\hat{U}$, which can be appreciated by first considering the spectrum of $\hat{U}^2$. The operator $\hat{U}^2$ is block diagonal in the sublattice basis, i.e., $$\begin{aligned}
  \hat{U}^2 = \sum_{j=e,o} \sum_{l=e,o} \hat{\Pi}_{j,l} \hat{U}^2
  \hat{\Pi}_{j,l}.\end{aligned}$$ Take an eigenstate ${|{\Psi}\rangle}$ of $\hat{U}^2$ on the $(e,e)$ sublattice, $$\begin{aligned}
  \hat{U}^2 {|{\Psi}\rangle} &= e^{-2i\varepsilon} {|{\Psi}\rangle},
  \end{aligned}$$ with $0\le \varepsilon < \pi$. First, note that $\hat{U}{|{\Psi}\rangle}$ is on the $(o,o)$ sublattice, and is an eigenstate of $\hat{U}^2$ with the same eigenvalue $e^{-2i\varepsilon}$. Second, we can use ${|{\Psi}\rangle}$ to generate two eigenstates of $\hat{U}$, since $$\begin{aligned}
  \label{eq:eigstates_from_U2}
  \hat{U} &
  \left( {|{\Psi}\rangle} \pm e^{i\varepsilon} \hat{U} {|{\Psi}\rangle} \right) =
  \pm e^{-i\varepsilon}\left( {|{\Psi}\rangle}
  \pm e^{i\varepsilon} \hat{U} {|{\Psi}\rangle} \right).\end{aligned}$$ These relations also hold for an eigenstate ${|{\Psi'}\rangle}$ on the $(e,o)$ sublattice, with $\hat{U}{|{\Psi'}\rangle}$ on the $(o,e)$ sublattice. The results of the previous paragraph can be rephrased[@asboth2015edge] as *sublattice symmetry* of $\hat U$, represented by a unitary and Hermitian sublattice operator, $$\begin{aligned}
  \hat{\Gamma} \hat{U} \hat{\Gamma} &= -\hat{U}, \,\,\text{with}& 
  \hat{\Gamma} &= \hat{\Pi}_{e,e} + \hat{\Pi}_{e,o}
  - \hat{\Pi}_{o,o} - \hat{\Pi}_{o,e}. \end{aligned}$$ Every eigenstate ${|{\Phi}\rangle}$ of the walk has a sublattice symmetric partner, $$\begin{aligned}
  \label{eq:sublattice_partners}
  \hat{U}{|{\Phi}\rangle}&=e^{-i\varepsilon} {|{\Phi}\rangle} \quad \Rightarrow \quad
  \hat{U}\hat{\Gamma}{|{\Phi}\rangle}&=e^{-i(\varepsilon +\pi)} \hat{\Gamma} {|{\Phi}\rangle},\end{aligned}$$ which is related to Eq.  by ${|{\Phi}\rangle}=({|{\Psi}\rangle}
+ e^{i\varepsilon}\hat{U}{|{\Psi}\rangle})/\sqrt{2}$.

The sublattice symmetry we defined above is *not* the same as the sublattice symmetry familiar from solid state physics, which is also known as chiral symmetry. Chiral symmetry states $\hat\Gamma
\hat{H}\hat\Gamma = -\hat{H}$, linking eigenstates of a Hamiltonian at energy $\varepsilon$ to eigenstates at $-\varepsilon$. The quantum walk we consider here has chiral symmetry only if the parameters are finetuned. For chiral symmetry, we need $\phi=0$, and $\alpha_2(x,y)=\beta_1(x,y)$ and $\alpha_1(x,y) = \beta_2(x,y)$, moreover, either $\theta_1(x,y)=\theta_2(x,y) +n \pi$ (with chiral symmetry by $\hat{\sigma}_x$), or $\theta_1(x,y)=-\theta_2(x,y)+n \pi$ (with chiral symmetry by $\hat{\sigma}_y$).

Parameters $\alpha$ and $\beta$ represent a vector potential
------------------------------------------------------------

The parameters $\alpha_{1,2}(x,y)$ and $\beta_{1,2}(x,y)$ of the coin rotations can be understood to represent a vector potential[@yalccinkaya2015two; @arnault2016quantum; @sajid2018creating; @cedzich2019quantum; @Mallick:2019hvd]. This is best seen by writing the time-step operator in a different timeframe [@asboth2013bulk], which amounts to a similarity transformation on the original time-step operator of Eq. : $$\begin{aligned}
  \hat{U}_2 &=  e^{-i \alpha_1 \hat{\sigma}_z} \hat{U} e^{i \alpha_1 \hat{\sigma}_z}  = e^{-i\phi} \hat{S}_y'  e^{-i\theta_2\hat{\sigma}_y}
  \hat{S}_x'e^{-i\theta_1\hat{\sigma}_y};\\
 \hat{S}_y' &= e^{-i\alpha_1\hat{\sigma}_z}\hat{S}_y e^{-i\beta_2
   \hat{\sigma}_z};\quad  \hat{S}_x' =
 e^{-i\alpha_2\hat{\sigma}_z}\hat{S}_x e^{-i\beta_1 \hat{\sigma}_z}.
 \label{eq:shift_gauge}\end{aligned}$$

First, if $\alpha_{1,2}$ and $\beta_{1,2}$ are independent of position, then their effect on the quantum walk is just a gauge transformation. To prove this, we rewrite Eq.  as $$\begin{aligned}
 \hat{S}_x' &=
         e^{-i(\beta_1+\alpha_2)\hat{x}} \hat{S}_x
         e^{i(\beta_1+\alpha_2) \hat{x}},\end{aligned}$$ where $\hat{x}$ = $\sum_x x {|{x}\rangle}{\langle{x}|}$ is the position operator corresponding to position coordinate $x$. Thus $\hat{S}'_x$ is a unitary transformed version of $\hat{S}_x$, and the transforming operator commutes with rotations, phase operations, and $\hat{S}_y$ as well. Similar statements hold for $\hat{S}'_y$. Thus, in this case, the timestep operator $\hat{U}_2$ of a quantum walk is unitary equivalent to the timestep operator of the same walk with $\alpha_{1,2}=\beta_{1,2}=0$.

For $\alpha_{1,2}(x,y)$ and $\beta_{1,2}(x,y)$ depending on position, we find a more complicated situation. For this more general case, it is helpful to think of the quantum walk as a nearest-neighbor hopping model on the square lattice with time-dependent hopping amplitudes[@asboth2015edge]. The angles $\alpha_{1,2}(x,y)$ and $\beta_{1,2}(x,y)$ can then be included in this hopping model by the use of Peierls phases: for a link from site $(x,y)$ to $(x\pm 1, y)$, the corresponding Peierls phase is $\pm(\beta_1(x,y) + \alpha_2(x\pm
1,y))$, and from $(x,y)$ to $(x,y\pm 1)$, it is $\pm(\beta_2(x,y) +
\alpha_1(x\pm 1,y))$. A caveat: these Peierls phases depend on the sublattice.

The quantum walk is topological
-------------------------------

The quantum walk we consider has topological phases and associated chiral edge states similar to the simpler two-dimensional quantum walk. That simpler quantum walk (i.e., the vector-potential-like parameters $\alpha=\beta=0$), without spatial disorder has been recognized to host robust edge states, even though the Chern numbers of both of its quasienergy bands vanish[@kitagawa2012topological]. The bulk topological invariant explaining the presence of the edge states is a winding number, the RLBL-invariant[@rudner2013anomalous; @asboth2015edge; @asboth2017spectral; @sajid2018creating] we will call $\nu$. Since spatially constant vector-potential-like parameters $\alpha,\beta$ can be gauged away, we can conclude that our quantum walk has a winding number that depends only on the values of $\theta_1$ and $\theta_2$. The invariant follows a harlequin pattern[@asboth2015edge], which can be put in a concise (although somewhat obfuscated) formula, $$\begin{aligned}
  \label{eq:clean_topo_number}
  \nu &= \text{sgn} \left[ \sin(\theta_1 - \theta_2)
    \sin(\theta_1+\theta_2) \right],\end{aligned}$$ with $\nu=0$ signaling the critical, gapless cases, where the topological invariant is not well defined. Since we would like to study disordered quantum walks, we need to find a way to calculate this topological invariant with all parameters of the quantum walk spatially dependent.

Topological invariant of the disordered quantum walk using scattering theory {#sec:topo_scattering}
============================================================================

We will compute the bulk topological invariant for disordered quantum walks by detecting the edge states, in a scattering setup, borrowing ideas from numerics on Floquet systems[@fulga2016scattering; @rodriguez2019topological; @liu2020anomalous]. Scattering theory has already been applied to one-dimensional quantum walks to obtain topological invariants even in the disordered case[@tarasinski2014scattering; @rakovszky2015localization], and such quantized signatures have even been measured experimentally[@barkhofen2017measuring]. We will first sketch the concept, and then give the concrete recipe for numerical implementation.

Conceptual setup and implementation {#subsec:scattering_setup}
-----------------------------------

To calculate the topological invariant, we need to (1) define a scattering setup for the quantum walk with semi-infinite leads attached, (2) implement periodic and open boundary conditions in the transverse direction, i.e., implement a “cut” connecting the two leads, (3) define the scattering matrix, and then (4) give a recipe for its calculation.

\(1) We start with a scattering setup. We take a rectangle, $1 \le x
\le L_x$ and $1 \le y \le L_y$ (with $L_y$ even), as the “system”, and semi-infinite extensions along $x$ as the leads. For later use, we define the projector $\hat{\Pi}_\text{sys}$ to the system as, $$\begin{aligned}
  \hat\Pi_\text{sys} &= \sum_{x=1}^{L_x}\sum_{y=1}^{L_y}\sum_{s=\pm
  1} {|{x,y,s}\rangle}{\langle{x,y,s}|}.\end{aligned}$$ The left lead is at $x\le 0$ and $1\le y \le L_y$, while the right lead is $L_x<x$ and $1\le y \le L_y $. In the leads, the quantum walk is simplified, with all coin parameters set to 0, and the $y$-shift omitted. Thus in the leads, the quantum walk just propagates the particle right if its spin is up (+1), and left if its spin is down (-1). We use periodic boundary conditions along $y$, i.e., in the $\hat{S}_y$ shift, Eq. , $y+s$ should be replaced by $(y+s-1) \text{ mod } L_y +1$.

\(2) To realize edge states, we need to define a cut region connecting the two leads, a quantum walk equivalent of open boundary conditions. We do this in the simplest possible way, while not breaking sublattice symmetry[@asboth2015edge]: the cut is a row of sites, $y=y_\text{cut}$ (or more rows of sites can be used to the same effect; we used two rows in the simulation, see black sites in Fig. \[fig:scatter\_pic\]), where the angle parameters are set so that a walker impinging on the cut is certainly reflected back from it. We have two different choices here,

\[eq:def\_cuts\] $$\begin{aligned}
  \text{cut $A$:}\quad \theta_1(x,y_\text{cut})=0,&\quad
  \theta_2(x,y_\text{cut})=\frac{\pi}{2};\\
  \text{cut $B$:}\quad \theta_1(x,y_\text{cut})=\frac{\pi}{2},&\quad
  \theta_2(y=y_\text{cut})=0.\end{aligned}$$

To be specific we set all other coin parameters in the cut to $\alpha_j=\beta_j= 0$, although giving them any other value would not have any important effect. These two choices of coin parameters both ensure reflection – constituting a “bulk” with all flat bands[@kitagawa2012topological] – but have different topological invariant: $\nu_A=-1$ for the first case and $\nu_B = +1$ for the second case. Thus, if the bulk of the scattering region has a topological invariant $\nu$ that is different from $\nu_\text{cut}$, with “cut” being A or B, we expect chiral edge states in the bulk propagating in opposite directions directly above and below the cut, ${\left|{\nu-\nu_\text{cut}}\right|}$ of them in both directions. (3) To define the scattering matrix we need to construct scattering eigenstates of the quantum walk, and isolate their reflected/transmitted parts. This is done in the same way as for the one-dimensional quantum walk[@tarasinski2014scattering], and so we only give the definition of the process here, and refer the readers to Ref.  for details. At quasienergy $\varepsilon$ (i.e., eigenvalue $e^{-i\varepsilon}$), there are $L_y$ scattering states, eigenstates of the quantum walk , originating from a particle incident from the left at position $y=n$ (with $1 \le n \le L_y$), which read, $$\begin{aligned}
  {|{\Psi_{n,\varepsilon}}\rangle} &= \sum_{t=-\infty}^\infty
    e^{i\varepsilon t} \hat{U}^t  {|{0,n,+1}\rangle}.\end{aligned}$$ The reflection (transmission) matrix element $\mathfrak{r}_{mn}$ ($\mathfrak{t}_{mn}$) is the probability amplitude of the part of the scattering state that is in the left (right) lead, propagating towards the left (right), at $y=m$. These can be obtained as $$\begin{aligned}
 \mathfrak{r}_{mn}(\varepsilon) &= {\langle{0,m,-1}|{\Psi_{n,\varepsilon}}\rangle};\\
 \mathfrak{t}_{mn}(\varepsilon) &=  {\langle{L_x+1,m,+1}|{\Psi_{n,\varepsilon}}\rangle} .\end{aligned}$$ The eigenvalues of the transmission matrix, $\mathfrak{t}(\varepsilon)^\dagger \mathfrak{t}(\varepsilon)$, are the transmission eigenvalues. In the Landauer-Büttiker formalism[@nazarov2009quantum], quantum transport takes place via independent channels, and there are $L_y$ of these in this system, since there are $L_y$ right-propagating modes in the lead at each quasienergy. The degree to which each channel is open for transport is quantified by the corresponding transmission eigenvalue[@nazarov2009quantum]. The total transmission at a given quasienergy is therefore the trace of the transmission matrix, $$\begin{aligned}
  T(\varepsilon) &= \sum_{mn} {\left|{\mathfrak{t}_{mn}(\varepsilon)}\right|}^2.
  \label{eq:total_transmission_def}\end{aligned}$$ (4) Having sketched the conceptual definition of the scattering matrix, we need to circumvent the need for semi-infinite leads, as in the one-dimensional case[@tarasinski2014scattering]. First, note that the leads are nonreflective, and thus the part of the wavefunction exiting the system can safely be projected out at the end of every timestep. For a quantum walk on an $L_x\times L_y$ rectangle, we can thus realize the (relevant part of the) $x$-leads by a single column of sites at $x=0, y=1,\ldots,L_y$. We now take periodic boundary conditions along $x$ (as well as $y$), i.e., in the $\hat{S}_x$ shift, Eq. , replace $x+s$ by $x+s \text{ mod } (L_x+1)$. However, at the end of every timestep, we first read out the contents of the extra column $x=0$ – with $s=+1$ being the transmitted part and $s=-1$ the reflected part of the wavefunction – and then erase it.

Summarizing all this, the recipe for the transmission matrix reads, $$\begin{aligned}
  \mathfrak{t}_{mn}(\varepsilon) &= 
           \lim_{t_{\text{max}}\to\infty} 
  \sum_{t=0}^{t_\text{max}} e^{i\varepsilon t} \mathfrak{t}_{mn}(t);\\
  \mathfrak{t}_{mn}(t+1) &= {\langle{0,m,+1}|}
  \,\hat{U} \,(\hat{\Pi}_\text{sys} \hat{U})^t \,{|{0,n,+1}\rangle}.\end{aligned}$$ From this, the total transmission at any quasienergy $\varepsilon$ can be calculated using Eq. .

Topological invariants from transmission
----------------------------------------

We can infer the bulk topological invariant via from the presence and number of topologically protected edge states at any quasienergy $\varepsilon$ via the calculation of the transmission with and without cut. First, by calculating the transmission without cut, we can check if the system is insulating, i.e., $$\begin{aligned}
    \text{no cut}:&\quad
  T(\varepsilon) \underset{L\to \infty}{\longrightarrow} 0,\end{aligned}$$ where $L$ gives the scale of the system size, i.e., $L_x/L$ and $L_y/L$ constant. This requirement is a prerequisite for topologically protected edge states. Then, we calculate the transmission with either type of cut, which will only be due to the edge states above or below the cut (whichever carries current from left to right), $$\begin{aligned}
    \text{with cut A/B}:&\quad
  T_{A/B}(\varepsilon) \underset{L\to \infty}{\longrightarrow} {\left|{\nu_{A/B}-\nu}\right|},\end{aligned}$$ where we use $T_A$ to denote total transmission with cut $A$, and $T_B$ for cut $B$. We can thus infer the bulk topological invariant, $$\begin{aligned}
  \label{eq:nu_from_T}
  \nu &= \text{sgn}(T_A - T_B) \frac{T_A+T_B}{2},\end{aligned}$$ where the dependence on quasienergy $\varepsilon$ was suppressed for readability, and the sign function is defined to be $\text{sgn}(x)=+1$ if $0<x$, $-1$ if $x<0$, and $0$ if $x=0$. Note that because of the sublattice symmetry of the quantum walk, for an input from a mode with even $n$, the transmission will only be to a mode with odd/even $m$ (and have nonvanishing amplitudes for even/odd $t$) for odd/even $L_x$.

\

Topological invariants of quantum walks with maximal phase disorder {#subsec:top_inv_disorder}
-------------------------------------------------------------------

With maximal phase disorder – i.e., $\phi$ in Eq.  chosen for each site randomly and uniformly from the interval $[-\pi,\pi)$ – averaging quantities over quasienergy amounts to some disorder averaging. This is because phase disorder randomizes the quasienergy of each eigenstate, and therefore all quantities should be on average quasienergy independent. This goes for the spectrum of transmission eigenvalues, the total transmission, even the shapes of the edge states’ wavefunctions. The quasienergy-averaged value of the total transmission, as per Eq. , thus reads, $$\begin{aligned}
  \label{eq:average_total_transmission}
  \overline{T} &=\lim_{t_{\text{max}}\to\infty}
    \sum_{m=1}^{L_y}  \sum_{n=1}^{L_y}  \sum_{t=0}^{t_\text{max}}
       {{\left|{\mathfrak{t}_{mn}(t)}\right|}^2}.\end{aligned}$$

We illustrate how the method of obtaining topological invariants from total transmission with and without cuts works by numerical results, for the two-dimensional quantum walk with position-independent $\theta$, $\alpha$ and $\beta$ parameters, but completely disordered $\phi$. We expect from previous work[@edge2015localization] that the topological invariant should be given by Eq. . We calculated transmission for $t_\text{max}=2000$ timesteps, on a system size of $L_x = 39$ and $L_y = 60$, considering only even input channels (odd channels have identical contribution). The results, shown in Fig. \[fig:topo\_trans\], largely confirm our expectations. We see anomalous Floquet-Anderson insulator phases characterized by low transmission in the case without cuts, and to a good approximation quantized transmission in the case with cuts (Total transmission of 1 for even input channels is shown, same results for odd input channels were obtained), topological invariants matching up with Eqs.  and . These insulating phases are separated by lines of critical states, where transmission is high both with and without cuts. We defer a more detailed analysis of the finite-size scaling of the transmission to the next Section.

Disorder in phase and magnetic parameters {#sec:phase_disorder}
=========================================

Before turning to the two-dimensional quantum walk with all parameters random, we revisit the problem of a quantum walk with fixed $\theta$ parameters and fully random $\alpha, \beta$ or $\phi$. When disorder is only in the phase $\phi$, this reduces to the case of the split-step quantum walk, studied in Ref. , since position independent angle parameters $\alpha$ and $\beta$ can always be gauged away. There, it was found that in the simple split-step walk, with angle parameters $\theta_1, \theta_2$ fixed to generic values, phase disorder leads to Anderson localization. On the other hand, when $\theta_1, \theta_2$ is tuned to a topological phase transition, phase disorder leads to a diffusive spread of the wavefunction – a consequence of the critical nature of the system.

The numerical results in this Section complement those of Ref.  by (1) a more in-depth analysis of the way the wavefunction of an initially localized particle spreads, (2) different numerical tools, the analysis of level repulsion statistics and of the finite-size scaling of transmission. We will also investigate the effect of having disorder not in $\phi$, but in the vector-potential-like parameters $\alpha$ and $\beta$. This is very similar to having onsite phase disorder but in a spin-dependent way[@edge2015localization]. Before showing our numerical results, we discuss a technical detail that makes the simulation more efficient: a rotated basis.

Rotated basis
-------------

For the time evolution of the quantum walk with fixed $\theta_{1,2}$, we used a rotated basis, with square-shaped unit cells containing two sites each from different sublattices, as shown in Fig. \[fig:basis\_tran\]. Because of the sublattice symmetry, for a walk started from one sublattice, at any time during its time evolution, its wavefunction will only have support on at most one site per unit cell. Thus it is enough to keep track of the integer valued unit cell indices $n_+$ and $n_-$, and discard the sublattice information. On the $(e,e)$ and $(o,o)$ sublattices, these unit cell indices, and the corresponding coordinates $x_\pm$, are related to the integer valued site coordinates $x,y$, by $$\begin{aligned}
  n_\pm &= \frac{y\pm x}{2};&
  x_\pm &= \frac{y\pm x}{\sqrt{2}}.\end{aligned}$$ When rewriting the timestep operator of the quantum walk in this rotated basis, we need to start with the shifts defined in Eq. . For the quantum walk started from the $(e,e)$ or $(o,o)$ sublattice (a square in Fig. \[fig:basis\_tran\]), we need

$$\begin{aligned}
  \hat{S}_{1\square} &=  \sum_{n_+,n_-} \Big(
  {|{n_+, n_- - 1, +1}\rangle} {\langle{n_+,n_-,+1}|} \nonumber \\
  \quad& \quad \quad + {|{n_+-1, n_-, -1}\rangle} {\langle{n_+,n_-,-1}|}
  \Big);\\
  \hat{S}_{2\square} &=  \sum_{n_+,n_-} \Big(
  {|{n_++1, n_-+1, +1}\rangle} {\langle{n_+,n_-,+1}|} \nonumber \\
  \quad & \quad \quad + {|{n_+, n_-, -1}\rangle} {\langle{n_+,n_-,-1}|}
  \Big),
  \end{aligned}$$

\[eq:S\_diag\_def\]

and obtain the timestep operator as $$\begin{aligned}
  \label{eq:quantum_walk_ee}
  \hat{U}_\square &= \hat{F} \hat{S}_{2\square} \hat{R}_2
  \hat{S}_{1\square} \hat{R}_1,\end{aligned}$$ where the parameters of the operators $\hat{F}$, $\hat{R}_j$ have to be chosen to match the position of the walker.

For the sake of completeness, a quantum walk started from the $(o,e)$ or $(e,o)$ sublattices (circle in Fig. \[fig:basis\_tran\]) has modified shifts,

$$\begin{aligned}
  \hat{S}_{1\circ} &=  \sum_{n_+,n_-} \Big(
  {|{n_++1, n_-, +1}\rangle} {\langle{n_+,n_-,+1}|} \nonumber \\
  \quad& \quad \quad + {|{n_+, n_-+1, -1}\rangle} {\langle{n_+,n_-,-1}|}
  \Big);\\
  \hat{S}_{2\circ} &=  \sum_{n_+,n_-} \Big(
  {|{n_+, n_-, +1}\rangle} {\langle{n_+,n_-,+1}|} \nonumber \\
  \quad & \quad \quad + {|{n_+-1, n_--1, -1}\rangle} {\langle{n_+,n_-,-1}|}
  \Big),
  \end{aligned}$$

\[eq:Scirc\_diag\_def\]

and timestep operator, $$\begin{aligned}
  \label{eq:quantum_walk_eo}
  \hat{U}_\circ &= \hat{F} \hat{S}_{2\circ} \hat{R}_2
  \hat{S}_{1\circ} \hat{R}_1.\end{aligned}$$

With the rotated basis, a factor of two is gained in efficiency for the representation of the wavefunction, i.e., an array of size $L_+\times L_-$ is enough for a simulated area of $L_+ \times L_-$ unit cells with $2 L_+ L_-$ sites. Besides, we will show later that during time evolution, the disorder-averaged position distribution is generically anisotropic, and elongated either along the diagonal or along the anti-diagonal direction (in the original, $x-y$ basis), depending on the values of $\theta_1$ and $\theta_2$. Thus, to minimize finite-size effects, it is more efficient to use rectangular rather than square shaped simulation area, elongated along the direction in which the position distribution is elongated, i.e., rotated by $45^\circ$. This is straightforward to implement in the rotated basis.

Time evolution of wavefunction
------------------------------

A direct test of Anderson localization is simulation of the time evolution of the wavefunction for a particle started from a single site. In case of Anderson localization, the root-mean-squared displacement (distance from the origin, square root of the variance of position) remains bounded, and the long-time limiting form of the disorder-averaged probability distribution falls off exponentially with distance from the origin. In contrast, for a critical system, diffusive spread of the wavefunction is expected, with position variance increasing linearly with time, and the disorder-averaged probability distribution approaching a Gaussian shape. Both signatures for both cases have been observed for the split-step quantum walk with phase disorder[@edge2015localization], although for a limited range of parameters, with $\theta_1+\theta_2=\pi/2$. As we show below, that departing from this limitation on the rotation angles $\theta$ changes the dynamics qualitatively, introducing anisotropy.

For the numerical simulation, we found it advantageous to work with absorbing boundary conditions. Absorbing boundary conditions, as already used for the scattering calculation in Sect. \[subsec:scattering\_setup\], simply means setting the value of the wavefunction at the boundary to 0 at the end of every timestep. This allows us to track the finite-size error, or the “leaving probability” directly by the norm of the wavefunction, $$\begin{aligned}
  \label{eq:pleave_def}
  p_\text{leave} = 1 - \sum\limits_{n_+ = 1}^{L_+} \sum\limits_{n_- =
    1}^{L_-} \sum\limits_{s = \pm 1}|\psi_{n_+, n_-, s}
  (t_\text{max})|^2.\end{aligned}$$ Note that the absorbing boundary conditions as defined here are a “quick-and-dirty” way to emulate a finite segment of an infinite plane: in fact, a small part of the wavefunction is reflected back from the absorbing boundaries. Such reflections are numerical artifacts, which we can estimate by monitoring the error $p_\text{leave}$.

We have noticed that the disorder-averaged shape of the position distribution as the wavefunction of the quantum walk spreads is anisotropic, with contour lines having elliptical shapes, ellipses elongated either along the diagonal ($x+y$) or along the antidiagonal ($x-y$) direction. Typical examples are shown in Fig. \[fig:fig1\_evolve\]. The direction of elongation depends on the parameters $\theta_j$, following a harlequin pattern, which can be put in a concise (although somewhat obfuscated) formula: the dominant direction is $$\begin{aligned}
  x \pm y, \quad \text{with }\,\, \pm &= \text{sgn} \left[ \cos(\theta_1 -
    \theta_2) \cos(\theta_1+\theta_2) \right],
  \label{eq:dominant_axis}\end{aligned}$$ with an isotropic position distribution if $\cos(\theta_1 - \theta_2)
\cos(\theta_1+\theta_2)=0$. We do not have a complete understanding of why the elliptical contour lines are not rotated by any angle other than $\pm 45^\circ$. There are, however, two special cases that are straightforward to check by considering two consecutive timesteps (as in [@asboth2012symmetries] for a one-dimensional quantum walk). If either $\theta_1$ or $\theta_2$ is set to $n \pi$, the quantum walk is one-dimensional, spreading only along the diagonal. If either $\theta_1$ or $\theta_2$ is set to $(n+1/2) \pi$, it spreads only along the anti-diagonal (hence if one of them is 0, the other $\pi/2$, no spreading at all). This is consistent with Eq. .

We now turn to the rate at which the position distribution spreads and its shape in the long-time limit, to find signatures of Anderson localization vs. diffusion. For generic cases, i.e., when $\theta_1 \neq \pm\theta_2+n\pi$, with $n\in\mathbb{N}$, with maximal phase disorder, we find signatures of Anderson localization. The envelope of the disorder-averaged position distribution in the long-time limit can be fitted very well with $$\begin{aligned}
  p(x_+,x_-,t) &\propto
  \exp\left[- \sqrt{ \left(\frac{x_+}{\zeta_+(t)}\right)^2
    + \left(\frac{x_-}{\zeta_-(t)}\right)^2 }\, \right],
  \label{eq:ansatz_localized}\end{aligned}$$ where $\zeta_+$ and $\zeta_-$ increase slowly with time (we expect them to saturate in the long-time limit); these denote the localization lengths along the diagonal and anti-diagonal directions, with the variance of the position given by ${\langle{x^2+y^2}\rangle} \approx
3(\zeta_+^2+\zeta_-^2)$. By “envelope”, we mean that on every second site the wavefunction is 0 because of sublattice symmetry. Eq.  is an exponentially localized form, with contours that have elliptic shapes, tilted by $45^\circ$. We show an example in Fig. \[fig:fig1\_evolve\] (a), where $\theta_1=0.2\pi$, $\theta_2=0.4\pi$, and the fitted values of the localization lengths are $\zeta_+=3.7$ and $\zeta_-=6.7$ (with maximal disorder in the magnetic parameters as well, $\zeta_+=3.5$ and $\zeta_-=6.8$). For the critical cases, i.e., when $\theta_1 = \pm\theta_2+n\pi$, the wavefunction spreads out in a diffusive way. We find that a good fit for the envelope is $$\begin{aligned}
  p(x_+,x_-,t) \propto
  \exp\left[- \frac{x_+^2}{4D_+ t}  - \frac{x_-^2}{4D_- t} \right],
  \label{eq:ansatz_diffusive}\end{aligned}$$ with $D_+$ and $D_-$ denoting the diffusion coefficients along the diagonal and antidiagonal directions, with the variance of the position given by ${\langle{x^2+y^2}\rangle} \approx 2(D_+ + D_- )t$. We show an example in Fig. \[fig:fig1\_evolve\] (b), where $\theta_1=\theta_2=0.2\pi$, and the fitted values of the diffusion coefficients are $D_+=1.1$ and $D_-=0.31$ (unaffected by maximal disorder in the magnetic parameters).

![Diagonal (D, along $x_-=0$, only even $x_+$) and antidiagonal (AD, along $x_+=0$, only even $x_-$) cuts of the position distributions of Fig. \[fig:fig1\_evolve\] (semi-log plot). Phase disorder (suffix $\phi$) and disorder in magnetic parameter (suffix $\alpha,\beta$) give almost indistinguishable results. Continuous lines show the fitted theory curves - they are almost completely covered by the numerical data points, indicating excellent fit, except for the close vicinity of the origin. The inset shows the time dependence of the square root of the variance (log-log plot), which scales diffusively for the critical case, and appears to approach saturation for the generic case.[]{data-label="fig:diagonal_ad_cuts"}](position_distribution_fitting){width="\columnwidth"}

We show more details on diagonal and antidiagonal cuts of the disorder-averaged position distributions in Fig. \[fig:diagonal\_ad\_cuts\]. We took these cuts from the simulation runs represented in Fig. \[fig:fig1\_evolve\]. We note that the parameters $\zeta_\pm$ and $D_\pm$ in Eqs.  and were obtained by fitting the analytical curves to these cuts. The fits, as already seen in Fig. \[fig:fig1\_evolve\], are quite good, except for the diffusive case in the vicinity of the origin, where we observe a spike. We also plot in Fig. \[fig:diagonal\_ad\_cuts\] the numerical results for cases with the maximal disorder taken in $\alpha$ and $\beta$, rather than in $\phi$; the effects of these two types of disorder are the same, consistent with previous results[@edge2015localization]. In the inset we show the measured time dependence of the root-mean-squared width. For the critical case, $\theta_1=\theta_2=0.2\pi$, this is a linear function on the log-log plot, with $\sqrt{{\langle{x^2+y^2}\rangle}} = 1.5
t^{0.5}$, i.e., diffusive scaling. For the generic case, $\theta_1=0.2\pi$, $\theta_2=0.4\pi$, the variance grows slower, consistent with the expectation that it would eventually saturate.

Level spacing statistics
------------------------

A frequently used tool to characterize Anderson localization/criticality is the level spacing statistics[@sepehrinia2011numerical]. As often applied in Hamiltonian systems, first an energy $E$ is fixed, and then from each disorder realization the gap around $E$ is taken, i.e., $\delta_j = \varepsilon_{+,j}-\varepsilon_{-,j}$, with $\varepsilon_{\pm,j}$ denoting the first level above/below $E$ in the $j$th disorder realization. The ensemble of normalized level spacings $s_j$ is then defined as $$\begin{aligned}
  s_j &= \frac{\delta_j}{{\langle{\delta}\rangle}},&\end{aligned}$$ with ${\langle{\ldots}\rangle}$ denoting the disorder average.

In an Anderson localized system, eigenstates separated by large distances cannot be coupled by local perturbations, therefore their energies (in our case, quasienergies) are essentially independent. Thus the normalized level spacing has exponential probability distribution (in this context also called Poissonian), $$\begin{aligned}
\label{eq:levelspacing_localized}
  p(s) &= p_\text{loc}(s) = e^{-s},\end{aligned}$$ with $p(s)$ denoting the probability density. In a critical system, on the other hand, the eigenstates have extended wavefunctions, which can be coupled by local perturbations and hybridize, and therefore we expect level repulsion to occur. The normalized level spacings in this case are expected to follow a Wigner–Dyson distribution, which for the case relevant for us (so-called Gaussian Unitary Ensemble) reads, $$\begin{aligned}
\label{eq:levelspacing_gue}
  p(s) &= p_\text{GUE}(s) =
  \frac{32}{\pi ^{2}} s^{2} e^{{-{\frac {4}{\pi }}s^{2}}}.\end{aligned}$$

When we compute level spacing statistics for quantum walks, we need to pay attention to the sublattice structure (see Sec. \[subsec:sublattice\]). The timestep operator $\hat{U}$ of a walk describes two independent quantum walks, taking place on different sublattices ($x+y=$ even vs $x+y=$ odd). Thus $\hat{U}$ has two sets of energy levels whose level spacing distributions should be calculated separately (just as in the case of Hamiltonians with unitary symmetries): there can be no level repulsion between levels from different sets. To account for this, we start from the spectrum of $\hat{U}^2$, to obtain that of $\hat{U}$, as in Sec. \[subsec:sublattice\]. For a quantum walk on $2N$ sites (with $N$ even, and with proper boundary conditions), we need to diagonalize the $N \times N$ blocks of $\hat{U}^2$ on the $(e,e)$ sublattice (denoted by $\square$ in Fig.  \[fig:basis\_tran\]), and on the $(e,o)$ sublattice (denoted by $\circ$). From the corresponding eigenvalues of $\hat{U}^2$, namely, $e^{2i\varepsilon^\square_j}$ and $e^{2i\varepsilon^\circ_j}$, with $1\le j < N$, and $0\le \varepsilon^{\square/\circ}_j < \pi$, we obtain the $2N$ independent quasienergy level spacings of the spectrum of $\hat{U}$ as $$\begin{aligned}
  \label{eq:delta_circ_def}
  \delta^{\square}_j &= (\varepsilon^{\square}_{(j+1) \text{ mod } N} -
  \varepsilon^{\square}_{j}) \mod \pi,\end{aligned}$$ with $j=1,\ldots,N$, and with analogous definitions for $\delta_j^\circ$.

Having maximal phase disorder, i.e., $\phi$ of Eq.  equidistributed in the interval $[-\pi,\pi)$, gives us a huge boost in numerical efficiency for obtaining the level statistics. As already discussed in Sect. \[subsec:top\_inv\_disorder\], we can treat all of the quasienergies on the same footing, and for a quantum walk on $2L_+ L_- = 2N$ sites obtain $2N$ values of normalized level spacing as $$\begin{aligned}
  s_j &= \frac{N}\pi \delta_j^\square;&
  s_{j+2N} &= \frac{N}\pi \delta_j^\circ,\end{aligned}$$ with the level spacings $\delta^\square_j$ and $\delta^\circ_j$ for $j=1,\ldots,N$ defined in Eq. . Thus, we obtain a level spacing ensemble of size $2N$ by only diagonalizing two unitary matrices – a boost in numerical efficiency. On the downside, to obtain all the level spacings we fully diagonalize these large unitaries of size $N\times N$, where $N=L_+ L_-$, with a numerical cost of $O(N^3)$. This constrains us to system sizes of the order of $N\approx 20000$, which turns out to be large enough for the cases we considered here. We could go beyond that by repeatedly sampling the spectrum of $\hat{U}^2$ using iterative algorithms that obtain only part of the spectrum, such as the Arnoldi method[@lehoucq1998arpack] (which works applies ideas used in the Lanczos algorithm to more general matrices).

We show numerically obtained level spacing distributions, and the theoretical expectations for them, in Fig. \[fig:levelspacings\_1\]. Here we used the same examples as in Fig. \[fig:diagonal\_ad\_cuts\]: a generic quantum walk with $\theta_1 = 0.2 \pi, \theta_2 = 0.4\pi$, and a critical quantum walk with $\theta_1 = \theta_2 = 0.2\pi$ – in both cases, results collected from a single sublattice are shown. For both cases we considered two types of maximal disorder: disorder in $\phi$ and $\alpha,\beta$ fixed, or the other way round. For the critical case, we see excellent agreement between the level spacing distributions and the Wigner-Dyson distribution of Eq. . For the generic, Anderson localized case, we observe a very good agreement with the exponential distribution of Eq. , with some level repulsion showing up at the $s\ll 1$ end – states that are almost degenerate – which we attribute to finite size effects. Using the diagonal basis for these calculations, allowed us to adapt the shape of the simulated area to the expected shape of the wavefunctions, and thus minimize finite-size effects.

Scaling of transmission {#subsec:scaling_transmission}
-----------------------

The calculation of the transmission matrix of a disordered two-dimensional quantum walk gives us yet another numerical tool to differentiate between Anderson localization and diffusive spread. Scaling up the system size while keeping the system shape constant, total transmission across an insulator should decrease exponentially, while for the diffusive case (e.g., a metal), we expect a total transmission that is roughly constant[@nazarov2009quantum].

For the finite-size scaling of the transmission, we calculated quasienergy-and disorder-averaged total transmissions for three different system sizes with the same shape of $L_y/L_x=3/2$ (together with the extra column of sites for the leads), for different values of $\theta_1$ and $\theta_2$. We kept $\theta_1+\theta_2=0.6\pi$, and used values of $\theta_1$ from $\pi/10$ to $\pi/2$, thus tuning the system across a topological phase transition which takes place at $\theta_1=\theta_2=0.3 \pi$. For each value of $\theta_1,\theta_2$, we calculated transmission for three different system sizes, $L_x=19, L_y=30$; $L_x=39, L_y=60$; and $L_x=59, L_y=90$. We used Eq. , for the numerics, with the number of timesteps $t_\text{max}$ chosen such that the incident walker from any input lead will have left the system with probability above $99\%$, namely, $t_\text{max} = 1000$ for $L_x=19,
L_y=30$; $t_\text{max} = 2000$ for $L_x=39, L_y=60$; and $t_\text{max}
= 3000$ for $L_x=59, L_y=90$. This criterion for choosing $t_\text{max}$ can be written using Eq.  as $p_\text{leave}<0.01$, with the conditional wavefunction $ {|{\Psi_n(t+1)}\rangle} = \hat{U}
\left(\hat{\Pi}_\text{sys}\hat{U}\right)^t{|{0,n,+1}\rangle}$. Our numerical results for the quasienergy-averaged total transmission, shown in Fig. \[fig:transmission\_phasetransition\], confirm that tuning $\theta_1$ and $\theta_2$ indeed drives the quantum walk with complete phase disorder across a quantum phase transition between different localized phases. We observe that for $\theta_1\neq\theta_2$, the total transmission decreases exponentially (see semilogarithmic plot in the inset). For $\theta_1=\theta_2=0.3\pi$, on the other hand, the transmission is unchanged as the system size is scaled up, a numerical signature of diffusive transport.

Haar random quantum walk is diffusive {#sec:haar_disorder}
=====================================

We finally turn to the question raised in the title: is the completely disordered two-dimensional quantum walk Anderson localized, or does it spread out to infinity? Thus, we take a two-dimensional quantum walk as in Eq. , with both rotation operators $\hat{R}_j$ chosen randomly from $U(2)$ according to the Haar measure[@medrazzi2006generate_random], i.e., Haar random operators. We can realize Haar random coins with the operators of Eqs.  and , by taking their parameters from properly defined distributions[@ozols2009generate]. The parameters $\alpha_j(x,y), \beta_j(x,y)$, with $j=1,2$, and $\phi(x,y)$ have to be uncorrelated random variables, uniform in the interval $[-\pi,\pi)$. The parameters $\theta_j(x,y)$ need to be generated as $$\begin{aligned}
  \label{eq:haar_random}
 \theta_j(x,y) &= \arcsin\left(\sqrt{\zeta_j(x,y)}\right)\end{aligned}$$ with $\zeta_j(x,y)$ uncorrelated uniform random in the interval $[0,1]$.

In Fig. \[fig:haar\_evolve\] the position distribution is shown after 10000 timesteps on a system of size $1001\times 1001$. Even without disorder averaging, the distribution is quite smooth, and roughly isotropic. It corresponds to a Gaussian, can be fitted quite well with the diffusive ansatz of Eq. , with diffusion coefficients $D_+=D_-=0.54$. In Fig. \[fig:haar\_cuts\], we show the fit with the diffusive curve of the cross-sectional cut of the disorder-averaged position distribution (from 200 random realizations, after 1500 timesteps) of the Haar random quantum walk, It is only near the origin that the Gaussian fit is not very good; as shown in the inset, we find here a pronounced peak, just as with the phase disordered quantum walks of Fig. \[fig:diagonal\_ad\_cuts\]. The simulations of Fig. \[fig:haar\_cuts\] were run on systems of $501 \times 501$ unit cells in the diagonal basis, with absorbing boundary condition, error $p_\text{leave}<10^{-15}$.

We also have numerical evidence – shown in Fig. \[fig:haar\_cuts\] – that coherence plays almost no role in the way the Haar random quantum walk spreads. We show on the plot the position distributions of two classicalized variants of the quantum walk, after the same number of timesteps, on a same system size. The first classicalized variant is a time-dependent Haar random quantum walk, obtained by generating new disorder realizations of the rotation matrices for every timestep (average of 200 random realizations of the walk is shown). The second classicalized variant was obtained with a single disorder realization, but with all coherence omitted from the quantum walk. Here, we replaced the unitary timestep operator $\hat U$ of Eq.  – more precisely, of Eq.  – with the corresponding stochastic operator, i.e., replaced all complex phase factors by 1, and the parameters $\pm\sin\theta_j(x,y)$ and $\cos\theta_j(x,y)$ by $\sin^2\theta_j(x,y)$ and $\cos^2\theta_j(x,y)$, respectively. For this second classicalized variant, the positive-valued function $\Psi(x,y)$ is interpreted as probability instead of probability amplitude. Here no disorder averaging was needed, as a single random realization already gave results with negligible statistical fluctuations. Figure \[fig:haar\_levelspacing\] shows the level spacing distribution, which matches very well the Wigner-Dyson distribution of Eq.. This indicates the presence of extended states – or, in other words, the absence of Anderson localization. The distribution was obtained for a single disorder realization, on a $160\times 160$ lattice, on the diagonal basis, with periodic boundary conditions. The quasienergy-averaged total transmission, shown in Fig. \[fig:Haar\_transmission\] is roughly unchanged as the system size is increased by a factor of 4, again indicating diffusive transmission. We used Eq.  for the transmission calculations, with parameters similar to those of Fig. \[fig:transmission\_phasetransition\], i.e., four different system sizes, with integration times $t_\text{max}$ matched for error $p_\text{leave}<0.01$, namely, $L_x=19, L_y=30, t_\text{max} = 1000$; $L_x=39, L_y=60, t_\text{max} = 2000$; $L_x=59, L_y=90, t_\text{max} = 3000$; $L_x=79, L_y=120, t_\text{max} = 4000$.

The diffusive spread of the Haar random two-dimensional quantum walk is due to disorder-induced topological criticality. Parameters $\alpha$ and $\beta$ have no impact on the topological phase, which is determined solely by $\theta_1$ and $\theta_2$. Although we are not sampling $\theta_1$ and $\theta_2$ uniformly, as in Ref. , we are using a distribution, using Eq. , shown in Fig. \[fig:haar\_theta\_probs\], which has the same probability density for $(\theta_1,\theta_2)=(\theta_A,\theta_B)$ as $(\theta_1,\theta_2)=(\theta_B,\theta_A)$, for any $\theta_A,\theta_B$. Thus, the qualitative argument used from Ref. [@edge2015localization], related to network models of topological phase transitions, still applies.

Conclusion, discussion {#sec:conclusion}
======================

We have studied Anderson localization and topological delocalization in two-dimensional split-step quantum walks with two internal (coin) states, and complete phase disorder. We first reviewed the known case of phase disorder on walks with real coin operators, which results in an anomalous Floquet-Anderson insulator. Here we used the numerical tools of wavefunction spread, level spacing distribution, and scaling of transmission. We simulated wavefunction spread in a diagonal basis, better adapted to the observed nonisotropic shape of the position distributions (elongated along the diagonal or antidiagonal direction in the $x-y$ basis). Moreover, we calculated the topological invariants for these disordered quantum walks using scattering theory, thus substantiating the topological explanation of the delocalization given in Ref. . We then have shown how the numerical tools and analytical arguments carry over to the case of completely disordered two-dimensional quantum walks, i.e., with $U(2)$ coin operators taken randomly and uniformly according to the Haar measure. We have found that this maximal disorder does not lead to Anderson localization, but results in a diffusive spread of the quantum walk. The absence of localization is explained by the observation that Haar random disorder tunes the system to a critical state, between different anomalous Floquet-Anderson phases.

Complete phase disorder was crucial in boosting the efficiency of the numerical tools we used, by ensuring that – statistically speaking – all quasienergies are equivalent. We could thus interpret wavefunction spread in a straightforward way, and use quasienergy averaging as a way of disorder averaging. We have found that this holds if complete disorder is in the phase $\phi$ or in the magnetic parameters $\alpha$ and $\beta$ – as could be expected from earlier work[@edge2015localization].

It would be interesting to consider whether the topological delocalization occurs in split-step quantum walks in higher dimensions. For the one-dimensional case, Haar random coins result in Anderson localization, which has even been rigorously proven[@ahlbrecht2011disordered]. There in the absence of any symmetries (Cartan class A), no topological invariant exists, and thus generic disorder destroys topological phases rather than drive the system to criticality. The same holds in every odd dimension[@hasan2010colloquium], and thus we expect Anderson localization for Haar random split-step quantum walks in any odd dimensions. However, in even dimensions we have the possibility of disorder-induced topological delocalization, and it would be interesting to find if this occurs, e.g., for four-dimensional split-step quantum walks with two internal states.

We also wonder how sensitive our conclusions are to the number of internal states of the quantum walk, which we set to two. Two internal states is the smallest number with which a discrete-time quantum walk can be constructed, but it requires the split-step construction for a quantum walks in two or more dimensions. Moving to higher number of internal states represents a challenge for the description of topological phases, because of the higher number of coin parameters. Moreover, with a larger internal coin space, it is quite possible that the localization length would be significantly larger, making it harder to observe numerical signatures of Anderson localization.

Our work also points to some open problems concerning a more complete picture of the localization-delocalization transition in disordered two-dimensional quantum walks. First, we lack the precise condition for Anderson localization of quantum walks for more general types of disorder – we have such a precise condition for the one-dimensional case[@rakovszky2015localization]). We have made some first steps in this direction for the numerical investigation of binary disorder, the results are in Appendix \[sec:binary\_disorder\]. Second, we don’t have an analytical understanding of why we obtain diffusion, rather than sub- or superdiffusion in the critical case. For both of these questions, mapping to network models of topological transitions[@delplace2017phase], or other theoretical tools[@woo2019quantum], could be used.

This work was partially supported by the Institute for Basic Science, Project Code IBS-R024-D1. A.M. thanks the Institute for Solid State Physics and Optics, Wigner Research Centre for hospitality and financial support. J.K.A. acknowledges support from the National Research, Development and Innovation Fund of Hungary within the Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017-00001), FK124723 and K124351.

Localization under binary disorder in rotation angles {#sec:binary_disorder}
=====================================================

In the main text we have observed delocalization of the system when the rotation angles are randomized, i.e., $(\theta_1, \theta_2)$ are position-dependent, picked from continuous random distributions. Thus, each $(x,y)$ has (typically) a different value of $\theta_1$ and $\theta_2$. Now we ask what happens with *binary disorder*, i.e., when the angles $(\theta_1, \theta_2)$ are position-dependent, but picked from the simplest discrete distribution: a probabilistic mixture (parameter $p_A$) of two sets of values, $$\begin{aligned}
  \theta_1(x,y) &= \theta_{1A};&
  \theta_2(x,y) &= \theta_{2A},\quad
  \text{prob.}~p_A~; \nonumber\\
  \theta_1(x,y) &= \theta_{1B};&
  \theta_2(x,y) &= \theta_{2B}, \quad
  \text{prob.}~(1-p_A), 
  \label{eq:binary_disorder_def}\end{aligned}$$ while the other coin parameters $\alpha_{1,2}$, $\beta_{1,2}$ and $\phi$ are uniformly random in range $[-\pi, \pi)$. Set $A$ corresponds to a winding number of $\nu_A=+1$.

   

For a numerical look into binary disorder, we varied the parameter sets $A$ and $B$ as follows. We fixed parameter set $A$, $$\begin{aligned}
  \label{eq:binary_theta_A}
 \theta_{1A} &= \frac{5\pi}{8};& \theta_{2A} &= - \frac{\pi}{8}.\end{aligned}$$ The winding number of set $A$, as per Eq.  – which holds for the case with or without phase disorder – is $+1$. We varied the $\theta_1$ and $\theta_2$ parameters of set $B$ together, such that $$\begin{aligned}
  \label{eq:binary_B_A}
\theta_{2B} - \theta_{2A} = \theta_{1A} - \theta_{1B}\end{aligned}$$ is always respected. Thus, when $\pi/4 <
\theta_{2B} < 3\pi/4$ , the winding number of set $B$ is $-1$, otherwise it is $+1$.

We show the numerical results on diffusive/localized behaviour as functions of the difference of the $\theta_2$ parameters and the probability $p_A$ of the set $A$ in Figs. \[fig:binary\_transmission\] and \[fig:binary\_rms\]. Figure \[fig:binary\_transmission\] shows increased transmission at the phase transition, accompanied by a change in the quantized transmission in the case with a cut (of type $B$), as in Fig. \[fig:topo\_trans\]; Fig. \[fig:binary\_rms\] shows the increased spreading of the walk in the critical case. When we have no values from set $A$, i.e., $p_A=0$, – bottom line of plots, – we see the expected signatures of the phase transition between anomalous Floquet-Anderson localized phases. Mixing in sites with coins in set $A$, i.e., increasing $p_A$, we find that the size of the phase with the winding number $-1$ shrinks. The critical value of $p_A$ needed for the winding number of $A$ to dominate the mixture is largest when $B$ is in a flat-band limit, $\theta_{1B}=0$ and $\theta_{2B} = \pi/2$. This critical $p_A$ is $1/2$ when $\theta_{1B}=\theta_{1A}-\pi/2$ and $\theta_{2B}=\theta_{2A}+\pi/2$; this is the case when the two sets $A$ and $B$ describe disorder-averaged dynamics that map unto each other under a mirror reflection. We remark that for these numerical results we used disorder in $\phi,
\alpha, \beta$, both in the system, the cuts, and the leads, and the $xy$ basis instead of the diagonal basis.
