Driven quantum dynamics: will it blend?

TL;DR

This paper demonstrates how driven many-body quantum systems can generate Haar-random unitaries, with convergence properties analyzed through Bethe-Ansatz and mean-field techniques, revealing potential for physical randomness generation.

Contribution

It introduces a method to produce Haar-uniform random operations in driven quantum systems using control and integrable models, and analyzes convergence times and spectral gaps.

Findings

Any fully controllable system converges to a unitary q-design over time.

The spectral gap of the driven spin chain's Liouvillean is independent of q.

Bethe-Ansatz shows the gap's independence from q, suggesting a universal property.

Abstract

Randomness is an essential tool in many disciplines of modern sciences, such as cryptography, black hole physics, random matrix theory and Monte Carlo sampling. In quantum systems, random operations can be obtained via random circuits thanks to so-called q-designs, and play a central role in the fast scrambling conjecture for black holes. Here we consider a more physically motivated way of generating random evolutions by exploiting the many-body dynamics of a quantum system driven with stochastic external pulses. We combine techniques from quantum control, open quantum systems and exactly solvable models (via the Bethe-Ansatz) to generate Haar-uniform random operations in driven many-body systems. We show that any fully controllable system converges to a unitary q-design in the long-time limit. Moreover, we study the convergence time of a driven spin chain by mapping its random evolution into a semigroup with an integrable Liouvillean and finding its gap. Remarkably, we find via Bethe-Ansatz techniques that the gap is independent of q. We use mean-field techniques to argue that this property may be typical for other controllable systems, although we explicitly construct counter-examples via symmetry breaking arguments to show that this is not always the case. Our findings open up new physical methods to transform classical randomness into quantum randomness, via a combination of quantum many-body dynamics and random driving.

Equations

$$\hat{H}(t)=H+g(t)V,$$

$$g(t)=\sum _{k=1}^K{A}_k\cos ({\omega }_{k}t+{\phi }_k),$$

$${\Vert {\mathbb{E}}_U\left[U^{\otimes q}(\cdot )U^{\otimes q}{}^{†}\right]-\int U^{\otimes q}(\cdot )U^{\otimes q}{}^{†}{\mu }_{\mathrm{Haar}}(dU)\Vert }_{\diamond }<ϵ,$$

$$e({\mu }_U,q)={\Vert {\mathbb{E}}_U\left[U^{\otimes q,q}\right]-\int U^{\otimes q,q}{\mu }_{\mathrm{Haar}}(dU)\Vert }_{\infty }<ϵ,$$

$${\mathbb{E}}_U\left[U^{\otimes q}\rho U^{\otimes q}{}^{†}\right]$$

$${\mathcal{L}}^q\rho$$

$${\mathcal{L}}_q=-\mathrm{i}{\mathring{\mathrm{H}}}^{\oplus \mathrm{q}}-\frac{\sigma }{2}({\mathring{\mathrm{V}}}^{\oplus \mathrm{q}}{)}^2,$$

$$\mathrm{Tr}[P_{\sigma }^{†}{\rho }_0]=\mathrm{Tr}[P_{\sigma }^{†}{\rho }_{\infty }]=\sum _{\pi \in S_q}{\rho }_{\pi }\mathrm{Tr}[P_{\sigma }^{†}{P}_{\pi }],$$

$${\rho }_{\infty }=\underset{t\to \infty }{\lim }{e}^{t{\mathcal{L}}_q}{\rho }_0=\sum _{\pi ,\sigma }{\left(M^{-1}\right)}_{\pi \sigma }\mathrm{Tr}\left[P_{\sigma }^{†}{\rho }_0\right]P_{\pi },$$

$$\underset{t\to \infty }{\lim }{\mathrm{e}}^{\mathrm{t}{\mathcal{L}}^{\mathrm{q}}}\rho =\int \mathrm{dU}{\mathrm{U}}^{\otimes \mathrm{q}}\rho {\mathrm{U}}^{\otimes \mathrm{q}†},$$

$${\mathcal{L}}_q$$

$${\mathcal{L}}_k^{\mathrm{loc}}$$

$${\mathcal{L}}_1=\sum _i{\lambda }_i{\Pi }_{(i)},$$

$${\Pi }_j^{(i)}={\Pi }_{(0)}^{\otimes (j-1)}\otimes {\Pi }_{(i)}\otimes {\Pi }_{(0)}^{(q-j)},$$

$$||e^{t{\mathcal{L}}^q}-{\mathcal{U}}_{\infty }|{|}_{1\to 1}\le d^{2q}||e^{t{\mathcal{L}}^q}-{\mathcal{U}}_{\infty }|{|}_{2\to 2},$$

$$||e^{t\mathcal{L}}-{\mathcal{U}}_{\infty }|{|}_{2\to 2}$$

$$||e^{t{\mathcal{L}}^q}-{\mathcal{U}}_{\infty }|{|}_{2\to 2}\le \sum _{j:{\lambda }_j\ne 0}{e}^{t\mathrm{Re}[{\lambda }_j]}\sqrt{||R_j||||L_j||},$$

$$||e^{t{\mathcal{L}}^q}-{\mathcal{U}}_{\infty }|{|}_{2\to 2}\le e^{-t{\lambda }^{\ast }}{d}^{2q}\kappa (\mathcal{S}{)}^2,$$

$$\underset{\rho }{\sup }||e^{t{\mathcal{L}}^q}(\rho )-{\mathcal{U}}_{\infty }(\rho )||\le e^{-t{\lambda }^{\ast }}{d}^{4q}\kappa (\mathcal{S}{)}^2.$$

$$||e^{t{\mathcal{L}}^q}-{\mathcal{U}}_{\infty }|{|}_{1\to 1}$$

$$\ge$$

$$||e^{t{\mathcal{L}}^q}-{\mathcal{U}}_{\infty }|{|}_{1\to 1}\ge e^{O(q)}{e}^{-t{\lambda }^{\ast }}.$$

$${\mathcal{L}}_q$$

$$-\frac{\sigma }{2}\sum _{\alpha \beta \gamma \delta }{V}_{\alpha \beta }{V}_{\gamma \delta }(B_{\alpha \beta }^{\uparrow }-B_{\beta \alpha }^{\downarrow })(B_{\gamma \delta }^{\uparrow }-B_{\delta \gamma }^{\downarrow }).$$

$$\sum _{\alpha u}{a}_{\alpha u}^{†}{a}_{\alpha u}=\sum _{\alpha }{B}_{\alpha \alpha }=q\text{openone}.$$

$${\mathcal{L}}_q$$

$$-\frac{\sigma }{2}\sum _{\alpha \beta uv}{V}_{\alpha \alpha }{V}_{\beta \beta }(n_{\alpha u\uparrow }-n_{\alpha u\downarrow })(n_{\beta v\uparrow }-n_{\beta v\downarrow }),$$

$${\mathcal{L}}_q=$$

$${\mathcal{L}}_q$$

$$Q{H}^T+HQ=Q{V}^T+VQ=0.$$

$$|{\psi }_Q⟩=\sum _{\tilde{\alpha }\tilde{\beta }}(Q\otimes Q^{\ast }{)}_{\tilde{\alpha }\tilde{\beta }}{a}_{\tilde{\alpha }}^{†}{a}_{\tilde{\beta }}^{†}|0⟩.$$

Full text

I Introduction

Randomness generating quantum operations play a central role in our understanding of very various physical phenomena Guhr et al. (1998). Recently, with the development of quantum information processing, random operations have found new applications, not only as a theoretical tool, but also in practical protocols. Indeed, they are used in quantum cryptography Hayden et al. (2004), quantum process tomography Bendersky et al. (2008), fidelity estimation Dankert et al. (2009), quantum communication and entanglement sharing Harrow et al. (2004); Abeyesinghe et al. (2009); Hastings (2009), quantum data-hiding DiVincenzo et al. (2002); Hayden et al. (2004); Piani et al. (2014) and entanglement generation Oliveira et al. (2007); Žnidarič (2008); Hamma et al. (2012); Zanardi (2014). Because of their crucial importance, several procedures have been developed to generate either truly random or pseudo-random operations via random quantum circuits Brandao et al. (2016); Brown and Viola (2010); Emerson et al. (2003); Dankert et al. (2009); Gross et al. (2007); Harrow and Low (2009); Turner and Markham (2016); Alexander et al. (2016). However, from the physical point of view, these protocols often have a complexity comparable with universal quantum computation, being based on the application of a sufficiently large set of quantum gates. Here, on the other hand, we consider a more physically inspired approach, based on quantum control, where the quantum system is controlled by random classical pulses.

Quantum control is an established research field at the overlap of control theory and quantum mechanics. Essentially it provides a framework to steer a quantum system through Hilbert space by applying time-dependent fields. Controllability is a powerful algebraic tool to fully characterise when any possible unitary evolution in the system’s Hilbert space can be obtained from the Schrödinger equation with a suitable choice of time-dependent fields. The central question of this paper is what happens when we apply random fields to a controllable system. We will show, under some conditions, that after a suitably long *mixing * time the corresponding random unitary evolutions of the system converge to a uniformly random set, as measured by the Haar measure. Therefore, one of the central result of this paper is that driving a controllable quantum systems with stochastic control pulses offers a natural approach to generate random unitary operations with physical processes.

Within this picture, the estimation of the mixing time is the crucial theoretical aspect. We use several tools from the theory of open quantum systems and many-body physics, such as low-energy effective Liouvilleans, mean-field techniques and the Bethe-Ansatz, to find an accurate estimation of the mixing time in several situations. In particular, we focus on a one-dimensional system with edge control due to the availability of analytical tools, as well as the intuitive interpretation available in such a system with Lieb-Robinson bounds and spin waves. This particular case is also motivated by the current experimental capabilities in integrated photonic circuits Perez-Leija et al. (2013); Pitsios et al. (2016), where different stochastic control pulses can be simulated by changing the spatial extent of the waveguides via electrically tuned on-chip heaters Carolan et al. (2015). In those systems a major recent result has been the experimental measurement of boson sampling Broome et al. (2013); Spring et al. (2013); Crespi et al. (2013), a problem which is believed to be hard to simulate classically. Random unitary operations and higher dimensional systems are required in boson sampling to have a convincing demonstration of quantum computational supremacy Aaronson and Arkhipov (2011). Pseudo-random operations in those experiments are currently obtained via a finely tuned network of several beam splitters and phase-shifters. The different approach presented here is based on the simpler implementation of noisy quantum walks and, therefore, can offer an advantage to perform boson sampling experiments on larger scale.

A further motivation for this paper comes from quantum control itself. The algebraic tools developed in quantum control are typically not able to provide an estimation of the control time needed to reach a given target operation. In view of practical applications, this is a big handicap, because noise will always limit the total time available to an experimenter. It is therefore of interest to find estimates of such times. The analytical expressions for the mixing time obtained in this paper provide also an easily computable upper bound for the control time. Indeed, by definition, after the mixing time the system has already explored all possible unitary evolutions with stochastic control pulses. This implies that, apart from measure zero sets, at this time any evolution is achievable with a suitable choice of the control field.

Finally, another motivation for the present work is for the problem of fast scrambling of quantum information. The problem was first identified in the setting of black hole physics Hayden and Preskill (2007); Sekino and Susskind (2008), where it was conjectured that black holes start evaporating information when most localized microscopic degrees of freedom become inaccessible without measuring a constant fraction of the whole system. Unfortunately, identifying mechanisms for fast scrambling has been challenging, and providing tools to rigorously analyze scrambling times even more so. Moreover, explicit constructions of fast scramblers Lashkari et al. (2013) are not directly inspired by physical models. Here we describe a physically motivated process that could lead to new insights in the design and analysis of fast scrambling models.

The paper is organized as follows: in section II we show how to obtain Haar-uniform unitary evolutions (i.e. a unitary design) via quantum control techniques. We will focus on $q$-design, not only for its applications in quantum information, but also to quantify the distance with the target uniform distribution. We will consider Markovian stochastic control pulses and introduce some general techniques for the estimation of the mixing time. In section III we map the problem of unitary design to a general many-body problem, studying its mean-field solution and discussing the limitations of the latter approach via symmetry breaking arguments. In section IV we focus on a specific one-dimensional model controlled at one of its boundaries. We show that this model in certain limits can be mapped to an exactly solvable model and we study its analytic solution via Bethe-Ansatz techniques. A central result of this section is that the mixing time for this particular model is independent of the number of copies $q$. Intuitively the $q$-independence implies that pseudo-random unitaries obtained with random control pulses approximate all the moments of the Haar distribution with the same accuracy. These predictions are then corroborated with numerical simulations. In Section V we show other applications for boson sampling, the decay of correlations in spin chains, and for the estimation of the control time. Conclusions and perspectives are written in section VI.

II Unitary designs via quantum control

Physical quantum systems are modeled via a Hamiltonian operator $H$, which describes the interactions between the components of the system. When external control is applied to the system, its evolution is represented by a time-dependent Hamiltonian

[TABLE]

