Almost Markovian processes from closed dynamics

It is common, when dealing with quantum processes involving a subsystem of a much larger composite closed system, to treat them as effectively memory-less (Markovian). While open systems theory tells us that non-Markovian processes should be the norm, the ubiquity of Markovian processes is undeniable. Here, without resorting to the Born-Markov assumption of weak coupling or making any approximations, we formally prove that processes are close to Markovian ones, when the subsystem is sufficiently small compared to the remainder of the composite, with a probability that tends to unity exponentially in the size of the latter. We also show that, for a fixed global system size, it may not be possible to neglect non-Markovian effects when the process is allowed to continue for long enough. However, detecting non-Markovianity for such processes would usually require non-trivial entangling resources. Our results have foundational importance, as they give birth to almost Markovian processes from composite closed dynamics, and to obtain them we introduce a new notion of equilibration that is far stronger than the conventional one and show that this stronger equilibration is attained.


Introduction
The quest towards understanding how thermodynamics emerges from quantum theory has seen a great deal of recent progress [1,2]. Most prominently, the conundrum of how to recover irreversible phenomena that obey the second law of thermodynamics, starting from reversible and recurrent Schrödinger dynamics, has been rectified by considering an analogous dynamical equilibrium. This is known generically as equilibration on average; it implies that time-dependent quantum properties evolve towards a certain fixed equilibrium value, and stay close to it for most times 1 .
Equilibration on average has been widely studied Pedro Figueroa-Romero: pedro.figueroaromero@monash.edu 1 Even stronger (but usually the one of interest) is the notion of equilibration in a finite time interval. The corresponding statement for a physical system replaces the term equilibration with thermalisation and the equilibrium state with a thermal state. For detail see e.g. [1] −→ for small non-degenerate subsystems of larger closed systems, which are usually taken to be (quasi) isolated, and for particular models (see e.g. [3]). In a typical setup, the closed dynamics of a system and environment, as illustrated in Fig. 1(a), is considered. Remarkably, it can be shown that the expectation value for any observable O on the system alone lies in the neighbourhood of its time-averaged value for almost all times [ Fig. 1(b)]. Yet, such equilibration results, in full generality, remain unsolved [1]; in particular, it is unknown whether they extend to correlations between observables at different times.
To see this, consider the case depicted in Fig. 1(c), where the subsystem is interrogated (by measuring, or otherwise manipulating it) multiple times as it evolves with the remainder of the closed system, which acts as an environment. In such a scenario, the flow of information between interventions must be taken into Figure 2: Markovian dynamics. A quantum circuit corresponding to the Markov approximation, which is mathematically equivalent to starting off with an uncorrelated ρ (1) E ⊗ ρ S state and then artificially refreshing the environment (discarding and replacing with a new state) at each step (cf. Fig. 1c) account; in particular, there may be intrinsic temporal correlations mediated by the environment, in other words memory built up within the dynamics itself, independently of the external interventions, which are the hallmark of so-called non-Markovian dynamics. Equilibration on average says that the process erases the information contained in the initial state of the system; however, there may still be non-Markovian memory of the initial state encoded in the temporal correlations between observables. It is therefore an open question whether this information is truly scrambled by the process and therefore inaccessible to the system at long times. Such memoryless dynamics are anticipated for the subparts of closed manybody dynamics [4]. This scrambling of information in correlations would represent a stronger notion of equilibration than the usual one, which only considers statistics of individual observables at a given time. Moreover, the absence of memory effects would generically imply the usual kind of equilibration.
Searching for this stronger notion of equilibration is not physically unmotivated. While, strictly speaking, non-Markovian dynamics is the norm for subsystems of a closed system evolving with a time-independent Hamiltonian [5,6], Markov approximations are commonly made throughout physics (and elsewhere) to good effect. However, the standard approach to making this approximation in open quantum systems theory involves a series of assumptions about weak interactions with an environment sufficiently large and ergodic that the memory of past interactions gets practically lost [7,8]; this is mathematically equivalent (through an analogue of Stinespring's dilation theorem) to continually refreshing (discarding and replacing) the environment's state [9], i.e., artificially throwing away information from the environment (see Fig. 2). In reality, this is not usually how a closed system evolves, leading to the question: In the general case, how does dissipative Markovian dynamics emerge from closed Schrödinger evolution? That this puzzle is similar to the ones concerning equilibration and thermalisation should not be too surprising, as thermalisation is often achieved by assuming a Markov process [10].
Here, we give an initial answer to this deep foundational question. We show that, for the class of evolu-tions given by the Haar measure (typically involving an interaction of the subsystem with a significant portion of the remainder), above a critical time scale the dynamics of a subsystem is exponentially close to a Markov process as a function of the relative size of the subsystem compared with the whole. Moreover, our result does not assume weak coupling between the system and the environment, and we do not employ any approximations. Specifically, we use averages over the unitary group together with the concentration of measure phenomenon [11][12][13][14][15] to determine when a subsystem dynamics might typically be close to Markovian. Our approach is broadly similar to that of Refs. [14,16] both in spirit and in terms of the mathematical tools employed. In other words, our main result gives a mechanism for equilibration starting from closed dynamics without employing any approximations.
Our results are possible thanks to recent efforts towards the understanding of general quantum processes, culminating in the development of the process tensor framework [17][18][19], which allows for the complete and compact description of an open system's evolution under the influence of a series of external manipulations or measurements [20,21].
Crucially, the process tensor framework leads to an unambiguous distinction between Markovian and non-Markovian dynamics, analogous to the usual one for classical stochastic processes. In this sense, it subsumes other approaches to quantum non-Markovianity (when compared with respect to the same process), including those based on increasing trace-distance distinguishability [22]. It allows for a precise quantification of memory effects (including across multiple time points) [5], that we employ here as a natural setting to explore temporal correlations and the role that they play in a generalized, or strong, notion of equilibration.
2 Quantum processes with memory: the process tensor We begin with the motivation of witnessing equilibration effects through a sequence of measurements on the system. We consider a closed environmentsystem (ES) with finite dimensional Hilbert space The system S is acted on with an initial measurement operation M 0 followed by a series of subsequent measurement operations {M i } k−1 i=1 at a sequence of fixed times, between which it evolves unitarily with the environment E; these operations are represented by trace non-increasing completely positive (CP) maps, and they include information about how the system is affected by the measurement procedure. Our results also hold for more general operations than just measurements, including unitary transformations or any other experimental manipulations of the system, but we focus on the former in order to be concrete. This scenario is depicted in Fig. 1(c). Formally, a kstep process tensor is a completely positive causal 2 map T k:0 from the set {M i } k−1 i=0 to output states, ]. It must be stressed, however, that although experimental interventions play an important role in the reconstruction of a process through a generalized quantum process tomography, they do not determine its nature, i.e., the process tensor itself accounts precisely for the intrinsic dynamics occurring independently of the control of the experimenter. Moreover, like any other CP map, the process tensor can be expressed as a (many-body) quantum state through the Choi-Jamiołkowski isomorphism [21,23,24]. For a k-step process tensor, the corresponding (normalized) Choi state Υ k:0 can be obtained by swapping the system with one half of a maximally entangled state 3 and the system is swapped with an A subsystem at each time. With this convention, the Choi state of the process tensor can be seen to take the form where all identity operators act on the ancillary systems, U i is the ES unitary operator that evolves the joint system following the control operation at time i and the S i are swap operators between system S and the ancillary system A i . The corresponding quantum circuit is depicted in Fig. 3. Figure 3: Choi state of a process tensor. A quantum circuit diagram for the Choi state of a k-step process tensor under joint ES unitary evolution: the final state Υ is a many-body quantum state capturing all the properties of the process. 2 Meaning that the output of a process tensor at time k does not depend on deterministic inputs at later times K > k. 3 The dimensions of each A i , B i must be fixed to d S to make the swap operator well defined.
Throughout the paper, we write Υ for the Choi state of a k-step process unless otherwise specified.

