Driven quantum dynamics: will it blend?

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.


I. INTRODUCTION
Randomness generating quantum operations play a central role in our understanding of very various physical phenomena [1]. 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 [2], quantum process tomography [3], fidelity estimation [4], quantum communication and entanglement sharing [5][6][7], quantum data-hiding [2,8,9] and entanglement generation [10][11][12][13]. Because of their crucial importance, several procedures have been developed to generate either truly random or pseudo-random operations via random quantum circuits [4,[14][15][16][17][18][19][20]. 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 [21,22], where different stochastic control pulses can be simulated by changing the spatial extent of the waveguides via electrically tuned on-chip heaters [23]. In those systems a major recent result has been the experimental measurement of boson sampling [24][25][26], 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 [27]. 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 [28,29], 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 [30] 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Ĥ 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 × d Hermitian matrices while g(t) is a scalar function depending on time t. For multiple pulsesĤ(t) = H + 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 = T exp −i T 0Ĥ (s) ds , where 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 ∈ SU(d) it is possible to find a control profile g(t) such that U = T exp −i T 0Ĥ (s) ds 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 [31]  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 [32]. 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 where the amplitudes A k , the phases ϕ k , and possibly even the frequencies ω k are random variables. We use the notation E[·] 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.
A. 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 [4,[14][15][16][17][18], 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 [33,34]. 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 timeindependent 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 [35]. 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. [36]). 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 for suitably small , where · refers to the diamond norm, E U denotes an average over some given distribution of unitaries µ U and µ 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 . A related notion [18] is that of quantum expanders, which are defined by where X ⊗q,q = X ⊗q ⊗(X ⊗q ) * . Eq. (4) can be regarded as the vertorised version of Eq. (3): given an operator X = i j X i j |i j|, its vectorized form is |X = i j X i j |i j . 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 [18]. It follows from the definition that |AX = A ⊗ 1 1|X and |XA = 1 1 ⊗ A T |X . Therefore, X ⊗q,q is the vectorization of the superoperator ρ → X ⊗q ρX ⊗q † . 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 [14], e(µ, q) ≤ 2q W(µ, µ Haar ), for all measures µ, being W the Wasserstein distance [37] , 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, 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 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 1, in view of the central limit theorem. (ii) We assume also that g(t) is harmonic, namely that 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 E[g(t)] = 0. In view of these assumptions, exploiting the results of [38,39], 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 T c(T s) σ 2 δ(s) where δ is the Dirac delta function and σ is a constant. In the long-time limit, t > T H , V , one finds then that where and X ⊕q =X⊕X⊕ . . . , being ⊕ the Kronecker sum X⊕Y=X⊗1 1+1 1⊗Y. Therefore, with these three approximations, the long-time dynamics of the stochastic process is Markovian and described by the above Lindblad equation [40,41], where the operator L q is called Liouvillean. Similarly to what happens with the replica trick in statistical physics [42], the average over the noise effectively couples the initially uncoupled copies. Sometimes we will use the more convenient vectorised form of the above equation whereX = X ⊗ 1 1 − 1 1 ⊗ X T is the vectorization of the commutator [X, ·]. If t → ∞ then E U U ⊗q ρU ⊗q † converges to one of the steady states of the Liouvillean L q . In the following section we prove that the steady state manifold of 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 → ∞. 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.

B. Steady state of the Liouvillean evolution
We start by describing the steady state of 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 [43]. Given an orthonormal basis {M µ } of the steady state space, equipped with the standard Hilbert-Schmidt product, there exists a dual operator set {J µ } such that L q † J µ = 0, where L q † is the Liouvillean operator (7) after the substitution H → −H. The latter substitution does not change the dynamical algebra, so algebraic considerations based on controllability hold also for L q † . From the conserved quantities J µ and their dual operators M µ one finds the steady state as ρ ∞ = µ M µ Tr(J µ ρ 0 ) where ρ 0 is the initial state [43]. Since the system is controllable, repeated commutators of H ⊕q and V ⊕q give rise to the algebra su(d) ⊕q . Therefore, because of the Schur-Weyl duality [44], the only operators that commute with both H ⊕q and V ⊕q , and more generally with Eq. (5), are index permutation operators. Let S q be the group of permutations of the set 1, . . . , q and let P σ , σ ∈ S q be the operator which permutes the index of the tensor copy H ⊗q , namely the operator that maps ψ i 1 ,i 2 ,··· ,i n to ψ σ(i 1 ),σ(i 2 ),··· ,σ(i n ) for each set of indices i j . It is simple to show that P π P σ =P πσ 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, L q (P ρ ) = L q † (P ρ ) = 0, so ρ ∞ = σ ρ σ P σ . However, since the operators P σ are not orthonormal, one has where in the first equality holds because P σ is a conserved quantity. By inverting the above equation we find that where M σπ = Tr[P † σ P π ]. It has been shown in Ref. [45] that M σπ = d l(σ −1 π) where l(σ) is the number of cycles in the cycle decomposition of σ. The dimensionality of the steady state manifold is then given by the matrix rank of M. One finds that the steady state degeneracy is ∼ e O(q) . The right-hand side of (10) is exactly equal to the integration over the Haar measure (see e.g. Proposition 3 in [45]). Therefore, we have shown that lim t→∞ e tL q ρ = dU U ⊗q ρU ⊗q † , 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 a uniform integration over the Haar measure.

C. 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 L q from Eq. (8) into local terms L loc k acting only on the k-th copy, and a non-local interaction. Indeed, whereH k ,V k , and accordingly L loc k , act only on the k-th copy. Therefore each L loc k for different k is equivalent to a singlecopy Liouvillean L 1 . We assume that the operator L 1 is diagonalizable (with right and left eigenvectors) and call its eivenvalue decomposition, where the eigenvalues λ j are ordered with decreasing real part (starting from zero) and Π j are the corresponding eigenprojections. The operators are then eigenprojections of L q , with eigenvalue λ i . To show this, we note indeed that Π (i) j is proportional to the vectorization of the identity operator in each copy, aside from the j-th one, since Π (0) is the projection onto the steady state and, accordingly, Π (0) (X) = ρ ∞ Tr[X], which is proportional to the identity operator. Therefore,V l Π (i) j (X)] = 0 for all X), as long as l j. On the other hand, for l = j, it isV kVl Π (i) j = 0, since by construction k j. This shows that (15) is a projector on the eigenspace of L q with eigenvalue λ i . Moreover, from the operators (15) one can also construct the eigenstates of L q that act on the irreducible representations of the symmetric group -indeed since the permutation operators P σ commute with the Liouvillean, then P σ (Π (i) j )P † σ is an eigenprojection of L q for all σ. In summary, the eigenstate of L 1 with the lowest gap can be used to construct some exact eigenstates of 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 ρ Π (i) j is also an eigenvector with eigenvalue λ j of L q .