where $g(t)$ is an external control pulse and $V$ is an operator. If $d$ is the dimension of the Hilbert space, then $H$ and $V$ are $d\times d$ Hermitian matrices while $g(t)$ is a scalar function depending on time $t$. For multiple pulses $\hat{H}(t)=H+\sum_{i}g_{i}(t)V_{i}$. After some time $T$, the combined action of the natural interactions and the external pulses is a unitary operation $U=\mathcal{T}\exp\left({-}\rm i\int_{0}^{T}\hat{H}(s)\,ds\right)$, where $\mathcal{T}$ represents the time order operator. In general, the amount of different unitary operations $U$ that can be obtained from the dynamics of the system is limited. However, if the system is fully controllable, then any operation can be obtained with a suitable engineering of the control pulse. In other terms, given any $U\in{\rm SU}(d)$ it is possible to find a control profile $g(t)$ such that $U=\mathcal{T}\exp\left({-}\rm i\int_{0}^{T}\hat{H}(s)\,ds\right)$ where the control time $T$ depends on the target unitary $U$. There are many powerful theorems to test controllability. In general a system described by the Hamiltonian as in Eq. (1) is controllable D’Alessandro (2007) if $H,V$ and their nested commutators $[A,[B,[C,\dots]]]$ (where $\{A,B,C,\dots\}{\in}\{H,V\}$) generate the Lie algebra of SU(d). Although the algebraic conditions for controllability are well known, it is still an open problem in quantum control to estimate the control time $T$, given also the knowledge of the target gate $U$ and the operators $H$ and $V$. For fully controllable systems there exists a minimal control time, generally unknown, such that all target gates can be obtained exactly at that time Jurdjevic and Sussmann (1972). For small dimensional systems, analytic bounds of such universal control time may be found in terms of quantum speed limits or Cartan decompositions of spin systems. In high dimensional system, such tools become intractable. If the system is drift-free ($H_{0}=0$), control times are trivial or only determined by energy bounds on the time-dependent fields. We are instead interested in systems where the controls need to work together with a drift to achieve full control (so-called weak controllability). In such a case, the timescale is bounded by the dynamics of the drift and provides insights into the many-body physics triggered by it.

We now consider the control pulse as a stochastic process, namely where a certain profile $g(t)$ can be applied to the system with a probability $p_{g(t)}$, and study the distribution of the resulting unitary operations. Such a random pulse can be obtained, for example, by considering the Fourier expansion of the control signal

[TABLE]

where the amplitudes $A_{k}$, the phases $\varphi_{k}$, and possibly even the frequencies $\omega_{k}$ are random variables. We use the notation $\mathbb{E}[\cdot]$ to denote the average over those random variables. Repeating the experiment with many random signals one obtains a distribution of unitary matrices, where each matrix $U$ is obtained with probability $p_{U}$. Random unitary operations play a central part in many quantum information protocols. A pivotal role in many applications is played by the uniform distribution, called also Haar distribution, which is invariant under the action of the unitary group itself. In the following sections we study when, and how rapidly, the distribution $p_{U}$ converges to the Haar-uniform distribution.

II.1 Comparing random evolutions: unitary $q$-design

Obtaining truly uniform random unitaries is a very hard task, and normally one observes pseudo-uniform distributions which approximate the uniform (Haar) measure up to some errors. Pseudo-uniform distributions can be obtained with random quantum circuits Brandao et al. (2016); Brown and Viola (2010); Emerson et al. (2003); Dankert et al. (2009); Gross et al. (2007); Harrow and Low (2009), but these circuits typically require many different gates that make demanding the implementation in physical systems. Recently, alternative protocols based on physically inspired time-dependent Hamiltonians have been proposed Nakata et al. (2017); Onorati et al. (2016). Nonetheless, these approaches still require that all the interactions inside the system should change in time, an assumption that currently is beyond reach in many experimental platforms. Here, on the other hand, we focus on a general scheme which occurs in most quantum systems, namely when the natural and time-independent interaction $H$ experienced by the system is paired with an external control, as in Eq. (1).

There are many ways of comparing the distance between two quantum processes. When dealing with randomness generating processes, it is often convenient and relevant to work with approximate $q$-designs Brandao et al. (2012). A unitary $q$-design is a distribution of unitaries, possibly discrete, that gives the same expectations of the Haar distribution for polynomial functions of degree at most $q$ (see e.g. Low (2010)). It is often inaccessible experimentally to distinguish between truely random processes and approximate $q$-designs. Formally, approximate $q$-designs are defined by the requirement that

[TABLE]

for suitably small $\epsilon$, where $\|\cdot\|_{\diamond}$ refers to the diamond norm, $\mathbb{E}_{U}$ denotes an average over some given distribution of unitaries $\mu_{U}$ and $\mu_{\rm Haar}(dU)$ is the Haar measure. This is the most stringent distinguishability measure between quantum processes, and guarantees that no single (global) measurement on the system and a possible ancilla can distinguish between the two processes with probability larger than $\epsilon$. A related notion Harrow and Low (2009) is that of quantum expanders, which are defined by

[TABLE]

where $X^{\otimes q,q}=X^{\otimes q}\otimes(X^{\otimes q})^{*}$. Eq. (4) can be regarded as the vertorised version of Eq. (3): given an operator $X=\sum_{ij}X_{ij}|{i}\rangle\langle{j}|$, its vectorized form is $|X\rangle\rangle=\sum_{ij}X_{ij}|{ij}\rangle$. However it is striclty weaker, and the separation between the two bounds can be exponential in the system size. However, Eq. (4) is often much easier to work with in practice Harrow and Low (2009). It follows from the definition that $|AX\rangle\rangle=A\otimes\openone|X\rangle\rangle$ and $|XA\rangle\rangle=\openone\otimes A^{T}|X\rangle\rangle$. Therefore, $X^{\otimes q,q}$ is the vectorization of the superoperator $\rho\mapsto X^{\otimes q}\rho X^{\otimes q}{}^{\dagger}$. Quantum expanders and $q$-design compare probability distributions of unitary matrices by comparing the “moments” of the distribution, namely random processes that depend polynomially on the random variable. Two close distributions of unitary matrices have similar moments, as shown in Brandao et al. (2016), $e(\mu,q)\leq 2q\,\mathcal{W}(\mu,\mu_{\rm Haar})$, for all measures $\mu$, being $\mathcal{W}$ the Wasserstein distance Oliveira (2009) $\mathcal{W}(\mu_{1},\mu_{2})=\sup_{f}\left|\int f(U)\,[\mu_{1}(dU)-\mu_{2}(dU)]~{}\right|$, where $f$ is a 1-Lipschitz function, and $U$ is a unitary matrix. The Wasserstein distance is a measure between classical probability distributions, and hence one can use a number of classical Markov chain mixing tricks to bound it. However, we will use not be using it, as we instead use tools from condensed matter physics to bound the mixing time.

In the quantum control setting, $\mathbb{E}_{U}$ in Eqs. (3) and (4) is the average over many unitary operations obtained after the application of random pulses up to a certain time $T$. Therefore

[TABLE]

To simplify the theoretical description of this problem we make two assumptions. (i) We assume that the stochastic process $g(t)$ is Gaussian. This is a reasonable approximation in many-cases and can be obtained e.g. via Eq. (2) when $K\gg 1$, in view of the central limit theorem. (ii) We assume also that $g(t)$ is harmonic, namely that $\mathbb{E}[g(t+s)g(t)]=c(s)$ is independent of $t$. Moreover, without loss of generality, the harmonic process can be chosen such that $\mathbb{E}[g(t)]=0$. In view of these assumptions, exploiting the results of Ishizaki and Fleming (2009); Banchi et al. (2013), in appendix A we find a closed form expression for Eq.(5). That expression can be drastically simplified if we assume that the correlation time is finite and there exists a suitably large $T$ such that $Tc(Ts)\simeq\frac{\sigma}{2}\delta(s)$ where $\delta$ is the Dirac delta function and $\sigma$ is a constant. In the long-time limit, $t>T\gg\|H\|,\|V\|$, one finds then that

[TABLE]

where

[TABLE]

and $X^{\oplus q}{=}X{\oplus}X{\oplus}\dots$, being $\oplus$ the Kronecker sum $X{\oplus}Y{=}X{\otimes}\openone{+}\openone{\otimes}Y$. Therefore, with these three approximations, the long-time dynamics of the stochastic process is Markovian and described by the above Lindblad equation Lindblad (1976); Gorini and Kossakowski (1976), where the operator $\mathcal{L}^{q}$ is called Liouvillean. Similarly to what happens with the replica trick in statistical physics Mézard et al. (1987), the average over the noise effectively couples the initially uncoupled copies. Sometimes we will use the more convenient vectorised form of the above equation

[TABLE]

where $\mathring{X}=X\otimes\openone-\openone\otimes X^{T}$ is the vectorization of the commutator $[X,\cdot]$. If $t\to\infty$ then $\mathbb{E}_{U}\left[U^{\otimes q}\rho U^{\otimes q}{}^{\dagger}\right]$ converges to one of the steady states of the Liouvillean $\mathcal{L}^{q}$.

In the following section we prove that the steady state manifold of $\mathcal{L}^{q}$ coincides with the state space after averaging over the Haar measure, namely that all the moments of the random unitary evolution converge to the averages over the uniform distribution for $t\to\infty$. Moreover, we will study the mixing time via the gap of the Liouvillean and show that, in several cases, the latter is independent on $q$. Physically this is important, because it implies that all the moments converge (in 2-norm) at the same time, as given by the inverse of the Liouvillean gap, and that, accordingly, we can use the latter to estimate the mixing time of the random unitary evolutions.

II.2 Steady state of the Liouvillean evolution

We start by describing the steady state of $\mathcal{L}^{q}$. In general, the dimensionality of the steady state set is in one-to-one relation with the conserved quantities of the Lindbladian evolution Albert and Jiang (2014). Given an orthonormal basis $\{M_{\mu}\}$ of the steady state space, equipped with the standard Hilbert-Schmidt product, there exists a dual operator set $\{J_{\mu}\}$ such that $\mathcal{L}^{q\dagger}J_{\mu}=0$, where $\mathcal{L}^{q\dagger}$ is the Liouvillean operator (7) after the substitution $H\to-H$. The latter substitution does not change the dynamical algebra, so algebraic considerations based on controllability hold also for $\mathcal{L}^{q\dagger}$. From the conserved quantities $J_{\mu}$ and their dual operators $M_{\mu}$ one finds the steady state as $\rho_{\infty}=\sum_{\mu}M_{\mu}\operatorname*{Tr}(J_{\mu}\rho_{0})$ where $\rho_{0}$ is the initial state Albert and Jiang (2014). Since the system is controllable, repeated commutators of $H^{\oplus q}$ and $V^{\oplus q}$ give rise to the algebra su$(d)^{\oplus q}$. Therefore, because of the Schur-Weyl duality Goodman and Wallach (2000), the only operators that commute with both $H^{\oplus q}$ and $V^{\oplus q}$, and more generally with Eq. (5), are index permutation operators. Let $S_{q}$ be the group of permutations of the set $1,\dots,q$ and let $P_{\sigma}$, $\sigma\in S_{q}$ be the operator which permutes the index of the tensor copy $\mathcal{H}^{\otimes q}$, namely the operator that maps $\psi_{i_{1},i_{2},\cdots,i_{n}}$ to $\psi_{\sigma(i_{1}),\sigma(i_{2}),\cdots,\sigma(i_{n})}$ for each set of indices $i_{j}$. It is simple to show that $P_{\pi}P_{\sigma}{=}P_{\pi\sigma}$ and that these operators form a unitary representation of the permutation group $S_{q}$. The index permutation operators are the only conserved quantities of the Liouvillean, $\mathcal{L}^{q}(P_{\rho})=\mathcal{L}^{q\dagger}(P_{\rho})=0$, so $\rho_{\infty}=\sum_{\sigma}\rho_{\sigma}P_{\sigma}$. However, since the operators $P_{\sigma}$ are not orthonormal, one has

[TABLE]

where in the first equality holds because $P_{\sigma}$ is a conserved quantity. By inverting the above equation we find that

[TABLE]

where $M_{\sigma\pi}=\operatorname*{Tr}[P_{\sigma}^{\dagger}P_{\pi}]$. It has been shown in Ref. Brandão et al. (2012) that $M_{\sigma\pi}=d^{l(\sigma^{-1}\pi)}$ where $l(\sigma)$ is the number of cycles in the cycle decomposition of $\sigma$. The dimensionality of the steady state manifold is then given by the matrix rank of $M$. One finds that the steady state degeneracy is $\sim e^{\mathcal{O}(q)}$. The right-hand side of (10) is exactly equal to the integration over the Haar measure (see e.g. Proposition 3 in Brandão et al. (2012)). Therefore, we have shown that

[TABLE]

namely that the infinite time-evolution of the system under the Liouvillean (7) is equivalent to an integration over the Haar measure.

In summary, we have shown that by driving a controllable system with random control pulses Eq. (2), where the stochastic process is Gaussian, harmonic and has a finite correlation time, then the resulting average evolution of the quantum system converges for $t\to\infty$ to a uniform integration over the Haar measure.

II.3 Construction of excited states

Certain excited states of the Liouvillean (8) can be built up directly from the excitations of the individual quantum systems. It is convenient to separate $\mathcal{L}_{q}$ from Eq. (8) into local terms ${\mathcal{L}}_{k}^{\rm loc}$ acting only on the $k$-th copy, and a non-local interaction. Indeed,

[TABLE]

where $\mathring{H}_{k}$, $\mathring{V}_{k}$, and accordingly $\mathcal{L}^{\rm loc}_{k}$, act only on the $k$-th copy. Therefore each $\mathcal{L}^{\rm loc}_{k}$ for different $k$ is equivalent to a single-copy Liouvillean $\mathcal{L}_{1}$. We assume that the operator $\mathcal{L}_{1}$ is diagonalizable (with right and left eigenvectors) and call

[TABLE]

its eivenvalue decomposition, where the eigenvalues $\lambda_{j}$ are ordered with decreasing real part (starting from zero) and $\Pi_{j}$ are the corresponding eigenprojections. The operators

[TABLE]