An unambiguous measure of non-Markovianity
Among other important properties, the process tensor leads to a well-defined Markov criterion, from which it is possible to construct a family of operationally meaningful measures of non-Markovianity, many of which can be stated simply as distances between a process tensor's Choi state Υ and that of a corresponding Markovian one Υ (M) . The latter must take the form of a tensor product of quantum maps E i:i−1 connecting adjacent pairs of time steps: [25]. In [5], a such measure of non-Markovianity was first introduced through the relative entropy between a process and its corresponding Markovian one, R(Υ Υ (M) ) ≡ tr[Υ(log Υ − log Υ (M) )], minimized over the latter. Here, however, in analogy with other studies on equilibration, we choose the measure of non-Markovianity defined in terms of the trace distance D as where X 1 ≡ tr √ XX † is the trace norm (or Schatten 1-norm); in particular, this is related to relative entropy through the so-called quantum Pinsker inequality, R(Υ Υ (M) ) ≥ 2 D 2 (Υ, Υ (M) ). The trace distance is a natural choice, as it represents an important distinguishability measure and has played a central role in the discovery of several major results on typicality [14] and equilibration [1,16,26] (or lack thereof, e.g. [27]). It has an experimental importance for comparing quantum processes [28] and is naturally related to other important distinguishability measures such as the diamond norm or the fidelity. However, it is important to point out that the above measure is fundamentally different (and more powerful) than other measures of non-Markovianity that have been proposed in recent years, including one that employs trace distance [22]. The measure we use here captures all non-Markovian features across multiple time steps and provides an operational meaning for non-Markovianity. Equipped with this measure, we are in a position to determine what value it takes for a typical process, answering the question posed at the beginning of this manuscript.

