Markovianization with approximate unitary designs

Memoryless processes are ubiquitous in nature, in contrast with the mathematics of open systems theory, which states that non-Markovian processes should be the norm. This discrepancy is usually addressed by subjectively making the environment forgetful. Here we prove that there are physical non-Markovian processes that with high probability look highly Markovian for all orders of correlations; we call this phenomenon Markovianization. Formally, we show that when a quantum process has dynamics given by an approximate unitary design, a large deviation bound on the size of non-Markovian memory is implied. We exemplify our result employing an efficient construction of an approximate unitary circuit design using two-qubit interactions only, showing how seemingly simple systems can speedily become forgetful. Conversely, since the process is closed, it should be possible to detect the underlying non-Markovian effects. However, for these processes, observing non-Markovian signatures would require highly entangling resources and hence be a difficult task. A question of foundational importance is ‘why is nature forgetful?’, playing an important role in our understanding of thermodynamics. Here, the authors study a class of quantum processes, called approximate unitary designs, to show that these processes are highly forgetful - i.e. Markovian.

A foundational question of modern physics is to understand the origins of irreversibility 1 . In particular, to determine whether fundamental laws, which are fully reversible, are consistent with phenomena like equilibration and thermalization. The dynamical version of this conundrum concerns the emergence of forgetful processes from isolated ones. In quantum mechanics, an isolated process is unitary, and cannot lose information; past behavior in one part of the system will always be remembered, eventually returning to influence the future.
However, there are many ways in which nature manifests forgetful processes, where a system's evolution is determined with a seeming disregard to its previous interactions with its surroundings. For example, a carbon atom does not typically remember its past and behaves like any other carbon atom. Such processes are not isolated, and the general intuition is that the dynamics of a system, in contact with a large environment, can be approximately described as memoryless 2 . Yet, formal derivations of memory-less quantum processes require several assumptions about the coupling strength with the environment, the timescales of dynamical correlations, and an infinitedimensional reservoir. For finite-sized environments, this can only be achieved exactly by continually refreshing (discarding and replacing) the environment's state, i.e., artificially throwing away information from the environment. The problem this poses is akin to the one made by the Fundamental Postulate of Statistical Mechanics 1 , which a-priori sets the probabilities of a closed system to be in any of its accessible microstates as equal.
Thus the foundational question remains open: can forgetful processes arise from isolated processes without any artificial discarding of information? Because forgetful processes are often called Markovian, we refer to the mechanism for forgetting as Markovianization, in the same spirit as the terms equilibration and thermalization 1,[3][4][5][6][7][8] . Indeed, Markovianization is likely to come about through mechanisms intimately related to these other processes. For instance, dissipative Markov processes have fixed points to which the system relaxes; this is a mechanism for equilibration, and also possibly for thermalization. We have previously argued for the emergence of Markovianization for mathematically typical processes, using averages with respect to the Haar measure 9 ; however, such processes are far from physically typical 1 .
In this paper, we identify a class of isolated physical processes which approximately Markovianize in a strong sense, where even the multi-time quantum correlations vanish. To do so, we employ large deviation bounds for approximate unitary designs derived by R. Low 10 , and apply them to the process tensor formalism [11][12][13] , which describes quantum stochastic processes. We show that, similar to the way in which quantum states thermalize, quantum processes can Markovianize in the sense that they can converge to a class of typical processes, satisfying a meaningful large deviation principle whenever they are undergone within a large environment and under complex enough-but not necessarily fully randomdynamics. As a proof of principle, we employ a recent efficient construction of approximate unitary designs with quantum circuits 14 to illustrate how a dilute gas would quickly Markovianize. These results directly impose bounds on complexity and timescales for standard master equations employed in the theory of open systems. Finally, we discuss possible extensions of our results to many-body systems with time-independent Hamiltonians. Our results are timely given the ever-increasing interest and relevance in determining the breakdown of the Markovian approximation in modern experiments [15][16][17][18] .