are then eigenprojections of $\mathcal{L}^{q}$, with eigenvalue $\lambda_{i}$. To show this, we note indeed that $\Pi^{(i)}_{j}$ is proportional to the vectorization of the identity operator in each copy, aside from the $j$-th one, since $\Pi_{(0)}$ is the projection onto the steady state and, accordingly, $\Pi_{(0)}(X)=\rho_{\infty}\operatorname*{Tr}[X]$, which is proportional to the identity operator. Therefore, $\mathring{V}_{l}\,\Pi^{(i)}_{j}=0$ (because $[V_{l},\Pi^{(i)}_{j}(X)]=0$ for all $X$), as long as $l\neq j$. On the other hand, for $l=j$, it is $\mathring{V}_{k}\,\mathring{V}_{l}\,\Pi^{(i)}_{j}=0$, since by construction $k\neq j$. This shows that (15) is a projector on the eigenspace of $\mathcal{L}_{q}$ with eigenvalue $\lambda_{i}$. Moreover, from the operators (15) one can also construct the eigenstates of $\mathcal{L}_{q}$ that act on the irreducible representations of the symmetric group – indeed since the permutation operators $P_{\sigma}$ commute with the Liouvillean, then $P_{\sigma}(\Pi^{(i)}_{j})P_{\sigma}^{\dagger}$ is an eigenprojection of $\mathcal{L}_{q}$ for all $\sigma$.

In summary, the eigenstate of $\mathcal{L}_{1}$ with the lowest gap can be used to construct some exact eigenstates of $\mathcal{L}_{q}$, although it remains to be shown that they have the smallest gap. These eigenvalues have degeneracy at least as large as the ground state degeneracy, since $P_{\rho}\Pi_{j}^{(i)}$ is also an eigenvector with eigenvalue $\lambda_{j}$ of $\mathcal{L}^{q}$.

II.4 Convergence time

Given the results of the previous section, we want to know how rapidly the semigroup converges to the uniform distribution Eq. (11). In Appendix B, we provide a brief introduction to the convergence theory of dynamical semigroups, and argue that when the generator is not reversible (detailed balance), the convergence is governed by the singular value gap of the channels rather than the spectral gap of the generator. In general we want to bound the trace norm, but it will be more convenient to analyze the $2\rightarrow 2$ norm:

[TABLE]

where $\mathcal{U}_{\infty}=\lim_{t\to\infty}e^{t\mathcal{L}^{q}}$ and $d$ is the dimension of the local Hilbert space. Let $s_{j}(t)$ be the singular values of $e^{t\mathcal{L}}$, ordered from largest to smallest. The largest has magnitude one. Then the singular values of $(e^{t\mathcal{L}}-\mathcal{U}_{\infty})$ are strictly smaller than one, and

[TABLE]

If the Liouvillian were reversible, then the singular values $s_{j}(t)$ would be given by $e^{t\lambda_{j}}$, where $\lambda_{j}$ are the eigenvalues of $\mathcal{L}$. Unfortunately the semigroups that we will be working with are not Hermitian. Nonetheless, from Eq. (64), we find that the $2\rightarrow 2$ norm can be bounded in terms of the eigenvalues and eigenvectors of $\mathcal{L}^{q}$ as

[TABLE]

where $\lambda_{j}$ are the eigenvalues of $\mathcal{L}^{q}$, and $R_{j},L_{j}$ are its right and left eigenvectors, satisfying ${\rm tr}[L^{\dagger}_{j}R_{k}]=\delta_{jk}$.

In general it is very difficult to bound Eq. (18), since the norms of the eigenvectors can be very large, and it is often difficult to get good bounds on the spectrum. Nonetheless, in Appendices B, C and D, we study both the weak $(\sigma=\epsilon\rightarrow 0)$ and strong $(\sigma=\epsilon^{-1}\rightarrow\infty)$ coupling limits, and show the following properties: (i) the spectral gap is $O(\epsilon)$, both in the strong and weak coupling limits – for strong driving, the decrease of the gap for larger $\sigma$ is consistent with the general occurrence in open systems Zanardi et al. (2016); (ii) the eigenvectors satisfy $|R_{j}\rangle=\mathcal{S}|\Phi_{j}\rangle$ and $|L_{j}\rangle=\mathcal{S}^{{\dagger},-1}|\Phi_{j}\rangle$, for some invertible matrix $\mathcal{S}$ and an orthonormal basis $|\Phi_{j}\rangle$. The condition number of $\mathcal{S}$ is $\kappa(\mathcal{S})\equiv||\mathcal{S}||~{}||\mathcal{S}^{-1}||$ and satisfies $\kappa(\mathcal{S})=O(1+\epsilon)$. Moreover, in sections III.2 and IV we will discuss some cases where the Liouvillean gap is independent on $q$. Models whose mixing time is independent on $q$ have been obtained also in Onorati et al. (2016), at the expense of more stringent requirements on the fluctuating terms of the Hamiltonian.

We then get that

[TABLE]

where $\lambda^{*}$ is the eigenvalue with the smallest non-zero real part and $\kappa(\mathcal{S})=\mathcal{O}(1+\epsilon)$. In terms of the trace norm, we then get that

[TABLE]

In the weak or strong coupling limits, the condition number will be of order one yielding a mixing time of $T^{*}\sim 4q\log(d)/\lambda^{*}$. We lost a lot in two steps of the bound, both times involving a term of order $d^{2q}$. In certain cases, this is overly pessimistic. For instances, for a tensor product of $n$ semigroups, the mixing time is $T^{*}\sim\log(n)T^{*}_{1}$, where $T^{*}_{1}$ is the mixing time of a single subsystem Kastoryano et al. (2012). We might ask whether the mixing time of Eq. (7) is also of the order $T^{*}\sim\log(q)T^{*}_{1}$, with $T^{*}_{1}=O(1/\lambda^{*})$?

We can see that this is not the case from the following argument:

[TABLE]

since the lower bound is saturated when $\mathcal{S}=\openone$, and we have isolated the subspace with eigenvalue $\lambda^{*}$. Now, in Section II.3 we have argued that if the gap of $\mathcal{L}^{q}$ is the same as the gap of $\mathcal{L}^{1}$, then we can construct the eigenvectors with minimal non-zero eigenvalue of $\mathcal{L}^{q}$ from those of $\mathcal{L}^{1}$. In particular, the size of this subspace is at least as large as the size of the ground state subspace. But we know that the ground state subspace has dimension $d_{0}\geq e^{O(q)}$. Hence the first excited subspace does as well. Then,

[TABLE]

Thus the mixing time is at least $T^{*}\sim O(q/\lambda^{*})$, even in the weak coupling limit.

Finally, we comment on the distinction between the singular value gap of $e^{t\mathcal{L}}$ and the eigenvalue gap of $\mathcal{L}$. We know that as $t\rightarrow\infty$, the singular value gap $s^{*}(t)$, namely the largest singular value $s_{j}(t)\neq 1$, converges to $e^{t\lambda^{*}}$, however it is not clear how rapidly this occurs. This will be discussed in the numerical studies of Sec. IV where we will show that, both in the strong and weak coupling limits, the difference between the spectral gap and the singular value gap vanishes on a time scale much smaller than $1/\lambda^{*}$.

III Many-body theory of unitary design

In the previous section we have argued that bounding the spectral gap of the dynamical semigroup is in many relevant cases sufficient to get good estimates on the mixing time of the process. Here we will study such a gap by introducing a general mapping from a control Liouvillian to a non-Hermitian many-body Hamiltonian, and then study its mean field solution. The mean field approach has been already successfully applied Brown and Viola (2010) to estimate the convergence time of permutationally invariant random quantum circuits, where at each step a gate from a universal set is applied to a random pair of qubits. Moreover, in Sec. IV we will analyze an integrable example via Bethe-Ansatz techniques, from whose solution it appears that the eigenstates with smallest gap are constructed from the steady states by changing the internal state of a single unpaired particle. This fact shares several similarities with what happens in bosonic condensates, and in particular with their mean field solution Blaizot and Ripka (1986). Motivated by these two examples, it is natural to apply the mean field analysis to generic Hamiltonian evolutions with random pulses. However, although the predictions of the mean field solution are consistent with several numerical simulations, we will clarify that this approach cannot be general by constructing explicit counterexamples via symmetry breaking arguments.

III.1 Mapping to a non-Hermitian many body Hamiltonian

A powerful method for estimating the spectral gap of the Liouvillean is to map Eq. (8) to a many-body problem, and then use powerful techniques developed in condensed matter systems to obtain the spectrum. In order to find this mapping we introduce a basis $b_{\alpha\beta}=|{\alpha}\rangle\langle{\beta}|$, $\alpha,\beta=1,\dots,d$ and call $B_{\alpha\beta}=b_{\alpha\beta}^{\oplus q}$. These operators satisfy the SU(d) commutation relation, $[B_{\alpha\beta},B_{\gamma\delta}]=B_{\alpha\delta}\delta_{\beta\gamma}-\delta_{\alpha\delta}B_{\beta\gamma}$ and therefore define a reducible representation of SU(d). Moreover, $\mathring{X}^{\oplus q}=X^{\oplus q}\otimes\openone-\openone\otimes X^{\oplus q}{}^{T}=\sum_{\alpha\beta}\left(X_{\alpha\beta}B^{\uparrow}_{\alpha\beta}-(X^{T})_{\alpha\beta}B^{\downarrow}_{\alpha\beta}\right)$ where we set $B^{\uparrow}_{\alpha\beta}=B_{\alpha\beta}\otimes\openone$ and $B^{\downarrow}_{\alpha\beta}=\openone\otimes B_{\alpha\beta}$. Hence, the Liouvillean can be written as

[TABLE]

The form (24) is a convenient starting point because it depends only on the original $d\times d$ operators introduced in (1), while the complicated action into the $q$-copy Hilbert space is transferred into the basis operators $B$.

The operators $B$ form a reducible representation of SU(d) and can be decomposed in terms of irreducible operators that act on different invariant subspaces of the original $(\mathbb{C}^{d})^{\otimes q}$ Hilbert space. Indeed, because of the Schur-Weyl duality, every irreducible representation of $(\mathbb{C}^{d})^{\otimes q}$ is decomposed as $(\mathbb{C}^{d})^{\otimes q}=\otimes_{\lambda}\mathcal{P}^{\lambda}\otimes U^{\lambda}$ where $\mathcal{P}^{\lambda}$ is an irreducible representation of the symmetric group $S_{q}$ and $\mathcal{U}^{\lambda}$ an irreducible representation of SU(d). A convenient expression for the fully-symmetric and fully-anti-symmetric subspaces is given by Rowe et al. (2012) $B_{\alpha\beta}=a_{\alpha}^{\dagger}a_{\beta}$, where $a_{\alpha}$ and $a_{\alpha}^{\dagger}$ are either bosonic or fermionic creation and annihilation operators. Moreover, even a generic (though reducible) representation can be constructed from either bosonic or fermionic annihilation operators by adding an extra index and writing $B_{\alpha\beta}=\sum_{u}a_{\alpha u}^{\dagger}a_{\beta u}$. From the definition of $B$ one realizes that in this generic representation there are exactly $q$ particles since

[TABLE]

For convenience, we also perform the calculation in the basis where $V$ is diagonal. Therefore, Eq.(24) becomes

[TABLE]

where $n_{x}=a_{x}^{\dagger}a_{x}$. Thanks to this general representation, the many-body Liouvillean has been mapped to a many-particle Hubbard-like problem (27) where the hopping part is anti-Hermitian. The original dependence on $q$ is mapped to the number of particles, namely to the constraint (25) that there are exactly $q$ particles in the “spin-up” and “spin-down” states, $\sum_{\alpha u}n_{\alpha u\uparrow}=\sum_{\alpha u}n_{\alpha u\downarrow}=q\,\openone$.

III.2 Mean-field approach

We consider here the decomposition (12) where each $\mathcal{L}^{\rm loc}_{k}$ for different $k$ is equivalent to a single-copy Liouvillean $\mathcal{L}_{1}$. From the above decomposition it is clear that if the gap of $\mathcal{L}^{\rm loc}=\sum_{k}\mathcal{L}^{\rm loc}_{k}$ equals the gap of $\mathcal{L}_{q}$ then the Liouvillean gap $\lambda^{*}$ is independent on $q$.

Extending the treatment of Section III.1, we define a local basis of operators $\tilde{B}_{\tilde{\alpha}\tilde{\beta}}=B^{\uparrow}_{\alpha_{\uparrow}\beta_{\uparrow}}\delta_{\alpha_{\downarrow}\beta_{\downarrow}}+\delta_{\alpha_{\uparrow}\beta_{\uparrow}}B^{\downarrow}_{\alpha_{\downarrow}\beta_{\downarrow}}$ where $\tilde{\alpha}=(\alpha_{\uparrow},\alpha_{\downarrow})$, and similarly for $\tilde{\beta}$, are multi-indices running from 1 to $d^{2}$. Therefore we can write the decomposition Eq. (12) as

[TABLE]

and, writing $\tilde{B}_{{\tilde{\alpha}}{\tilde{\beta}}}=a_{\tilde{\alpha}}^{\dagger}a_{\tilde{\beta}}$ with bosonic operators, then

[TABLE]