Typicality of Markovian processes
In order to quantify non-Markovianity, as defined in Eq. (2), for a typical process, we first need a concrete notion of typicality. There are many ways in which one could sample from the space of quantum processes, but as we are particularly interested in the typicality (or absence thereof) of Markovian processes independently of ad-hoc assumptions like weak coupling, we would like our measure to assign nonvanishing probabilities to mathematically generic unitary dynamics on the closed ES system.
Here, we achieve this by sampling the evolution from the unitarily invariant probability measure, the so-called Haar measure. This has the additional advantage of allowing us to average by means of random matrix theory techniques [15,[29][30][31][32] and leads to the relatively straightforward application of concentration of measure results [11][12][13]. Specifically, we use the Haar measure to sample two distinct types of ES evolution: 1. Random interaction: All U i independently chosen, where U i represents the closed (unitary) ES evolution between the application of measurement operation M i−1 and M i , as depicted in Fig. 1(c); the entire set enters into the process tensor through Eq. (1). In the first case, the global system will quickly explore its entire (pure) state space for any initial state. The second case corresponds more closely to what one might expect for a truly closed system, where the Hamiltonian remains the same throughout the process. These correspond to two extremes; more generally, the dynamics from step to step may be related but not identical.
We highlight that under this sampling method, and thus this notion of typicality, the interaction or information flow between all parts of the whole system plus environment is relevant, i.e. no parts of the environment dimension are superfluous. We will now show that this kind of dynamics is almost always close to Markovian.

Main result
Our main result states that, for a randomly sampled k-step quantum process Υ undergone by a d Sdimensional subsystem of a larger d E d S -dimensional composite, the probability for the non-Markovianity N to exceed a function of k, d S and d E , that becomes very small in the large d E limit, itself becomes small in that limit. Precisely, we prove that, for an arbitrary > 0, where with c = 1/4 for a constant interaction process and c = (k + 1)/4 for a random interaction process. The function is an upper bound on the expected non-Markovianity E[N ] whose details depend on the way in which processes are sampled. The proof is presented in full in Appendix E, and it can be outlined as follows: we first bound the difference between N for two different processes in terms of the distance between the unitaries used to generate them; Levy's lemma [11,13,14] then states that the fraction of processes with non-Markovianity more than away from the expectation value (which we upper-bound by B k ) is bounded by a concentration function. By considering the geometry of the spaces of unitaries we are sampling from, we are able to use the exponential function appearing on the right hand side of Eq. (3). This implies that our result holds particularly for this class of evolutions where all parts of system and environment effectively interact, in contrast with some commonly considered open systems models.
Our result is meaningful when both B k + and e −C 2 are small; the latter is fulfilled in the small subsystem or large environment limit, which in our setting means . The size of B k depends entirely on the size of the average purity of the process tensor E[tr(Υ 2 )], which we have computed analytically for the case of Haar-randomly sampled evolutions (which equivalently could be sampled from (2k + 2)-unitary designs [33,34], as we will discuss below) in Appendix F.
The purity tr(Υ 2 ) is a quantifier of the mixedness (uniformity of eigenvalues) of a positive operator. When computed on reduced states of bipartite systems it can also serve as a quantifier of entanglement. In Appendix D, we analytically compute the expected purity of the Choi state of a process tensor E[Υ 2 ] in both the constant and the random interaction pictures, which can be directly translated as a quantifier for noisiness of the quantum process itself and entanglement between system and environment. The expressions take the form for the random interaction picture and in the constant interaction case, where S n is the symmetric group over n elements, Wg is known as the Weingarten function [30], τ is a product of entries of ρ depending on permutations τ and ∆ is a product, scaling with k, of monomials in d E and d S depending on permutations σ and τ .

Limiting cases
In both the constant and random interaction cases, B k is a well-behaved rational function of d E , d S and k. Eq. (6) takes a non-trivial form mainly because of the Wg function (which is intrinsic to the Haarunitary averaging in the constant interaction picture). However, due to results found in [30,31], we can still study analytically the behaviour of the bound B k for both cases in the two following limits.
In the large environment limit, which corresponds to the purity of the maximally mixed state. Note that the averaging occurs after computing the purity; that is, independently of which case is considered. Every process sampled will be indistinguishable from the maximally noisy (and hence Markovian) one, with probability that tends to one as d E is increased: Furthermore, from Eq. (5) it is seen that it does so at The other interesting limiting case is the one where the ES dimension is fixed, but the number of time steps is taken to be very large. The resulting process encodes all high order correlation functions between observables over a long period of time. Again for both cases, the expected purity in this limit goes as which corresponds to the maximally mixed purity of the environment. This implies that the Choi state of the full ES unitary process is maximally entangled between S and E. In this limit, we get correspondingly meaning that nothing can be said about typical non-Markovianity; a typical process in this limit could be highly non-Markovian. Indeed, we expect this to be the case, since the finite-dimensional ES space will have a finite recurrence time. Both our main result in Eq. (3), and the bound on non-Markovianity, Eq. (4), generalize well-known results for quantum states, i.e., when k = 0. In this case, our main result reduces to the usual one for the typicality of maximally mixed states (or bipartite maximally entangled states), with C(d E , d S ) = d E d S /4 (see e.g. [35]). As we also detail in Appendix F, when we take the Haar measure, in either case, we recover the average purity where ρ S ≡ Υ 0:0 = tr E (U ρ U † ), as found in [15,[36][37][38][39][40][41]. This then leads to which goes as d S /4d E when d E d S , and is also a standard result [35].

Numerical examples
To support our results, we sampled Choi states Υ numerically in the random interaction case and computed their corresponding average non-Markovianity . Process tensors were generated by sampling Haar random unitaries according to [42].
as a function of environment dimension d E for a fixed system dimension d S = 2, obtaining the behaviour shown in Fig. 5. For constant interaction the numerical results are practically indistinguishable from those in the random case, but as mentioned, the analytical bound B k is much harder to compute exactly. This suggests that either a simpler bound exists or that it might be possible to simplify the one we have obtained. As expected, our numerical results fall within the bound B k (d E , d S = 2) and they behave similarly; we notice that the bound in general seems to be somewhat loose, and become loosest when d E d 2k+1

S
, implying that non-Markovianity might be hard to detect even when not strictly in the large environment limit. However, it does saturate rapidly as k increases.
So far, our results are valid for process tensors constructed with Haar random unitaries at k evenly spaced steps; we are effectively considering a strong interaction between system and environment which rapidly scrambles quantum information in both [33]. The mechanism for equilibration is precisely that of dephasing [1], or effectively, the scrambling of information on the initial state of the system. This suggests that, even when timescales will differ with the type of evolution considered, most physical evolutions fall within our result, with e.g. a weaker behavior in k.
Our results could potentially be extended to their analogues sampled from unitary n-designs, approximate or otherwise [43], that reproduce the Haar distribution up to the n-th moment. As mentioned above, our expression for the upper bound in Eq. (4) is identical for samples from a 2k + 2-designs, and, if the non-Markovianity N could be approximated by a polynomial function of U , Eq. (3) would also hold for n-designs with large enough n [44]. This type of distribution has recently been shown to be approximated by a wide class of physical evolutions [45]. In particular, the unitary circuits involved in complex quantum computations typically look Haar random [34]; if it holds in this case, our result would therefore suggest that the behaviour of components of a quantum computer would not typically exhibit non-Markovian memory (at least due to the computation itself). Furthermore, even when it's known that ensembles {e −iHt } t≥0 for time-independent Hamiltonians H cannot strictly become Haar random [33], certain complex Hamiltonian systems can be treated as random up to a certain Haar-moment via unitary designs [45].
These issues should be investigated further, however, we will now show that our results still hold at a coarse-grained level, where the intermediate dynamics corresponds to products of Haar random unitaries, which are not themselves Haar random.

Observing non-Markovianity
The choice of unitaries going into the process tensor in the previous sections (in our case drawn from the Haar measure), dictates a time scale for the system, up to a freely chosen energy scale. However, we can straightforwardly construct process tensors on a longer, coarse-grained time scale by simply allowing the system to evolve (with an identity operation) between some subset of time steps. In fact, process tensors at all time scales should be related in this way to an underlying process tensor with an infinite number of steps [46].
To see that our main result directly applies to any process tensor that can be obtained through coarse graining, we again consider the definition of our non-Markovianity measure, given in Eq. (2). Consider the coarser grained process tensor Υ k: are replaced by identity operations (i.e., the system is simply left to evolve).
since the set of allowed Υ (M) c strictly contains the allowed Υ (M) at the finer-grained level. This renders the new process less distinguishable from a Markovian one, i.e., coarse-graining can only make processes more Markovian.
The physical intuition behind this result is that the amount of information which can be encoded in a coarse grained process tensor is strictly less than that in its parent process tensor. In fact, this is a key feature of non-Markovian memory: the memory should decrease under coarse graining. On the other hand, due to the same reasoning, we cannot say anything about finer-grained dynamics. One approach to tackling this issue would be to choose a different sampling procedure which explicitly takes scales into account. However, we leave this for future work.
There is another important limitation for observing non-Markovianity. The operational interpretation of the trace distance, discussed in the previous section, implies that observing non-Markovianity requires applying a measurement that is an eigen-projector operator of Υ − Υ (M) . The optimal measurement will, in general, be entangled across all time steps. In practice, this is hard to achieve and typically one considers a sequence of local measurements m ∈ M. In general, for an any set of measurements M we can define a restricted measure of non-Markovianity detectable with that set: , which means that the detectable non-Markovianity will be smaller. This is akin to the eigenstate thermalization hypothesis [47], where all eigenstates of a physical Hamiltonian look uniformally distributed with respect to most 'physically reasonable' observables. In our setting, this means that looking for non-Markovianity with observables that are local in time -i.e., 'physically reasonable' -we find almost no temporal correlations.
The locality constraint, along with monotonicity of non-Markovianity under coarse graining, have further important consequences for a broad class of open systems studies where master equations are employed [48]. Since master equations usually only account for two-point correlations with local measurements, they will be insensitive to most of the temporal correlations being accounted for by our measure, leading to an even greater likelihood for their descriptions to be Markovian. We now discuss the broader implications of our results.

Discussions and Conclusions
Our results imply that it is fundamentally hard to observe non-Markovianity in a typical process and thus go some way to explaining the overwhelming success of Markovian theories.
Specifically, we have shown that, even when interacting strongly with the wider composite system, a subsystem will typically undergo highly Markovian dynamics when the rest of the system has a sufficiently large dimension, and that the probability to be significantly non-Markovian vanishes with the latter. Our main result formalizes the notion that in the large environment limit a quantum process, taken uniformly at random, will be almost Markovian with very high probability. This corroborates the common understanding of the Born-Markov approximation [7,8], but, crucially, we make no assumptions about weak coupling between E and S. Instead, in the Haar random interactions we consider, every part of the sys-tem typically interacts significantly with every part of E. This is in contrast to many open systems models, even those with superficially infinite dimensional baths, where the effective dimension of the environment is relatively small [49]; it can always be bounded by a function of time scales in the system-environment Hamiltonian [50], which could be encoded in a bath spectral density. Our result is also more general than the scenario usually considered, since it accounts for interventions and thus the flow of information between S and E across multiple times.
While it may still be possible to observe non-Markovian behaviour at a time scale that is smaller than the fundamental time scale set by the chosen unitaries, Eq. (13) tells us that any coarse grained process will remain concentrated around the Markovian ones in the large environment limit. Otherwise, for larger and larger systems, one needs an ever increasing number of time steps, corresponding to higher order correlations, in order to increase the probability of witnessing non-Markovianity. However, even in this case, from the discussion in the previous section, we know that the measurement on this large number of times steps will be temporally entangled, which may also be difficult to achieve.
We have introduced a stronger notion of equilibration than the conventional one [1], so that strong equilibration implies weak equilibration. Our main result shows that it is possible to attain strong equilibration, where the information of the initial state of the system does not exist in multi-time correlations, much less expectation values at different times. Our result can be interpreted as a mechanism for equilibration and thermalisation. That is, we have shown that the dynamics of the system alone is almost Markovian, and such processes have well defined fixed points [51]. Even if we do not characterise the equilibrium state here, e.g. as done for quantum states under randomized local Hamiltonians [52], we may conclude that the state of the system approaches this fixed point and the expectation value for any observable will lie in the neighbourhood of the time-averaged value [as depicted in Fig. 1(b)] for most times. However, one point of departure between the two results is our reliance on Haar sampling.
More precisely, we have considered the case where the joint unitary evolution between each time step is drawn uniformly at random according to the Haar measure, for the two limiting cases of the random unitary operators at each step being the same or independent. This is, however, different from simply sampling the Choi state of the process Haar randomly. It is also not equivalent to picking Haar random states for the system at every time, as one would expect if the results of Ref. [14] were to hold independently at all times.
Nevertheless, tackling the issues we have pointed out, a notion of equilibration in general processes can potentially be made precise and other issues such as equilibration time scales can also be readily explored; the approach that we have presented constitutes a first significant step towards achieving this goal.

A The process tensor
The Choi state of a process tensor with initial state ρ is given by with where all identities are in the total ancillary system and the U i are ES unitary operators at step i (here the initial unitary U 0 is taken simply to randomize the initial state according to the Haar measure), and with S αβ = 1 E ⊗ |α β|. Notice that S ασ = S † σα , S ab S † cd = δ bd S ac and tr(S ab ) = d E δ ab . Also, Υ ∈ S(H S k i=1 H AiBi ), i.e. the resulting Choi process tensor state belongs to the S-ancillary system, which has dimension d 2k+1 S .

B Derivation of an upper-bound on non-Markovianity
As a first approach, given the complexity of computing (averages over) the Markovian Choi state, we may upper bound this distance by the one with respect to the maximally mixed state (the noisiest Markovian process possible) We may further bound this by considering the following cases separately.
where | · | denotes the standard absolute value, so using the inequality X 1 ≤ dim(X) X 2 for a square matrix X (which can be derived from the Cauchy-Schwarz inequality applied to the eigenvalues of X), where X 2 = tr(XX † ) is the Schatten 2-norm, Furthermore, applying Jensen's inequality (in particular for the square-root, E[ √ X] ≤ E[X]) with expectation over evolution (either in the ergodic or time-independent case), this gives

S
: This case is a small subsystem limit for most k.
Directly applying X 1 ≤ dim(X) X 2 as before, Similarly, taking the average over evolution, by means of Jensen's inequality, This proves that the function which can also be seen as a so-called 1-fold twirl and be computed by the Schur-Weyl duality [29,33]. When we deal with averages over the Haar measure with several unitaries we will denote by E U i and E U the ergodic and time-independent cases, respectively.

C.2 Schur-Weyl duality and the moments of the unitary group
The Schur-Weyl duality allows to compute integrals of twirled maps over the Haar measure given by where X ∈ H ⊗k ; in particular, the case U † AU XU † BU can be taken to a twirled form and be evaluated as in e.g. in [29], which we here use to compute the average purity of the Choi (ergodic) process tensor. In particular, Φ (k) commutes with all V ⊗k such that V ∈ U(d), and by such property, the Schur-Weyl duality assures that also where the sum is over the symmetric group S k on k symbols, P τ is a permutation operator defined by P τ |x 1 , . . . , x k = |x τ (1) , . . . , x τ (k) and f τ is a linear function of X [33]. The k = 1 case is well known [15] and in general for any operator acting on C d results in a completely depolarizing channel, A related result, which can equivalently be used to compute integrals over the unitary group for polynomials in the unitaries with the Haar measure, and which we also make use of in this work, is the one for the k-moments of the d-dimensional unitary group [30], where U ij is the ij component of U ∈ U(d) and Wg is known as the Weingarten function (here specifically for the unitary group) [31,53]. The Weingarten function can be defined in different ways, the most common being in terms of a sum over partitions (or, equivalently, Young tableaux ) of k and the so-called characters of the symmetric group; it is a fairly complicated function to evaluate explicitly and we refer to [29,30] for the details, in any case it is simply a rational function in the dimension argument d (particular cases are often given in the literature for small k, see e.g. [30,33]). An alternative is to perform numerical calculations, as a computer package in the Mathematica software was developed in [32] that allows to numerically evaluate Eq. (28) and the Weingarten function (a subsequent version for Maple is also now available in [54]).
Here we obtain the expected states E(Υ) and expected purities E[tr(Υ 2 )] over the Haar measure of Choi process tensor states Υ as defined by Eq. (14). We compute these averages for both the ergodic E U i and timeindependent E U cases. The time-independent averages are written in terms of the Wg function; we express them explicitly for the simplest cases and we show that in the small subsystem limit they reduce to the corresponding ergodic averages.
The asymptotic behavior of Eq. (28) boils down to that of the Wg function and in [30] it has been shown that as a refinement of a result in [31], where #σ is the number of cycles of the permutation σ counting also fixed points (assignments from an element to itself, σ(x) = x).

D.1 Ergodic average process Choi state
The ergodic average Choi state of a k step process is given by the maximally mixed state, This can be seen by a simple guess, as there are k + 1 independent integrals of a 1-fold twirl kind, and both the initial state and the S operators have trace one. This is also an intuitive result, here corresponding to the noisiest process possible. The detail of the calculation is as follows.
From the definition of the Choi process tensor state (14), in the ergodic case, where the sum is over repeated indices (Greek letters), hence by repeatedly evaluating the 1-fold twirl (27), with d = d E d S (also, the unitary groups are implicitly of dimension d E d S ), we get where the Haar measure symbols are implied but omitted just to economize notation.

D.2 Time-independent average process Choi state
The average Choi process tensor for k steps in the time-independent case can be given in terms of a set {|s of S system bases for i = 0, . . . , k as with implicit sum over all repeated basis (s ( ) i ) indices, where here S k+1 is the symmetric group on {0, . . . , k}, and with the definitions where {|e , with i = 0, . . . , k, also summed over all its elements, is a set of environment bases and the ∆ term is simply a monomial in d E with degree determined by σ and τ .
The case k = 0 recovers E U (Υ 0:0 ) = 1 S /d S as expected, as no process occurs.
In the case of the superchannel k = 1, we get Eq. (41), which is quite different from the maximally mixed state that arises in the ergodic case, although it converges to it as d E → ∞. We notice that tr[E U (Υ 1:0 )] = 1 as expected and one may also verify the purity of this average state to be given by Eq. (42) which is in fact close to that of the maximally mixed state for all d E , d S and 1/d S ≤ tr(ρ 2 S ) ≤ 1, and converges to it when The detail of the calculation is as follows. By definition, We now consider directly decomposing the unitaries as U = U ab |a b| and U † = U * a b |b a | -notice that the a and b labels refer to the whole ES space-, introducing an E basis in the S operators as S ab = |ea eb|, and then evaluating the k + 1 moments of U(d E d S ) by means of Eq. (28), where here S k+1 denotes the symmetric group on {0, . . . , k}, and which taking i → ς and j → ς to recover each E and S part explicitly, turns into and thus as stated by Eq.(33).

D.2.2 Small subsystem limit
When looking at the limit d E → ∞, the term that does not vanish is the one with σ, τ = 1 k+1 , i.e. with both permutations being identities, as these generate the most numerator powers in d E via the δ terms (correspondingly, ∆ k,σ,τ in Eq.(35)) when summed over all e i 's, i.e.
all other permutations will vanish because of the d E powers in the denominator generated by the Wg function, the least powers produced by it are those when στ −1 = 1 k+1 because #1 n = #[(1)(2) · · · (n)] = n, i.e. identity produces the greatest number of cycles, being the number of fixed points. Finally 0 ς 0 0 ς 0 |ρ| 0 ς 0 = tr ρ = 1, and thus matching the average over an ergodic process.
Below we obtain the Lipschitz constant for the trace distance between a Choi process tensor and the maximally mixed state as a function of unitary operators. We consider both the time-independent (U i = U j = U for all steps i = j) and the ergodic case (U i = U j for all steps i = j).

E.1.2 Levy's lemma
Let M be a manifold with metric δ M and probability measure µ M , and let f : M → R be an η-Lipschitz function, then for all ε > 0, where α M is known as a concentration function for M , and is defined as (an upper bound on) the measure of the set of points in M more than a distance x from the minimal-boundary volume enclosing half of M . See [11,55].

E.2 Main proof
We . We highlight that the evolution for Υ is sampled from the Haar measure, or equivalently one may replace P with something like P Υ∼µ U i .

E.2.1 The Lipschitz constants
We take τ : Consider for now different labelings for the unitary at each time step. By the triangle inequality, where it is clear by context that U stands for U ⊗ 1. Now we use Lemma 1 of [56], which states that for two unitaries A, B and any σ. For the simplest case with k = 1, this gives and doing similarly, for the i-th step, Thus bounding iteratively expression (48), it follows that finally giving On the other hand, for the ergodic case, let U(d) = U(d) × · · · × U(d) k+1 times be the k + 1 Cartesian product space of d-dimensional unitary groups, then we define ζ : Similarly as before, we now have and we may let the metric on U(d) be the 2-product metric δ U defined [57] by which then satisfies and thus we conclude that so the (bound on) Lipschitz constants coincide with η ≤

E.2.2 The concentration function
We now make use of a result related to the Gromov-Bishop inequality (see e.g. Theorem 7 in [55]) stating that if Ric(M ) ≥ Ric(Σ n (R)) = n−1 R 2 for an n-dimensional manifold M , where Σ n (R) is the n-dimensional sphere of radius R and Ric(X) is the infimum of diagonal elements of the Ricci curvature tensor on X [55], then the respective concentration functions satisfy For the time-independent case, the corresponding manifold is the group manifold U of U(d) (where here d = d E d S ), which is diffeomorphic to SU(d) × Σ 1 (1) and thus has Ric(U ) = d/2 and dim(U ) = d 2 . Then it follows that Ric(U ) ≥ Ric(Σ d 2 (r)) if r 2 ≥ 2(d 2 − 1)/d and hence, taking the minimal case, For the ergodic case, the corresponding manifold is the group manifold U of the k + 1 Cartesian product space U(d) × · · · × U(d). In this case Ric(U) = (k + 1)d/2 and dim(U) = (k + 1)d 2 . Then it follows that and hence, taking the minimal case, The factor dη −2 behaves as O(d E d −2k+1 S ), so both concentration functions will be small whenever d E d 2k S .

F Average purity
The purity of a quantum state is defined by the trace of its square and has an interpretation and relevance of its own in terms of correlations and/or indeterminacy of the state [41]. In terms of the Choi process tensor, it can offer insights about correlations and entanglement developed along a given process, as well as the noisiness generated. Aside from a direct interpretation of the purity, we compute its average over the Haar measure for the process tensor in order to obtain explicit consequences on the average non-Markovianity.

F.1 Ergodic case
The average purity of the Choi process tensor for k-steps in the ergodic case is given by We first notice that the simplest case k = 0, i.e. with no process occurring, reproduces the now well known [15,[36][37][38][39][40][41]. The average purity Eq. (5) decreases polynomially in d E and in the small subsystem limit satisfies thus leading to the purity of the maximally mixed state at a rate O(1/d E ). It is also manifestly decreasing in k and in the long time limit satisfies becoming independent of the dimension of the system. This limit corresponds strictly to a Choi process tensor with an infinite number of ancillary systems of dimension d S (making the purity bounded by zero from below) or more intuitively to a process after a large number of steps. The detail of the calculation is as follows.
With Θ = ρ ⊗ Ψ ⊗k , we first (following the approach in section 2 of [58]) write the trace as where Γ E (·) ≡ 1 E ⊗ tr E (·) is a CP map and can be expressed as an operator-sum, i.e. for any X acting on where the identities are on the ancillary system and K ij = |i j| ⊗ 1 S are the (ES system) Kraus operators with {|i } d E i=1 a given basis for E. And so we need to compute the integrals and in fact these are really in the ES part only, which has the form ω (ρ) where we dropped the Haar measure symbols just to economize notation. Summation over repeated indices x, y and i, j on the Kraus operators is implied in this expression and the ancillary part takes the form We may then evaluate the integrals from U 0 to U k with the result where sum is implied over e i 's, and so on, only relying on Eq. (67). First, let us state that Before carrying on, let us notice that the case ω k:k is special because here one evaluates the terms on the Kraus operators tr(K ij K † ij ) and tr(K ij ) tr(K † ij ) summed over their indices, which leaves Notice that we can get rid of the Kronecker deltas easily in the full Ω k:0 when summing over Greek indices, as δ αθ δ βλ give rise to maximally mixed states in the ancillas while terms δ αβ δ λθ give rise to identities, and there are only deltas of this kind. Thus we may simply assign δ αθ δ βλ → 1/d S and δ αβ δ λθ → 1 when plugging the corresponding expressions in Ω 2:0 6 . Furthermore, a direct consequence of this is that Ω k:0 will be a linear combination of 2 k+1 outer products between 1, ρ and Ψ, implying that (as ρ is pure, Ψ is idempotent with trace one and the trace of an outer product is the product of traces) the average purity E U i [tr(Υ 2 k:0 )] will be a sum of the scalar terms in ω (ρ) k:0 together with an extra overall factor of d − 1. Let us denote by the terms that appear in expression (70) after having taken δ αθ δ βλ → 1/d S and δ αβ δ λθ → 1. Doing similarly for the ω k:k case, we find that where the terms with a C factor are taken as zero for k = 1 and as C for k = 2. We may then simplify this expression by means of the geometric series to the right hand side of Eq. (5).

F.2 Time-independent case
We make use of Eq. (28) with the dimension argument of the Weingarten function taken implicitly as d E d S .
The average purity of the Choi process tensor for k steps in a time-independent process can be expressed with implicit sum over repeated indices, as where the notation in ρ is as in Eq. (34), and where now the symmetric group is on {0, . . . , 2k + 1}; we also make the definitions in Eq. (84) for∆ (d E ) ,∆ (d S ) , which are monomials in d E and d S correspondingly with degree determined by σ, τ . The number of elements in the symmetric group scale as (2k + 2)! and expressions like Eq. (6) are generally non-trivial (see also e.g. [41]); however, we can gain valuable insights from it by looking at their asymptotics. Overall, numerical sampling suggest that this average purity remains quite close to that of the ergodic case. For k = 0, with no process, we may also verify that we recover E U [tr(Υ 2 0:0 )] = (d E + d S )/(d E d S + 1) as expected [15,[36][37][38][39][40][41].
In the small subsystem limit, as detailed below in F.2.1, we obtain corresponding to the purity in the ergodic case in this limit. Similarly, in the long time limit, we obtain in F.2.2, and the average purity in this limit also coincides with that of ergodic processes. The detail of the calculation is as follows.
We may write the integral in the ES part, as in Eq.(65) now with same unitary (throughout we take sums over repeated indices and omit the Haar measure symbol), We decompose the unitaries in the whole ES space as U = U ab |a b| and U † = U * a b |b a | (i.e. the a, b labels here refer to the whole ES space) and we enumerate the labels of the ones to the left of ρ in Eq. (79) from 0 to k (priming, , the adjoint ones) and the remaining from k + 1 to 2k + 1, we also do this in increasing order for the adjoint unitaries and decreasing order for the original unitaries so that they match the order of the S operators, i.e. the unitary components will originally appear as U * i 0 j 0 · · · U * i k j k U i k j k · · · U i0j0 U * i k+1 j k+1 · · · U * i 2k+1 j 2k+1 U i 2k+1 j 2k+1 · · · U i k+1 j k+1 , and we then rearrange them to the form of Eq. (28) just keeping track of the correct order in the operator part, for which we introduce sets of E-basis into each S operator as S ab = |ea eb| labeled by e for the ones to the left of ρ in Eq. (79), priming the ones between adjoint unitaries, and by the corresponding ones to the right, also priming the ones between unitaries.