D. 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 → 2 norm: where U ∞ = lim t→∞ e tL q and d is the dimension of the local Hilbert space. Let s j (t) be the singular values of e tL , ordered from largest to smallest. The largest has magnitude one. Then the singular values of (e tL − U ∞ ) are strictly smaller than one, and If the Liouvillian were reversible, then the singular values s j (t) would be given by e tλ j , where λ j are the eigenvalues of L.
Unfortunately the semigroups that we will be working with are not Hermitian. Nonetheless, from Eq. (B5), we find that the 2 → 2 norm can be bounded in terms of the eigenvalues and eigenvectors of L q as where λ j are the eigenvalues of L q , and R j , L j are its right and left eigenvectors, satisfying tr[L † j R k ] = δ 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 (σ = → 0) and strong (σ = −1 → ∞) coupling limits, and show the following properties: (i) the spectral gap is O( ), both in the strong and weak coupling limits -for strong driving, the decrease of the gap for larger σ is consistent with the general occurrence in open systems [46]; (ii) the eigenvectors satisfy |R j = S|Φ j and |L j = S †,−1 |Φ j , for some invertible matrix S and an orthonormal basis |Φ j . The condition number of S is κ(S) ≡ ||S|| ||S −1 || and satisfies κ(S) = O(1 + ). Moreover, in sections III B 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 [34], at the expense of more stringent requirements on the fluctuating terms of the Hamiltonian.
We then get that where λ * is the eigenvalue with the smallest non-zero real part and κ(S) = O(1 + ). In terms of the trace norm, we then get that In the weak or strong coupling limits, the condition number will be of order one yielding a mixing time of T * ∼ 4q log(d)/λ * . 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 * ∼ log(n)T * 1 , where T * 1 is the mixing time of a single subsystem [47]. We might ask whether the mixing time of Eq. (7) is also of the order T * ∼ log(q)T * 1 , with T * We can see that this is not the case from the following argument: since the lower bound is saturated when S = 1 1, and we have isolated the subspace with eigenvalue λ * . Now, in Section II C we have argued that if the gap of L q is the same as the gap of L 1 , then we can construct the eigenvectors with minimal nonzero eigenvalue of L q from those of 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 ≥ e O(q) . Hence the first excited subspace does as well. Then, Thus the mixing time is at least T * ∼ O(q/λ * ), even in the weak coupling limit. Finally, we comment on the distinction between the singular value gap of e tL and the eigenvalue gap of L. We know that as t → ∞, the singular value gap s * (t), namely the largest singular value s j (t) 1, converges to e tλ * , 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/λ * .

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 [15] 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 [48]. 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.
A. 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 αβ = |α β|, α, β = 1, . . . , d and call Hence, the Liouvillean can be written as The form (24) is a convenient starting point because it depends only on the original d×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 (C d ) ⊗q Hilbert space. Indeed, because of the Schur-Weyl duality, every irreducible representation of (C d ) ⊗q is decomposed as (C d ) ⊗q = ⊗ λ P λ ⊗ U λ where P λ is an irreducible representation of the symmetric group S q and U λ an irreducible representation of SU(d). A convenient expression for the fullysymmetric and fully-anti-symmetric subspaces is given by [49] where a α and a † α 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 αβ = u a † αu a βu . From the definition of B one realizes that in this generic representation there are exactly q particles since For convenience, we also perform the calculation in the basis where V is diagonal. Therefore, Eq.(24) becomes where n x = a † x 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, αu n αu↑ = αu n αu↓ = q 1 1.