Results
Quantum stochastic processes. A classical stochastic process on a discrete set of times is the joint probability distribution of a time-ordered random variable, Pðx k ; ; x 0 Þ. A process is said to have finite memory whenever the state of the system at a given time is only conditionally dependent on its previous m states: Pðx k jx kÀ1 ; ; x 0 Þ ¼ Pðx k jx kÀ1 ; ; x kÀm Þ. Here, m is the Markov order; when m = 1 the process is called Markovian, and when m = 0 the process is called random. Finite memory processes, and in particular Markov processes, have garnered significant attention in the sciences for two principal reasons. First, the complexity of a process grows with the Markov order and thus it is easier to work with finite memory processes. Second, many physical processes tend to be well approximated by those with finite memory.
Generalizations of Markov processes and Markov order to the quantum realm have been plagued with technical difficulties 19 , which have their origin in the fundamentally invasive nature of quantum measurement. However, recently, a generalized and unambiguous characterization of quantum stochastic processes within the process tensor framework 11,20 has paved the way to alleviating these difficulties. The success of this framework lies in generalizing the notion of time-ordered events in the quantum realm.
Consider a system-environment composite SE of dimension d SE = d S d E with an initial state ρ (0) that undergoes a evolution U 0 . An intervention A 0 is then made on the system S alone, followed by evolution U 1 . For concreteness, onward we will consider U i ≠ U j . Then a second intervention A 1 on S alone. This continues until a final intervention A k is performed following U k . A quantum event x i at the i th time step corresponds to an outcome of the corresponding intervention, and is represented by a completely positive (CP) map A x i ðÁÞ : In other words, an intervention is the action of an instrument J ¼ fA is a completely positive trace preserving (CPTP) map. This is depicted schematically in Fig. 1. In general, the evolution U is allowed to be a CPTP map on SE. In this paper, however, we are interested in an isolated SE, where the Us are unitary transformations: UðÁÞ :¼ UðÁÞU y , with U a unitary operator.
The probability to observe a sequence of quantum events is given by This can be rewritten, clearly separating the influence of the environment from that of the interventions, in a multi-time generalization of the Born rule 21-23 : where T denotes transpose, Λ :¼ A x 0 Á Á Á A x k , and the effects on the system due to interaction with the environment have been isolated in the so-called process tensor ϒ. We have depicted ϒ and Λ in Fig. 1a as the red and green comb-like regions, respectively. A circuit depiction of the same process ϒ, along with the instruments Λ is given in Fig. 1b. Maps like the process tensor are abstract objects with many different representations 12 . In this manuscript, for convenience, we work with the Choi state representation 12,24 of the process tensor, shown in Eq. (10) of the Methods section. The process tensor ϒ is a complete representation of the stochastic quantum process, containing all accessible multi-time correlations [25][26][27][28] . Similarly, the tensor Λ contains all of the details of the instruments and their outcomes. This tensor, in general, is also a quantum comb, where the bond represents information fed forward through an ancillary system. Finally, the process tensor can be formally shown to be the quantum generalization of a classical stochastic process, satisfying a generalized extension theorem with consistency conditions for a family of joint probabilities to guarantee the existence of an underlying continuous quantum stochastic process 13 , and reducing to classical stochastic process in the correct limit 29,30 .
Measuring non-Markovianity. The convenience of using the Choi state ϒ is that it translates temporal correlations between timesteps into spatial correlations. Furthermore, as detailed in the Methods section on the process tensor, ϒ can be efficiently described when written as a matrix product operator 11,31 , whose bond dimension represents the dimension of a quantum environment that could mediate the non-Markovian correlations. In particular, when the bond dimension is one, the process is Markovian. Specifically, a process ϒ (M) is Markovian if and only if it has the form with E j:i a CPTP map on the system connecting the i th to the i + 1 th time 12,20 . This quantum Markov condition in Eq. (2) allows for a precise quantification of memory effects; it is fully consistent with the classical Markov condition, and contains all of the popular witnesses of quantum non-Markovianity 19 . Importantly, it allows for operationally meaningful measures of non-Markovianity: for instance, the relative entropy of the process tensor with respect to its marginals, which happen to be the closest Markovian process tensor, i.e. N S :¼ min ϒ ðMÞ Sðϒ k ϒ ðMÞ Þ, quantifies the probability of mistaking ϒ and ϒ (M) , which decreases in the number of realizations of the process n as expðÀnN S Þ. For the current considerations, a natural choice is the so-called diamond norm. Just as trace distance is a natural metric for differentiating two quantum states, in the sense of having a clear operational definition, the natural distance for differentiating two quantum channels is the diamond norm, which allows for the use of additional ancillas 32 . We are interested in optimally differentiating between a non-Markovian process from a Markovian one, which leads to the multi-time diamond distance: where 27,33 , with the supremum over i ≥ 1 and a set of CP maps fO i g. This definition generalizes the diamond norm for quantum channel distinguishability 34 (also called cb-norm 35 or completely bounded trace norm 24 ), reducing to it for a single step process tensor, and similarly being interpreted as the optimal probability to discriminate a process from the closest Markovian one in a single shot, given any set of measurements, which can be made together with an ancilla.
Vanishing non-Markovianity in Eq. (3) would imply that the process must have the form of Eq. (2). The derivations of such processes make ad-hoc assumptions such as artificially refreshing the environment between time-steps (i.e., assumption of an infinite bath) that render approximations such as Born-Markov. Classical processes additionally require randomness injection by hand for stochasticity. Here, we show that a class of underlying quantum mechanisms lead to the emergence of Markovianity without ad-hoc assumptions. Namely, We show that the above measure of non-Markovianity in Eq. (3) vanishes as the global SE dynamics becomes more complex. This is entirely analogous to entanglement being the underlying mechanism explaining the emergence of statistical mechanics from quantum dynamics alone and accounting for the artificial postulate of equal a-priori probabilities 3 . Quantum processes and the process tensor. a A k-step quantum process ϒ on system S alone is due to the time evolution of an initial systemenvironment (SE) state ρ (0) with distinct unitary transformations U i with i = 0, 1, …, k. In between each pair of unitaries, an external operation (e.g., a measurement) A i for i = 0, 1, …, k is applied; this can also be described by a tensor Λ. b An n-qubit SE-system ( 0 j i depicting a single qubit) with two-qubit gate interactions (depicted by vertical lines between squares) only: a subsystem qubit is probed at the i th step through A i . While the standard approach towards typicality or equilibrium properties concerns the whole SE dynamics and/or a single measurement on system S as in Standard Statistical Mechanics, we show that complex-not necessarily uniformly random-dynamics within large environments will be highly Markovian with high probability. Markovianization with unitary designs. The generic form of open quantum dynamics is non-Markovian, but, despite this, it is often very well approximated by simpler Markovian dynamics. How this memorylessness emerges is not dissimilar to questions, regarding the emergence of thermodynamic behavior, which have pervaded quantum mechanics since its conception. Indeed, it can be shown that canonical quantum states are typical [36][37][38][39] , and we now know that the fundamental postulate of equal a-priori probabilities of statistical mechanics can be traced back to the entanglement between subsystems and their environment 3 . It turns out that, very similarly, if we sample a generic quantum process occurring in a large finite environment at random, it will be almost Markovian with very high probability 9 .
This sampling procedure can be formalized through the socalled Haar probability measure, μ h , over the d-dimensional unitary group UðdÞ, which is the unique (up to a multiplicative constant) measure with the property that, if U 2 UðdÞ is distributed according to the Haar measure, then so is any composition UV or VU, with a fixed V 2 UðdÞ. It can be normalized to one, so as to constitute a legitimate probability measure 40 . The Haar measure allows one to swiftly obtain statistical properties of uniformly distributed quantities [40][41][42][43][44][45] and, furthermore, to prove concentration of measure results [46][47][48] ; these somewhat surprisingly imply that, when drawn from the right distribution, certain quantities will become overwhelmingly likely to be close to another fixed quantity as the Hilbert space dimension is increased. Henceforth, we write U~μ h to refer to U as distributed according to the Haar measure and, similarly, we use P h and E h to denote probabilities and expectations with respect to the Haar measure.
The Result by Modi et al. on Markovian Typicality 9 , which is reproduced in detail in the Methods section on, gives a mathematically sound result of concentration of measure around Markovian processes. However, it assumes a Haar-distributed uniform sampling of unitary dynamics, and we know that nature seldom behaves randomly 49,50 . The dynamics of a vast number of physically relevant models can be approximated as Markovian 51 , so can we say that these also satisfy a concentration of measure with respect to Markovianity?
In some circumstances, sets of physical processes can approximate some of the statistical features of the Haar measure 1,52-54 ; for example, consider the toy model depicted in Fig. 2, comprising a dilute gas of n particles evolving autonomously in a closed box. The gas particles interact with each other in one of two ways as they randomly move inside the box. Following and intervening on a special impurity particle, taken to be the system, this model can be approximately thought to be described by a circuit such as the one in Fig. 1b. The simplicity of this system suggests that it can only uniformly randomize after a large number of random two-qubit interactions, progressively resembling genuine Haar random dynamics.
One possible way to quantify this progressive resemblance of the Haar measure is given by the concept of unitary designs. In general an ϵ-approximate t-design, which we denote μ t ϵ , can be defined through for a suitable metric ∥ ⋅ ∥, where UðÁÞ :¼ UðÁÞU y and VðÁÞ :¼ VðÁÞV y are unitary maps with U; V 2 UðdÞ. Here, as above, the notation E Ω indicates the expectation value with respect to a given probability measure μ Ω , i.e. V $ μ t ϵ and U~μ h . That is, μ t ϵ approximates the Haar measure up to the t th moment with a small error ϵ. In the case we are interested in, the unitary maps will correspond to SE unitaries, as depicted in Fig. 1(a), according to the either the Haar measure or a unitary design. We also do not assume anything about the parameter t other than it is a positive non-zero integer. Notice what this would mean for a model similar to that of Fig. 2: as individual random two-body interactions of each kind accumulate, what we expect is for the dynamics to start scrambling their information across the whole gas in the box, progressively becoming more complex and uniformly random 55 . Unitary designs give us this finite quantification of the approximation to uniform Haar randomness and, in this case, it can give us a precise way to account for the progressive emergence of complexity from seemingly simple individual twobody interactions.
Unitary designs for t = 2, 3 have been widely studied 56-65 and efficient constructions are known for larger values of t 14,59,66 . The latter are of particular relevance, precisely as designs for large t, i.e., those with a higher complexity 55 , are expected to satisfy tighter large deviation bounds, approaching concentration of measure as the level and quality of the design increases.
Such large deviation bounds over approximate unitary designs were derived in a general form by R. Low 10 for a polynomial function satisfying a concentration of measure bound, and we now use them to demonstrate the phenomenon of Markovianization for corresponding classes of processes.