We assume that $\mathcal{L}_{q}^{\rm loc}$ is diagonalizable (with left and right eigenvectors) as $(\mathcal{L}_{q}^{\rm loc})_{{\tilde{\alpha}}{\tilde{\beta}}}=\sum_{j}Z_{{\tilde{\alpha}}j}\,\lambda_{j}\,Z^{-1}_{j{\tilde{\beta}}}$ for a non-singular matrix $Z$, where $j=0$ corresponds to the steady state. Then we define new bosonic operators via the non-unitary Bogoliubov transformation $\tilde{a}^{\prime}_{i}=\sum_{\tilde{\alpha}}Z_{{\tilde{\alpha}}i}a_{\tilde{\alpha}}^{\dagger}$, $\tilde{a}_{i}=\sum_{\tilde{\alpha}}(Z^{-1})_{i{\tilde{\alpha}}}a_{\tilde{\alpha}}$. These operators still satisfy the canonical commutation relations [$\tilde{a}_{i},\tilde{a}_{j}^{\prime}]={\tilde{\delta}}_{ij}$, though $\tilde{a}_{i}^{\prime}\neq\tilde{a}_{i}^{\dagger}$. As shown in Appendix G, in this language, the steady state of the many-body Liouvillean (28) is therefore the boson “condensate” $|{\Omega}\rangle=\frac{(\tilde{a}_{0}^{\prime})^{q}}{\sqrt{q}!}|{0}\rangle$ where $|{0}\rangle$ is the bosonic vacuum. Elementary excitations with respect to this state can be constructed with a Bogoliubov (mean-field) approach by defining a variational wave-function $|{\psi}\rangle=\sum_{j}\psi_{j}\frac{(\tilde{a}_{0}^{\prime})^{q-1}}{\sqrt{(q-1)!}}a_{j}|{0}\rangle$, for $j\neq 0$ and optimising over the amplitudes $\psi_{j}$. These states are motivated by the analytic solution of the integrable model considered in Section IV, where the excited states with minimal gap have a single quasi-particle excitation. Although mean-field techniques have been highly studied mostly for Hamiltonian systems Blaizot and Ripka (1986), they can be extended also to non-normal operators Laestadius and Kvaal (2017) where left and right eigenvectors form a bi-orthonormal basis. Within this variational formalism we show in Appendix G that the four-body interaction in (28) does not alter the eigenstates, which are therefore exactly given by the bare single-particle eigenstates $|{\Omega^{\rm exc.}_{j}}\rangle=\frac{(\tilde{a}_{0}^{\prime})^{q-1}}{\sqrt{(q-1)!}}a_{j}|{0}\rangle$ with exact eigenvalue $\lambda_{j}$, for any $q$. This shows that the eigenvalues, at least in the low-energy subspace, are not “renormalized” for larger values of $q$. The obtained states $|{\Omega^{\rm exc.}_{j}}\rangle$ are indeed the symmetric combination of (15), which, as shown before, are an exact eigenstate of $\mathcal{L}_{q}$. Within this simple mean-field treatment there are no other eigenvalues with a smaller gap than $\min_{j}|\Re[\lambda_{j}]|$. Therefore, the final outcome of the mean field treatment is that, at least for fully symmetric states, the Liouvillean gap is constant as a function of $q$.

III.3 Counterexample to the mean-field treatment

The mean field treatment of the previous section, based on single particle excitations, predicts that the Liovillean gap is independent on $q$, as long as the mean field approach is accurate. Also the rigorous Bethe-Ansatz treatment of Section IV, valid for a particular integrable model, will show that the Liouvillean gap is independent on $q$, by explicitly showing that the states with minimal gap are made by unpaired particles. That rigorous treatment thus justifies the mean-field approach, at least for that particular model. However, here we show that the predictions of the mean-field theory cannot be general by finding a counterexample where a state with two bounded particles (hence appearing for $q\geq 2$) may have a lower gap.

We construct this counterexample via symmetry arguments. Clearly in the fully controllable case $H$ and $V$ must not share a symmetry – otherwise only symmetric unitaries can be obtained – but this lack of common symmetries is not sufficient. Indeed, generically, in tensor copies there may be other non-trivial symmetries but, because of the Schur-Weyl duality, in the fully controllable case only the permutation symmetries can remain. Suppose now that our system is not controllable because there exists an operator $\tilde{X}$, different from a permutation operator, such that $[H^{\oplus p},\tilde{X}]=[V^{\oplus p},\tilde{X}]=0$ and that the solutions of $[H^{\oplus q},X]=[V^{\oplus q},X]=0$ for $q<p$ are only permutation operators. In this case, Eq. (11) would be valid for $q<p$, but not when $p=q$, as the symmetry $\tilde{X}$ introduces an extra steady state. Then, suppose that we restore full-controllability by adding a small $\mathcal{O}(\epsilon)$ term in either $H$ or $V$ such that the operator $\tilde{X}$ is not a symmetry anymore (we say that the symmetry $\tilde{X}$ is explicitly broken). This splits the extra steady state into an eigenvector with small $\mathcal{O}(\epsilon)$ eigenvalue which, for small enough $\epsilon$ can be smaller than the gap, obtained when $q<p$. If this counterexample can be constructed, then the gap for $q<p$ may be different from the gap at $q=p$. Below we show that this construction is indeed possible already with $p=2$ and that these extra eigenstates correspond to bound particles in the many-body framework.

As shown in Refs. Zeier and Schulte-Herbrüggen (2011); Zimborás et al. (2015) a rather surprising necessary and sufficient condition for controllability is that there are exactly two independent solutions of the equations $[H^{\oplus 2},X]=[V^{\oplus 2},X]=0$. Nontheless, a simpler necessary condition (though not sufficient Zeier and Schulte-Herbrüggen (2011)) is the absense of non-zero solutions to the set of equations

[TABLE]

Taking the complex conjugate of Eq.(29) we find that $Q$ satisfies $Q^{*}H+H^{T}Q^{*}=Q^{*}V+V^{T}Q^{*}=0$, as $H$ and $V$ are Hermitian. Because of this, $QQ^{*}$ commutes with both $H$ and $V$ and, owing to the Schur’s lemma, $QQ^{*}$ is proportional to the identity. Refs. Obata (1958) proved that $QQ^{*}=1$ when $Q$ is symmetric and $QQ^{*}=-1$ when $Q$ is anti-symmetric. If there are non-zero solutions of (29), then the system is not controllable and there are extra steady states such as the bosonic paired state for $q=2$

[TABLE]

Indeed, for both $Q$ symmetric and anti-symmetric $Q\otimes Q^{*}$ is symmetric, thus justifying the bosonic approach. The proof can be readily obtained from (28), indeed for both $X=H,V$

[TABLE]

so because of Eq.(29) we find $\mathring{H}_{{\tilde{\gamma}}{\tilde{\delta}}}a^{\dagger}_{\tilde{\gamma}}a_{\tilde{\delta}}|{\psi_{Q}}\rangle=\mathring{V}_{{\tilde{\gamma}}{\tilde{\delta}}}a^{\dagger}_{\tilde{\gamma}}a_{\tilde{\delta}}|{\psi_{Q}}\rangle=0$, namely $\mathcal{L}_{2}|{\psi_{Q}}\rangle=0$. Hence, the extra symmetry $Q$ introduces a pairing between bosons in the steady state, which is expressed by Eq. (30) – note that it is indeed a pairing because $[Q\otimes Q^{*}]_{{\tilde{\alpha}}{\tilde{\beta}}}\neq Q_{\tilde{\alpha}}Q^{*}_{\tilde{\beta}}$ since $Q$ is a matrix.

As discussed before, we can restore controllability by explicitly breaking the symmetry (29) with small terms: $QH^{T}+HQ=\epsilon_{H}$, $QV^{T}+VQ=\epsilon_{V}$ where at least one between $\epsilon_{V}$ or $\epsilon_{H}$ has to be non-zero, otherwise the system is not controllable. In this case $|{\psi_{Q}}\rangle$ is not a steady state but, within first order perturbation theory, can be used to create a state with eigenvalue ${\tilde{\delta}}=\mathcal{O}(\epsilon_{V},\epsilon_{H})$. In particular, one can construct specific examples where $\epsilon_{V}$ and $\epsilon_{H}$ are much smaller than the gap $\lambda^{*}$ of $\mathcal{L}_{1}$ so that ${\tilde{\delta}}<\lambda^{*}$. Therefore, exploiting these broken symmetries we can construct counterexamples where the gap changes as a function of $q$. The simplest example is a two spin system with $H=(\sigma_{1}^{x}\sigma_{2}^{x}+\sigma_{1}^{y}\sigma_{2}^{y}+\sigma_{1}^{x})+\epsilon\sigma_{1}^{z}\sigma_{2}^{z}$ and $V=\sigma_{1}^{y}$, where $\sigma^{\alpha}_{j}$ are the Pauli matrices acting on the spin $j$. For instance, for $\epsilon=0.1$ the gap of $\mathcal{L}_{1}$ is $\approx 0.45$ while the gap of $\mathcal{L}_{2}$ is $\approx 0.05$.

In spite of this counterexample, we have observed that in most numerical examples, performed for small values of $d$ and $q$ with a random choice of $H$ and $V$, the Liouvillean gap is constant as a function of $q$. This allows us to conjecture that “typically”, namely for most choices of $H$ and $V$, the Liouvillean has a constant gap, as predicted by the mean-field approach. Since in Eq. (12) each copy interacts with all the others, this conjecture is supported by the well-known validity (see e.g. Mézard et al. (1987)) of the mean-field solution in long-range models.

IV The controllable quantum walk

We focus on a specific model that is of experimental interest, namely a single-particle hopping in a one-dimensional lattice; see Fig. 1. This framework can describe different physical systems, such as a spin impurity in a spin chain, a single electronic excitation in quantum dot arrays and a photon traveling in a one-dimensional photonic chip. The resulting quantum walk can be modeled via the Hamiltonian

[TABLE]

where $|{n}\rangle$ represents the state in which the walker is in position $n$, and $L$ is the length of the chain. This Hamiltonian has found numerous applications in quantum transport problems and remote entanglement generation in spin chains Nikolopoulos et al. (2014); Banchi et al. (2011a); Burgarth and Bose (2005); Banchi et al. (2011b).

Moreover, we consider a local control field on a single site of the chain, namely the $c$-th site, which is modeled by Hamiltonian term $g(t)V$, where $V{=}|{c}\rangle\langle{c}|$ and $g(t)$ is a time dependent control profile. One can show that the chain is controllable provided that $c$ and $L{+}1$ are co-prime numbers Wang et al. (2012); Burgarth et al. (2013). For simplicity, in the following we set $c=1$. The above hopping Hamiltonian with local control can be realized in many physical systems; for example, in reconfigurable photonic chips Perez-Leija et al. (2013); Pitsios et al. (2016), where the different control pulses can be obtained by electrically tuned on-chip heaters Carolan et al. (2015).

In the following we evaluate the Liouvillean gap for all possible values of $q$ in the strong driving limit, namely when $\sigma\gg 1$. The opposite weak driving limit is discussed in appendix C for the single-particle $q=1$ case. We start by considering two important cases, namely the fully symmetric and fully anti-symmetric representation where $B_{\alpha\beta}=a^{\dagger}_{\alpha}a_{\beta}$ for either bosonic or fermionic degrees of freedom. We then extend our analysis to the general case.

IV.1 Gap analysis: fully-symmetric representation

We consider first the fully symmetric representation where $B_{\alpha\beta}=a^{\dagger}_{\alpha}a_{\beta}$ so one can omit the index $u$ from the equations of Section III.1. Plugging the operators $H$ and $V$ of the controllable chain into Eq. (27) one finds the following Liouvillean

[TABLE]

To diagonalize the above operator we assume that $\sigma\gg 1$ and we study the “low-energy” effective dynamics. In that limit the dissipative part $\sigma\mathcal{D}=\frac{\sigma}{2}(n^{\uparrow}_{1}-n^{\downarrow}_{1})(n^{\uparrow}_{1}-n^{\downarrow}_{1})$ has either eigenvalue 0 or $\sigma\gg 1$. With a perturbative approach, discussed in Appendix D, we decouple the latter “high-energy” subspace and obtain an effective Liouvillean acting in the low-energy sector. From a first order expansion as a function of $\sigma^{-1}$ the effective Liouvillean is given by

[TABLE]

where $g_{k}=\frac{2}{L}\sin^{2}\left(\frac{\pi k}{L}\right)$, $\tilde{a}_{k\updownarrow}=\sum_{\alpha=1}^{L-1}\frac{2}{L}\sin^{2}\left(\frac{\pi k\alpha}{L}\right)a_{\alpha{+}1,\updownarrow}$ and $a_{1\updownarrow}\equiv\tilde{a}_{0\updownarrow}$. We call now $K_{i}^{+}=\tilde{a}_{i\uparrow}^{\dagger}\tilde{a}_{i\downarrow}^{\dagger}$, $K_{i}^{-}=(K_{i}^{+})^{\dagger}$ and $K_{i}^{z}=(\tilde{n}_{i\uparrow}+\tilde{n}_{i\downarrow}+1)/2$ and note that these operators satisfy the SU(1,1) commutation relations

[TABLE]

With these definitions we find then

[TABLE]

where $K_{i}\cdot K_{j}\equiv-(K_{i}^{+}K_{j}^{-}+K_{i}^{-}K_{j}^{+})/2+K_{i}^{z}K_{j}^{z}$ is the SU(1,1) invariant product, namely the analogue of the Heisenberg interaction. The model (36) is a SU(1,1) Gaudin model Gaudin (1976), which is known to be exactly solvable with the Bethe-Ansatz approach. We explicitly diagonalize it in the appendix E by applying Richardson’s method Richardson (1968). We find that the eigenvalues of the Liouvillean $\hat{\mathcal{L}_{q}}$ are

[TABLE]

where the non-negative integers $n_{k}$ parametrize the number of unpaired particles in mode $k$ (see the discussion in Appendix E) and the $E_{\alpha}$ are either zero or the solution of the non-linear set of equations

[TABLE]

where $E_{\alpha}=1/\omega_{\alpha}$. From that expression it is clear that the steady state corresponds to $E_{\alpha}=0$ and $n_{k}=0$, for each $\alpha$ and $k$. Solutions to the above equations are known to be related with the roots of Heine-Stieltjes polynomials (see e.g. Szegö (1959)). By exploiting this relationship, one finds that all the solutions $\omega_{\alpha}$ of (38) are real, different from each other, and different from the poles of (38). Moreover, $g_{k}=g_{L-k}$ so the sum in (38) can be restricted to the first half where $g_{k}<g_{k+1}$. The roots of the Heine-Stieltjes polynomials have also the important property that they lie inside the intervals $2g_{k+1}^{-1}<\omega_{\alpha}<2g_{k}^{-1}$ for some $k$, so that $2E_{\alpha}>\min_{k}g_{k}=g_{1}$. This constraint allows us to find rigorously the gap of the Liouvillean $\hat{\mathcal{L}_{q}}$. Indeed, thanks to the latter inequality, the paired states have a larger gap than the unpaired ones, so we can focus only on the solutions where $E_{\alpha}=0$. The minimum gap is then obtained when $n_{1}=n_{L-1}=1$ and $n_{k}=0$ otherwise. This is an allowed state (for $L>2$) as it satisfies all the constraints and provides the gap