B. Mean-field approach
We consider here the decomposition (12) where each L loc k for different k is equivalent to a single-copy Liouvillean L 1 . From the above decomposition it is clear that if the gap of L loc = k L loc k equals the gap of L q then the Liouvillean gap λ * is independent on q.
Extending the treatment of Section III A, we define a local basis of operatorsBαβ , and similarly forβ, are multi-indices running from 1 to d 2 . Therefore we can write the decomposition Eq. (12) as and, writingBαβ = a † α aβ with bosonic operators, then We assume that L loc q is diagonalizable (with left and right eigenvectors) as (L loc q )αβ = j Zα j λ j Z −1 jβ 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ã i = α Zα i a † α ,ã i = α (Z −1 ) iα aα. These operators still satisfy the canonical commutation relations [ã i ,ã j ] =δ i j , thoughã i ã † i . As shown in Appendix G, in this language, the steady state of the many-body Liouvillean (28) is therefore the boson "condensate" |Ω = where |0 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 |ψ = j ψ j , for j 0 and optimising over the amplitudes ψ 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 quasiparticle excitation. Although mean-field techniques have been highly studied mostly for Hamiltonian systems [48], they can be extended also to non-normal operators [50] where left and right eigenvectors form a bi-orthonormal basis. Within this variational formalism we show in Appendix G that the fourbody interaction in (28) does not alter the eigenstates, which are therefore exactly given by the bare single-particle eigenstates |Ω exc.
with exact eigenvalue λ 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 |Ω exc.
j are indeed the symmetric combination of (15), which, as shown before, are an exact eigenstate of L q . Within this simple mean-field treatment there are no other eigenvalues with a smaller gap than min j | [λ 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.
C. 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 ≥ 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 nontrivial 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 operatorX, different from a permutation operator, such that [H ⊕p ,X] = [V ⊕p ,X] = 0 and that the solutions of [H ⊕q , X] = [V ⊕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 symmetryX introduces an extra steady state. Then, suppose that we restore full-controllability by adding a small O( ) term in either H or V such that the operatorX is not a symmetry anymore (we say that the symmetryX is explicitly broken). This splits the extra steady state into an eigenvector with small O( ) eigenvalue which, for small enough 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. [51,52] a rather surprising necessary and sufficient condition for controllability is that there are exactly two independent solutions of the equations [H ⊕2 , X] = [V ⊕2 , X] = 0. Nontheless, a simpler necessary condition (though not sufficient [51]) is the absense of non-zero solutions to the set of equations 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. [53] proved that QQ * = 1 when Q is symmetric and QQ * = −1 when Q is anti-symmetric. If there are nonzero solutions of (29), then the system is not controllable and there are extra steady states such as the bosonic paired state for q = 2 Indeed, for both Q symmetric and anti-symmetric Q ⊗ Q * is symmetric, thus justifying the bosonic approach. The proof can be readily obtained from (28), indeed for both X = H, V so because of Eq.(29) we findHγ˜δa † γ a˜δ|ψ Q =Vγ˜δa † γ a˜δ|ψ Q = 0, namely L 2 |ψ Q = 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 ⊗ Q * ]αβ QαQ * β since Q is a matrix. As discussed before, we can restore controllability by explicitly breaking the symmetry (29) with small terms: QH T + HQ = H , QV T +V Q = V where at least one between V or H has to be non-zero, otherwise the system is not controllable. In this case |ψ Q is not a steady state but, within first order perturbation theory, can be used to create a state with eigenvaluẽ δ = O( V , H ). In particular, one can construct specific examples where V and H are much smaller than the gap λ * of L 1 so thatδ < λ * . 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 = (σ x 1 σ x 2 + σ y 1 σ y 2 + σ x 1 ) + σ z 1 σ z 2 and V = σ y 1 , where σ α j are the Pauli matrices acting on the spin j. For instance, for = 0.1 the gap of L 1 is ≈ 0.45 while the gap of L 2 is ≈ 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. [42]) 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 where |n 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 [54][55][56][57]. 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 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 [58,59]. 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 [21,22], where the different control pulses can be obtained by electrically tuned on-chip heaters [23].
In the following we evaluate the Liouvillean gap for all possible values of q in the strong driving limit, namely when σ 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 αβ = a † α a β for either bosonic or fermionic degrees of freedom. We then extend our analysis to the general case.