Theorem 1
Given a k-step process ϒ on a d S dimensional subsystem, generated from global unitary d SE dimensional SE dynamics distributed according to an ϵ-approximate unitary t-design μ t ϵ , the likelihood that its non-Markovianity exceeds any δ > 0 is bounded as where B is defined as Fig. 2 A toy model analogous to a system with dynamics given by an approximate unitary design with two kinds of two-qubit interactions only. An impurity particle (teal) immersed in a gas of n E particles (arrows depicting direction of motion) within a closed box, where all particles interact in pairs in one of two ways (dashed circles) at random, can be similarly described by an approximate unitary design. The result of Theorem 1 ensures that for a large enough n E and number of interactions, most processes analogous to this one with approximate unitary designs will be almost Markovian.
for any m ∈ (0, t/4] and where C is defined in Eq. (14) and B an upper bound on the expected norm-1 non- The proof is displayed in full in the Methods section. The overall strategy is as done by R. Low 10 : a bound on the moments E t ϵ ½N m ♦ is given in terms of B, C and η, followed by Markov's inequality. The quantity η is related to the ϵ-approximate unitary t-design μ t ϵ through for any m > 0 and corresponds to the sum of the moduli of the coefficients of N 2 2 . We explicitly determine a bound on this quantity within the proof of Theorem 1 in the Methods section, which is the one we take as definition in Eq. (7).
The choice of 0 < m≤t/4 can be made to optimize the right-handside of the inequality, which ideally should be small whenever δ is. The term d 3ð2kþ1Þ S =δ 2 arises from bounding N ♦ and Markov's inequality, while the three summands within square brackets will be small provided i) C is large, ii) B is small and iii) the unitary design sufficiently small ϵ and large t is well-approximate and high enough. For conditions i) and ii), we require a fixed k such that d E ) d 2kþ1 S : this implies B % 0, so that ignoring subleading terms, we require for a meaningful bound, as detailed in the Methods section on Convergence towards Markovianity.
Overall, the bound in Eq. (5) approaches concentration whenever d E is large relative to d S and k, together with large enough t, as shown in Fig. 3. Generally, it can be seen by inspection that the scaling in these cases will be polynomially vanishing in d E , exponentially vanishing in t (upon appropriate choice of parameter m), and becomes loose, polynomially in d S and exponentially in k. Therefore, the vast majority of processes sampled from such a t-design are indistinguishable from Markovian ones in this limit. This can be intuitively understood as that for processes of small subsystems in large environments (d E ) d 2kþ1 S ) undergoing complex enough dynamics (large enough t) will look almost Markovian with high probability if the system is probed not too many times (small k). We will now show how these processes can be modeled in terms of random circuits.
Markovianization by circuit design. While no explicit sets forming unitary t-designs for t ≥ 4 are known to date, several efficient constructions generating approximate unitary designs by quantum circuits are known. Using these constructions we can highlight the physical implications of the theorem above. We begin by discussing the details of one such construction. As suggested in Fig. 1b, this construction only requires simple twoqubit interactions and, under certain conditions, yields an approximate unitary design, from which we can use Eq. (5) in our main Theorem to verify that Markovianization emerges.
We focus specifically on Result 2 by Winter et al. 14 , reproduced in the Methods section on efficient unitary designs, where a circuit with interactions mediated by two-qubit diagonal gates with three random parameters is introduced. The intuition behind such construction is that repeated alternate applications of these diagonal gates quickly randomizes the system. Notice that this idea now fully captures the gas scenario depicted in Fig. 2, where we only have two types of random two-body interactions repeatedly occurring, and we focus on one of the particles of the gas. The detail of this construction is reproduced in the Methods section on efficient circuit unitary designs.
We can illustrate this idea in Fig. 4, where we depict an n-qubit SE composite with k interventions on one of the qubits, with the interactions within the circuit being only between pairs of qubits and of only two kinds; these form blocks of unitaries between each time-step i that we label W ' i , where ℓ is related to the amount of two-qubit interactions as explicitly defined in Eq. (42). The main Result 2 of by Winter et al. 14 states that for an n-qubit system, when t is of order ffiffiffi n p , a circuit W ' yields an ϵapproximate unitary t-design if ℓ ≥ t − log 2 (ϵ)/n, up to leading order in n and t.
Furthermore, of great relevance in this result is the fact that almost all 2-qubit gates in each repetition of W ' can be applied simultaneously because they commute 64,67 . Therefore, if W ' yields an approximate unitary design as above, the order of the non-commuting gate depth D, defined by Winter et al. 67 as the circuit depth when each commuting part of the circuit is counted as a single part, will coincide with the bound on the order of the number of repetitions ℓ. That is, the non-commuting gate depth asymptotes to We can now think of the system from the toy model of Fig. 2 as given by a spin locally interacting with a large, n E -qubit environment, via a random time-independent Hamiltonian, with Eq. (5) statistically predicting under which conditions memory effects can be neglected. Notice that this is only a physical picture evoked by the W ' circuits rather than exactly being the model described by it. In Fig. 5 we take such a system for a single qubit and demand a bound B≤0.01 on the probability P t ϵ ½N ♦ ≥ 0:1 for a k = 2 timestep process; with this, we plot the scaling of the noncommuting gate depth D required to achieve an ϵ = 10 −12 approximate unitary t-design using W ' circuits for different values of 2 ≤ t ≤ 10. While the number of 2-qubit gates is on the order of 10 4 , the number of repetitions ℓ is at most 12 for an approximate 10-design and stays mostly constant as the number of environment qubits increases.
This construction naturally accommodates the cartoon example in Fig. 2. As long as the two interactions in the example together generate the necessary level of complexity, Markovianization will emerge. This shows, in principle, how simple dynamics described by approximate unitary designs can Markovianize under the right conditions. Moreover, taking the physical interpretation of a qubit locally interacting through two-qubit diagonal unitaries with a large environment, it also hints at how , on P t ϵ ½N ! 0:1, the probability P t ϵ over an ϵ-approximate t-design for the non-Markovianity N to exceed δ = 0.1, against log 2 (d E ), where d E is environment dimension, for a subsystem qubit undergoing a joint closed approximate unitary design interaction between a given number of interventions k. We fix an ϵ = 10 −12 approximate unitary t-design for different values 2 ≤ t ≤ 10 and fixed values of timesteps k, optimizing m for each case. macroscopic systems can display Markovianization of small subsystem dynamics in circuits requiring just a small gate depth. Furthermore, for macroscopic systems with coarse observables, the same Markovianization behavior would remain resilient to a much larger number of interventions.