[TABLE]

This gap is exact in the strong driving limit, can be achieved already at $q=1$ and is the same for all higher values of $q$, as we have shown that there are no smaller non-zero eigenvalues. Therefore, we proved here explicitly that in the strong driving limit the gap is independent on the number of copies $q$. In the following sections we extend this result, which up to now is restricted to the fully-symmetric representation, to show that (39) is indeed the gap, irrespective of the chosen representation.

IV.2 Gap analysis: anti-symmetric representations

We first consider another particular case, namely the fully anti-symmetric representation, that will be used as a basis for the general solution discussed in the next section. We start from (27) and we write $B^{\alpha}_{ij}=a_{i\alpha}^{\dagger}a_{j\alpha}$ with fermionic creation and annihilation operators. Repeating the effective Liouvillean description of the previous section we find

[TABLE]

where $S_{0}\cdot S_{k}=\sum_{\alpha=x,y,z}S_{0}^{\alpha}S_{k}^{\alpha}$ refers to the SU(2)-invariant, product, namely the spin Heisenberg interaction, $S^{\pm}_{j}=S_{j}^{x}\pm iS_{j}^{y}$, and where we have defined $S_{j}^{-}=\tilde{a}_{j\uparrow}\tilde{a}_{j\downarrow}$, $S_{j}^{+}=(S_{j}^{-})^{\dagger}$ and $S_{j}^{z}=(\tilde{a}_{j\uparrow}^{\dagger}\tilde{a}_{j\uparrow}+\tilde{a}_{j\downarrow}^{\dagger}\tilde{a}_{j\downarrow}-1)/2$. It is simple to verify that the above operators satisfy the SU(2) commutation relations on the same site, and commute on different sites, so that Eq. (40) is equivalent to the central spin model first studied by Gaudin Gaudin (1976). The diagonalization of the Gaudin Heisenberg Hamiltonian proceeds along the same lines of the SU(1,1) one. There are two main differences: (i) the different sign in (40) and (36) and (ii) becayse of the Pauli exclusion principle the number of particles $n_{k}$ per mode $k$ is limited to either 0 or 1. We find then that the eigenvalues are given by Eq. (37), where the non-zero energies $E_{\alpha}$ are the solutions of

[TABLE]

However, because of the different sign in (41), we cannot relate the solutions of (41) to the roots of the Heine-Stieltjes polynomials, so we cannot bound the gap using the argument of the fully symmetric case. Nonetheless, in the next section we consider a more general technique, valid for all the representations, where such a bound can be obtained using physical arguments borrowed from classical electrostatics.

IV.3 General gap analysis

As we have discussed in Section III.1, a general representation of the SU(L) algebra can can be obtained via extended creation and annihilation operators Rowe et al. (2012), namely $B_{\alpha\beta}=\sum_{u}\tilde{a}^{\dagger}_{\alpha u}\tilde{a}_{\beta u}$ for either bosonic or fermionic operators. We use the fermionic representation for convenience, since our derivation uses the particle-hole symmetry that is a non-unitary operation in bosonic systems (see e.g. Blaizot and Ripka (1986)). Because of the Pauli exclusion principle, in order to satisfy the constrain $\sum_{\alpha}B_{\alpha\alpha}=q$, the auxiliary index $u$ has to run from 1 to $q$. Performing the same perturbative approach of Appendix D, valid in the strong driving limit $\sigma\gg 1$, one finds that the effective Liouvillean $\hat{\mathcal{L}_{q}}$ can be written in the diagonal basis of the Hamiltonian as

[TABLE]

The above Hermitian operator corresponds to the purely dissipative Liouvillean

[TABLE]

where $\tilde{V}_{k}=|{\omega_{k}}\rangle\langle{\omega_{0}}|$, being $|{\omega_{k}}\rangle=\sum_{j=1}^{L-1}\frac{2}{L}\sin\left(\frac{\pi jk}{L}\right)^{2}|{j{+}1}\rangle$ and $|{\omega_{0}}\rangle=|{1}\rangle$. One can check that the operators $\tilde{V}_{k}$ and their Hermitian conjugate form a controllable set, so the steady state of the effective Liouvillean coincides with the original one. We now perform two transformations. The first one is the Jordan-Wigner transformation to obtain proper fermionic degrees of freedom, namely where creation/annihilation operators with different indices $\uparrow$ and $\downarrow$ anti-commute. The second-one is a particle-hole transformation in the spin-down sector. These transformations are implemented together by defining $W=\prod_{ju}e^{i\tilde{a}^{\dagger}_{ju\uparrow}\tilde{a}_{ju\uparrow}}$ and setting $a_{ju\uparrow}=\tilde{a}_{ju\uparrow}$ and $a_{ju\downarrow}=W\tilde{a}^{\dagger}_{ju\uparrow}$. Eq. (42) then becomes

[TABLE]

where $X^{(j)}_{\alpha\beta}=\left(a^{\dagger}_{j\alpha}a_{j\beta}-a_{j\beta}a^{\dagger}_{j\alpha}\right)/2$ and the Greek letters refer to the multi-index composed by the auxiliary index and the “effective spin” index, i.e. $\alpha=(us)$ where $u=1,\dots,q$ and $s=\{\uparrow,\downarrow\}$. The traceless operators $X^{(j)}_{\alpha\beta}$ satisfy the SU(2q)⊕L commutation relations,

[TABLE]

so that Eq.(43) represents a SU(2q) version of the Gaudin model. Indeed, Eq.(43) is invariant under the Bogoliubov transofmation $a_{j\alpha}\to\sum_{\beta}U_{\alpha,\beta}\,a_{j\beta}$, where $U$ is a unitary $(2q)\times(2q)$ matrix. SU(2q) has $(2q)^{2}-1$ generators, so one operator in (44) is dependent on the others. This is shown by the equation $[\sum_{\alpha}X^{(j)}_{\alpha\alpha},X^{(k)}_{\beta\gamma}]=0$ for each $\beta$ and $\gamma$. Going back to the original representation, namely performing back the particle-hole transformation, one finds that

[TABLE]

The Gaudin-like model (43) has been solved for different algebras (namely not only the SU(1,1) and SU(2) cases discussed before) in Refs. Ushveridze (1994); Falceto and Gawȩdzki (1997), while the duality between the different models that can be obtained by exploiting the auxiliary indices has different ramifications in mathematical physics (see e.g. Mukhin et al. (2008) and references therein), especially due to its connections with the Knizhnik-Zamolodchikov equation Mukhin et al. (2008); Feigin et al. (1994). In Appendix F we exploit the general solution Ushveridze (1994); Falceto and Gawȩdzki (1997) of the Gaudin model (43), valid when the operators $X$ define any semi-simple Lie algebra, to obtain the eigenvalues of the Liouvillean (43) when the SU(2q) operators are defined via the fermionic representation (45). As in the fully-symmetric and fully-antisymmetric case discussed in the previous sections, the eigenvalues of $\hat{\mathcal{L}}_{q}$ are parametrized by non-negative integers $n_{\uparrow j}$ and $n_{\downarrow j}$, and are given by

[TABLE]

where $\omega_{j,\alpha}$ for $j=1,\dots,2q-1$ are the solutions of

[TABLE]

being $z_{0}=0$, $\mu^{0}_{j}=\delta_{qj}$, and, for $k>0$, $z_{k}=2g_{k}^{-1}$ and $\mu^{k}_{j}=\delta_{j,q}(1-\delta_{n_{\downarrow k}>0}-\delta_{n_{\uparrow k}>0})+\delta_{j,q+n_{\uparrow k}}+\delta_{j,q-n_{\downarrow k}}$. In (47) we set $\omega_{0,\beta}=\omega_{2q,\beta}\to-\infty$ namely, in other terms, for $j=1$ or $j=2q-1$ one of the two fractions in the second line is zero.

Owing to the similarity between Eqs. (46) and (37), if we can show that the solutions of (47) satisfy the inequality $2\omega^{-1}_{q,\alpha}>g_{k}$ for each $\alpha$ and $k$, then we can straightforwardly apply the reasoning of Section IV.1 to prove that the gap is indeed given by Eq. (39) for any representation. However, the sign difference between Eqs. (47) and (38) prevents us from using the theory of Heine-Stieltjes polynomials to prove that inequality, as we did in Section IV.1. Here we use a different approach, used also in Ref. Ushveridze (1994) for a different purpose, which is based on mapping the mathematical equations (47) to an electrostatic problem, and then use our classical physics intuition. Following Ref. Ushveridze (1994) we define the two-dimensional vector $\vec{\omega}_{j\alpha}$ whose real components are the real and imaginary part of $\omega_{j\alpha}$ and interpret those vectors as the positions of some particles with index $\alpha$ and species $j=1,\dots,2q-1$. The equations (47) can then be interpreted as the conditions for an extremum of the function $\mathcal{W}(\{\omega\})$ defined as

[TABLE]

where $\vec{z}_{k}=(z_{k},0)$ and the Cartan matrix $C_{ij}$ has non-zero components only on the diagonal, where $C_{ii}=2$, and for $|i-j|=1$, where $C_{ij}=-1$.

This shows that the problem of finding a solution to the system of equations (47) is equivalent to the problem of finding the equilibrium positions of a set of particles in a two-dimensional plane interacting via the logarithmic potential (48). That potential is analogous to the electrostatic potential since the Coulomb interaction in 2D is logarithmic. Particles of the same species repel each other, while particles with nearest-neighbour species attract each other. Finding the equilibrium positions of those particles is in general quite complicated, although the problem can be solved explicitly in the thermodynamic limit Roman et al. (2002). At first sight one may think that the problem has no solutions since the potential (49) is unstable. However, because of the Z2 symmetry ($\Im[\omega_{j,\alpha}]\to-\Im[\omega_{j,\alpha}]$), due to the fact that the $z_{k}$’s are reals, all the forces on the real line are longitudinal. This property allows us to seek for solutions of Eq. (47) in the class of real numbers Ushveridze (1994). On the real line, the problem becomes stable and one-dimensional. An example of this effective one-dimensional potential is shown on Fig. 2 where one can see the two unbounded regions for $\omega<\min_{k}z_{k}$ and for $\omega>\max_{k}z_{k}$, where no solutions can exist. Therefore, this electrostatic analogy shows that the only stable solutions with finite $\omega_{i\alpha}$ can be found only between poles of $\mathcal{V}_{i}(\omega)$, or, in other terms, that the solutions of the non-linear set of equations (47) satisfy the constraint $\min_{k}z_{k}<\omega_{j\alpha}<\max_{k}z_{k}$, i.e. $2\omega_{j\alpha}^{-1}>\min_{k}g_{k}$. This, together with the discussion of Section IV.1, shows that Eq.(39) is indeed the gap of the Liouvillean $\hat{\mathcal{L}}_{q}$ in the strong-driving limit.

IV.4 Numerical results for the controllable chain

In the previous sections we have done an extensive theoretical analysis to show that, in a chain controlled on one boundary, the Liouvillean gap in the strong-driving limit is constant as a function of $q$ and scales as $\propto L^{-3}$ as a function of the length $L$ of the chain – this scaling is consistent with what has been obtained in spin chains with boundary dissipation Žnidarič (2015). The scaling $\propto L^{-3}$ is obtained also in the weak driving limit discussed in Appendix C, though that analysis is valid only for $q=1$. Nontheless, in all our numerical experiments obtained for small values of $L$ and $q$ we found that the gap is constant as a function of $q$ over the whole range of $\sigma$.

In Fig. 3 we study the Liouvillean gap and show that the theoretical predictions of the strong and weak driving limits are very accurate in their respective limit of validity. Moreover, we found that the accuracy of the strong driving limit is not affected by the length of the chain. This is shown indeed in the inset Fig. 3 where one observes an almost constant behaviour as a function of $L$.

In Fig. 4, on the other hand, we show that the Liouvillean gap scales as $L^{-3}$ for different values of $\sigma$. This scaling has been predicted in the strong and weak driving limits by Eqs. (39) and (76). However, Fig. 4 shows that such scaling is valid also for $\sigma\approx 2$ where neither the strong nor the weak coupling limit holds (compare e.g. the values of Fig. 4 and Fig. 3).

In Fig. 5 we study the relationship between the Liouvillean gap and the gap $s^{*}(t)$ in the singular values of $e^{t\mathcal{L}_{q}}$ which is a good estimate of the convergence time (see section II.4). As expected, both in the strong and weak coupling limit the $s^{*}(t)$ converges to $e^{-\lambda^{*}t}$ much earlier than mixing time-scales. Therefore, in these regimes, one finds that the convergence time is basically $1/\lambda^{*}$. On the other hand, for $\sigma=2$ the matching between $e^{-\lambda^{*}t}$ and $s^{*}(t)$ only happens at longer times. Therefore, as expected from the analysis of Section II.4, in this regime there is a correction to the mixing time due to the norm of the left and right eigenvectors. Nonetheless, similarly to the Liovillean gap, our numerical simulations for small values of $L$ and $q$ show that also the singular value gap is independent on $q$ over the whole range of $\sigma$. Therefore, we argue that it may be a general feature of this model that the resulting convergence time is independent on $q$.

Finally we consider a stochastic simulation of the evolution of a controllable chain with random fields: we generate several random driving functions (2) and, for each function, we calculate the corresponding unitary evolution and then study the statistics of the generated unitary matrices. To test whether the resulting distribution approximates the Haar measure we decompose each unitary into the $L^{2}$ angles introduced in Ref. Spengler et al. (2012). Using a simple reparametrization of these angles one can write the Haar measure as