A. Gap analysis: fully-symmetric representation
We consider first the fully symmetric representation where B αβ = a † α a β so one can omit the index u from the equations of Section III A. Plugging the operators H and V of the control-lable chain into Eq. (27) one finds the following Liouvillean To diagonalize the above operator we assume that σ 1 and we study the "low-energy" effective dynamics. In that limit the dissipative part σD = σ 2 (n ↑ 1 − n ↓ 1 )(n ↑ 1 − n ↓ 1 ) has either eigenvalue 0 or σ 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 σ −1 the effective Liouvillean is given bŷ where g k = 2 L sin 2 πk L ,ã k = L−1 α=1 2 L sin 2 πkα L a α+1, and a 1 ≡ã 0 . We call now and K z i = (ñ i↑ +ñ i↓ + 1)/2 and note that these operators satisfy the SU(1,1) commutation relations With these definitions we find then where K i · K j ≡ −(K + i K − j + K − i K + j )/2 + K z i K z j is the SU(1,1) invariant product, namely the analogue of the Heisenberg interaction. The model (36) is a SU(1,1) Gaudin model [60], 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 [61]. We find that the eigenvalues of the LiouvilleanL q are 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 α are either zero or the solution of the non-linear set of equations where E α = 1/ω α . From that expression it is clear that the steady state corresponds to E α = 0 and n k = 0, for each α and k. Solutions to the above equations are known to be related with the roots of Heine-Stieltjes polynomials (see e.g. [62]). By exploiting this relationship, one finds that all the solutions ω α 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 −1 k+1 < ω α < 2g −1 k for some k, so that 2E α > min k g k = g 1 . This constraint allows us to find rigorously the gap of the LiouvilleanL 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 α = 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 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.