Discussion
We have shown how physical quantum processes Markovianize, i.e., forget the past, for a class of physically motivated systems that can finitely can approximate random ones. Forgetfulness is indeed a common feature of the world around us, and one that is crucial for doing science. Without forgetfulness, repeatability would be impossible. After all, if each carbon atom remembered its own past then it will be unique and there would be no sense in classifying atoms and molecules. Beyond these foundational considerations, our results have direct consequences for the study of open systems using standard tools, such as master equations and dynamical maps. The latter of which can be seen as a family of one-step process tensors (with initial SE correlations a minimum of two steps must be considered 16,68 ). Specifically, our results, for the case of k ≤ 2, can be used to estimate the time scale, using gate depth as a proxy, on which an approximate unitary design's open dynamics can be described (with high probability) with a truncated memory kernel 2,69,70 , or even a Markovian master equation.
Conversely, for larger k, our results would have implications for approximations made in computing higher order correlation functions, such as the quantum regression theorem 71 . These higher order approximations are independent of those at the level of dynamical maps, which can, e.g., be divisible, even when the process is non-Markovian 72 . This is reflected in the loosening behavior of the bound in Eq. (5) as the number of timesteps increases, which can be interpreted as a growing potential for temporal correlations to become relevant when more information about the process is accessible.
This breadth of applicability is in contrast with the results of Modi et al. 9 , where it was shown that quantum processes satisfy a concentration of measure with respect to Haar measure around Markovian ones, which has two main drawbacks: first, as stated above, Haar random interactions do not exist in nature and hence the relevance of the result is limited. Second, the rate of Markovianization is far too strong. Almost all processes, sampled according to the Haar measure, will simply look random, i.e., Markov order m = 0 even for a large k. This, unlike our current result, misses almost all interesting physical dynamical processes. While the behaviour of our large deviation bound is polynomial, rather than exponential, thus not exhibiting concentration per-se, we have nevertheless exemplified how, with modestly large environments and relatively simple interactions, almost Markovian processes can come about with high probability. Physical macroscopic environments will be far larger than the scale shown in Figs. 3 and 5.
Despite the fundamental relevance of our result, it is well known that typicality arguments can have limited reach. For instance, the exotic Hamiltonians, introduced by Gemmer et al. 73 , which lead to strange relaxation, may not Markovianize even though the SE process is highly complex with a large E. There is also still significant scope for further addressing physical aspects, such as the question of whether, and how, a time-independent Hamiltonian can give rise to an approximate unitary design 14 , the relevant time scales of Markovianization, or the potential role of different approaches to pseudo-randomness such as that by Kastoryano et al. 74 , where it is shown that driven quantum systems can converge rapidly to the uniform distribution. Furthermore, a renewed wave of interest in thermalization has come along with the socalled Eigenstate Thermalization Hypothesis (ETH), which is a stronger and seemingly more fundamental condition on thermalization [75][76][77][78][79][80] , and we would thus expect a deep connection in the sense of ETH between Markovianization and thermalization to be forthcoming. In any case, it is clear that many physical systems Markovianize at some scale, and it only remains to discover how. j i is a single qubit), the unitaries W ' , composed of ℓ alternate repetitions only two distinct types of random interactions (depicted by diamonds and squares joined by the interacting qubits), and defined by Eq. (42), generate an ϵ-approximate unitary t-design whenever ℓ ≥ t − log 2 (ϵ)/n, as shown bt Winter et al. 14 . This can be thought as stemming from repeated alternate applications of random 2-qubit gates diagonal in only two Pauli bases. A qubit probed with a set of operations fA i g on a system undergoing ϵ-approximate unitary t-design dynamics W ' on a large environment will Markovianize for small design error ϵ and large complexity t as specified in the main text.  14 for a 2-step process on a single qubit to Markovianize with respect to environment size. Scaling of the noncommuting gate depth D, given by the minimum amount of alternate repetitions ℓ of the two kinds of random two-qubit diagonal gates within the unitary W ' , plotted against the environment qubits n E = log 2 (d E ), to generate an ϵ = 10 −12 approximate unitary t-design for 2 ≤ t ≤ 10. This is such that for a single-qubit system undergoing a process with k = 2 timesteps, the probability P t ϵ for the non-Markovianity N exceeding 0.1 is less or equal than 0.01, i.e. P t ϵ ½N ! 0:1 B 0:01.