[TABLE]

namely as a uniform distriution of the angles $\varphi_{j}$ in the range $[0,2\pi]$. Therefore, testing whether the resulting distribution approximates a Haar measure is equivalent to testing whether the angles $\varphi_{j}$ are distributed as a multinomial uniform distribution. In Fig. 6 we do a simple test to verify the distribution of the angles $\varphi_{j}$: we divide the interval $[0,2\pi]$ into 25 bins and plot, as a 3D histrogram, the matrix whose elements $(i,j)$ are the number of times that the angle $\varphi_{i}$ is found in the $j$-th bin. As Fig. 6 shows, the distribution of the unitary matrices is far from uniform both in the noncontrollable case and in the controllable case after a short time (upper panel). Nonetheless, in spite of the finite number of samples, after a long time ($t\approx 55$) in the controllable case the angles’ distribution is almost flat (lower panel), thus showing that the resulting unitary matrices are approximately distributed according to the Haar measure.

V Other applications

V.1 Multi-point correlation functions

Here we discuss some direct applications, beyond $q$-design, of the main findings of our paper. In boson sampling experiments the output probability is proportional to $|{\rm per}(\tilde{U})|^{2}$, being ${\rm per}(\tilde{U})$ the matrix permanent of the $q\times q$ matrix $\tilde{U}$, where $\tilde{U}$ is built from some columns and rows of a $L\times L$ Haar-uniform matrix $U$ Broome et al. (2013); Spring et al. (2013); Crespi et al. (2013). Therefore

[TABLE]

where $\sigma,\sigma^{\prime}$ are permutations in the symmetric group $S_{q}$, $\mathcal{K}_{\rm b.s.}$ is a suitable index contraction operator and $U^{\otimes q,q}=U^{\otimes q}\otimes(U^{\otimes q})^{*}$ as in Eq. (4).

A similar expression arises in the evaluation of multi-point correlation functions in quasi-free particle-preserving bosonic and fermionic models. If $U$ is the $L\times L$ one particle evolution matrix from time [math] to time $t$ and $a_{j}(t)=\sum_{k}U_{jk}a_{k}(0)$, then because of the Wick’s theorem

[TABLE]

where $\mathcal{K}_{\rm m.p}$ depends on the initial two-point correlation functions $\langle a_{i}^{\dagger}(0)a_{j}(0)\rangle$. Expressions like (53) arise also in XY spin chains, which can be mapped to a quasi-free fermionic model via the Jordan-Wigner transformation Lieb et al. (1961). For instance, the driven XY model

[TABLE]

can be mapped, in the single-particle subspace, to the driven quantum walk of Section IV. Calling $U$ the resulting single-particle evolution, then in any subspace long-range spin operators $S^{\alpha}_{i}S^{\beta}_{j}$, for $\alpha,\beta\in\{x,y\}$ can be written as a combination of fermion strings as in (53) where $q=|i-j|$ for $i\neq j$. Therefore, with a suitable $\mathcal{K}_{\rm XY}$ that depends on the initial correlations, one can write the dynamical long-range correlations between spin operators in an XY chain as

[TABLE]

for $\alpha,\beta\in\{x,y\}$. Similarly, $\langle S^{z}_{i}(t)S^{z}_{j}(t)\rangle=\operatorname*{Tr}\left[U^{\otimes 2,2}\mathcal{K}^{zz}_{\rm XY}\right]$.

In all the above cases we can bound the convergence of the random dynamics to the values expected from the Haar distribution. Indeed, for any $\mathcal{K}$

[TABLE]

where we used (4). Thanks to the analysis of Section II.4, and since the gap (39) for the controllable quantum walk is independent on $q$, one can then bound the expected errors in all the above cases. For boson sampling experiments, this shows how the error depends on the number $q$ of bosons, while for XY spin chains, it shows how the error decays as a function of the distance $q$ between spins.

V.2 Estimation of the control time

We show here that the mixing time, which is easy to compute especially for $q=1$, can give an estimation of the control time. Fixing $H$ and $V$, for how long does one have to drive the system in order to achieve a generic target gate? If after the time $T^{*}_{ex}$ the random evolutions are Haar-randomly distributed, then the control time to obtain a certain gate $U$ satisfies $T_{c}(U)<T^{*}_{\rm ex}$. However, for approximate $q$-design, $T^{*}$ provides only a rate of convergence, rather than a sharp bound. This results into an error, which may also be due to the fact that the target gate $U$ is not achievable yet at time $T^{*}$. However, after a time $\tau T^{*}$ this error probability exponentially decreases as a function of $\tau$. We can thus regard $T^{*}$ as an estimation for $T_{c}$. An estimation of the mixing time $T^{*}$ can be easily obtained for any choice of $H$ and $V$ via the inverse of the gap $\lambda^{*}$, which depends on $\sigma$ (see e.g. Fig. 3). Since $T_{c}$ does not involve any specific properties (amplitudes, frequencies) of the pulse, one has to compare it with $T^{*}_{\rm min}=\min_{\sigma}T^{*}(\sigma)\simeq T^{*}(\sigma\simeq 2.5)\approx 0.055\,L^{3}$.

In order to estimate $T_{c}$ we perform a numerical experiment with the QuTip quantum control package Johansson et al. (2012). We consider the model (31) and, for each length $L=10,\dots,20$, we generate a Haar-random unitary $U$ and find the time $T_{c}$ as the minimal time for which the program converges. We find that $T_{c}$ obtained in this way scales as $T_{c}\approx 0.069\,L^{3}$. This shows two remarkable facts: (i) the values of $T_{c}$ and $T^{*}_{\rm min}$ are very close for $L=10,\dots,20$; (ii) both $T_{c}$ and $T^{*}_{\rm min}$ exhibit the same scaling with the length $L$, so it is expected that this close relationship is maintained also for larger $L$. In view of our findings, one can find an empirical upper bound on $T_{c}$ as $3T^{*}_{\rm min}/2$.

VI Conclusions and perspectives

In this paper we have studied the quantum dynamics resulting from a stochastic driving of quantum many-body systems, and we have answered the following questions: when, and how rapidly, the dynamics of a driven quantum system is equivalent to to a fully uniform random evolution, namely under unitaries sampled from the Haar measure. The first major finding is that, when the system is fully controllable and the stochastic signal has finite correlation time, then its random dynamics converges to the Haar distribution in the “long time” limit. The second major result is about the estimation of the driving time $T^{*}$: this is done by studying the deviations from the Haar distribution using the framework of approximate $q$-design, and using second-quantization to map the problem into the estimation of the mixing time in an open quantum many-body Liouvillean with $2q$ virtual particles.

We have performed a thorough analysis of the Markovian limit (e.g. white noise) using tools from the theory of dynamical semigroups, and we found upper bounds on $T^{*}$ in terms of the gap of the Liouvillean operator. We studied the mean field solution of the resulting many-body model, which predicts a constant Liouvillean gap as a function of $q$, and we have shown its limitations via symmetry breaking arguments. Nonetheless, we found that the mean-field predictions are correct in a wide variety of different numerical studies, obtained with random choices of $H$ and $V$, and match with the analytic solution of a particular model, namely a one dimensional system with strong control on one of its boundaries. The latter analytic solution has been obtained by mapping the effective Liovillean to an exactly solvable model, and then using Bethe-Ansatz techniques to explicitly show that the excited states with smallest gap are built from unpaired quasi-particles, as in the mean field treatment. We have then corroborated our predictions with numerical simulations, putting strong evidence that the considered one-dimensional model provides a quantum expander with a constant mixing time as a function of $q$. Therefore, our results show that certain driven physical systems can provide a significant advantages over random quantum circuits where the mixing time increases polynomially as a function of $q$ Brandao et al. (2012).

The results presented in this paper have many applications. The first one, already discussed, is a physically motivated approach to generate pseudo-uniform random unitary operations, which have many applications in quantum information processing protocols. The one-dimensional system that is extensively analyzed in this paper is motivated by the recent experiments with integrated photonic circuits Perez-Leija et al. (2013); Pitsios et al. (2016), where random unitary operations have been used in the first small-scale experimental observations of boson sampling Broome et al. (2013); Spring et al. (2013); Crespi et al. (2013). The results presented in this paper enable the implementation of random operations in integrated photonic chips that, being based on noisy quantum walks rather than carefully designed multi-mode beam splitters and phase shifters, are much simpler to fabricate for a larger number of modes. Therefore, our results provide a new avenue to prove quantum supremacy in boson sampling experiments.

Moreover, we have considered other applications, such as the dynamics of correlation functions in an XY spin chain, and the estimation of the control time $T_{c}$, one of the major open problems for quantum control. Given a target unitary $U$ and the physical interactions described by $H_{0}$ and $V$, how can we chose $T_{c}$ such that $U$ is achievable by driving the system for a time $T_{c}$? With numerical experiments, performed on $L$-site chains, we found that both $T_{c}$ and $T_{*}$ are very close for $L=10,\dots,20$, and both scale as $L^{3}$. Hence, the mixing time $T_{*}$ under random signals provides an easily computable estimation of $T_{c}$, for any $H_{0}$ and $V$.

Finally, there are several applications in quantum many-body physics, where the interplay between quantum many-body effects and noise is currently a subject of intensive study in many area, such as spin glass Mézard et al. (1987), the fast scrambling of quantum information Hayden and Preskill (2007); Sekino and Susskind (2008), and many-body localization Pal and Huse (2010); Ponte et al. (2015). The explicit one dimensional model discussed in Section IV is a single-particle model, where many-body physics arises due to unitary $q$ design, which introduces $2q$ virtual particles. An interesting future perspective is the study of random driving in physical interacting many-body systems (e.g. interacting spin systems and/or cold atoms optical lattices). In fact, the competition between physical many-body effects, and those arising from the unitary design, may give rise to novel states of matters and phase transitions Medvedyeva et al. (2016); Žnidarič (2015); Banchi et al. (2014); Diehl et al. (2011); Prosen and Pižorn (2008), produce large amount of entanglement Nahum et al. (2017), and give new insights into the process of thermalization and equilibration Eisert et al. (2015). Haar-random quantum states are known to have, typically, an extensive amount of entanglement Hayden et al. (2006). Since we have shown that any controllable quantum system converges to a maximally mixing dynamics, the real time dynamics will be very hard to simulate numerically in the many-body settings, because of the large amount of entanglement involved. Nonetheless, the controllability requirement provides a sufficient algebraic method to infer, a priori, whether a randomly driven condensed matter system is expected to produce lots of entanglement in the long time limit.

Acknowledgements.

The authors thank S. Bose, E. Compagno, F. Falceto, J. Links, A. Marcus, L. Maccone, S. Maniscalco, E. Mukhin, S. Paesani, R. Santagati, S. Severini, A. Werner, R. Zeier, and Z. Zimborás for interesting discussions. L.B. has received funding for this research from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement No. 308253 PACOMANEDIA. MJK was supported by the Villum foundation. D.B. acknowledges support from the EPSRC grant No. EP/M01634X/1.

Appendix A Gaussian harmonic pulses

To simplify the theoretical description, in this section we consider only $q=1$ and call $\mathcal{E}_{t}$ the quantum channel resulting from the average evolution of the quantum system

[TABLE]

Extensions to higher values of $q$ is straightforward. As described in section II, we now make two assumptions, namely that $g(t)$ is Gaussian and harmonic, where $\mathbb{E}[g(t+s)g(t)]=c(s)$ is independent on $t$ and $\mathbb{E}[g(t)]=0$. In view of these assumptions, we can simplify (57) by expanding the exponentials into the Dyson series, then using the Wick’s theorem to decompose the expectation values and finally resumming the series. The result in the interaction picture is then Ishizaki and Fleming (2009); Banchi et al. (2013)

[TABLE]

where $(I)$ refers to the interaction picture with respect to $H$. If the correlation time is finite then there exist a suitably large $T$ such that $Tc(Ts)\simeq\frac{\sigma}{2}\delta(s)$ where $\delta$ is the Dirac delta function ans $\sigma$ is a constant. In the long-time limit one finds that

[TABLE]

when $t>T$, namely in the Schrödinger picture

[TABLE]

Appendix B Semigroup convergence times

There exist several measures to estimate convergence of a semigroup of completely positive trace preserving (CPTP) maps. The one with the most natural operational interpretation is trace norm convergence, as it reflects the likelihood that the time evolved state can be distinguished from the stationary state at a given time $t$.

[TABLE]

where $T_{\infty}=\lim_{t\rightarrow\infty}e^{t\mathcal{L}}$, and $\epsilon(t)$ is the distinguishability error. A less stringent convergence requirement is to ask whether $e^{t\mathcal{L}^{q}}$ is an expander for a given value of $t$. Then, we want to estimate

[TABLE]

where a hat indicates that the CPTP maps are represented as channels (see Ref. Wolf for more details on the representation of channels). Trace norm convergence and ’spectral convergence’ are related, by noting that

[TABLE]

where $d$ is the dimension of the Hilbert space, and recalling that $||e^{t\mathcal{L}}-T_{\infty}||_{1\rightarrow 1}=\sup_{\rho}||e^{t\mathcal{L}}(\rho)-T_{\infty}(\rho)||_{1}$.

In order to estimate the above norms it is important to recall the spectral properties of quantum dynamical semigroups. The spectrum of a Liouvillian $\mathcal{L}$ has non-positive real part, and there always exists at least one eigenvalue of magnitude zero, corresponding to a stationary state of the semigroup: $\mathcal{L}(\rho)=0$. The rest of the spectrum comes in complex conjugate pairs. The Liouvillian is called unital if it annihilates the identity $\mathcal{L}(\openone)=0$. The Liouvillian in Eq. (7) has this property. A unital Liouvillian is called reversible if $\hat{\mathcal{L}}=\hat{\mathcal{L}}^{\dagger}$, in which case its spectrum is real. Unfortunately, Eq. (7) is not reversible. Convergence of a non-reversible semigroup is governed by the singular values of $e^{t\mathcal{L}}$ rather than its eigenvalues. The singular spectrum of $e^{t\mathcal{L}}$ is equal to the spectrum of $\sqrt{e^{t\hat{\mathcal{L}}}e^{t\hat{\mathcal{L}}^{\dagger}}}$.