B. 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 α i j = a † iα a jα with fermionic creation and annihilation operators. Repeating the effective Liouvillean description of the previous section we find where S 0 · S k = α=x,y,z S α 0 S α k refers to the SU(2)-invariant, product, namely the spin Heisenberg interaction, S ± j = S x j ± iS y j , and where we have defined S − j =ã j↑ã j↓ , S + j = (S − j ) † and S z j = (ã † j↑ã j↑ +ã † j↓ã j↓ − 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 [60]. 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 α are the solutions of k g k (n k − 1) 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.

C. General gap analysis
As we have discussed in Section III A, a general representation of the SU(L) algebra can can be obtained via extended creation and annihilation operators [49], namely B αβ = uã † αuãβ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 nonunitary operation in bosonic systems (see e.g. [48]). Because of the Pauli exclusion principle, in order to satisfy the constrain α B αα = 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 σ 1, one finds that the effective LiouvilleanL q can be written in the diagonal basis of the Hamiltonian aŝ The above Hermitian operator corresponds to the purely dissipative Liouvillean whereṼ k = |ω k ω 0 |, being |ω k = L−1 j=1 2 L sin π jk L 2 | j+1 and |ω 0 = |1 . One can check that the operatorsṼ 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 ↑ and ↓ anti-commute. The secondone is a particle-hole transformation in the spin-down sector. These transformations are implemented together by defining W = ju e iã † ju↑ã ju↑ and setting a ju↑ =ã ju↑ and a ju↓ = Wã † ju↑ . Eq. (42) then becomeŝ where X ( j) αβ = a † jα a jβ − a jβ a † jα /2 and the Greek letters refer to the multi-index composed by the auxiliary index and the "effective spin" index, i.e. α = (us) where u = 1, . . . , q and s = {↑, ↓}. The traceless operators X ( j) αβ satisfy the SU(2q) ⊕L commutation relations, so that Eq.(43) represents a SU(2q) version of the Gaudin model. Indeed, Eq. (43) is invariant under the Bogoliubov transofmation a jα → β U α,β a jβ , where U is a unitary (2q) × (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 [ α X ( j) αα , X (k) βγ ] = 0 for each β and γ. Going back to the original representation, namely performing back the particle-hole transformation, one finds that 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. [63,64], 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. [65] and references therein), especially due to its connections with the Knizhnik-Zamolodchikov equation [65,66]. In Appendix F we exploit the general solution [63,64] 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 ofL q are parametrized by non-negative integers n ↑ j and n ↓ j , and are given by where ω j,α for j = 1, . . . , 2q − 1 are the solutions of being z 0 = 0, µ 0 j = δ q j , and, for k > 0, z k = 2g −1 k and µ k j = δ j,q (1 − δ n ↓k >0 − δ n ↑k >0 ) + δ j,q+n ↑k + δ j,q−n ↓k . In (47) we set ω 0,β = ω 2q,β → −∞ 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ω −1 q,α > g k for each α and k, then we can straightforwardly (49) with three different values of z k and µ k i = 1.
apply the reasoning of Section IV A 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 A. Here we use a different approach, used also in Ref. [63] 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. [63] we define the two-dimensional vector ω jα whose real components are the real and imaginary part of ω jα and interpret those vectors as the positions of some particles with index α and species j = 1, . . . , 2q − 1. The equations (47) can then be interpreted as the conditions for an extremum of the function W({ω}) defined as where z k = (z k , 0) and the Cartan matrix C i j has non-zero components only on the diagonal, where C ii = 2, and for |i − j| = 1, where C i j = −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 [67]. At first sight one may think that the problem has no solutions since the potential (49)  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 ω < min k z k and for ω > max k z k , where no solutions can exist. Therefore, this electrostatic analogy shows that the only stable solutions with finite ω iα can be found only between poles of V i (ω), or, in other terms, that the solutions of the non-linear set of equations (47) satisfy the constraint min k z k < ω jα < max k z k , i.e. 2ω −1 jα > min k g k . This, together with the discussion of Section IV A, shows that Eq.(39) is indeed the gap of the Liouvil-leanL q in the strong-driving limit.

D. 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 ∝ 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 [68]. The scaling ∝ 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 σ. 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 σ. This scaling has been predicted in the strong and weak driving limits by Eqs. (39) and (C12). However, Fig. 4  shows that such scaling is valid also for σ ≈ 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 tL q which is a good estimate of the convergence time (see section II D). As expected, both in the strong and weak coupling limit the s * (t) converges to e −λ * t much earlier than mixing time-scales. Therefore, in these regimes, one finds that the convergence time is basically 1/λ * . On the other hand, for σ = 2 the matching between e −λ * t and s * (t) only happens at longer times. Therefore, as expected from the analysis of Section II D, 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 σ. 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. [69]. Using a simple reparametrization of these angles one can write the Haar measure as namely as a uniform distriution of the angles ϕ j in the range [0, 2π]. Therefore, testing whether the resulting distribution approximates a Haar measure is equivalent to testing whether the angles ϕ j are distributed as a multinomial uniform distribution. In Fig. 6 we do a simple test to verify the distribution of the angles ϕ j : we divide the interval [0, 2π] into 25 bins and plot, as a 3D histrogram, the matrix whose elements (i, j) are the number of times that the angle ϕ 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 ≈ 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.

A. 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 |per(Ũ)| 2 , being per(Ũ) the matrix permanent of the q × q matrixŨ, whereŨ is built from some columns and rows of a L × L Haar-uniform matrix U [24][25][26]. Therefore where σ, σ are permutations in the symmetric group S q , K b.s. is a suitable index contraction operator and U ⊗q,q = U ⊗q ⊗ (U ⊗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 × L one particle evolution matrix from time 0 to time t and a j (t) = k U jk a k (0), then because of the Wick's theorem where K m.p depends on the initial two-point correlation functions a † i (0)a j (0) . Expressions like (53) arise also in XY spin chains, which can be mapped to a quasi-free fermionic model via the Jordan-Wigner transformation [70]. For instance, the driven XY model can be mapped, in the single-particle subspace, to the driven quantum walk of Section IV. Calling U the resulting singleparticle evolution, then in any subspace long-range spin operators S α i S β j , for α, β ∈ {x, y} can be written as a combination of fermion strings as in (53) where q = |i− j| for i j. Therefore, with a suitable K XY that depends on the initial correlations, one can write the dynamical long-range correlations between spin operators in an XY chain as for α, β ∈ {x, y}. Similarly, S z i (t)S z j (t) = Tr U ⊗2,2 K zz XY . 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 K where we used (4). Thanks to the analysis of Section II D, 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.

B. 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 * 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 τT * this error probability exponentially decreases as a function of τ. 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 λ * , which depends on σ (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 * min = min σ T * (σ) T * (σ 2.5) ≈ 0.055 L 3 .
In order to estimate T c we perform a numerical experiment with the QuTip quantum control package [71]. We consider the model (31) and, for each length L = 10, . . . , 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 ≈ 0.069 L 3 . This shows two remarkable facts: (i) the values of T c and T * min are very close for L = 10, . . . , 20; (ii) both T c and T * 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 * 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 dynam-ical 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 [35].
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 [21,22], where random unitary operations have been used in the first smallscale experimental observations of boson sampling [24][25][26]. 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, . . . , 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 manybody physics, where the interplay between quantum manybody effects and noise is currently a subject of intensive study in many area, such as spin glass [42], the fast scrambling of quantum information [28,29], and many-body localization [72,73]. 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 [68,[74][75][76][77], produce large amount of entanglement [78], and give new insights into the process of thermalization and equilibration [79]. Haar-random quantum states are known to have, typically, an extensive amount of entanglement [80]. 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.
To simplify the theoretical description, in this section we consider only q = 1 and call E t the quantum channel resulting from the average evolution of the quantum system 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 E[g(t+s)g(t)] = c(s) is independent on t and E[g(t)] = 0. In view of these assumptions, we can simplify (A1) 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 [38,39] 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 T c(T s) σ 2 δ(s) where δ is the Dirac delta function ans σ is a constant. In the long-time limit one finds that 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.
where T ∞ = lim t→∞ e tL , and (t) is the distinguishability error. A less stringent convergence requirement is to ask whether e tL q is an expander for a given value of t. Then, we want to estimate where a hat indicates that the CPTP maps are represented as channels (see Ref. [81] for more details on the representation of channels). Trace norm convergence and 'spectral convergence' are related, by noting that where d is the dimension of the Hilbert space, and recalling that ||e tL − T ∞ || 1→1 = sup ρ ||e tL (ρ) − T ∞ (ρ)|| 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 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: L(ρ) = 0. The rest of the spectrum comes in complex conjugate pairs. The Liouvillian is called unital if it annihilates the identity L(1 1) = 0. The Liouvillian in Eq. (7) has this property. A unital Liouvillian is called reversible ifL =L † , 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 tL rather than its eigenvalues. The singular spectrum of e tL is equal to the spectrum of e tL e tL † . It is not difficult to see that the 2 → 2 norm is related to the singular spectrum. Let s j (t) be the singular values of e tL , ordered from largest to smallest. The largest has magnitude one. We know that asymptotically s * (t) = e tRe[λ * ] , where now λ j are the eigenvalues of L written in decreasing (real part) order, and λ * is the gap of L; i.e. the smallest (in magnitude) non-zero real part of any eigenvalue of L. To see this, note that, assuming it has no Jordan blocks, the Liouvillian can be written in its spectral decomposition as where R j , L j are a bi-orthonormal basis of operators: i.e. tr[L † j R k ] = δ 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, Hence, for very large t, the convergence is governed by the gap, and s * (t) → e tλ * . In principle we do not know at what scale e −tλ * | R j |R j | | ψ|L j L j |ψ |. 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.
Case 3: We note that ω i +ω¯i = 0 whereī = L−i+1. Therefore, if l =ī and k =j the resonance condition is achieved. To avoid double counting with case 1 we write l =ī, k =j, i j, i j so where we use the fact that V i j = V ji = V¯i j = V jī . All the other elements are zero. All the non-zero elements of D RWA are discusse in Case 1,2,3. Since most of the terms are zero, it is quite easy to find the eigenvalues of D RWA . We call those eigenvalues |S = i j S i j |i j . 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 = i j S i j |i j be an off-diagonal state, then the eigenvalue equation D RWA |S o = λ|S o written as kl|D RWA |S o = λS kl for k l is when l k and Therefore, for each pair k, l Eq (C10) is a 2 × 2 matrix eigenvalue problem whose minimum (in absolute value) eigenvalue is On the other hand When L 1 we can neglect the O(L −2 ) correction, and since V ll is minimized for l = 1 we find that the gap is We now show that the other "diagonal" eigenvalues |S d = i S ii |ii have a larger gap. Writing the eigenvalue equation we find − σ 2 (2V ki δ ki S kk − 2V 2 ki S ii ) = λS kk , namely we have to find the eigenvalues of the matrix R ik = σ(V ik δ ik −V 2 ik ). Calling V d = σdiagV and a i = 2 √ σ L+1 sin 2 k i then R = −V d + a T a. Using the matrix determinant lemma in the eigenvalue equation we find The first term in the above equation gives the solutions λ = −V ll = − 2σ L+1 sin 2 k l which have a higher gap. On the other hand, the second term in (C13) provides the equation where in the last equation we use the identity 1 a+b = 1 a − b a(a+b) . Since l sin 2 k l = (L + 1)/2 we are left with the equation (C14) A solution to that equation is clearly λ = 0, namely the steady state. On the other hand all the other solutions must satisfy λ < − 2σ 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 |λ| > 2σ L+1 sin 2 k 1 > gap. This concludes the proof that the gap is given by (C12).

Appendix D: Strong driving limit
We focus here in the derivation of the effective Liouvillean (33). Let us define then P as the projector onto the low-energy (eigenvalue zero) subspace of D = 1 2 (n ↑ 1 − n ↓ 1 )(n ↑ 1 − n ↓ 1 ). This space is generated by all the states such that n ↑ 1 = n ↓ 1 . We set also Q = 1 1 − P and call H the Hamiltonian part such that L q = −iH − σD. We call then also X PP = PXP, with similar definitions for X PQ , X QP , X QQ . We can therfore write L q in the block form where σ H , D and where we used the fact that PD = DP = 0. The low-energy eigenvalues can then be obtained using the determinant identity det A B C D = det(D) det(A − BD −1 C) -see also [46,83] for a related approach. Indeed, using a first order expansion for σ → ∞ it is simple to see that the small eigenvalues are the eigenvalues of the effective operator The above effective operator can be obtained also with a (possibly non-unitary) similarity transformation e S D to decouple the "low-energy" and "high-energy" subspaces. Namely one can find S D such that One finds that (D3) is valid up to the first order in σ −1 , with L eff q given by (D2), by choosing such that where S 2, * = H PP H PQ D −2 QQ − H PQ D −1 QQ H QQ D −1 QQ . Note that iS 1 is a Hermitian operator, unlike iS 2 .
We now obtain the effective operator explicitly. Since P commutes with all the operators acting on all but the first sites, one realizes that H PQ and H PQ are only composed by the projections of a † 1 a 2 and their complex conjugate. Moreover, where the |mn is a short-hand notation for (a † 1↑ ) m (a † 1↓ ) n √ m!n! |0 . Similarly we find Pa 1↑ Q = n 1 √ n 1 |n 1 , n 1 n 1 +1, n 1 | , Pa 1↓ Q = n 1 n 1 + 1|n 1 , n 1 n 1 , n 1 +1| .
Since in H QP the up/down states on the first site differ only for one paritcle it is D −1 QQ H QP = 2H QP . Hence the effective operator is given by −iH PP − 2 σ H PQ H QP . This can be computed from Pa † 1 Qa 1 P = n 1 P (D8) Pa 1 Qa † 1 P = (n 1 + 1)P (D9) and their Hermitian conjugate (all the other terms are zero). Moreover, n 1↑ P = n 1↓ P. We find then H PQ H QP = − 2(a † 1↑ a † 1↓ a 2↑ a 2↓ + h.c.) + n 1↑ (n 2↓ + 1) (D11) + n 1↓ (n 2↑ + 1) + n 2↓ (n 1↑ + 1) + n 2↑ (n 1↓ + 1) = − 2(a † 1↑ a † 1↓ a 2↑ a 2↓ + h.c.) − 1+ + (n 1↑ + n 1↓ + 1)(n 2↑ + n 2↓ + 1) (D12) 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 L eff q = −iH PP − 2 σ H PQ H QP and 2/σ is small. We note that H PP = L α=2 (a † α↑ a α+1,↑ − a † α↓ a α+1,↓ + h.c.). The above operator can be diagonalized with a Bogoliobov transformation: defining the operatorsã k = L−1 α=1 2 L sin 2 πkα L a α+1, we find that H PP = L k=1 2 cos kπ L (ñ k↑ −ñ k↓ ). Because of this particular form, the rotating wave approximation in (D12) corresponds to expanding the operators a 2 into the diagonal basisã k , neglecting the "oscillating" off-diagonal terms. In other terms, we can write whereL q is the Hermitian Liouvillean in the rotating wave approximation shown in Eq. (33), where O(L q ) = O(σ −1 ), while L osc. q , of order O(σ 0 ), is composed by the oscillating terms that are neglected in the long-time limit. In particular, from (C4) one finds that RWA holds for t O(L 2 ). This approximation is therefore consistent with the results of section IV, where one finds a Liouvillean gap O(L −3 ) that provides a lower bound to the convergence time t > O(L 3 ). However, while the eigenvalues depend only on the Hermitian operator L q , the eigenvectors depend on the oscillating terms via (C6). By mixing (C5) with (D3) we find then that the eigenvalues with small O(σ −1 ) real part have right eigenvectors given by We perform explicitly the diagonalization of the Richardson-Gaudin model (36) in the bosonic representation discussed in Section IV A, where K − i =ã i↑ãi↓ , K + i = (K − i ) † and K z i = (ñ i↑ +ñ i↓ + 1)/2. We start by defining a trial eigenstate |Ω ν with no pairing, namely such that These equations force the constraints ν i = (n i↑ + n i↓ + 1)/2 , n i↑ n i↓ = 0 , namely there cannot be in the same site both up-particles and down-particles. Moreover, ν 0 ≡ 1/2 because the model has been obtained by projecting the Liouvillean into the states where n 0↑ = n 0↓ . The eigenvalue of state |Ω ν is thuŝ Since there are extra constraints, q = L−1 i=0 n i↑ = L−1 i=0 n i↓ , for a given set of allowed "quantum numbers" v k the number N of paired particles satisfies i (2v i − 1) + 2N = 2q, namely By defining the ansatz one sees that Moreover, We now first consider the N = 1 case and impose the eigenvalue equationL q |ψ = λ|ψ where we define λ = E 0 − namely that in the single-excitation subspace the variational Liouvillean is already diagonal. This shows that the eigenval-ues, at least in the low-energy subspace, are not "renormalized" for larger values of q.