Methods
The process tensor. The Choi state representation of the process tensor is given by where each Ψ is a maximally entangled state on an ancillary space of dimension d 2 S , and where is a unitary operator acting on the whole SE together with the 2k ancillas. All identities act on the ancillary system, the U i are SE unitary operators at step i, and S i is a SWAP operator between system S and half of the i th ancillary space at the ith time-step of the process. The definition in Eq. (10) is a generalization of the standard Choi state for quantum channels, as given by the Choi-Jamiołkowski isomorphism (CJI) 24 . The CJI for quantum channels establishes a one to one correspondence with a quantum state on a larger Hilbert space, given as the action of the channel onto half a maximally entangled state. The standard definition uses unnormalized maximally entangled states; however, here we are concerned with distinguishability of Choi states through the diamond norm in Eq. (3) and the Schatten norms in Eq. (IV B), so we avoid a normalization factor in these by normalizing the Choi states by definition. A discussion in full depth about the process tensor, its different representations and its properties and relevance is given by Modi et al. 12 .
As stated in the main text, ϒ can be efficiently described when written as a matrix product operator 11,31 . A matrix product operator (MPO) gets its name from the representation of an n-body operatorÔ aŝ O ¼ ∑ fp;qg O q 1 q n p 1 p n p 1 p n q 1 q n , where the coefficients can be represented as a product of matrices, O q 1 q n p 1 p n ¼ tr½M Á Á Á M p n q n n . In particular, a matrix product density operator is a MPO with M i known as the bond dimension. For the process tensor ϒ, these matrices generically correspond to M and s ð0Þ are subsystem S basis vectors and U i is an SE unitary at timestep i 11 . This means the bond dimension of ϒ is d E , which in practice should be much smaller, given that only part of the environment interacts with the system at any given time.
A non-ambiguous measure of non-Markovianity. As with any distinguishability measure, the non-Markovianity metric of Eq. (3) is not unique, and we choose the diamond norm for its mentioned operational significance. However, more generally, for any Schatten p-norm k Xk p :¼ trðjXj p Þ 1 p , a similar quantity can be defined N p :¼ 1 2 min ϒ ðMÞ k ϒ À ϒ ðMÞ k p , as done with p = 1 in by Modi et al. 9 , whenever ϒ is normalized such that tr½ϒ ¼ tr½ϒ ðMÞ ¼ 1. Then, we have the hierarchy N 1 ≥ N 2 ≥ , induced by that of the Schatten norms. As the black diamond norm is generally difficult to compute exactly, a particularly useful relation to Eq. (3) is d À2kÀ1 S N ♦ ≤ N 1 ≤ N ♦ , in the sense that once any Schatten norm is known, the black diamond norm is automatically bounded.
Nevertheless, we highlight that, in general, any distinguishability measure N between a process ϒ and the closest Markovian one ϒ (M) will capture all non-Markovian features across multiple time steps, i.e., all multi-time phenomena and memory effects 20 . This is in contrast to other measures of non-Markovianity, e.g. trace-distance based measure 19 and other based on divisibility 81 , that have been proposed in recent years. In particular, all other measures relying on completely positive divisibility are only able to account for temporal correlations across at most three time-steps and are not sufficient to enforce the multi-time Markov condition 82 . This is even true in the classical case. Concretely, there are explicit examples of multi-time non-Markovian processes that are shown to be completely positive divisible processes, thus also deemed to be Markovian by the trace-distance based measure 20,82 . On the other hand, if a process satisfies the multi-time Markov condition, then it will be completely positive divisible.
In other words, the multi-time Markov condition is a stronger one that contains Markov conditions based on completely positive divisibility. This is why we consider the multi-time Markov condition in this manuscript.
Markovian typicality. In general, we say that a function f from a metric space S with metric Δ S and probability measure μ σ , to the real numbers, satisfies a concentration of measure around its mean if, for any point x 2 S and any δ > 0, where as done in the remainder of this manuscript, P σ and E σ explicitly refer to the probability and expectation with x~μ σ , and where L > 0 is the so-called Lipschitz constant of f, which can be determined according to jf ðxÞ À f ðyÞj ≤ L Δ S ðx; yÞ for any two points x; y 2 S. Whenever L is small, intuitively this implies that f varies slowly in such space. Finally, the function α σ is called a concentration rate; it generally must be vanishing in increasing δ in order for (12) to constitute concentration of measure, and it intuitively tells us how strong such concentration is.
Particularly well-known is the example of concentration of measure in the hypersphere of a high dimension, where for all functions that do not change too rapidly, i.e. with a small Lipschitz constant L, the function evaluated on a point picked uniformly at random will be close to its mean value with high probability, i.e. specifically α σ decays exponentially with − δ 2 . This is also known as Levy's lemma 46 and it has, remarkably, also been used by Winter et al. 3 to show that the fundamental theorem of statistical mechanics arises from entanglement.
Similarly, Modi et al. 9 showed that quantum processes satisfy a concentration of measure around Markovian ones, explaining the emergence of Markovianity without a-priori assumptions. In particular, there, the trace distance N 1 was used as a measure of non-Markovianity, which strictly speaking gives the distinguishability between explicitly constructed Choi states of corresponding process tensors and has no operational meaning; however, we can use the relation d À2kÀ1 S N ♦ ≤ N 1 ≤ N ♦ to relate this to the stricter notion of non-Markovianity defined in terms of the diamond norm in Eq. (3). This implies that the main result by Modi et al. 9 , where all SE unitaries of Eq. (10) were randomly sampled according to the Haar measure, can be written equivalently as where is the Lipschitz constant of N 1 , and is an upper bound on E h ½N 1 , the expected non-Markovianity over the Haar measure, with x : , and the expected purity, i.e., the noisiness of the process ϒ, over the Haar measure.
Holding everything else constant, the bound B ≥ E h ½N 1 satisfies so that the expected non-Markovianity vanishes in the d E → ∞ limit and becomes loosest in the k → ∞ limit case. The significance of Eq. (15) is thus that quantum processes with not too many interventions in high dimensional environments will look to be almost Markovian with high probability. This means that, even when processes generically carry temporal correlations, these are typically low, explaining the emergence of Markovian processes without ad-hoc assumptions such as the Born-Markov approximation of weak coupling 51 .
Unitary designs. The result in Eq. (15) assumes that the dynamics are Haar distributed; however, implementing a Haar random unitary requires an exponential number of two-qubit gates and random bits 83 , thus Haar random dynamics cannot be obtained efficiently in a physical setting.
An exact unitary t-design is defined 10 as a probability measure μ t on UðdÞ such that for all positive s ≤ t, and all d s × d s complex matrices X, As per the definition in Eq. (18), a unitary t-design reproduces up to the t th moment over the uniform distribution given by the Haar measure. In particular, μ t can consist of a finite ensemble fV i ; p i g N i¼1 of unitaries V i and probabilities pi , as is now common in applications such as so-called randomized benchmarking of error rates in quantum gates 60,62 .
Moreover, this definition can be relaxed by letting a unitary design approximate the Haar measure with a small error ϵ. In this manuscript we specifically employ the definition by R. Low 10 for unitary designs. It uses the fact that the definition of an exact t-design, μ t , can be written in terms of a balanced monomial Θ of degree less or equal to t in the elements of the unitaries U. A balanced monomial of degree t is a monomial in the unitary elements with precisely t conjugated and t unconjugated elements: for example, U ab U cd U Ã ef U Ã hg is a balanced monomial of degree 2. Thus, writing Eq. (18) in terms of matrix elements, this can be seen to be equivalent to requiring E t ½ΘðVÞ ¼ E h ½ΘðUÞ for all monomials Θ of degree s≤t. Similarly, for an ϵ-approximate t-design we adopt the definition by R. Low 10 with Eq. (4) implying for monomials Θ of degree s≤t. From now on, we will focus on the more general approximate designs. We will see below that the degree ϵ to which the distribution of the unitary dynamics on μ t ϵ differs from an exact design for given t depends on the complexity of the model.
Large deviation bounds for t-designs. The general idea for the main result by R. Low 10 (similarly applied before by Horodecki et al. 66 ) is that given a μ t ϵ distribution as an ϵ-approximate unitary t-design and a concentration result for a polynomial X of degree p, then one can compute the last term f t ϵ in with m≤t/2p, which will generally have a dependence f t ϵ ¼ f t ϵ ðϵ; t; X Þ. Using Markov's inequality which is the form of the main large-deviation bound. Specifically, the results that we employ are the following, proved R. Low 10 .