It is not difficult to see that the $2\rightarrow 2$ norm is related to the singular spectrum. Let $s_{j}(t)$ be the singular values of $e^{t\mathcal{L}}$, ordered from largest to smallest. The largest has magnitude one. We know that asymptotically $s^{*}(t)=e^{t{\rm Re}[\lambda^{*}]}$, where now $\lambda_{j}$ are the eigenvalues of $\mathcal{L}$ written in decreasing (real part) order, and $\lambda^{*}$ is the gap of $\mathcal{L}$; i.e. the smallest (in magnitude) non-zero real part of any eigenvalue of $\mathcal{L}$. To see this, note that, assuming it has no Jordan blocks, the Liouvillian can be written in its spectral decomposition as

[TABLE]

where $R_{j},L_{j}$ are a bi-orthonormal basis of operators: i.e. ${\rm tr}[L^{\dagger}_{j}R_{k}]=\delta_{jk}$. Importantly, the norm of any given $L_{j},R_{j}$ can be large, which prevents us from getting any rigorous (universal) bounds between the singular values and the eigenvalues. Then,

[TABLE]

Hence, for very large $t$, the convergence is governed by the gap, and $s^{*}(t)\rightarrow e^{t\lambda^{*}}$. In principle we do not know at what scale $e^{-t\lambda^{*}}\gg|\langle R_{j}|R_{j}\rangle|~{}|\langle\psi|L_{j}\rangle\langle L_{j}|\psi\rangle|$.

We argue in the main text, that for the specific model of a controllable quantum walk, the prefactors do not contribute to the asymptotics in the weak or strong coupling limits.

Appendix C Weak driving limit

A convenient approximation for the long-time dynamics in the weak coupling limit $\sigma\ll 1$ is the rotating wave approximation (RWA) Fleming et al. (2010). We consider the case $q=1$ and assume that $V$ is a matrix of real numbers and call $\mathcal{D}=-\sigma\mathring{V}^{2}/2$ the dissipative part in (8). Going to the interaction picture with respect to the Hamiltonian part one finds that $\dot{\rho}_{I}(t)=\mathcal{D}_{I}(t)\rho_{I}(t)$ where in the eigenbasis of $H=\sum_{j}\omega_{j}|{\omega_{j}}\rangle\langle{\omega_{j}}|$ it is

[TABLE]

where $\omega_{ij}=\omega_{i}-\omega_{j}$. The rotating wave approximation consists in neglecting all the terms where $\omega_{ij}\neq\omega_{kl}$, because for large $t$ they are highly oscillating and average out:

[TABLE]

This approximation is expected to hold when

[TABLE]

RWA is related to degenerate perturbation theory. Indeed, the unperturbed ($\sigma=0$) eigenvalues of (8) are given by $|{\Phi^{(0)}_{ij}}\rangle=\sum_{kl}\delta_{\omega_{ij},\omega_{kl}}\alpha_{kl}^{ij}|{\omega_{k}\omega_{l}}\rangle$ with eigenvalue $-i\omega_{ij}$. From degenerate first order perturbation theory we know that, for small $\sigma$, the eigenvalues of Eq. (8) are obtained by diagonalising $\mathcal{D}^{\rm RWA}$, which is block diagonal where each block acts on different degenerate subspaces. The eigenvectors of $\mathcal{D}^{\rm RWA}$ provide the matrices $\alpha^{ij}_{kl}$. Note that since $\mathcal{D}^{\rm RWA}$ is Hermitian, the states $|{\Phi_{ij}^{(0)}}\rangle$ form an orthonormal basis which depends both on $H$ (from the basis $|{\omega_{k}}\rangle$) and $V$ (via the diagonalization of $\mathcal{D}^{\rm RWA}$). Moreover, the real eigenvalues $\Delta_{ij}$ of (67) provide the first order correction to the eigevectors of Eq. (8) that, to the first order in $\sigma^{-1}$ are $-i\omega_{ij}+\Delta_{ij}$. The Liovillean gap is given by the minimum non-zero value of $-\Delta_{ij}$. Similarly one finds the correction to the (right) eigenvector

[TABLE]

where

[TABLE]

Since $\mathcal{S}^{\rm RWA}$ is a Hermitian operator the new vectors in (69) do not form an orthonormal basis.

We now focus on the the the chain discussed in Section IV where $\omega_{k_{j}}=2\cos k_{j}$, $k_{j}=\pi j/(L+1)$, $V^{2}=V$ and we call $W_{ij}=\langle{\omega_{i}}|V|{\omega_{j}}\rangle=\frac{2}{L+1}\sin k_{i}\sin k_{j}$. To simplify the equations we use the compact notation $|{i}\rangle\equiv|{\omega_{i}}\rangle$ and we use $c=1$, namely we assume that the controlled site is the first one. We note that the resonance condition $\omega_{i}-\omega_{j}=\omega_{k}-\omega_{l}$ is achieved in three different cases:

Case 1: $i=k$ and $j=l$

[TABLE]

Case 2: $i=j\neq k=l$

[TABLE]

Case 3: We note that $\omega_{i}+\omega_{\bar{i}}=0$ where $\bar{i}=L-i+1$. Therefore, if $l=\bar{i}$ and $k=\bar{j}$ the resonance condition is achieved. To avoid double counting with case 1 we write $l=\bar{i}$, $k=\bar{j}$, $i\neq j$, $i\neq\bar{j}$ so

[TABLE]

where we use the fact that $V_{ij}=V_{ji}=V_{\bar{i}j}=V_{j\bar{i}}$. All the other elements are zero.

All the non-zero elements of $\mathcal{D}^{\rm RWA}$ are discusse in Case 1,2,3. Since most of the terms are zero, it is quite easy to find the eigenvalues of $\mathcal{D}^{\rm RWA}$. We call those eigenvalues $|{S}\rangle=\sum_{ij}S_{ij}|{ij}\rangle$. From the cases 1 and 3, one can see that the off-diagonal states where $S_{ii}=0$ are decoupled from the diagonal ones. Therefore we consider these two cases separately. Let $|{S_{o}}\rangle=\sum_{i\neq j}S_{ij}|{ij}\rangle$ be an off-diagonal state, then the eigenvalue equation $\mathcal{D}^{\rm RWA}|{S^{o}}\rangle=\lambda|{S^{o}}\rangle$ written as $\langle{kl}|\mathcal{D}^{\rm RWA}|{S^{o}}\rangle=\lambda S_{kl}$ for $k\neq l$ is

[TABLE]

when $l\neq\bar{k}$ and

[TABLE]

Therefore, for each pair $k,l$ Eq (74) is a $2\times 2$ matrix eigenvalue problem whose minimum (in absolute value) eigenvalue is

[TABLE]

On the other hand

[TABLE]

When $L\gg 1$ we can neglect the $\mathcal{O}(L^{-2})$ correction, and since $V_{ll}$ is minimized for $l=1$ we find that the gap is

[TABLE]

We now show that the other “diagonal” eigenvalues $|{S^{d}}\rangle=\sum_{i}S_{ii}|{ii}\rangle$ have a larger gap. Writing the eigenvalue equation we find $-\frac{\sigma}{2}(2V_{ki}\delta_{ki}S_{kk}-2V_{ki}^{2}S_{ii})=\lambda S_{kk}$, namely we have to find the eigenvalues of the matrix $R_{ik}=\sigma(V_{ik}\delta_{ik}-V_{ik}^{2})$. Calling $V^{d}=\sigma{\rm diag}V$ and $a_{i}=\frac{2\sqrt{\sigma}}{L+1}\sin^{2}k_{i}$ then $R=-V^{d}+a^{T}a$. Using the matrix determinant lemma in the eigenvalue equation we find

[TABLE]

The first term in the above equation gives the solutions $\lambda=-V_{ll}=-\frac{2\sigma}{L+1}\sin^{2}k_{l}$ which have a higher gap. On the other hand, the second term in (77) provides the equation

[TABLE]

where in the last equation we use the identity $\frac{1}{a+b}=\frac{1}{a}-\frac{b}{a(a+b)}$. Since $\sum_{l}\sin^{2}k_{l}=(L+1)/2$ we are left with the equation

[TABLE]

A solution to that equation is clearly $\lambda=0$, namely the steady state. On the other hand all the other solutions must satisfy $\lambda<-\frac{2\sigma}{L+1}\sin^{2}k_{l}$ for some $l$ because otherwise all the elements in the sum are positive and there is clearly no solution. Therefore all the solutions must satisfy $|\lambda|>\frac{2\sigma}{L+1}\sin^{2}k_{1}>{\rm gap}$. This concludes the proof that the gap is given by (76).

Appendix D Strong driving limit

We focus here in the derivation of the effective Liouvillean (33). Let us define then $\mathcal{P}$ as the projector onto the low-energy (eigenvalue zero) subspace of $\mathcal{D}=\frac{1}{2}(n^{\uparrow}_{1}-n^{\downarrow}_{1})(n^{\uparrow}_{1}-n^{\downarrow}_{1})$. This space is generated by all the states such that $n_{1}^{\uparrow}=n_{1}^{\downarrow}$. We set also $\mathcal{Q}=\openone-\mathcal{P}$ and call $\mathcal{H}$ the Hamiltonian part such that $\mathcal{L}_{q}=-i\mathcal{H}-\sigma\mathcal{D}$. We call then also $X_{PP}=\mathcal{P}X\mathcal{P}$, with similar definitions for $X_{PQ}$, $X_{QP}$, $X_{QQ}$. We can therfore write $\mathcal{L}_{q}$ in the block form

[TABLE]

where $\sigma\gg\|\mathcal{H}\|,\|\mathcal{D}\|$ and where we used the fact that $\mathcal{P}\mathcal{D}=\mathcal{D}\mathcal{P}=0$. The low-energy eigenvalues can then be obtained using the determinant identity $\det{\begin{pmatrix}A&B\\ C&D\end{pmatrix}}=\det(D)\det(A-BD^{-1}C)$ – see also Zanardi et al. (2016); Zanardi and Venuti (2014) for a related approach. Indeed, using a first order expansion for $\sigma\to\infty$ it is simple to see that the small eigenvalues are the eigenvalues of the effective operator

[TABLE]

The above effective operator can be obtained also with a (possibly non-unitary) similarity transformation $e^{\mathcal{S}_{\rm D}}$ to decouple the “low-energy” and “high-energy” subspaces. Namely one can find $\mathcal{S}_{\rm D}$ such that

[TABLE]

One finds that (81) is valid up to the first order in $\sigma^{-1}$, with $\mathcal{L}_{q}^{\rm eff}$ given by (80), by choosing

[TABLE]

such that

[TABLE]

where $\mathcal{S}_{2,*}=\mathcal{H}_{PP}\mathcal{H}_{PQ}\mathcal{D}_{QQ}^{-2}-\mathcal{H}_{PQ}\mathcal{D}_{QQ}^{-1}\mathcal{H}_{QQ}\mathcal{D}_{QQ}^{-1}$. Note that $i\mathcal{S}_{1}$ is a Hermitian operator, unlike $i\mathcal{S}_{2}$.

We now obtain the effective operator explicitly. Since $\mathcal{P}$ commutes with all the operators acting on all but the first sites, one realizes that $\mathcal{H}_{PQ}$ and $\mathcal{H}_{PQ}$ are only composed by the projections of $a_{1\updownarrow}^{\dagger}a_{2\updownarrow}$ and their complex conjugate. Moreover,

[TABLE]

where the $|{mn}\rangle$ is a short-hand notation for $\frac{(a_{1\uparrow}^{\dagger})^{m}(a_{1\downarrow}^{\dagger})^{n}}{\sqrt{m!n!}}|{0}\rangle$. Similarly we find

[TABLE]

Since in $\mathcal{H}_{QP}$ the up/down states on the first site differ only for one paritcle it is $\mathcal{D}_{QQ}^{-1}\mathcal{H}_{QP}=2\mathcal{H}_{QP}$. Hence the effective operator is given by $-i\mathcal{H}_{PP}-\frac{2}{\sigma}\mathcal{H}_{PQ}\mathcal{H}_{QP}$. This can be computed from

[TABLE]

and their Hermitian conjugate (all the other terms are zero). Moreover, $n_{1\uparrow}\mathcal{P}=n_{1\downarrow}\mathcal{P}$. We find then

[TABLE]

In order to make further analytical progresses we also use the rotating wave approximation which is consistent with the perturbative treatment (see Appendix C) since $\mathcal{L}_{q}^{\rm eff}=-i\mathcal{H}_{PP}-\frac{2}{\sigma}\mathcal{H}_{PQ}\mathcal{H}_{QP}$ and $2/\sigma$ is small. We note that $\mathcal{H}_{PP}=\sum_{\alpha=2}^{L}(a_{\alpha\uparrow}^{\dagger}a_{\alpha{+}1,\uparrow}-a_{\alpha\downarrow}^{\dagger}a_{\alpha{+}1,\downarrow}+{\rm h.c.})$. The above operator can be diagonalized with a Bogoliobov transformation: defining the operators $\tilde{a}_{k\updownarrow}=\sum_{\alpha=1}^{L-1}\frac{2}{L}\sin^{2}\left(\frac{\pi k\alpha}{L}\right)a_{\alpha{+}1,\updownarrow}$ we find that $\mathcal{H}_{PP}=\sum_{k=1}^{L}2\cos\left(\frac{k\pi}{L}\right)(\tilde{n}_{k\uparrow}-\tilde{n}_{k\downarrow})$. Because of this particular form, the rotating wave approximation in (90) corresponds to expanding the operators $a_{2\updownarrow}$ into the diagonal basis $\tilde{a}_{k\updownarrow}$, neglecting the “oscillating” off-diagonal terms. In other terms, we can write