Theorem 2
(Large deviation bounds for t-designs by R. Low 10 ) Let X be a polynomial of degree T. Let f ðUÞ ¼ ∑ i α i Θ s i ðUÞ where Θ s i ðUÞ are monomials and let α(f) = ∑ i |α i |. Suppose that f has probability concentration and let μ t ϵ , be an ϵ-approximate unitary t-design, then for any integer m with 2mT≤t.
This is the most general result providing a large-deviations bound on approximate unitary designs, where ζ can be any quantity, in particular the expectation of f. The main idea from this result (similarly applied before by Horodecki et al. 66 ) is that given a μ t ϵ distribution as an ϵ-approximate unitary tdesign and a concentration result for a polynomial f of degree T, then one can compute where m≤t/2T. Using Markov's inequality we have which is the form of the main large deviations bound in Eq. (23). More precisely, the other two main results that come along with the proof of Theorem 2 by R. Low 10 , and allowing to compute the right hand-side of Eq. (25) are the following.

Lemma 3
(3.4 of by R. Low 10 ) Let X be a polynomial of degree T and ζ any constant. Let f ðUÞ ¼ ∑ i α i Θ s i ðUÞ where Θ s i ðUÞ are monomials and let α(f) = ∑ i |α i |. Then for an integer m such that 2mT≤t and μ t ϵ an ϵ-approximate unitary t-design, Lemma 4 (5.2 of by R. Low 10 ) Let X be any non-negative random variable with probability concentration where γ ≥ 0, then for any m > 0.
So, in essence, given these results, we determine the right-hand sides of Eq. (26) and Eq. (28) through the measure of non-Markovianity in Eq. (3) and all the other relevant quantities in such terms.