[TABLE]

where $\hat{\mathcal{L}}_{q}$ is the Hermitian Liouvillean in the rotating wave approximation shown in Eq. (33), where $\mathcal{O}(\hat{\mathcal{L}}_{q})=\mathcal{O}(\sigma^{-1})$, while $\mathcal{L}_{q}^{\rm osc.}$, of order $\mathcal{O}(\sigma^{0})$, is composed by the oscillating terms that are neglected in the long-time limit. In particular, from (68) one finds that RWA holds for $t\gg\mathcal{O}(L^{2})$. This approximation is therefore consistent with the results of section IV, where one finds a Liouvillean gap $\mathcal{O}(L^{-3})$ that provides a lower bound to the convergence time $t>\mathcal{O}(L^{3})$. However, while the eigenvalues depend only on the Hermitian operator $\hat{\mathcal{L}}_{q}$, the eigenvectors depend on the oscillating terms via (70). By mixing (69) with (81) we find then that the eigenvalues with small $\mathcal{O}(\sigma^{-1})$ real part have right eigenvectors given by

[TABLE]

where $|{\Phi_{i}}\rangle$ form an orthonormal basis (depdent on both $H$ and $V$), $\mathcal{S}\approx\mathcal{S}_{\rm RWA}+\mathcal{S}_{\rm D}+[\mathcal{S}_{\rm RWA},\mathcal{S}_{\rm D}]/2=\mathcal{O}(\sigma^{-1})$, but $e^{\mathcal{S}}{}^{\dagger}\neq e^{-\mathcal{S}}$. The corresponding left eigenvectors are then $\langle{\Phi_{i}}|e^{-\mathcal{S}}$.

Appendix E Diagonalization of the Richardson-Gaudin model

We perform explicitly the diagonalization of the Richardson-Gaudin model (36) in the bosonic representation discussed in Section IV.1, where $K_{i}^{-}=\tilde{a}_{i\uparrow}\tilde{a}_{i\downarrow}$, $K_{i}^{+}=(K_{i}^{-})^{\dagger}$ and $K_{i}^{z}=(\tilde{n}_{i\uparrow}+\tilde{n}_{i\downarrow}+1)/2$. We start by defining a trial eigenstate $|{\Omega_{\nu}}\rangle$ with no pairing, namely such that

[TABLE]

These equations force the constraints

[TABLE]

namely there cannot be in the same site both up-particles and down-particles. Moreover, $\nu_{0}\equiv 1/2$ because the model has been obtained by projecting the Liouvillean into the states where $n_{0\uparrow}=n_{0\downarrow}$. The eigenvalue of state $|{\Omega_{\nu}}\rangle$ is thus

[TABLE]

Since there are extra constraints, $q=\sum_{i=0}^{L-1}n_{i\uparrow}=\sum_{i=0}^{L-1}n_{i\downarrow}$, for a given set of allowed “quantum numbers” $v_{k}$ the number $N$ of paired particles satisfies $\sum_{i}(2v_{i}-1)+2N=2q$, namely

[TABLE]

By defining the ansatz

[TABLE]

one sees that

[TABLE]

Moreover,

[TABLE]

We now first consider the $N=1$ case and impose the eigenvalue equation $\hat{\mathcal{L}_{q}}|{\psi}\rangle=\lambda|{\psi}\rangle$ where we define $\lambda=E_{0}-\frac{8}{\sigma}\sum_{\alpha}E_{\alpha}$. The eigenvalue equation becomes then

[TABLE]

From that equation we get the relationship

[TABLE]

namely

[TABLE]

By using the last equation we find

[TABLE]

where $M_{\alpha\beta}=\frac{E_{\alpha}E_{\beta}}{E_{\alpha}-E_{\beta}}\frac{u_{0\beta}}{\nu_{0}}$. Using all the above results the eigenvalue equation becomes

[TABLE]

By evaluating $Z_{\alpha}=0$ one gets the equations

[TABLE]

for $\alpha=1,\dots,N$ where $N$ is given by (97). Clearly, $E_{\alpha}=0$ is a solution, while the solutions different from zero are found by solving the equation

[TABLE]

where we used the fact that $\nu_{0}=1/2$. In conclusion, the eigenvalues of the Liouvillean $\hat{\mathcal{L}_{q}}$ are

[TABLE]

where $n_{k}=2\nu_{k}-1$ and the $E_{\alpha}$ are either zero or the solution of (107). From that expression it is clear that the steady state corresponds to $E_{\alpha}=0$ and $n_{k}=0$. The eigenvalues for larger values of $q$ are given by all the previous solutions with smaller $q$ (this can be seen by adding some $E_{\beta}=0$ for larger values of $N$) together with new solutions due to the larger values of $N$ and the larger set of allowed configurations for $n_{k}$.

Appendix F Solution of the SU(2q)-invariant Gaudin model

We describe here the algebraic approach to general Gaudin models and then apply it to our general fermionic representation introduced in Section IV.3. We fix a basis $h_{\alpha}^{(j)}$ of the Cartan subalgebra acting on the $j$-th copy formed by the diagonal operators $X^{(j)}_{\alpha\alpha}$. A state $|{\Omega_{\nu}}\rangle$ which is a simultaneous eigenvector of all the operators $h_{\alpha}^{(j)}$ is called a weight vector. We write $h^{(j)}_{\alpha}|{\Omega_{\nu}}\rangle=\nu_{j\alpha}|{\Omega_{\nu}}\rangle$ where $\nu^{j}_{\alpha}$ is called weight. On the other hand, in the Cartan-Weyl basis the eigenvalue $\chi$ of the adjoint transformation, namely $[h^{(j)}_{\alpha},e_{\chi}]=\chi_{j\alpha}e_{\chi}$, for a given $e_{\chi}$ in the representation, is called a root. Because of Eq.(44) a root can only have eigenvalue -1,0,1. If one fixes an ordering $c_{\alpha}>c_{\alpha+1}$ and writes $h^{(j)}=\sum_{\alpha}c_{\alpha}h^{(j)}_{\alpha}$ then the eigen-operators of $h^{(j)}$ with positive eigenvalue are called the “raising operators”. They correspond to $X^{(j)}_{\alpha\beta}$ for any $\alpha<\beta$. A highest weight vector is a weight vector $|{\Omega_{\nu}}\rangle$ such that all the other vectors in an irreducible representation can be obtained from $|{\Omega_{\nu}}\rangle$ via some lowering operators. As such, a highest weight vector is annihilated by all the raising operators. We call $\chi^{j}$ the simple roots of the algebra, and we fix an inner product between roots $(\chi^{j},\chi^{k})=\sum_{\alpha}\chi^{j}_{\alpha}\chi^{k}_{\alpha}$, and write $|\chi^{j}|^{2}=(\chi^{j},\chi^{j})$. The matrix $C_{jk}=\frac{2(\chi^{j},\chi^{k})}{|\chi^{j}|^{2}}$ is called the Cartan matrix. We call also $F_{jk}=(C^{-1})_{jk}\frac{2}{|\chi^{k}|^{2}}$. Moreover, we call $z_{0}=0$ and $z_{k}=2g_{k}^{-1}$ for $k>0$.

Thanks to the above definitions, and owing to the results of Refs. Ushveridze (1994); Falceto and Gawȩdzki (1997), we can write the eigenvalues of the Gaudin model (43) as

[TABLE]

where $\mu^{j}_{k}$ are the eigenvalues of the Chevalley operators $H^{(j)}_{k}=\frac{2}{|\chi^{k}|^{2}}\sum_{\alpha}\chi^{k}_{\alpha}h^{(j)}_{\alpha}$, namely $H^{(j)}_{k}|{\Omega_{\nu}}\rangle=\mu^{j}_{k}|{\Omega_{\nu}}\rangle$ and so $\mu^{j}_{k}=\frac{2}{|\chi^{k}|^{2}}\sum_{\alpha}\chi^{k}_{\alpha}\nu^{j}_{k}$, and where the Bethe roots satisfy the equations

[TABLE]

The above expressions for the eigenvalues hold whenever the operators $X^{(j)}_{\alpha\beta}$ define any semi-simple Lie algebra. In the particular case discussed in Section IV.3 those operators define a SU(2q)-invariant Gaudin model, in a specific multi-fermion representation. For SU(2q) the simple roots are $\chi^{j}_{\alpha}=\delta_{\alpha j}-\delta_{\alpha+1,j}$ so $|\chi^{j}|^{2}=2$ and $C_{ij}=2\delta_{ij}-(\delta_{i,j-1}+\delta_{i,j+1})$, where $i,j=1,\dots 2q{-}1$. Therefore, $F_{ij}=\sum_{\ell}\frac{2}{2q}\frac{\sin(\pi i\ell/2q)\sin(\pi j\ell/2q)}{2-2\cos(\ell\pi/2q)}$ and the Chevalley operators are given by $H^{(j)}_{\alpha}=X^{(j)}_{\alpha,\alpha}-X^{(j)}_{\alpha+1,\alpha+1}$. We fix the ordering $\{({\downarrow},1),({\downarrow},2)\dots,({\uparrow},1),({\uparrow},2),\dots\}$ so that

[TABLE]

Because of the above equations, the raising operators are given by $a^{\dagger}_{j,i,\uparrow}a_{j,k,\uparrow}$ with $i>k$, by $a^{\dagger}_{j,i,\downarrow}a_{j,k,\downarrow}$ with $i<k$, and by $a_{j,i,\sigma}a_{j,k,\sigma}$. Therefore, the highest weight vectors may contain in the same mode $j$ either spin-$\uparrow$ particles or spin-$\downarrow$ particles, but not both. The only possible highest weight states are then either $\prod_{i}^{n_{j\uparrow}}a^{\dagger}_{j,i,\uparrow}|{0}\rangle$ or $\prod_{i}^{n_{j\downarrow}}a^{\dagger}_{j,q-i+1,\downarrow}|{0}\rangle$. These states are parametrized by the numbers $n_{j\uparrow}$ and $n_{j\downarrow}$ that satisfy $n_{j\uparrow}n_{j\downarrow}=0$. Therefore

[TABLE]

Moreover, $n_{0\uparrow}=n_{0\downarrow}$ so $\mu^{0}_{j}=\delta_{j,q}$. By explicit calculation for $j\leq q$ one finds $F_{qj}=F_{q,2q-j}=j/2$. Therefore, (108) becomes

[TABLE]

where $\omega_{j,\alpha}$ are the solutions of Eq. (109), namely of Eq.  (47).

Appendix G Explicit mean-field analysis

In this section we perform explicitly the mean-field calculations discussed in section III.2, and we closely follow the notation of that section. We remind that Eq. (28) can be written as

[TABLE]

where the $\lambda$’s are ordered with decreasing (negative) real part, $\lambda_{0}=0$, $\mathring{V}=Z\tilde{V}Z^{-1}$ and we remind that the new bosonic creation operators are obtained via the non-unitary Bogoliubov transformation $\tilde{a}^{\prime}_{i}=\sum_{\alpha}Z_{\alpha i}a_{\alpha}^{\dagger}$, $\tilde{a}_{i}=\sum_{\alpha}(Z^{-1})_{i\alpha}a_{\alpha}$. The steady state is therefore the boson condensate $|{\Omega}\rangle=\frac{(\tilde{a}_{0}^{\prime})^{q}}{\sqrt{q}!}|{0}\rangle$ where $|{0}\rangle$ is the bosonic vacuum. Indeed, clearly this state is annihilated by the quadratic term. To see that even the second one annihilates it is important to remind that $S^{-1}_{0\alpha}$ is the right eigenvector of the steady state (corresponding to the steady state) and the corresponding left eigenvalue $S_{\alpha 0}$ is the identity operator. Therefore, $\tilde{V}_{i0}=\sum_{\alpha\beta}S^{-1}_{j\alpha}\mathring{V}_{\alpha\beta}S_{\beta 0}=0$ since $\sum_{\beta}\mathring{V}_{\alpha\beta}S_{\beta 0}$ is a vectorization of the expression $[V,\openone]$. Similarly $\tilde{V}_{0i}=0$. To study the elementary excitations with respect to this state, one can use Bogoliubov (mean field) approach starting from the variational states $|{\psi}\rangle=\sum_{j}\psi_{j}\frac{(\tilde{a}_{0}^{\prime})^{q-1}}{\sqrt{(q-1)!}}a^{\prime}_{j}|{0}\rangle$, for $j\neq 0$ and the corresponding $\langle{\psi^{\prime}}|=\sum_{j}\psi^{\prime}_{j}\langle{0}|\frac{(\tilde{a}_{0})^{q-1}}{\sqrt{(q-1)!}}a_{j}$, where $\sum_{j}\psi^{\prime}_{j}\psi_{j}=1$. The variational Liouvillean then becomes

[TABLE]

which, similarly to the Rayleigh-Ritz method, has to satisfy $\partial_{\psi}\mathcal{L}^{V}=\partial_{\psi^{\prime}}\mathcal{L}^{V}=0$ with the constraint $\sum_{j}\psi_{j}^{\prime}\psi_{j}=1$ (see e.g. Laestadius and Kvaal (2017)). However, because $\tilde{V}_{i0}=\tilde{V}_{0i}=0$ for all $i$, one can restrict the sum in the above equation to the values $i,j,k,l>0$, but because there is only one particle in $|{\psi}\rangle$ in the states $i>0$ one finds that

[TABLE]

namely that in the single-excitation subspace the variational Liouvillean is already diagonal. This shows that the eigenvalues, at least in the low-energy subspace, are not “renormalized” for larger values of $q$.