Proof of Theorem 1
A bound on the Haar moments of N 2 . Let us start by noticing that ∥X∥ 1 ≥∥X∥ 2 , so a concentration for N 1 given by (here 4 times the one defined in Eq. (14) in the main text), and B is defined in Eq. (15), also implies so that in turn Lemma 4 through Eq. (28) implies that for any m > 0.
A bound on the design moments of N 2 . For the case of all unitaries at each step being independently sampled, N 2 2 is a polynomial of degree p = 2 when the unitaries are all distinct (random interaction type). We can thus take N 2 2 and apply Lemma 3 for a unitary t-design μ t ϵ with t≥4m, which actually holds for real m > 0, as where η is the sum of the moduli of the coefficients of The proof of Lemma 3.4 by R. Low 10 requires m to be an integer through the multinomial theorem; in the notation of the cited paper, this can be relaxed to be a real number by applying the multinomial theorem for a real power: convergence requires an ordering such that jα t EM t j>2 1Àn jα tÀn EM tÀn j for each n = 1, …, t − 1 for both the approximate design and Haar expectations.
Let us explicitly write the process ϒ, defined in Eq. (10) in the main text, as a function of the set of unitaries U :¼ fU i g k i¼0 , i.e.
where here implicitly U ℓ stands for U ℓ ⊗ 1 2k−ancillas and the maximally entangled states Ψ are taken to be normalized. As the swaps between the system and the i th half ancillary system are given by S i ¼ ∑S αβ 1 β α h j i 1 where S αβ :¼ 1 E α j i β S , this can be written as Now, the standard approach to compute the sum of the moduli of the coefficients of a given polynomial is to evaluate on an argument (here a d SE × d SE matrix) full of ones (so that all single monomials equal to one) and take each summand to the corresponding modulus. We follow this approach, however, we first notice that the environment part in Eq. (34) is just a product of the environment parts of all unitaries and initial state. To see this, let U ¼ ∑U es e 0 s 0 es j i e 0 s 0 h j where e j i and s j i are E and S bases. Unitarity then implies ∑U ab es U ab ϵσ ¼ δ eϵ δ sσ , where the overline denotes complex conjugate, and so this means that tr E ½VS αβ UρU y S γδ V y ¼ ∑V es e 0 s 0 V eσ e 0 σ 0 U e 0 s 2 bs 0 2 U e 0 σ 2 bσ 0 2 ρ br bt ϕðSÞ where ϕ(S) stands for the system S part; for each b index the rest of the terms are summed over e; this generalizes similarly for any number of unitaries. This implies that at most d E terms need to be set to one and we can evaluate ϒ in a set of matrices J ¼ f1 E J S ; Á Á Á ; 1 E J S ; J E J S g with J a matrix with each element equal to one in the respective E or S systems: let ρ ¼ ∑ρ ese 0 s 0 es j i e 0 s 0 h j, then ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00629-w and hence (we now omit the subindex S on the J matrices for simplicity), where to obtain the second line we used the fact that J n = d n−1 J for positive integers n, here applied for n = 2, together with the trace over system S given by ∑ γ k Á γ k . This is similarly done to get the third line by ∑ δ k δ k ¼ 1 S , and taking the trace summing over γ kÀ1 , which can subsequently be done for all γ i and δ i . For the fourth line, the cyclicity of the trace was used, followed by an identity taken by summing up over α k , using J 2 = dJ, and taking the trace. This can be done through all remaining steps, giving the last equality. This, together with Eq. (32), implies that (now writing simply i, j for SE indices), where in the second line we used k Xk 2 1 ≤ d k Xk 2 2 for element-wise norms k Xk p p ¼ ð∑jx ij j p Þ and in the third line we used k ρk 2 2 ≤ 1.
Markov's inequality on N ♦ . As d À2k1À1 where in the third line we used Markov's inequality. This concludes the proof of Theorem 1.
Convergence towards Markovianity. We may first examine the third and penultimate lines leading to Eq. (38) for meaningful bounds P t ϵ ½N ♦ ≥ δ. The term d 3ð2kþ1Þ S =δ 2 arises from bounding the diamond norm and Markov's inequality; while δ is arbitrary, the d 3ð2kþ1Þ S could still be relevant when multiplied with E t ϵ N 2m 2 . This latter term will be small provided 1) C is large, 2) B is small and 3) the unitary design is approximate and high enough.
For 1) and 2), as detailed by Modi et al. 9 , we require a fixed k such that d E ) d 2kþ1

S
, together with m≤t/4. On the other hand if ϵ is non-zero, we require ϵ ( δ 2 2 The choice of real m is only restricted by 0 < m ≤ t/4, but otherwise is arbitrary. The right-hand side of Eq. (38) is not monotonic in m over all the remaining parameters, so it won't always be optimal for some fixed choice. One may thus optimize the choice of m numerically for each particular case. Efficient circuit unitary designs. As mentioned in the main text, we focus on Result 2 of of Winter et al. 14 . To begin with, an efficient approximation for a unitary design on a system composed of n-qubits is shown by Winter et al. 14 for a circuit labeled RDCðI 2 Þ, where the name stands for Random Diagonal Circuit, and refers to a circuit where I 2 ¼ fI i g is a set of subsets of qubit labels I i ⊂ {1, …, n}, such that |I i | = 2, i.e., at step i, I i picks a pair of qubits, to which a Pauli-Z-diagonal gate with three random parameters is applied. This construction can already be seen in the results of Winter et al. 64 as arising from only two types of random diagonal interactions, which can be simplified into a product of Z-diagonal ones.
A particular case which further simplifies things is then denoted by RDC ðtÞ disc ðI 2 Þ, where the subscript disc and the superscript t refer to discrete sets from which the parameters of the diagonal gates will be sampled, and which are determined by a given natural number t. Specifically, all gates in RDC ðtÞ disc ðI 2 Þ have the simplified form ðdiagf1; e iϕ 1 g diagf1; e iϕ 2 gÞ diagf1; 1; 1; e iϑ g; ð41Þ where diag denotes Pauli-Z diagonal, and with ϕ 1 , ϕ 2 chosen independently from the discrete set {2π m/(t + 1): m ∈ {0, …, t}} and ϑ chosen from {2π m/(⌊t/2⌋ + 1): m ∈ {0, …, ⌊t/2⌋}}. We emphasise that this is still a circuit with 2-qubit diagonal gates with only three random parameters each, and therein lies its simplicity. Now let H n = H ⊗n be n copies of the Hadamard gate, then the main Result 2 by Winter et al. 14 states that for an n-qubit system, when t is of order ffiffiffi n p , a circuit of the form yields an ϵ-approximate unitary t-design if up to leading order in n and t.
All the 2-qubit gates in each repetition of W ' , except those in H n , can be applied simultaneously because they commute 64,67 . Thus, as explained in the main text, if W ' yields an approximate unitary design, the order of the non-commuting gate depth D will coincide with the bound on the order of the number of repetitions ℓ.

Data availability
No datasets were generated or analyzed during the current study

Code availability
The code used in the analysis of the datasets is available from the corresponding authors on reasonable request.