<EMAIL_ADDRESS>
# Almost Markovian processes from closed dynamics

Pedro Figueroa–Romero  [ Kavan Modi  Felix A. Pollock  School of Physics and
Astronomy, Monash University, Clayton, Victoria 3800, Australia

(August 27, 2024)

###### Abstract

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.

Suggested keywords

##  1 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 111 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].

Equilibration on average has been widely studied 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 $\mathcal{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.

$\longrightarrow$$e^{-iH\tau}$$t=0$$t=\tau$$E$$S$$E$$S$(a)

$\langle\mathcal{O}\rangle_{\rho(t)}-\langle\mathcal{O}\rangle_{\omega}$$t$$0$(b)

$\cdots$$\rho$${E}$${S}$$\mathcal{M}_{0}$$U_{1}$$\mathcal{M}_{1}$$U_{2}$$\cdots$$\mathcal{M}_{{}_{k-1}}$$U_{k}$\$\rho^{(S)}_{k}$(c)

Figure 1: Weak and strong equilibration. (a) A notion of subsystem
equilibration: $S$ is in a nonequilibrium state at time $0$ and evolves
jointly with $E$ according to Schrödinger’s equation until it reaches an
equilibrium state at some time $\tau$; more precisely, in (b) the expectation
$\langle\,\cdot\,\rangle$ of an observable $\mathcal{O}$ on the state
$\rho(t)$ differs only slightly from the corresponding one on the time-
averaged state $\omega$ for most times, even if it eventually deviates greatly
from it. In (c) we illustrate the generalization we consider in this work,
where an initial joint state $\rho$ evolves unitarily between the application
of measurements $\mathcal{M}_{i}$ on the system $S$, rendering the final
reduced state $\rho_{k}^{(S)}$; the dotted box encloses the process tensor
which describes the process independently of the (experimentally) controllable
operations.

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 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 many-body
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].

$\rho_{{}_{E}}^{(1)}$$\rho_{{}_{S}}$${E}$${S}$$\mathcal{M}_{0}$$U_{1}$\$\scriptstyle{\rho_{{}_{E}}^{(2)}}$$\mathcal{M}_{1}$$U_{2}$$\cdots$\$\scriptstyle{\rho_{{}_{E}}^{(k)}}$$\mathcal{M}_{{}_{k-1}}$$U_{k}$\$\rho^{(S)}_{k}$
Figure 2: Markovian dynamics. A quantum circuit corresponding to the Markov
approximation, which is mathematically equivalent to starting off with an
uncorrelated $\rho_{{}_{E}}^{(1)}\otimes\rho_{{}_{S}}$ state and then
artificially refreshing the environment (discarding and replacing with a new
state) at each step (cf. Fig. 1c)

Here, we give an initial answer to this deep foundational question. We show
that, for the class of evolutions 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 environment-
system ($ES$) with finite dimensional Hilbert space
$\mathcal{H}_{ES}=\mathbb{C}^{d_{E}}\otimes\mathbb{C}^{d_{S}}$ initialised in
a pure joint state $\rho$. The system $S$ is acted on with an initial
measurement operation $\mathcal{M}_{0}$ followed by a series of subsequent
measurement operations $\\{\mathcal{M}_{i}\\}_{i=1}^{k-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 $k$-step
process tensor is a completely positive causal 222 Meaning that the output of
a process tensor at time $k$ does not depend on deterministic inputs at later
times $K>k$. map $\mathcal{T}_{k:0}$ from the set
$\\{\mathcal{M}_{i}\\}_{i=0}^{k-1}$ to output states,
$\rho^{\prime}_{k}=\mathcal{T}_{k:0}[\\{\mathcal{M}_{i}\\}_{i=0}^{k-1}]$.

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 [23, 24, 21]. For a $k$-step process
tensor, the corresponding (normalized) Choi state $\Upsilon_{k:0}$ can be
obtained by swapping the system with one half of a maximally entangled state
$\Psi_{A_{i}B_{i}}\equiv|\psi^{+}\rangle\\!\langle\psi^{+}|\equiv\sum_{j,\ell=1}^{d_{S}}|jj\rangle\\!\langle{\ell\ell}|/d_{S}$
at each time $i$ [25]. The $2k$ ancillary systems
$A_{1},B_{1},\cdots{A}_{k},B_{k}$ have Hilbert space dimension
$\dim(\mathcal{H}_{A})=\dim(\mathcal{H}_{B})=d_{S}$ 333 The dimensions of each
$A_{i}$, $B_{i}$ must be fixed to $d_{S}$ to make the swap operator well
defined. 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

$\displaystyle\begin{split}&\Upsilon_{k:0}=\operatorname{tr}_{E}[\,\mathcal{U}_{k:0}(\rho\otimes\Psi^{\otimes{k}})\,\mathcal{U}^{\dagger}_{k:0}],\\\
\mbox{with}\quad&\mathcal{U}_{k:0}\equiv(U_{k}\otimes\mathds{1})\mathcal{S}_{k}\cdots(U_{1}\otimes\mathds{1})\mathcal{S}_{1}(U_{0}\otimes\mathds{1}),\end{split}$
(1)

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 $\mathcal{S}_{i}$ are swap operators between system $S$
and the ancillary system $A_{i}$. The corresponding quantum circuit is
depicted in Fig. 3.

$\rho$${E}$${S}$$U_{0}$$U_{1}$$\cdots$$U_{k}$\$\cdots$$\Psi_{A_{1}B_{1}}$$\cdots$$\Psi_{A_{2}B_{2}}$$\vdots$$\cdots$$\vdots$$\Psi_{A_{k}B_{k}}$$\Upsilon_{k:0}$$\times$$\mathcal{S}_{1}$$\times$$\times$$\mathcal{S}_{2}$$\times$$\times$$\mathcal{S}_{k}$$\times$$A_{1}$$B_{1}$$A_{2}$$B_{2}$$A_{k}$$B_{k}$
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 $\Upsilon$ is a many-body quantum state capturing all the
properties of the process.

Throughout the paper, we write $\Upsilon$ for the Choi state of a $k$-step
process unless otherwise specified.

##  3 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 $\Upsilon$
and that of a corresponding Markovian one
$\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}$. The latter must take the form
of a tensor product of quantum maps $\mathcal{E}_{i:i-1}$ connecting adjacent
pairs of time steps:
$\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}=\bigotimes_{i=1}^{k}\mathcal{E}_{i:i-1}$
[25].

In [5], a such measure of non-Markovianity was first introduced through the
relative entropy between a process and its corresponding Markovian one,
$\mathcal{R}(\Upsilon\|\Upsilon^{\scriptscriptstyle{(\mathrm{M})}})\equiv\operatorname{tr}[\Upsilon(\log\Upsilon-\log\Upsilon^{\scriptscriptstyle{(\mathrm{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 $\mathcal{D}$ as

$\displaystyle\mathcal{N}\equiv\min_{\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}}\mathcal{D}(\Upsilon,\Upsilon^{\scriptscriptstyle{(\mathrm{M})}})\equiv\frac{1}{2}\min_{\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}}\|\Upsilon-\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}\|_{1}$
(2)

where $\|X\|_{1}\equiv\operatorname{tr}\sqrt{XX^{\dagger}}$ is the trace norm
(or Schatten 1-norm); in particular, this is related to relative entropy
through the so-called quantum Pinsker inequality,
$\mathcal{R}(\Upsilon\|\Upsilon^{\scriptscriptstyle{(\mathrm{M})}})\geq
2\,\mathcal{D}^{2}(\Upsilon,\Upsilon^{\scriptscriptstyle{(\mathrm{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, 26, 16] (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.

##  4 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 non-vanishing 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 [29, 15, 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. 1.

Random interaction: All $U_{i}$ independently chosen,

  2. 2.

Constant interaction: $U_{i}=U_{j},\quad\forall\leq{i}\neq{j}\leq{k}$,

where $U_{i}$ represents the closed (unitary) $ES$ evolution between the
application of measurement operation $\mathcal{M}_{i-1}$ and
$\mathcal{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.

###  4.1 Main result

Our main result states that, for a randomly sampled $k$-step quantum process
$\Upsilon$ undergone by a $d_{S}$-dimensional subsystem of a larger
$d_{E}d_{S}$-dimensional composite, the probability for the non-Markovianity
$\mathcal{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 $\epsilon>0$,

$\displaystyle\mathbb{P}[\mathcal{N}\geq\mathcal{B}_{k}(d_{E},d_{S})+\epsilon]\leq\mathrm{e}^{-\mathcal{C}(d_{E},d_{S})\epsilon^{2}},$
(3)

where
$\mathcal{C}(d_{E},d_{S})=c\,d_{E}d_{S}\left(\frac{d_{S}-1}{d_{S}^{k+1}-1}\right)^{2}$
with $c=1/4$ for a constant interaction process and $c=(k+1)/4$ for a random
interaction process. The function $\mathcal{B}_{k}(d_{E},d_{S})$ is an upper
bound on the expected non-Markovianity $\mathds{E}[\mathcal{N}]$ whose details
depend on the way in which processes are sampled.

$\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}$$\Upsilon$$\epsilon$$\mathcal{B}_{k}$$\mathbb{E}[\mathcal{N}]$$\color[rgb]{0,0,0.8}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.8}{\mathcal{N}}$(a)

$\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}$$\Upsilon$(b)

Figure 4: Concentration of Markovian processes in large environments. (a) A
geometric _cartoon_ of our main result in a space of process tensors: the
probability of the non-Markovianity $\mathcal{N}$ of deviating from
$\mathcal{B}_{k}$ by some $\epsilon>0$ decreases exponentially in
$\epsilon^{2}$. In (b), quantum processes $\Upsilon$ on large dimensional
environments (such that $d_{E}\gg{d}_{S}^{2k+1}$) concentrate around the
Markovian ones $\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}$.

The proof is presented in full in Appendix E, and it can be outlined as
follows: we first bound the difference between $\mathcal{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 $\epsilon$ away from the expectation value
(which we upper-bound by $\mathcal{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 $\mathcal{B}_{k}+\epsilon$ and
$\mathrm{e}^{-\mathcal{C}\epsilon^{2}}$ are small; the latter is fulfilled in
the _small subsystem_ or _large environment_ limit, which in our setting means
$d_{E}\gg{d}_{S}^{2k+1}$. We may also state the minimal value for $\epsilon$
that, assuming the large environment limit, renders both sides small, i.e.,
such that $\epsilon^{2}d_{E}\gg 1\gg\epsilon$; this is fulfilled for
$\epsilon=d_{E}^{-1/3}$ 444Similar to [14], here we look for an $x>0$ such
that $\epsilon=d_{E}^{-x}$ and $\epsilon^{2}d_{E}=d_{E}^{x}$. A geometrical
cartoon to illustrate the result is presented in Fig. 4.

This is the case because the upper bound
$\mathcal{B}_{k}\geq\mathbb{E}(\mathcal{N})$ is given by

$\displaystyle\mathcal{B}_{k}(d_{E},d_{S})\equiv\begin{cases}\frac{\sqrt{d_{E}\,\mathbb{E}[\operatorname{tr}(\Upsilon^{2})]-x}+y}{2}&\text{if}\quad
d_{E}<d_{S}^{2k+1}\\\\[7.22743pt]
\frac{\sqrt{d_{S}^{2k+1}\mathbb{E}[\operatorname{tr}(\Upsilon^{2})]-1}}{2}&\text{if}\quad
d_{E}\geq{d}_{S}^{2k+1},\end{cases}$ (4)

with $x\equiv\frac{d_{E}}{d_{S}^{2k+1}}\left(1+y\right)$, $y\equiv
1-\frac{d_{E}}{d_{S}^{2k+1}}$. The size of $\mathcal{B}_{k}$ depends entirely
on the size of the average purity of the process tensor
$\mathbb{E}[\operatorname{tr}(\Upsilon^{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 $\operatorname{tr}(\Upsilon^{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 $\mathbb{E}[\Upsilon^{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

$\displaystyle\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon^{2})]=\frac{d_{E}^{2}-1}{d_{E}(d_{E}d_{S}+1)}\left(\frac{d_{E}^{2}-1}{d_{E}^{2}d_{S}^{2}-1}\right)^{k}+\frac{1}{d_{E}},$
(5)

for the random interaction picture and

$\displaystyle\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon^{2})]=d_{S}^{-2k}\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\,\varrho\,_{\tau}\,\Delta_{k,\sigma,\tau}^{(d_{E},d_{S})},$
(6)

in the constant interaction case, where $\mathcal{S}_{n}$ is the symmetric
group over $n$ elements, $\mathrm{Wg}$ is known as the Weingarten function
[30], $\varrho_{\tau}$ is a product of entries of $\rho$ depending on
permutations $\tau$ and $\Delta$ is a product, scaling with $k$, of monomials
in $d_{E}$ and $d_{S}$ depending on permutations $\sigma$ and $\tau$.

###  4.2 Limiting cases

In both the constant and random interaction cases, $\mathcal{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 $\mathrm{Wg}$ function (which is
intrinsic to the Haar-unitary averaging in the constant interaction picture).
However, due to results found in [30, 31], we can still study analytically the
behaviour of the bound $\mathcal{B}_{k}$ for both cases in the two following
limits.

In the large environment limit,

$\displaystyle\lim_{d_{E}\to\infty}\mathds{E}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{S}^{2k+1}},$
(7)

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:

$\displaystyle\lim_{d_{E}\to\infty}\mathds{E}[\mathcal{N}]=\lim_{d_{E}\to\infty}\mathcal{B}_{k}(d_{E},d_{S})=0.$
(8)

Furthermore, from Eq. (5) it is seen that it does so at a rate
$\mathcal{O}(1/d_{E})$ in the random interaction case.

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

$\displaystyle\lim_{k\to\infty}\mathds{E}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{E}},$
(9)

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

$\displaystyle\lim_{k\to\infty}\mathds{E}[\mathcal{N}]\leq\lim_{k\to\infty}\mathcal{B}_{k}(d_{E},d_{S})=1,$
(10)

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
$\mathcal{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

$\mathds{E}[\operatorname{tr}(\rho_{S}^{2})]=\frac{d_{E}+d_{S}}{d_{E}d_{S}+1},$
(11)

where
$\rho_{S}\equiv\Upsilon_{0:0}=\operatorname{tr}_{E}(U\rho\,U^{\dagger})$, as
found in [36, 37, 38, 39, 15, 40, 41]. This then leads to

$\displaystyle\mathds{E}\left[\mathcal{D}\left(\rho_{S},\frac{\mathds{1}}{d_{S}}\right)\right]\leq\frac{1}{2}\sqrt{\frac{d_{S}^{2}-1}{d_{E}d_{S}+1}},$
(12)

which goes as $\sqrt{d_{S}/4d_{E}}$ when $d_{E}\gg{d}_{S}$, and is also a
standard result [35].

###  4.3 Numerical examples

To support our results, we sampled Choi states $\Upsilon$ numerically in the
random interaction case and computed their corresponding average non-
Markovianity $\mathbb{E}[\mathcal{N}]$ 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 $\mathcal{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 $\mathcal{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}\simeq d_{S}^{2k+1}$, 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.

Figure 5: Average non-Markovianity $\mathbb{E}_{{}_{U_{i}}}[\mathcal{N}]$ of a
random interaction process for a qubit in the environment dimension $d_{E}$ at
fixed time steps $k$. Discrete values are shown for numerical averages over
$\lfloor 40/k\rfloor$ randomly generated process tensors $\Upsilon$ at time
steps $k=1,2,3$ and with fixed $d_{S}=2$; error bars denote the standard
deviation due to sampling error. The lines above each set of points denote the
upper bound $\mathcal{B}_{k}(d_{E},2)$. Process tensors were generated by
sampling Haar random unitaries according to [42].

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 $\mathcal{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\geq 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.

##  5 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 $\Upsilon_{k:0/\\{i\in[0,k-1]\\}}\equiv\Upsilon_{c}$, where a subset of
control operations $\\{\mathcal{M}_{i}\\}_{i\in[0,k-1]}$ are replaced by
identity operations (i.e., the system is simply left to evolve). Letting
$\mathcal{N}_{\text{coarse}}\equiv\min_{\Upsilon_{c}^{\scriptscriptstyle{(\mathrm{M})}}}\mathcal{D}(\Upsilon_{c},\Upsilon_{c}^{\scriptscriptstyle{(\mathrm{M})}})$,
we have

$\displaystyle\mathcal{N}_{\text{coarse}}\leq\mathcal{N},$ (13)

since the set of allowed $\Upsilon_{c}^{\scriptscriptstyle{(\mathrm{M})}}$
strictly contains the allowed $\Upsilon^{\scriptscriptstyle{(\mathrm{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
$\Upsilon-\Upsilon^{\scriptscriptstyle{(\mathrm{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\in\mathbb{M}$. In general, for an any set of measurements
$\mathbb{M}$ we can define a restricted measure of non-Markovianity detectable
with that set:
$\mathcal{D}_{\mathbb{M}}(\Upsilon,\Upsilon^{\scriptscriptstyle{(\mathrm{M})}})\equiv\max_{m\in\mathbb{M}}\frac{1}{2}|\operatorname{tr}[m(\Upsilon-\Upsilon^{\scriptscriptstyle{(\mathrm{M})}})]|\leq\mathcal{D}(\Upsilon,\Upsilon^{\scriptscriptstyle{(\mathrm{M})}})$,
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.

##  6 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 system 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.

We are grateful to Simon Milz and Andrea Collevecchio for valuable
discussions. PFR is supported by the Monash Graduate Scholarship (MGS) and the
Monash International Postgraduate Research Scholarship (MIPRS). KM is
supported through Australian Research Council Future Fellowship FT160100073.

## References

  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, “Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems,” \bibfield journal \bibinfo journal Rep. Prog. Phys.\ \textbf \bibinfo volume 79,\ \bibinfo pages 056001 (\bibinfo year 2016).
  * Goold  _et al._ [2016] J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, “The role of quantum information in thermodynamics–a topical review,” \bibfield journal \bibinfo journal J. Phys. A: Math. Theor.\ \textbf \bibinfo volume 49,\ \bibinfo pages 143001 (\bibinfo year 2016).
  * García-Pintos  _et al._ [2017] L. P. García-Pintos, N. Linden, A. S. L. Malabarba, A. J. Short, and A. Winter, “Equilibration time scales of physically relevant observables,” \bibfield journal \bibinfo journal Phys. Rev. X\ \textbf \bibinfo volume 7,\ \bibinfo pages 031027 (\bibinfo year 2017).
  * Srednicki [1999] M. Srednicki, “The approach to thermal equilibrium in quantized chaotic systems,” \bibfield journal \bibinfo journal J. Phys. A: Math. Gen.\ \textbf \bibinfo volume 32,\ \bibinfo pages 1163 (\bibinfo year 1999).
  * Pollock  _et al._ [2018a] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, “Operational markov condition for quantum processes,” \bibfield journal \bibinfo journal Phys. Rev. Lett.\ \textbf \bibinfo volume 120,\ \bibinfo pages 040405 (\bibinfo year 2018\natexlaba).
  * Arenz  _et al._ [2015] C. Arenz, R. Hillier, M. Fraas, and D. Burgarth, “Distinguishing decoherence from alternative quantum theories by dynamical decoupling,” \bibfield journal \bibinfo journal Phys. Rev. A\ \textbf \bibinfo volume 92,\ \bibinfo pages 022102 (\bibinfo year 2015).
  * Schlosshauer [2007] M. A. Schlosshauer, \emph \bibinfo title Decoherence and the Quantum-To-Classical Transition (Springer Berlin Heidelberg, 2007).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, “Dynamics of non-Markovian open quantum systems,” \bibfield journal \bibinfo journal Rev. Mod. Phys.\ \textbf \bibinfo volume 89,\ \bibinfo pages 015001 (\bibinfo year 2017).
  * Ciccarello [2017] F. Ciccarello, “Collision models in quantum optics,” \bibfield journal \bibinfo journal Quantum Meas. Quantum Metrol.\ \textbf \bibinfo volume 4,\ \bibinfo pages 53 (\bibinfo year 2017).
  * Spohn and Lebowitz [2007] H. Spohn and J. L. Lebowitz, “Irreversible thermodynamics for quantum systems weakly coupled to thermal reservoirs,” in \emph \bibinfo booktitle Advances in Chemical Physics, edited by S. A. Rice (John Wiley & Sons, Ltd, 2007) pp. 109–142.
  * Ledoux [2001] M. Ledoux, \emph \bibinfo title The Concentration of Measure Phenomenon, Mathematical surveys and monographs (American Mathematical Society, 2001).
  * Milman and Schechtman [1986] V. D. Milman and G. Schechtman, \emph \bibinfo title Asymptotic Theory of Finite Dimensional Normed Spaces, Lecture Notes in Mathematics No. 1200 (Springer-Verlag, 1986).
  * Boucheron  _et al._ [2013] S. Boucheron, G. Lugosi, and P. Massart, \emph \bibinfo title Concentration Inequalities: A Nonasymptotic Theory of Independence (OUP Oxford, 2013).
  * Popescu  _et al._ [2006] S. Popescu, A. J. Short, and A. Winter, “Entanglement and the foundations of statistical mechanics,” \bibfield journal \bibinfo journal Nat. Phys.\ \textbf \bibinfo volume 2,\ \bibinfo pages 754 (\bibinfo year 2006).
  * Li [2013] Y. Li, “Seminar 6 of selected topics in mathematical physics: Quantum information theory,” (year 2013)NoStop
  * Masanes  _et al._ [2013] author author L. Masanes, author A. J. Roncaglia, and author A. Acín, title title Complexity of energy eigenstates as a mechanism for equilibration, 10.1103/PhysRevE.87.032137 journal journal Phys. Rev. E volume 87, pages 032137 (year 2013)NoStop
  * Chiribella  _et al._ [2008] author author G. Chiribella, author G. M. D’Ariano, and author P. Perinotti, title title Quantum circuit architecture, 10.1103/PhysRevLett.101.060401 journal journal Phys. Rev. Lett. volume 101, pages 060401 (year 2008)NoStop
  * Chiribella  _et al._ [2009] author author G. Chiribella, author G. M. D’Ariano, and author P. Perinotti, title title Theoretical framework for quantum networks, 10.1103/PhysRevA.80.022339 journal journal Phys. Rev. A volume 80, pages 022339 (year 2009)NoStop
  * Costa and Shrapnel [2016] author author F. Costa and author S. Shrapnel, title title Quantum causal modelling, 10.1088/1367-2630/18/6/063032 journal journal New J. Phys. volume 18, pages 063032 (year 2016)NoStop
  * Modi [2012] author author K. Modi, title title Operational approach to open dynamics and quantifying initial correlations, 10.1038/srep00581 journal journal Sci. Rep. volume 2, pages 581 (year 2012)NoStop
  * Pollock  _et al._ [2018b] author author F. A. Pollock, author C. Rodríguez-Rosario, author T. Frauenheim, author M. Paternostro, and author K. Modi, title title Non-Markovian quantum processes: Complete framework and efficient characterization, 10.1103/PhysRevA.97.012127 journal journal Phys. Rev. A volume 97, pages 012127 (year 2018b)NoStop
  * Breuer  _et al._ [2016] author author H.-P. Breuer, author E.-M. Laine, author J. Piilo, and author B. Vacchini, title title Colloquium: Non-Markovian dynamics in open quantum systems, 10.1103/RevModPhys.88.021002 journal journal Rev. Mod. Phys. volume 88, pages 021002 (year 2016)NoStop
  * Choi [1975] author author M.-D. Choi, title title Completely positive linear maps on complex matrices, 10.1016/0024-3795(75)90075-0 journal journal Linear Algebra Its Appl. volume 10, pages 285 (year 1975)NoStop
  * Jamiołkowski [1972] author author A. Jamiołkowski, title title Linear transformations which preserve trace and positive semidefiniteness of operators, 10.1016/0034-4877(72)90011-0 journal journal Rep. Math. Phys. volume 3, pages 275 (year 1972)NoStop
  * Milz  _et al._ [2017a] author author S. Milz, author F. A. Pollock, and author K. Modi, title title An introduction to operational quantum dynamics, 10.1142/S1230161217400169 journal journal Open Syst. Inf. Dyn. volume 24, pages 1740016 (year 2017a)NoStop
  * Short and Farrelly [2012] author author A. J. Short and author T. C. Farrelly, title title Quantum equilibration in finite time, 10.1088/1367-2630/14/1/013063 journal journal New J. Phys. volume 14, pages 013063 (year 2012)NoStop
  * Gogolin  _et al._ [2011] author author C. Gogolin, author M. P. Müller, and author J. Eisert, title title Absence of thermalization in nonintegrable systems, 10.1103/PhysRevLett.106.040401 journal journal Phys. Rev. Lett. volume 106, pages 040401 (year 2011)NoStop
  * Gilchrist  _et al._ [2005] author author A. Gilchrist, author N. K. Langford, and author M. A. Nielsen, title title Distance measures to compare real and ideal quantum processes, 10.1103/PhysRevA.71.062310 journal journal Phys. Rev. A volume 71, pages 062310 (year 2005)NoStop
  * Zhang [2014] author author L. Zhang, title title Matrix integrals over unitary groups: An application of Schur-Weyl duality, https://arxiv.org/abs/1408.3782 journal journal arXiv:1408.3782 (year 2014)NoStop
  * Gu [2013] author author Y. Gu,  _title Moments of Random Matrices and. Weingarten Functions_ , http://hdl.handle.net/1974/8241 Master’s thesis, school Queen’s University, address Ontario, Canada (year 2013)NoStop
  * Collins and Śniady [2006] author author B. Collins and author P. Śniady, title title Integration with respect to the Haar measure on unitary, orthogonal and symplectic group, 10.1007/s00220-006-1554-3 journal journal Commun. Math. Phys. volume 264, pages 773 (year 2006)NoStop
  * Puchała and Miszczak [2017] author author Z. Puchała and author J. A. Miszczak, title title Symbolic integration with respect to the Haar measure on the unitary groups, 10.1515/bpasts-2017-0003 journal journal Bull. Pol. Acad. Sci. Tech. Sci. volume 65, pages 21 (year 2017)NoStop
  * Roberts and Yoshida [2017] author author D. A. Roberts and author B. Yoshida, title title Chaos and complexity by design, 10.1007/JHEP04(2017)121 journal journal J. High Energy Phys. volume 2017, pages 121 (year 2017)NoStop
  * Brandão  _et al._ [2016] author author F. Brandão, author A. W. Harrow, and author M. Horodecki, title title Local random quantum circuits are approximate polynomial-designs, 10.1007/s00220-016-2706-8 journal journal Commun. Math. Phys. volume 346, pages 397 (year 2016)NoStop
  * Gläßle [2013] author author T. Gläßle, @noop title Seminar 8 of selected topics in mathematical physics: Quantum information theory, howpublished http://www.mpmueller.net/seminar/talk8.pdf (year 2013)NoStop
  * Lubkin [1978] author author E. Lubkin, title title Entropy of an n-system from its correlation with a k-reservoir, 10.1063/1.523763 journal journal J. Math. Phys. volume 19, pages 1028 (year 1978)NoStop
  * Page [1993] author author D. N. Page, title title Average entropy of a subsystem, 10.1103/PhysRevLett.71.1291 journal journal Phys. Rev. Lett. volume 71, pages 1291 (year 1993)NoStop
  * Lloyd and Pagels [1988] author author S. Lloyd and author H. Pagels, title title Complexity as thermodynamic depth, 10.1016/0003-4916(88)90094-2 journal journal Ann. Phys. New York volume 188, pages 186 (year 1988)NoStop
  * Scott and Caves [2003] author author A. J. Scott and author C. M. Caves, title title Entangling power of the quantum baker’s map, 10.1088/0305-4470/36/36/308 journal journal J. Phys. A: Math. Gen. volume 36, pages 9553 (year 2003)NoStop
  * Giraud [2007] author author O. Giraud, title title Purity distribution for bipartite random pure states, 10.1088/1751-8113/40/49/f03 journal journal J. Phys. A: Math. Theor. volume 40, pages F1053 (year 2007)NoStop
  * Pasquale  _et al._ [2012] author author A. D. Pasquale, author P. Facchi, author V. Giovannetti, author G. Parisi, author S. Pascazio, and author A. Scardicchio, title title Statistical distribution of the local purity in a large quantum system, 10.1088/1751-8113/45/1/015308 journal journal J. Phys. A: Math. Theor. volume 45, pages 015308 (year 2012)NoStop
  * Mezzadri [2007] author author F. Mezzadri, title title How to generate random matrices from the classical compact groups, https://www.ams.org/journals/notices/200705/fea-mezzadri-web.pdf journal journal Notices of the AMS volume 54, pages 592 (year 2007)NoStop
  * Cotler  _et al._ [2017] author author J. Cotler, author N. Hunter-Jones, author J. Liu, and author B. Yoshida, title title Chaos, complexity, and random matrices, 10.1007/JHEP11(2017)048 journal journal J. High Energy Phys. volume 2017, pages 48 (year 2017)NoStop
  * Low [2009] author author R. A. Low, title title Large deviation bounds for k-designs, 10.1098/rspa.2009.0232 journal journal Proc. R. Soc. A volume 465, pages 3289 (year 2009)NoStop
  * Nakata  _et al._ [2017] author author Y. Nakata, author C. Hirche, author M. Koashi, and author A. Winter, title title Efficient quantum pseudorandomness with nearly time-independent Hamiltonian dynamics, 10.1103/PhysRevX.7.021006 journal journal Phys. Rev. X volume 7, pages 021006 (year 2017)NoStop
  * Milz  _et al._ [2017b] author author S. Milz, author F. Sakuldee, author F. A. Pollock, and author K. Modi, title title Kolmogorov extension theorem for (quantum) causal modelling and general probabilistic theories, https://arxiv.org/abs/1712.02589 journal journal arXiv:1712.02589 (year 2017b)NoStop
  * Anza  _et al._ [2018] author author F. Anza, author C. Gogolin, and author M. Huber, title title Eigenstate thermalization for degenerate observables, 10.1103/PhysRevLett.120.150603 journal journal Phys. Rev. Lett. volume 120, pages 150603 (year 2018)NoStop
  * Pollock and Modi [2018] author author F. A. Pollock and author K. Modi, title title Tomographically reconstructed master equations for any open quantum dynamics, 10.22331/q-2018-07-11-76 journal journal Quantum volume 2, pages 76 (year 2018)NoStop
  * Tamascelli  _et al._ [2018] author author D. Tamascelli, author A. Smirne, author S. F. Huelga, and author M. B. Plenio, title title Nonperturbative treatment of non-Markovian dynamics of open quantum systems, 10.1103/PhysRevLett.120.030402 journal journal Phys. Rev. Lett. volume 120, pages 030402 (year 2018)NoStop
  * Luchnikov  _et al._ [2019] author author I. A. Luchnikov, author S. V. Vintskevich, author H. Ouerdane, and author S. N. Filippov, title title Simulation complexity of open quantum dynamics: Connection with tensor networks, 10.1103/PhysRevLett.122.160401 journal journal Phys. Rev. Lett. volume 122, pages 160401 (year 2019)NoStop
  * Arias  _et al._ [2002] author author A. Arias, author A. Gheondea, and author S. Gudder, title title Fixed points of quantum operations, 10.1063/1.1519669 journal journal J. Math. Phys. volume 43, pages 5872 (year 2002)NoStop
  * Cramer [2012] author author M. Cramer, title title Thermalization under randomized local Hamiltonians, 10.1088/1367-2630/14/5/053051 journal journal New J. Phys. volume 14, pages 053051 (year 2012)NoStop
  * Weingarten [1978] author author D. Weingarten, title title Asymptotic behavior of group integrals in the limit of infinite rank, 10.1063/1.523807 journal journal J. Math. Phys. volume 19, pages 999 (year 1978)NoStop
  * Ginory and Kim [2016] author author A. Ginory and author J. Kim, title title Weingarten calculus and the IntHaar package for integrals over compact matrix groups, https://arxiv.org/abs/1612.07641 journal journal arXiv:1612.07641 (year 2016)NoStop
  * Taranto  _et al._ [2018] author author P. Taranto, author K. Modi, and author F. A. Pollock, title title Emergence of a fluctuation relation for heat in nonequilibrium Landauer processes, 10.1103/PhysRevE.97.052111 journal journal Phys. Rev. E volume 97, pages 052111 (year 2018)NoStop
  * Epstein and Whaley [2017] author author J. M. Epstein and author K. B. Whaley, title title Quantum speed limits for quantum-information-processing tasks, 10.1103/PhysRevA.95.042314 journal journal Phys. Rev. A volume 95, pages 042314 (year 2017)NoStop
  * Deza and Deza [2009] author author M. Deza and author E. Deza, 10.1007/978-3-642-00234-2 _title Encyclopedia of Distances_ , Encyclopedia of Distances (publisher Springer Berlin Heidelberg, year 2009)NoStop
  * Zhang and Xiang [2017] author author L. Zhang and author H. Xiang, title title Average entropy of a subsystem over a global unitary orbit of a mixed bipartite state, 10.1007/s11128-017-1570-6 journal journal Quantum Inf. Process. volume 16, pages 112 (year 2017)NoStop

##  Appendix A The process tensor

The Choi state of a process tensor with initial state $\rho$ is given by

$\displaystyle\Upsilon_{k:0}=\operatorname{tr}_{E}[\,\mathcal{U}_{k:0}(\rho\otimes\Psi^{\otimes{k}})\,\mathcal{U}^{\dagger}_{k:0}],$
(14)

with

$\displaystyle\mathcal{U}_{k:0}\equiv(U_{k}\otimes\mathds{1})\mathcal{S}_{k}\cdots(U_{1}\otimes\mathds{1})\mathcal{S}_{1}(U_{0}\otimes\mathds{1}),$
(15)

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

$\displaystyle\mathcal{S}_{i}\equiv\sum_{\alpha,\beta}\mathfrak{S}_{\alpha\beta}\otimes\mathds{1}_{{}_{A_{1}B_{1}\cdots{A}_{i-1}B_{i-1}}}\otimes|\beta\rangle\\!\langle\alpha|\otimes\mathds{1}_{{}_{B_{i}A_{i+1}B_{i+1}\cdots{A}_{k}B_{k}}}$

with
$\mathfrak{S}_{\alpha\beta}=\mathds{1}_{E}\otimes|\alpha\rangle\\!\langle\beta|$.
Notice that
$\mathfrak{S}_{\alpha\sigma}=\mathfrak{S}_{\sigma\alpha}^{\dagger}$,
$\mathfrak{S}_{ab}\mathfrak{S}^{\dagger}_{cd}=\delta_{bd}\mathfrak{S}_{ac}$
and $\operatorname{tr}(\mathfrak{S}_{ab})=d_{E}\delta_{ab}$. Also,
$\Upsilon\in\mathcal{S}(\mathcal{H}_{S}\bigotimes_{i=1}^{k}\mathcal{H}_{A_{i}B_{i}})$,
i.e. the resulting Choi process tensor state belongs to the $S$-ancillary
system, which has dimension $d_{S}^{2k+1}$.

##  Appendix 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)

$\displaystyle\mathcal{N}\leq{D}(\Upsilon,d_{S}^{-(2k+1)}\mathds{1}).$ (16)

We may further bound this by considering the following cases separately.

1\. Case $d_{E}<d_{S}^{2k+1}$: We notice that555One may see this by looking at
the pure state
${|\Phi\rangle\\!\langle\Phi|\equiv\mathcal{U}(\Theta\otimes\Psi^{\otimes{k}})\mathcal{U}^{\dagger}}$
with Schmidt decomposition
${|\Phi\rangle=\sum_{i=1}^{n}\sqrt{\lambda_{i}}|e_{i}s_{i}\rangle}$ where
$n=\min(d_{E},d_{S}^{2k+1})$. $\mathrm{rank}(\Upsilon)\leq{d_{E}}$. Letting
$\gamma$ be the diagonal matrix of up to $d_{E}$ non-vanishing eigenvalues
$\lambda_{\gamma_{i}}$ of the Choi state, we may write

$\displaystyle\|\Upsilon-\frac{\mathds{1}}{d_{S}^{2k+1}}\|_{1}=\sum_{i=1}^{d_{E}}\left|\lambda_{\gamma_{i}}-\frac{1}{d_{S}^{2k+1}}\right|+\sum_{j=d_{E}+1}^{d_{S}^{2k+1}}\left|-\frac{1}{d_{S}^{2k+1}}\right|=\|\gamma-\frac{\mathds{1}_{E}}{d_{S}^{2k+1}}\|_{1}+1-\frac{d_{E}}{d_{S}^{2k+1}},$
(17)

where $|\cdot|$ denotes the standard absolute value, so using the inequality
${\|X\|_{1}\leq\sqrt{\mathrm{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}=\sqrt{\operatorname{tr}(XX^{\dagger})}$
is the Schatten 2-norm,

$\displaystyle\|\Upsilon-\frac{\mathds{1}}{d_{S}^{2k+1}}\|_{1}$
$\displaystyle\leq\sqrt{d_{E}}\|\gamma-\frac{\mathds{1}_{E}}{d_{S}^{2k+1}}\|_{2}+1-\frac{d_{E}}{d_{S}^{2k+1}}$
$\displaystyle=\sqrt{d_{E}\operatorname{tr}[\Upsilon^{2}]+\frac{d^{2}_{E}}{d_{S}^{4k+2}}-\frac{2d_{E}}{d_{S}^{2k+1}}}+1-\frac{d_{E}}{d_{S}^{2k+1}},$
(18)

Furthermore, applying Jensen’s inequality (in particular for the square-root,
$\mathbb{E}[\sqrt{X}]\leq\sqrt{\mathbb{E}[X]}$) with expectation over
evolution (either in the ergodic or time-independent case), this gives

$\displaystyle\mathbb{E}\left[\mathcal{N}\right]\leq\frac{1}{2}\left(\sqrt{d_{E}\mathbb{E}[\operatorname{tr}(\Upsilon^{2})]+\frac{d_{E}^{2}}{d_{S}^{4k+2}}-\frac{2d_{E}}{d_{S}^{2k+1}}}+1-\frac{d_{E}}{d_{S}^{2k+1}}\right).$
(19)

2\. Case $d_{E}\geq{d}_{S}^{2k+1}$: This case is a small subsystem limit for
most $k$. Directly applying ${\|X\|_{1}\leq\sqrt{\mathrm{dim}(X)}\|X\|_{2}}$
as before,

$\displaystyle\|\Upsilon-\frac{\mathds{1}}{d_{S}^{2k+1}}\|_{1}\leq\sqrt{d_{S}^{2k+1}}\|\Upsilon-\frac{\mathds{1}}{d_{S}^{2k+1}}\|_{2}=\sqrt{d_{S}^{2k+1}\operatorname{tr}[\Upsilon^{2}]-1},$
(20)

Similarly, taking the average over evolution, by means of Jensen’s inequality,

$\displaystyle\mathbb{E}\left[\mathcal{N}\right]$
$\displaystyle\leq\frac{1}{2}\sqrt{d_{S}^{2k+1}\mathbb{E}\left[\operatorname{tr}(\Upsilon^{2})\right]-1}.$
(21)

This proves that the function

$\displaystyle\mathcal{B}_{k}(d_{E},d_{S})\equiv\begin{cases}\frac{\sqrt{d_{E}\,\mathbb{E}[\operatorname{tr}(\Upsilon^{2})]-x}+y}{2}&\text{if}\quad
d_{E}<d_{S}^{2k+1}\\\\[7.22743pt]
\frac{\sqrt{d_{S}^{2k+1}\mathbb{E}[\operatorname{tr}(\Upsilon^{2})]-1}}{2}&\text{if}\quad
d_{E}\geq{d}_{S}^{2k+1},\end{cases}$ (22)

with $x\equiv\frac{d_{E}}{d_{S}^{2k+1}}\left(1+y\right)$, $y\equiv
1-\frac{d_{E}}{d_{S}^{2k+1}}$, as given in the main text, provides an upper
bound on the average non-Markovianity $\mathbb{E}[\mathcal{N}]$.

##  Appendix C Haar random averages

###  C.1 The Haar measure

The Haar measure is the unique probability measure $\mu$ on any of the
orthogonal, unitary or symplectic groups that is both left and right
invariant, for our purposes meaning that
$\int{f}(VU)\,d\mu(U)=\int{f}(UV)\,d\mu(U)=\int{f}(U)\,d\mu(U)$ for all
unitaries $V$, with normalization $\int{d}\mu(U)=1$. A random pure state can
thus be generated by taking $|\phi\rangle\in\mathcal{H}$ with
$d=\dim(\mathcal{H})$, sampling a unitary matrix $U\in\mathrm{U}(d)$ according
to the Haar measure and letting
$|\phi\rangle\\!\langle\phi|\equiv\phi\mapsto{U}\phi{U}^{\dagger}$. Its
average over the Haar measure is then given by

$\displaystyle\mathbb{E}_{{}_{H}}(\phi)\equiv\int\limits_{\mathrm{U}(d)}U\phi\,{U}^{\dagger}\,d\mu(U),$
(23)

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 $\mathbb{E}_{{}_{U_{i}}}$ and
$\mathbb{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

$\displaystyle\Phi^{(k)}(X)\mapsto\int\limits_{U(d)}{U}^{\otimes{k}}X({U}^{\otimes{k}})^{\dagger}\,d\mu(U),$
(24)

where $X\in\mathcal{H}^{\otimes{k}}$; in particular, the case
$U^{\dagger}{A}UX{U}^{\dagger}{B}U$ 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, $\Phi^{(k)}$ commutes
with all $V^{\otimes{k}}$ such that $V\in\mathrm{U}(d)$, and by such property,
the Schur-Weyl duality assures that also

$\displaystyle\Phi^{(k)}(X)=\sum_{\sigma\in\mathcal{S}_{k}}f_{\tau}(X)\mathcal{P}_{\tau},$
(25)

where the sum is over the symmetric group $\mathcal{S}_{k}$ on $k$ symbols,
$\mathcal{P}_{\tau}$ is a permutation operator defined by
$\mathcal{P}_{\tau}|x_{1},\ldots,x_{k}\rangle=|x_{\tau(1)},\ldots,x_{\tau(k)}\rangle$
and $f_{\tau}$ is a linear function of $X$ [33]. The $k=1$ case is well known
[15] and in general for any operator acting on $\mathds{C}^{d}$ results in a
completely depolarizing channel,

$\displaystyle\Phi^{(1)}(X)=\mathbb{E}({X})=\frac{\operatorname{tr}(X)}{d}\mathds{1}.$
(26)

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],

$\displaystyle\int\limits_{\mathrm{U}(d)}\prod_{\ell=1}^{k}U_{i_{\ell}{j}_{\ell}}U_{i^{\prime}_{\ell}{j}^{\prime}_{\ell}}^{*}\,d\mu(U)=\sum_{\sigma,\tau\in\mathcal{S}_{k}}\prod_{\ell=1}^{k}\delta_{i_{\ell}{i}^{\prime}_{\sigma(\ell)}}\delta_{j_{\ell}{j}^{\prime}_{\tau(\ell)}}\,\mathrm{Wg}(\tau\sigma^{-1},d),$
(27)

where $U_{ij}$ is the $ij$ component of $U\in\mathrm{U}(d)$ and $\mathrm{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 [30, 29]
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. (27) and the Weingarten function (a
subsequent version for Maple is also now available in [54]).

Here we obtain the expected states $\mathbb{E}(\Upsilon)$ and expected
purities $\mathbb{E}[\operatorname{tr}(\Upsilon^{2})]$ over the Haar measure
of Choi process tensor states $\Upsilon$ as defined by Eq. (14). We compute
these averages for both the ergodic $\mathbb{E}_{{}_{U_{i}}}$ and time-
independent $\mathbb{E}_{{}_{U}}$ cases. The time-independent averages are
written in terms of the $\mathrm{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. (27) boils down to that of the $\mathrm{Wg}$
function and in [30] it has been shown that

$\displaystyle\mathrm{Wg}(\sigma\in\mathcal{S}_{n},d)\sim\frac{1}{d^{2n-\\#\sigma}},\,\,\text{as}\,\,d\to\infty,$
(28)

as a refinement of a result in [31], where $\\#\sigma$ is the number of cycles
of the permutation $\sigma$ counting also _fixed points_ (assignments from an
element to itself, $\sigma(x)=x$).

##  Appendix D Average Processes

###  D.1 Ergodic average process Choi state

The ergodic average Choi state of a $k$ step process is given by the maximally
mixed state,

$\displaystyle\mathbb{E}_{{}_{U_{i}}}(\Upsilon_{k:0})=\frac{1}{d_{S}^{2k+1}}\mathds{1}_{SA_{1}B_{1}\cdots{A}_{k}B_{k}}.$
(29)

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
$\operatorname{\mathfrak{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,

$\displaystyle\Upsilon_{k:0}$
$\displaystyle=\operatorname{tr}_{E}\left[(U_{k}\otimes\mathds{1})\mathcal{S}_{k}\cdots(U_{1}\otimes\mathds{1})\mathcal{S}_{1}(U_{0}\otimes\mathds{1})(\rho\otimes\Psi^{\otimes{k}})(U_{0}^{\dagger}\otimes\mathds{1})\mathcal{S}_{1}^{\dagger}(U_{1}^{\dagger}\otimes\mathds{1})\cdots\mathcal{S}_{k}^{\dagger}(U_{k}^{\dagger}\otimes\mathds{1})\right]$
$\displaystyle=\operatorname{tr}_{E}\sum\left[{U}_{k}\operatorname{\mathfrak{S}}_{\alpha_{k}\beta_{k}}\cdots{U}_{1}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U_{0}\rho{U}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k}\delta_{k}}^{\dagger}{U}_{k}^{\dagger}\right]$
$\displaystyle\hskip
18.06749pt\otimes[(|\beta_{1}\rangle\\!\langle\alpha_{1}|\otimes\mathds{1}_{B_{1}})\Psi(|\gamma_{1}\rangle\\!\langle\delta_{1}|\otimes\mathds{1}_{B_{1}})\otimes\cdots\otimes(|\beta_{k}\rangle\\!\langle\alpha_{k}|\otimes\mathds{1}_{B_{k}})\Psi(|\gamma_{k}\rangle\\!\langle\delta_{k}|\otimes\mathds{1}_{B_{k}})],$
(30)

where the sum is over repeated indices (Greek letters), hence by repeatedly
evaluating the 1-fold twirl (26), with $d=d_{E}d_{S}$ (also, the unitary
groups are implicitly of dimension $d_{E}d_{S}$), we get

$\displaystyle\mathbb{E}_{{}_{U_{i}}}(\Upsilon_{k:0})=\operatorname{tr}_{E}\sum\int\limits_{U_{k}\cdots{U}_{0}}U_{k}\operatorname{\mathfrak{S}}_{\alpha_{k}\beta_{k}}\cdots{U}_{1}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U_{0}\rho{U}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k}\delta_{k}}^{\dagger}{U}_{k}^{\dagger}$
$\displaystyle\hskip
90.3375pt\otimes[(|\beta_{1}\rangle\\!\langle\alpha_{1}|\otimes\mathds{1}_{B_{1}})\Psi(|\gamma_{1}\rangle\\!\langle\delta_{1}|\otimes\mathds{1}_{B_{1}})\otimes\cdots\otimes(|\beta_{k}\rangle\\!\langle\alpha_{k}|\otimes\mathds{1}_{B_{k}})\Psi(|\gamma_{k}\rangle\\!\langle\delta_{k}|\otimes\mathds{1}_{B_{k}})]$
$\displaystyle=\frac{1}{d_{S}^{k}d}\operatorname{tr}_{E}\sum\int\limits_{U_{k-1}\cdots{U}_{0}}\langle{e}_{k}\beta_{k}|U_{k-1}\operatorname{\mathfrak{S}}_{\alpha_{k-1}\beta_{k-1}}\cdots{U}_{1}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U_{0}\rho{U}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k-1}\delta_{k-1}}^{\dagger}{U}_{k-1}^{\dagger}|e_{k}\delta_{k}\rangle$
$\displaystyle\hskip
90.3375pt\mathds{1}_{S}\otimes|\beta_{1}\alpha_{1}\cdots\beta_{k}\alpha_{k}\rangle\\!\langle\delta_{1}\gamma_{1}\cdots\delta_{k}\alpha_{k}|$
$\displaystyle=\frac{d_{E}}{d_{S}^{k}d^{2}}\operatorname{tr}_{E}\sum\int\limits_{U_{k-2}\cdots{U}_{0}}\langle{e}_{k-1}\beta_{k-1}|U_{k-2}\operatorname{\mathfrak{S}}_{\alpha_{k-2}\beta_{k-2}}\cdots{U}_{1}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U_{0}\rho{U}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k-2}\delta_{k-2}}^{\dagger}{U}_{k-2}^{\dagger}|e_{k-1}\delta_{k-1}\rangle$
$\displaystyle\hskip
90.3375pt\mathds{1}_{S}\otimes|\beta_{1}\alpha_{1}\cdots\beta_{k-1}\alpha_{k-1}\rangle\\!\langle\delta_{1}\gamma_{1}\cdots\delta_{k-1}\alpha_{k-1}|\otimes\mathds{1}_{A_{k}B_{k}}$
$\displaystyle\vdots$
$\displaystyle=\frac{d_{E}^{k+1}}{d_{S}^{k}d^{k+1}}\mathds{1}_{SA_{1}B_{1}\cdots{A}_{k}B_{k}}$
$\displaystyle=\frac{1}{d_{S}^{2k+1}}\mathds{1}_{SA_{1}B_{1}\cdots{A}_{k}B_{k}},$
(31)

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_{i}^{(\prime)}\rangle\\}_{s_{i}=1}^{d_{S}}$
of $S$ system bases for $i=0,\ldots,k$ as

$\displaystyle\mathbb{E}_{{}_{U}}(\Upsilon_{k:0})=\frac{1}{d_{S}^{k}}\sum_{\sigma,\tau\in\mathcal{S}_{k+1}}\rho_{{}_{\tau(0);0}}\mathrm{Wg}(\tau\sigma^{-1})\Delta_{k,\sigma,\tau}^{(d_{E})}|s_{\sigma(k)}\rangle\\!\langle{s}_{k}|\bigotimes_{j=1}^{k}|s_{\sigma(j-1)}s^{\prime}_{\tau(j)}\rangle\langle{s}_{j-1}s^{\prime}_{j}|,$
(32)

with implicit sum over all repeated basis ($s^{(\prime)}_{i}$) indices, where
here $\mathcal{S}_{k+1}$ is the symmetric group on $\\{0,\ldots,k\\}$, and
with the definitions

$\displaystyle\rho_{{}_{\tau(0);0}}$
$\displaystyle=\langle{e}^{\prime}_{\tau(0)}s^{\prime}_{\tau(0)}|\rho|e^{\prime}_{0}s^{\prime}_{0}\rangle,$
(33) $\displaystyle\Delta_{k,\sigma,\tau}^{(d_{E})}$
$\displaystyle=\delta_{e_{\sigma(k)}e_{k}}\prod_{\ell=1}^{k}\delta_{e_{\sigma(\ell-1)}e^{\prime}_{\tau(\ell)}}\delta_{e_{\ell-1}e^{\prime}_{\ell}},$
(34)

where $\\{|e^{(\prime)}_{i}\rangle\\}_{e_{i}=1}^{d_{E}}$, with $i=0,\ldots,k$,
also summed over all its elements, is a set of environment bases and the
$\Delta$ term is simply a monomial in $d_{E}$ with degree determined by
$\sigma$ and $\tau$.

The case $k=0$ recovers
$\mathbb{E}_{{}_{U}}(\Upsilon_{0:0})=\mathds{1}_{S}/d_{S}$ as expected, as no
process occurs.

In the case of the superchannel $k=1$, we get Eq. (40), which is quite
different from the maximally mixed state that arises in the ergodic case,
although it converges to it as $d_{E}\to\infty$. We notice that
$\operatorname{tr}[\mathbb{E}_{{}_{U}}(\Upsilon_{1:0})]=1$ as expected and one
may also verify the purity of this average state to be given by Eq. (41) which
is in fact _close_ to that of the maximally mixed state for all $d_{E}$,
$d_{S}$ and $1/d_{S}\leq\operatorname{tr}(\rho_{S}^{2})\leq 1$, and converges
to it when $d_{E}\to\infty$ at a $\mathcal{O}(1/d_{E})$ rate.

The detail of the calculation is as follows.

By definition,

$\displaystyle\mathbb{E}_{{}_{U}}(\Upsilon_{k:0})=\frac{\operatorname{tr}_{E}}{d_{S}^{k}}\sum\int\limits_{U}U\operatorname{\mathfrak{S}}_{\alpha_{k}\beta_{k}}\cdots{U}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U\rho{U}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k}\delta_{k}}^{\dagger}{U}^{\dagger}\,d\mu(U)\otimes|\beta_{1}\alpha_{1}\ldots\beta_{k}\alpha_{k}\rangle\\!\langle\delta_{1}\gamma_{1}\ldots\delta_{k}\gamma_{k}|.$
(35)

We now consider directly decomposing the unitaries as
$U=\sum{U}_{ab}|a\rangle\\!\langle{b}|$ and
$U^{\dagger}=\sum{U}^{*}_{a^{\prime}{b}^{\prime}}|b^{\prime}\rangle\\!\langle{a}^{\prime}|$
–notice that the $a$ and $b$ labels refer to the whole $ES$ space–,
introducing an $E$ basis in the $\operatorname{\mathfrak{S}}$ operators as
$\operatorname{\mathfrak{S}}_{ab}=\sum|ea\rangle\\!\langle{eb}|$, and then
evaluating the $k+1$ moments of $\mathrm{U}(d_{E}d_{S})$ by means of Eq. (27),

$\displaystyle\int\limits_{U}$ $\displaystyle
U\operatorname{\mathfrak{S}}_{\alpha_{k}\beta_{k}}\cdots{U}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U\rho{U}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k}\delta_{k}}^{\dagger}{U}^{\dagger}\,d\mu(U)$
$\displaystyle=\sum\int\limits_{U}U_{i_{0}j_{0}}\cdots{U}_{i_{k}j_{k}}U^{*}_{i^{\prime}_{0}j^{\prime}_{0}}\cdots{U}^{*}_{i^{\prime}_{k}j^{\prime}_{k}}\,d\mu(U)\,|i_{k}\rangle\\!\langle{j}_{0}|\rho|j^{\prime}_{0}\rangle\\!\langle{i}^{\prime}_{k}|\,\prod_{\ell=1}^{k}\delta_{j_{\ell}(e\alpha)_{\ell}}\delta_{(e\beta)_{\ell}{i}_{\ell-1}}\delta_{j^{\prime}_{\ell}(e^{\prime}\gamma)_{\ell}}\delta_{(e^{\prime}\delta)_{\ell}{i^{\prime}}_{\ell-1}}$
$\displaystyle=\sum\sum_{\sigma,\tau\in\mathcal{S}_{k+1}}\langle{j}_{0}|\rho|j^{\prime}_{0}\rangle\delta_{i_{0}i^{\prime}_{\sigma(0)}}\delta_{j_{0}j^{\prime}_{\tau(0)}}\prod_{\ell=1}^{k}\delta_{i_{\ell}{i}^{\prime}_{\sigma(\ell)}}\delta_{j_{\ell}{j}^{\prime}_{\tau(\ell)}}\delta_{j_{\ell}(e\alpha)_{\ell}}\delta_{(e\beta)_{\ell}{i}_{\ell-1}}\delta_{j^{\prime}_{\ell}(e^{\prime}\gamma)_{\ell}}\delta_{(e^{\prime}\delta)_{\ell}{i^{\prime}}_{\ell-1}}\,\mathrm{Wg}(\tau\sigma^{-1})\,|i_{k}\rangle\\!\langle{i}^{\prime}_{k}|$
$\displaystyle=\sum\sum_{\sigma,\tau\in\mathcal{S}_{k+1}}\langle{j}^{\prime}_{\tau(0)}|\rho|j^{\prime}_{0}\rangle\delta_{i_{k}i^{\prime}_{\sigma(k)}}\prod_{\ell=1}^{k}\delta_{(e\beta)_{\ell}{i}^{\prime}_{\sigma(\ell-1)}}\delta_{(e^{\prime}\delta)_{\ell}{i^{\prime}}_{\ell-1}}\delta_{(e\alpha)_{\ell}{j}^{\prime}_{\tau(\ell)}}\delta_{(e^{\prime}\gamma)_{\ell}{j}^{\prime}_{\ell}}\,\mathrm{Wg}(\tau\sigma^{-1})\,|i_{k}\rangle\\!\langle{i}^{\prime}_{k}|,$
(36)

where here $\mathcal{S}_{k+1}$ denotes the symmetric group on
$\\{0,\ldots,k\\}$, and which taking $i\to\epsilon\varsigma$ and
$j\to\epsilon^{\prime}\varsigma^{\prime}$ to recover each $E$ and $S$ part
explicitly, turns into

$\displaystyle\int\limits_{U}$ $\displaystyle
U\operatorname{\mathfrak{S}}_{\alpha_{k}\beta_{k}}\cdots{U}\operatorname{\mathfrak{S}}_{\alpha_{1}\beta_{1}}U\rho{U}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma_{1}\delta_{1}}^{\dagger}{U}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{\gamma_{k}\delta_{k}}^{\dagger}{U}^{\dagger}\,d\mu(U)$
$\displaystyle=\sum\sum_{\sigma,\tau\in\mathcal{S}_{k+1}}\langle\epsilon^{\prime}_{\tau(0)}\varsigma^{\prime}_{\tau(0)}|\rho|\epsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle\prod_{\ell=1}^{k}\delta_{\epsilon_{\sigma(\ell-1)}\epsilon^{\prime}_{\tau(\ell)}}\delta_{\epsilon_{\ell-1}\epsilon^{\prime}_{\ell}}\delta_{\beta_{\ell}\varsigma_{\sigma(\ell-1)}}\delta_{\delta_{\ell}\varsigma_{\ell-1}}\delta_{\alpha_{\ell}\varsigma^{\prime}_{\tau(\ell)}}\delta_{\gamma_{\ell}\varsigma^{\prime}_{\ell}}\,\mathrm{Wg}(\tau\sigma^{-1})\,|\epsilon_{\sigma(k)}\varsigma_{\sigma(k)}\rangle\\!\langle\epsilon_{k}\varsigma_{k}|,$
(37)

and thus

$\displaystyle\mathbb{E}_{{}_{U}}(\Upsilon_{k:0})=\frac{\sum}{d_{S}^{k}}$
$\displaystyle\sum_{\sigma,\tau\in\mathcal{S}_{k+1}}\langle\epsilon^{\prime}_{\tau(0)}\varsigma^{\prime}_{\tau(0)}|\rho|\epsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle\mathrm{Wg}(\tau\sigma^{-1})\delta_{\epsilon_{\sigma(k)}\epsilon_{k}}\prod_{\ell=1}^{k}\delta_{\epsilon_{\sigma(\ell-1)}\epsilon^{\prime}_{\tau(\ell)}}\delta_{\epsilon_{\ell-1}\epsilon^{\prime}_{\ell}}|\varsigma_{\sigma(k)}\varsigma_{\sigma(0)}\varsigma^{\prime}_{\tau(1)}\cdots\varsigma_{\sigma(k-1)}\varsigma^{\prime}_{\tau(k)}\rangle\langle\varsigma_{k}\varsigma_{0}\varsigma^{\prime}_{1}\cdots\varsigma_{k-1}\varsigma^{\prime}_{k}|,$
(38)

as stated by Eq.(32).

####  D.2.1 Superchannel case

For the superchannel case, $k=1$, we have $\mathcal{S}_{2}=\\{(0,1),(0)(1)\\}$
where the elements are permutations stated in _cycle_ notation, representing
the assignments $(0,1)=\begin{cases}0\to 1\\\ 1\to 0\end{cases}$ and
$(0)(1)=\bm{1}^{2}=\begin{cases}0\to 0\\\ 1\to 1\end{cases}$, then (we write
$\alpha,\beta,\gamma,\delta$ for
$\varsigma_{0},\varsigma_{1},\varsigma_{0}^{\prime},\varsigma_{1}^{\prime}$,
respectively to ease the notation)

$\displaystyle\mathbb{E}_{{}_{U}}(\Upsilon_{1:0})=\frac{\sum}{d_{S}}\bigg{\\{}$
$\displaystyle\left[\langle\delta|\rho_{S}|\gamma\rangle|\alpha\beta\gamma\rangle\\!\langle\beta\alpha\delta|+d_{E}^{2}|\beta\alpha\delta\rangle\\!\langle\beta\alpha\delta|\right]\mathrm{Wg}[\bm{1}^{2}]+d_{E}\left[\langle\delta|\rho_{S}|\gamma\rangle|\beta\alpha\gamma\rangle\\!\langle\beta\alpha\delta|+|\alpha\beta\delta\rangle\\!\langle\beta\alpha\delta|\right]\mathrm{Wg}[(0,1)]\bigg{\\}},$
(39)

so that with $\mathrm{Wg}[\bm{1}^{2},d]=\frac{1}{d^{2}-1}$ and
$\mathrm{Wg}[(0,1),d]=\frac{-1}{d(d^{2}-1)}$ [30] we get

$\displaystyle\mathbb{E}_{{}_{U}}(\Upsilon_{1:0})=\frac{1}{d_{E}^{2}d_{S}^{2}-1}\bigg{[}\frac{d_{E}^{2}}{d_{S}}\mathds{1}_{SAB}+\frac{\operatorname{\text{{swap}}}}{d_{S}}\otimes\rho_{{}_{S}}^{\intercal}-\frac{\operatorname{\text{{swap}}}}{d_{S}}\otimes\frac{\mathds{1}_{B}}{d_{S}}-\frac{\mathds{1}_{SA}}{d_{S}^{2}}\otimes\rho_{{}_{S}}^{\intercal}\bigg{]},$
(40)

where $\operatorname{\text{{swap}}}=\sum_{i,j}|ij\rangle\\!\langle{j}i|$ is
the usual swap operator, and hence for the corresponding purity one may verify
that

$\displaystyle\operatorname{tr}[\mathbb{E}_{{}_{U}}(\Upsilon_{1:0})^{2}]=\frac{2}{(d_{E}^{2}d_{S}^{2}-1)^{2}}\left[\frac{1}{d_{S}^{3}}+\operatorname{tr}(\rho_{S}^{2})\frac{d_{S}^{2}-d_{S}-1}{2d_{S}^{2}}-\frac{d_{E}^{2}}{d_{S}}+\frac{d_{E}^{4}d_{S}}{2}\right].$
(41)

####  D.2.2 Small subsystem limit

When looking at the limit $d_{E}\to\infty$, the term that does not vanish is
the one with $\sigma,\tau=\bm{1}^{k+1}$, i.e. with both permutations being
identities, as these generate the most numerator powers in $d_{E}$ via the
$\delta$ terms (correspondingly, $\Delta_{k,\sigma,\tau}^{(d_{E})}$ in
Eq.(34)) when summed over all $e_{i}$’s, i.e.

$\displaystyle\sum_{\begin{subarray}{c}e^{(\prime)}_{i}=1\\\
e^{\prime}\neq{e}^{\prime}_{0}\end{subarray}}^{d_{E}}\delta_{\epsilon_{k}\epsilon_{k}}\prod_{\ell=1}^{k}\delta_{\epsilon_{\ell-1}\epsilon^{\prime}_{\ell}}\delta_{\epsilon_{\ell-1}\epsilon^{\prime}_{\ell}}=d_{E}^{k+1},$
(42)

all other permutations will vanish because of the $d_{E}$ powers in the
denominator generated by the $\mathrm{Wg}$ function, the least powers produced
by it are those when $\sigma\tau^{-1}=\mathbf{1}^{k+1}$ because
$\\#\bm{1}^{n}=\\#[(1)(2)\cdots(n)]=n$, i.e. identity produces the greatest
number of cycles, being the number of fixed points. Finally
$\sum_{\epsilon^{\prime}_{0}\varsigma^{\prime}_{0}}\langle\epsilon^{\prime}_{0}\varsigma^{\prime}_{0}|\rho|\epsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle=\operatorname{tr}\rho=1$,
and thus

$\displaystyle\mathbb{E}_{{}_{U}}(\Upsilon_{k:0})\sim\frac{d_{E}^{k+1}}{d_{S}^{k}}\mathrm{Wg}(\bm{1}^{k+1})\,\mathds{1}_{SA_{1}B_{1}\ldots{A}_{k}B_{k}}=\frac{1}{d_{S}^{2k+1}}\,\mathds{1}_{SA_{1}B_{1}\ldots{A}_{k}B_{k}},\hskip
7.22743pt\text{when}\,\,d_{E}\to\infty,$ (43)

matching the average over an ergodic process.

##  Appendix E Proof of Main Result

###  E.1 Preamble: Concentration of Measure (Levy’s lemma)

####  E.1.1 Lipschitz functions

A function $f:X\to{Y}$ between metric spaces $(X,\delta_{{}_{X}})$ and
$(Y,\delta_{{}_{Y}})$ is $L$-Lipschitz if there is a real constant $\eta\geq
0$, known as the Lipschitz constant, such that

$\displaystyle\delta_{{}_{Y}}\left(f(x_{1}),f(x_{2})\right)\leq{L}\,\delta_{{}_{X}}(x_{1},x_{2}).$
(44)

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\neq{j}$) and the ergodic case ($U_{i}\neq{U}_{j}$ for all steps
$i\neq{j}$).

####  E.1.2 Levy’s lemma

Let $M$ be a manifold with metric $\delta_{{}_{M}}$ and probability measure
$\mu_{{}_{M}}$, and let $f:M\to\mathbb{R}$ be an $\eta$-Lipschitz function,
then for all $\varepsilon>0$,

$\displaystyle\mathbb{P}_{x\sim\mu_{{}_{M}}}[f(x)\geq\mathbb{E}_{x}(f)+\varepsilon]\leq\alpha_{{}_{M}}(\varepsilon/\eta),$
(45)

where $\alpha_{{}_{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 make use of Levy’s lemma and the upper-bound
$\mathbb{E}[\mathcal{N}]\leq\mathcal{B}_{k}(d_{E},d_{S})$ for the average non-
Markovianity over the unitary group with the Haar measure, to prove our main
result; specifically here we compute the corresponding Lipschitz constants and
concentration functions for the function
$D(\Upsilon,d_{S}^{-(2k+1)}\mathds{1})$. We highlight that the evolution for
$\Upsilon$ is sampled from the Haar measure, or equivalently one may replace
$\mathbb{P}$ with something like $\mathbb{P}_{\Upsilon\sim\mu_{{}_{U_{i}}}}$.

####  E.2.1 The Lipschitz constants

We take $\tau:\mathrm{U}(d_{E}d_{S})\to\mathbb{R}$ defined by
$\tau(U)=D(\Upsilon(U),\frac{\mathds{1}}{d_{S}^{2k+1}})=\frac{1}{2}\|\Upsilon(U)-\frac{\mathds{1}}{d_{S}^{2k+1}}\|_{1}$
for the time-independent case, where we explicitly mean
$\Upsilon(U)=\operatorname{tr}_{E}[\mathcal{U}_{k:0}\Theta\mathcal{U}_{k:0}^{\dagger}]$
as in definition (14) with $\Theta=\rho\otimes\Psi^{\otimes{k}}$ a pure state,
then for all $V\in\mathrm{U}(d_{E}d_{S})$,

$\displaystyle|\tau({U})-\tau({V})|$
$\displaystyle=\left|D(\Upsilon(U),\frac{\mathds{1}}{d_{S}^{2k+1}})-D(\Upsilon(V),\frac{\mathds{1}}{d_{S}^{2k+1}})\right|$
$\displaystyle\leq{D}(\Upsilon(U),\Upsilon(V))$
$\displaystyle\leq{D}(\mathcal{U}_{k:0}\Theta\mathcal{U}_{k:0}^{\dagger},\mathcal{V}_{k:0}\Theta\mathcal{V}_{k:0}^{\dagger}).$
(46)

Consider for now different labelings for the unitary at each time step. By the
triangle inequality,

$\displaystyle\|\mathcal{U}_{k:0}\Theta\mathcal{U}_{k:0}^{\dagger}-\mathcal{V}_{k:0}\Theta\mathcal{V}_{k:0}^{\dagger}\|_{1}$
$\displaystyle\leq\|\mathcal{U}_{k:0}\Theta\mathcal{U}_{k:0}^{\dagger}-U_{k}\mathcal{S}_{k}\mathcal{V}_{k-1:0}\Theta\mathcal{V}_{k-1:0}^{\dagger}\mathcal{S}_{k}^{\dagger}{U}_{k}^{\dagger}\|_{1}+\|U_{k}\mathcal{S}_{k}\mathcal{V}_{k-1:0}\Theta\mathcal{V}_{k-1:0}^{\dagger}\mathcal{S}_{k}^{\dagger}{U}_{k}^{\dagger}-\mathcal{V}_{k:0}\Theta\mathcal{V}_{k:0}^{\dagger}\|_{1}$
$\displaystyle=\|\mathcal{U}_{k-1:0}\Theta\mathcal{U}_{k-1:0}^{\dagger}-\mathcal{V}_{k-1:0}\Theta\mathcal{V}_{k-1:0}^{\dagger}\|_{1}+\|U_{k}\mathcal{S}_{k}\mathcal{V}_{k-1:0}\Theta\mathcal{V}_{k-1:0}^{\dagger}\mathcal{S}_{k}^{\dagger}{U}_{k}^{\dagger}-\mathcal{V}_{k:0}\Theta\mathcal{V}_{k:0}^{\dagger}\|_{1},$
(47)

where it is clear by context that $U$ stands for $U\otimes\mathds{1}$.

Now we use Lemma 1 of [56], which states that

$\displaystyle\|A\sigma{A}^{\dagger}-B\sigma{B}^{\dagger}\|_{1}\leq
2\|A-B\|_{2}$ (48)

for two unitaries $A$, $B$ and any $\sigma$. For the simplest case with $k=1$,
this gives

$\displaystyle\|U_{1}\mathcal{S}_{1}V_{0}\Theta{V}_{0}^{\dagger}\mathcal{S}_{1}^{\dagger}{U}_{1}^{\dagger}-V_{1}\mathcal{S}_{1}V_{0}\Theta{V}_{0}^{\dagger}\mathcal{S}_{1}^{\dagger}{V}_{1}^{\dagger}\|_{1}=\|\frac{\sum}{d_{S}}(U_{1}\operatorname{\mathfrak{S}}_{\alpha\beta}V_{0}\rho{V}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma\delta}U_{1}^{\dagger}-{V}_{1}\operatorname{\mathfrak{S}}_{\alpha\beta}V_{0}\rho{V}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma\delta}V_{1}^{\dagger})\otimes|\beta\alpha\rangle\\!\langle\gamma\delta|\|_{1}$
$\displaystyle\leq\frac{\sum_{\alpha\beta}}{d_{S}}\operatorname{tr}\sqrt{\sum_{\gamma\delta}(U_{1}\operatorname{\mathfrak{S}}_{\alpha\beta}V_{0}\rho{V}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma\delta}U_{1}^{\dagger}-{V}_{1}\operatorname{\mathfrak{S}}_{\alpha\beta}V_{0}\rho{V}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\gamma\delta}V_{1}^{\dagger})(U_{1}\operatorname{\mathfrak{S}}_{\beta\alpha}V_{0}\rho{V}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\delta\gamma}U_{1}^{\dagger}-{V}_{1}\operatorname{\mathfrak{S}}_{\beta\alpha}V_{0}\rho{V}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{\delta\gamma}V_{1}^{\dagger})}$
$\displaystyle\leq\frac{\sum_{\alpha\beta}}{d_{S}}\|U_{1}\operatorname{\mathfrak{S}}_{\alpha\beta}V_{0}\rho{V}_{0}^{\dagger}(\sum\operatorname{\mathfrak{S}}_{\gamma\delta})U_{1}^{\dagger}-{V}_{1}\operatorname{\mathfrak{S}}_{\alpha\beta}V_{0}\rho{V}_{0}^{\dagger}(\sum\operatorname{\mathfrak{S}}_{\gamma^{\prime}\delta^{\prime}})V_{1}^{\dagger}\|_{1}$
$\displaystyle\leq\frac{2}{d_{S}}\sum_{\alpha,\beta=1}^{d_{S}}\|U_{1}-V_{1}\|_{2}=2d_{S}\|U_{1}-V_{1}\|_{2},$
(49)

and doing similarly, for the $i$-th step,

$\displaystyle\|U_{i}\mathcal{S}_{i}\mathcal{V}_{i-1:0}\Theta\mathcal{V}_{i-1:0}^{\dagger}\mathcal{S}_{i}^{\dagger}{U}_{i}^{\dagger}-\mathcal{V}_{i:0}\Theta\mathcal{V}_{i:0}^{\dagger}\|_{1}$
$\displaystyle\leq 2d_{S}^{i}\|U_{i}-V_{i}\|.$ (50)

Thus bounding iteratively expression (47), it follows that

$\displaystyle\|\mathcal{U}_{k:0}\Theta\mathcal{U}_{k:0}^{\dagger}-\mathcal{V}_{k:0}\Theta\mathcal{V}_{k:0}^{\dagger}\|_{1}$
$\displaystyle\leq 2\sum_{\ell=0}^{k}d_{S}^{\ell}\|U_{\ell}-V_{\ell}\|_{2},$
(51)

finally giving

$\displaystyle|\tau(U)-\tau(V)|$
$\displaystyle\leq\left(\frac{d_{S}^{k+1}-1}{d_{S}-1}\right)\|U-V\|_{2}.$ (52)

On the other hand, for the ergodic case, let
$\mathbb{U}(d)=\underbrace{\mathrm{U}(d)\times\cdots\times\mathrm{U}(d)}_{k+1\,\text{times}}$
be the $k+1$ Cartesian product space of $d$-dimensional unitary groups, then
we define $\zeta:\mathbb{U}(d_{E}d_{S})\to\mathbb{R}$ by
$\zeta(\vec{U})=D(\Upsilon(\vec{U}),\frac{\mathds{1}}{d_{S}^{2k+1}})$ where
now $\vec{U}=(U_{0},\cdots,U_{k})$. Similarly as before, we now have

$\displaystyle|\zeta(\vec{U})-\zeta(\vec{V})|\leq\sum_{\ell=0}^{k}d_{S}^{\ell}\|U_{\ell}-V_{\ell}\|_{2},$
(53)

and we may let the metric on $\mathbb{U}(d)$ be the 2-product metric
$\delta_{\mathbb}{U}$ defined [57] by

$\displaystyle\delta_{\mathbb}{U}(\vec{x},\vec{y})=\sqrt{\sum_{\ell=0}^{k}\|x_{\ell}-y_{\ell}\|_{2}^{2}},$
(54)

which then satisfies

$\displaystyle\sum_{\ell=0}^{k}d_{S}^{\ell}\|U_{\ell}-V_{\ell}\|_{2}$
$\displaystyle\leq\sum_{\ell=0}^{k}d_{S}^{\ell}\sqrt{\sum_{\ell^{\prime}=0}^{k}\|U_{\ell^{\prime}}-V_{\ell^{\prime}}\|_{2}^{2}}=\left(\frac{d_{S}^{k+1}-1}{d_{S}-1}\right)\delta_{\mathbb}{U}(\vec{U},\vec{V}),$
(55)

and thus we conclude that

$\displaystyle|\zeta(\vec{U})-\zeta(\vec{V})|\leq\left(\frac{d_{S}^{k+1}-1}{d_{S}-1}\right)\delta_{\mathbb}{U}(\vec{U},\vec{V}),$
(56)

so the (bound on) Lipschitz constants coincide with
$\eta\leq\frac{d_{S}^{k+1}-1}{d_{S}-1}$, which essentially behaves as
$\mathcal{O}(d_{S}^{k})$.

####  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
$\mathrm{Ric}(M)\geq\mathrm{Ric}(\Sigma^{n}(R))=\frac{n-1}{R^{2}}$ for an
$n$-dimensional manifold $M$, where $\Sigma^{n}(R)$ is the $n$-dimensional
sphere of radius $R$ and $\mathrm{Ric}(X)$ is the infimum of diagonal elements
of the Ricci curvature tensor on $X$ [55], then the respective concentration
functions satisfy
$\alpha_{M}(x)\leq\alpha_{\Sigma^{n}(R)}(x)\leq\exp\left[-\frac{x^{2}(n-1)}{2R^{2}}\right]$
[11].

For the time-independent case, the corresponding manifold is the group
manifold $\mathscr{U}$ of $\mathrm{U}(d)$ (where here $d=d_{E}d_{S}$), which
is diffeomorphic to $\mathrm{SU}(d)\times\Sigma^{1}(1)$ and thus has
$\mathrm{Ric}(\mathscr{U})=d/2$ and $\dim(\mathscr{U})=d^{2}$. Then it follows
that $\mathrm{Ric}(\mathscr{U})\geq\mathrm{Ric}(\Sigma^{d^{2}}(r))$ if
$r^{2}\geq 2(d^{2}-1)/d$ and hence, taking the minimal case,

$\displaystyle\alpha_{{}_{\mathscr{U}}}(\varepsilon/\eta)\leq\alpha_{\Sigma^{d^{2}}(r)}(\varepsilon/\eta)\leq\exp\left(-\frac{\varepsilon^{2}d}{4\eta^{2}}\right).$
(57)

For the ergodic case, the corresponding manifold is the group manifold
$\mathbf{U}$ of the $k+1$ Cartesian product space
$\mathrm{U}(d)\times\cdots\times\mathrm{U}(d)$. In this case
$\mathrm{Ric}(\mathbf{U})=(k+1)d/2$ and $\dim(\mathbf{U})=(k+1)d^{2}$. Then it
follows that
$\mathrm{Ric}(\mathbf{U})\geq\mathrm{Ric}(\Sigma^{(k+1)d^{2}}(\mathcal{R}))$
if $\mathcal{R}^{2}\geq 2[(k+1)d^{2}-1]/[d(k+1)]$ and hence, taking the
minimal case,

$\displaystyle\alpha_{{}_{\mathbf{U}}}(\varepsilon/\eta)\leq\alpha_{\Sigma^{(k+1)d^{2}}(\mathcal{R})}(\varepsilon/\eta)\leq\exp\left(-\frac{\varepsilon^{2}(k+1)d}{4\eta^{2}}\right).$
(58)

The factor $d\eta^{-2}$ behaves as $\mathcal{O}(d_{E}d_{S}^{-2k+1})$, so both
concentration functions will be small whenever $d_{E}\gg{d}_{S}^{2k}$.

##  Appendix 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

$\displaystyle\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon^{2})]=\frac{d_{E}^{2}-1}{d_{E}(d_{E}d_{S}+1)}\left(\frac{d_{E}^{2}-1}{d_{E}^{2}d_{S}^{2}-1}\right)^{k}+\frac{1}{d_{E}}.$
(59)

We first notice that the simplest case $k=0$, i.e. with no process occurring,
reproduces the now well known result
$\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon_{0:0}^{2})]=(d_{E}+d_{S})/(d_{E}d_{S}+1)$
found in [36, 37, 38, 39, 15, 40, 41].

The average purity Eq. (5) decreases polynomially in $d_{E}$ and in the small
subsystem limit satisfies

$\displaystyle\lim_{d_{E}\to\infty}\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{S}^{2k+1}},$
(60)

thus leading to the purity of the maximally mixed state at a rate
$\mathcal{O}(1/d_{E})$. It is also manifestly decreasing in $k$ and in the
long time limit satisfies

$\displaystyle\lim_{k\to\infty}\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{E}},$
(61)

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 $\Theta=\rho\otimes\Psi^{\otimes{k}}$, we first (following the approach
in section 2 of [58]) write the trace as

$\displaystyle\operatorname{tr}(\Upsilon_{k:0}^{2})$
$\displaystyle=\operatorname{tr}\left[\left(\operatorname{tr}_{E}(\mathcal{U}_{k:0}\Theta\,\mathcal{U}^{\dagger}_{k:0})\right)^{2}\right]$
$\displaystyle=\operatorname{tr}\left[\operatorname{tr}_{E}\left((\mathds{1}_{E}\otimes\Upsilon_{k:0})\,\mathcal{U}_{k:0}\Theta\,\mathcal{U}^{\dagger}_{k:0}\right)\right]$
$\displaystyle=\operatorname{tr}\left[\Gamma_{{}_{E}}(\mathcal{U}_{k:0}\Theta\,\mathcal{U}^{\dagger}_{k:0})\,\mathcal{U}_{k:0}\Theta\,\mathcal{U}^{\dagger}_{k:0}\right]$
$\displaystyle=\operatorname{tr}\left[\Theta\,\mathcal{U}_{k:0}^{\dagger}\Gamma_{{}_{E}}(\mathcal{U}_{k:0}\Theta\,\mathcal{U}^{\dagger}_{k:0})\,\mathcal{U}_{k:0}\right],$
(62)

where
$\Gamma_{{}_{E}}(\cdot)\equiv\mathds{1}_{E}\otimes\operatorname{tr}_{E}(\cdot)$
is a CP map and can be expressed as an operator-sum, i.e. for any
$\mathcal{X}$ acting on $ESA_{1}B_{1}\cdots{A}_{k}B_{k}$ we have
$\Gamma_{{}_{E}}(\mathcal{X})=\sum_{i,j=1}^{d_{E}}(K_{ij}\otimes\mathds{1})\mathcal{X}(K_{ij}\otimes\mathds{1})^{\dagger}$
where the identities are on the ancillary system and
$K_{ij}=|i\rangle\\!\langle{j}|\otimes\mathds{1}_{S}$ are the ($ES$ system)
Kraus operators with $\\{|i\rangle\\}_{i=1}^{d_{E}}$ a given basis for $E$.
And so we need to compute the integrals

$\displaystyle\Omega_{k:0}\equiv\int\limits_{U_{k}\cdots\,U_{0}}\mathcal{U}_{k:0}^{\dagger}\Gamma_{E}(\mathcal{U}_{k:0}\Theta\,\mathcal{U}^{\dagger}_{k:0})\,\mathcal{U}_{k:0}\,d\mu(U_{0})\cdots\,d\mu(U_{k}),$
(63)

and in fact these are really in the $ES$ part only, which has the form

$\displaystyle\omega_{k:0}^{(\rho)}\equiv\int\limits_{U}U_{0}^{\dagger}\operatorname{\mathfrak{S}}_{x_{1}\alpha_{1}}U_{1}^{\dagger}\cdots{\operatorname{\mathfrak{S}}}_{x_{k}\alpha_{k}}U_{k}^{\dagger}{K}_{\tilde{\imath}\tilde{\jmath}}U_{k}\operatorname{\mathfrak{S}}_{\beta_{k}x_{k}}\cdots{U_{1}}\operatorname{\mathfrak{S}}_{\beta_{1}x_{1}}U_{0}\rho{U}_{0}^{\dagger}\operatorname{\mathfrak{S}}_{y_{1}\lambda_{1}}U_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{y_{k}\lambda_{k}}U_{k}^{\dagger}{K}^{\dagger}_{\tilde{\imath}\tilde{\jmath}}{U}_{k}\operatorname{\mathfrak{S}}_{\theta_{k}y_{k}}\cdots{U}_{1}\operatorname{\mathfrak{S}}_{\theta_{1}y_{1}}U_{0},$
(64)

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

$\displaystyle\Omega_{k:0}^{(AB)}$
$\displaystyle=(|\alpha_{1}\rangle\\!\langle\beta_{1}|\otimes\mathds{1}_{B_{1}}\otimes\cdots\otimes|\alpha_{k}\rangle\\!\langle\beta_{k}|\otimes\mathds{1}_{B_{k}})\Psi^{\otimes\,k}(|\lambda_{1}\rangle\\!\langle\theta_{1}|\otimes\mathds{1}_{B_{1}}\otimes\cdots\otimes|\lambda_{k}\rangle\\!\langle\theta_{k}|\otimes\mathds{1}_{B_{k}}).$
(65)

We may then evaluate the integrals from $U_{0}$ to $U_{k}$ with the result

$\displaystyle\int\limits_{\mathrm{U}(d)}U^{\dagger}{A}UX{U}^{\dagger}{B}U\,d\mu(U)=\frac{d\operatorname{tr}(AB)-\operatorname{tr}(A)\operatorname{tr}(B)}{d(d^{2}-1)}\operatorname{tr}(X)\mathds{1}+\frac{d\operatorname{tr}(A)\operatorname{tr}(B)-\operatorname{tr}(AB)}{d(d^{2}-1)}X,$
(66)

also, here $d=d_{E}d_{S}$, then we have

$\displaystyle\omega_{k:0}^{(\rho)}=\frac{1}{d(d^{2}-1)}$
$\displaystyle\bigg{(}\delta_{y_{1}x_{1}}\delta_{x_{1}y_{1}}\int\limits_{U_{k}\cdots{U}_{1}}\operatorname{tr}\left[\operatorname{\mathfrak{S}}_{\theta_{1}\alpha_{1}}U_{1}^{\dagger}\cdots{\operatorname{\mathfrak{S}}}_{x_{k}\alpha_{k}}U_{k}^{\dagger}{K}_{ij}U_{k}\operatorname{\mathfrak{S}}_{\beta_{k}x_{k}}\cdots{U}_{1}\operatorname{\mathfrak{S}}_{\beta_{1}\lambda_{1}}U_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{y_{k}\lambda_{k}}U_{k}^{\dagger}{K}^{\dagger}_{ij}{U}_{k}\operatorname{\mathfrak{S}}_{\theta_{k}y_{k}}\cdots{U}_{1}\right](d\mathds{1}-\rho)$
$\displaystyle+\delta_{x_{1}x_{1}}\delta_{y_{1}y_{1}}\int\limits_{U_{k}\cdots{U}_{1}}\operatorname{tr}\left[\operatorname{\mathfrak{S}}_{\beta_{1}\alpha_{1}}U_{1}^{\dagger}\cdots{\operatorname{\mathfrak{S}}}_{x_{k}\alpha_{k}}U_{k}^{\dagger}{K}_{ij}U_{k}\operatorname{\mathfrak{S}}_{\beta_{k}x_{k}}\cdots{U}_{1}\right]\operatorname{tr}\left[\operatorname{\mathfrak{S}}_{\theta_{1}\lambda_{1}}U_{1}^{\dagger}\cdots\operatorname{\mathfrak{S}}_{y_{k}\lambda_{k}}U_{k}^{\dagger}{K}^{\dagger}_{ij}{U}_{k}\operatorname{\mathfrak{S}}_{\theta_{k}y_{k}}\cdots{U}_{1}\right](d\rho-\mathds{1})\bigg{)}$
$\displaystyle=\frac{1}{d(d^{2}-1)}\bigg{[}d_{S}\langle{e}_{1}\alpha_{1}|\omega_{k:1}^{(\operatorname{\mathfrak{S}}_{\beta_{1}\lambda_{1}})}|e_{1}\theta_{1}\rangle(d\mathds{1}-\rho)+d_{S}^{2}\langle{e}_{1}\alpha_{1}|\omega_{k:1}^{(|e_{1}\beta_{1}\rangle\\!\langle{e}^{\prime}_{1}\lambda_{1}|)}|e^{\prime}_{1}\theta_{1}\rangle(d\rho-\mathds{1})\bigg{]},$
(67)

where sum is implied over $e_{i}$’s, and so on, only relying on Eq. (66).
First, let us state that

$\displaystyle\omega_{k:i-1}^{(X)}$
$\displaystyle=\frac{1}{d(d^{2}-1)}\bigg{[}d_{S}\langle{e}_{i}\alpha_{i}|\omega_{k:i}^{(\operatorname{\mathfrak{S}}_{\beta_{i}\lambda_{i}})}|e_{i}\theta_{i}\rangle(d\operatorname{tr}(X)\mathds{1}-X)+d_{S}^{2}\langle{e}_{i}\alpha_{i}|\omega_{k:i}^{(|e_{i}\beta_{i}\rangle\\!\langle{e}^{\prime}_{i}\lambda_{i}|)}|e^{\prime}_{i}\theta_{i}\rangle(dX-\operatorname{tr}(X)\mathds{1})\bigg{]},$
(68)

for $k>i-1\geq 0$. Then we get

$\displaystyle\omega_{k:0}^{(\rho)}=\frac{1}{d^{2}(d^{2}-1)^{2}}\bigg{\\{}$
$\displaystyle
d_{S}\langle{e}_{2}\alpha_{2}|\omega_{k:2}^{(\operatorname{\mathfrak{S}}_{\beta_{2}\lambda_{2}})}|e_{2}\theta_{2}\rangle\bigg{[}d_{S}(dd_{E}^{2}\delta_{\alpha_{1}\theta_{1}}\delta_{\beta_{1}\lambda_{1}}-d_{E}\delta_{\alpha_{1}\beta_{1}}\delta_{\lambda_{1}\theta_{1}})(d\mathds{1}-\rho)+d_{S}^{2}(dd_{E}\delta_{\alpha_{1}\theta_{1}}\delta_{\beta_{1}\lambda_{1}}-d_{E}^{2}\delta_{\alpha_{1}\beta_{1}}\delta_{\lambda_{1}\theta_{1}})(d\rho-\mathds{1})\bigg{]}$
$\displaystyle+d_{S}^{2}\langle{e}_{2}\alpha_{2}|\omega_{k:2}^{(|e_{2}\beta_{2}\rangle\\!\langle{e}_{2}^{\prime}\lambda_{2}|)}|e_{2}^{\prime}\theta_{2}\rangle\bigg{[}d_{S}(dd_{E}\delta_{\alpha_{1}\beta_{1}}\delta_{\lambda_{1}\theta_{1}}-d_{E}^{2}\delta_{\alpha_{1}\theta_{1}}\delta_{\beta_{1}\lambda_{1}})(d\mathds{1}-\rho)+d_{S}^{2}(dd_{E}^{2}\delta_{\alpha_{1}\beta_{1}}\delta_{\lambda_{1}\theta_{1}}-d_{E}\delta_{\alpha_{1}\theta_{1}}\delta_{\beta_{1}\lambda_{1}})(d\rho-\mathds{1})\bigg{]}\bigg{\\}}.$
(69)

Before carrying on, let us notice that the case $\omega_{k:k}$ is special
because here one evaluates the terms on the Kraus operators
$\operatorname{tr}(K_{ij}K_{ij}^{\dagger})$ and
$\operatorname{tr}(K_{ij})\operatorname{tr}(K_{ij}^{\dagger})$ summed over
their indices, which leaves

$\displaystyle\langle{e}_{k}\alpha_{k}|\omega_{k:k}^{(\operatorname{\mathfrak{S}}_{\lambda_{k}\beta_{k}})}|e_{k}\theta_{k}\rangle=\frac{1}{d(d^{2}-1)}[dd_{E}^{2}(dd_{E}\delta_{\alpha_{k}\theta_{k}}\delta_{\beta_{k}\lambda_{k}}-\delta_{\alpha_{k}\beta_{k}}\delta_{\lambda_{k}\theta_{k}})+d^{2}d_{E}(d_{S}\delta_{\alpha_{k}\beta_{k}}\delta_{\lambda_{k}\theta_{k}}-\delta_{\alpha_{k}\theta_{k}}\delta_{\beta_{k}\lambda_{k}})],$
(70)
$\displaystyle\langle{e}_{k}\alpha_{k}|\omega_{k:k}^{(|e_{k}\beta_{k}\rangle\\!\langle{e}^{\prime}_{k}\lambda_{k}|)}|e^{\prime}_{k}\theta_{k}\rangle=\frac{1}{d(d^{2}-1)}[dd_{E}^{3}(d_{S}\delta_{\alpha_{k}\theta_{k}}\delta_{\beta_{k}\lambda_{k}}-\delta_{\alpha_{k}\beta_{k}}\delta_{\lambda_{k}\theta_{k}})+d^{2}(dd_{E}\delta_{\alpha_{k}\beta_{k}}\delta_{\lambda_{k}\theta_{k}}-\delta_{\alpha_{k}\theta_{k}}\delta_{\beta_{k}\lambda_{k}})].$
(71)

Notice that we can get rid of the Kronecker deltas easily in the full
$\Omega_{k:0}$ when summing over Greek indices, as
$\delta_{\alpha\theta}\delta_{\beta\lambda}$ give rise to maximally mixed
states in the ancillas while terms
$\delta_{\alpha\beta}\delta_{\lambda\theta}$ give rise to identities, and
there are only deltas of this kind. Thus we may simply assign
$\delta_{\alpha\theta}\delta_{\beta\lambda}\to 1/d_{S}$ and
$\delta_{\alpha\beta}\delta_{\lambda\theta}\to 1$ when plugging the
corresponding expressions in $\Omega_{2:0}$666 In particular this implies that
the terms
$d_{S}\delta_{\alpha\theta}\delta_{\beta\lambda}-\delta_{\alpha\beta}\delta_{\lambda\theta}$
and
$dd_{E}\delta_{\alpha\beta}\delta_{\lambda\theta}-d_{E}^{2}\delta_{\alpha\theta}\delta_{\beta\lambda}$
won’t contribute to the average purity.. Furthermore, a direct consequence of
this is that $\Omega_{k:0}$ will be a linear combination of $2^{k+1}$ outer
products between $\mathds{1}$, $\rho$ and $\Psi$, implying that (as $\rho$ is
pure, $\Psi$ is idempotent with trace one and the trace of an outer product is
the product of traces) the average purity
$\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon^{2}_{k:0})]$ will be a sum
of the scalar terms in $\omega^{(\rho)}_{k:0}$ together with an extra overall
factor of $d-1$. Let us denote by

$\displaystyle A=d_{E}(d_{E}^{2}-1),\quad B=d_{E}^{2}-d_{E}^{2}=0,\quad
C=d_{E}^{2}(d_{S}-1/d_{S}),\quad D=d_{E}(dd_{E}-1/d_{S}),$ (72)

the terms that appear in expression (69) after having taken
$\delta_{\alpha\theta}\delta_{\beta\lambda}\to 1/d_{S}$ and
$\delta_{\alpha\beta}\delta_{\lambda\theta}\to 1$. Doing similarly for the
$\omega_{k:k}$ case,

$\displaystyle\mathcal{A}=dd_{E}^{2}(d_{E}^{2}+d_{S}^{2}-2),\quad\mathcal{B}=dd_{E}(d^{2}-1),$
(73)

we find that777 A way to deduce this is to sub-label each factor $A,B,\cdots$
by the $i$ index of the term
${\langle{e}_{i}\alpha_{i}|\omega^{X_{i}}|e_{i}^{(\prime)}\theta_{i}\rangle}$
that they come from and then substituting recursively in Eq. (69) as if one
were already evaluating the whole purity. The vanishing of $B$ simplifies this
process a great deal.

$\displaystyle\mathbb{E}_{{}_{U_{i}}}[\operatorname{tr}(\Upsilon_{k:0}^{2})]$
$\displaystyle=\frac{d_{S}^{k}(d-1)}{d^{k+1}(d^{2}-1)^{k+1}}\left[\mathcal{A}\,A^{k-1}+d_{S}\mathcal{B}\left(CA^{k-2}+C\sum_{i=1}^{k-2}d_{S}^{i}D^{i}A^{k-i-2}+d_{S}^{k-1}D^{k-1}\right)\right],$
(74)

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. (27) 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

$\displaystyle\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{S}^{2k}}\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\,\rho_{{}_{\tau(0);k+1}}\rho_{{}_{\tau(k+1);0}}\tilde{\Delta}_{k,\sigma,\tau}^{(d_{E})}\tilde{\tilde{\Delta}}_{k,\sigma,\tau}^{(d_{S})},$
(75)

where the notation in $\rho$ is as in Eq. (33), and where now the symmetric
group is on $\\{0,\ldots,2k+1\\}$; we also make the definitions in Eq. (83)
for $\tilde{\Delta}^{(d_{E})}$, $\tilde{\tilde{\Delta}}^{(d_{S})}$, which are
monomials in $d_{E}$ and $d_{S}$ correspondingly with degree determined by
$\sigma$, $\tau$. 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
$\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon_{0:0}^{2})]=(d_{E}+d_{S})/(d_{E}d_{S}+1)$
as expected [36, 37, 38, 39, 15, 40, 41].

In the small subsystem limit, as detailed below in F.2.1, we obtain

$\displaystyle\lim_{d_{E}\to\infty}\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{S}^{2k+1}},$
(76)

corresponding to the purity in the ergodic case in this limit. Similarly, in
the long time limit, we obtain in F.2.2,

$\displaystyle\lim_{k\to\infty}\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon^{2})]=\frac{1}{d_{E}},$
(77)

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.(64) now with same
unitary (throughout we take sums over repeated indices and omit the Haar
measure symbol),

$\displaystyle\omega_{k:0}^{(\rho)}\equiv\int\limits_{U}U^{\dagger}\operatorname{\mathfrak{S}}_{x_{1}\alpha_{1}}U^{\dagger}\cdots{\operatorname{\mathfrak{S}}}_{x_{k}\alpha_{k}}U^{\dagger}{K}_{\tilde{\imath}\tilde{\jmath}}U\operatorname{\mathfrak{S}}_{\beta_{k}x_{k}}\cdots{U}\operatorname{\mathfrak{S}}_{\beta_{1}x_{1}}U\rho{U}^{\dagger}\operatorname{\mathfrak{S}}_{y_{1}\lambda_{1}}U^{\dagger}\cdots\operatorname{\mathfrak{S}}_{y_{k}\lambda_{k}}U^{\dagger}{K}^{\dagger}_{\tilde{\imath}\tilde{\jmath}}{U}\operatorname{\mathfrak{S}}_{\theta_{k}y_{k}}\cdots{U}\operatorname{\mathfrak{S}}_{\theta_{1}y_{1}}U.$
(78)

We decompose the unitaries in the whole $ES$ space as
$U=\sum{U}_{ab}|a\rangle\\!\langle{b}|$ and
$U^{\dagger}=\sum{U}^{*}_{a^{\prime}{b}^{\prime}}|b^{\prime}\rangle\\!\langle{a}^{\prime}|$
(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 $\rho$ in Eq. (78) 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
$\operatorname{\mathfrak{S}}$ operators, i.e. the unitary components will
originally appear as
$U^{*}_{i^{\prime}_{0}j^{\prime}_{0}}\cdots{U}^{*}_{i^{\prime}_{k}j^{\prime}_{k}}U_{i_{k}j_{k}}\cdots{U_{i_{0}j_{0}}}U^{*}_{i^{\prime}_{k+1}j^{\prime}_{k+1}}\cdots{U}^{*}_{i^{\prime}_{2k+1}j^{\prime}_{2k+1}}U_{i^{\prime}_{2k+1}j^{\prime}_{2k+1}}\cdots{U_{i^{\prime}_{k+1}j^{\prime}_{k+1}}}$,
and we then rearrange them to the form of Eq. (27) just keeping track of the
correct order in the operator part, for which we introduce sets of $E$-basis
into each $\operatorname{\mathfrak{S}}$ operator as
$\operatorname{\mathfrak{S}}_{ab}=\sum|ea\rangle\\!\langle{eb}|$ labeled by
$e$ for the ones to the left of $\rho$ in Eq. (78), priming the ones between
adjoint unitaries, and by $\epsilon$ the corresponding ones to the right, also
priming the ones between unitaries.

This then leads to

$\displaystyle\omega_{k:0}^{(\rho)}$
$\displaystyle=\int\limits_{U}U_{i_{0}j_{0}}\cdots{U}_{i_{2k+1}j_{2k+1}}U^{*}_{i^{\prime}_{0}j^{\prime}_{0}}\cdots{U}^{*}_{i^{\prime}_{2k+1}j^{\prime}_{2k+1}}\,\langle{i}^{\prime}_{k}|K_{\tilde{\imath}\tilde{\jmath}}|i_{k}\rangle\\!\langle{j}_{0}|\rho|j^{\prime}_{k+1}\rangle\\!\langle{i}^{\prime}_{2k+1}|K^{\dagger}_{\tilde{\imath}\tilde{\jmath}}|i_{2k+1}\rangle|j_{0}^{\prime}\rangle\\!\langle{j}_{k+1}|$
$\displaystyle\hskip
144.54pt\prod_{\ell=1}^{k}\delta_{i_{\ell-1}(ex)_{\ell}}\delta_{j_{\ell}(e\beta)_{\ell}}\delta_{i_{\ell+k}(\epsilon{y})_{\ell}}\delta_{j_{\ell+k+1}(\epsilon\theta)_{\ell}}\delta_{i^{\prime}_{\ell-1}(e^{\prime}{x})_{\ell}}\delta_{j_{\ell}^{\prime}(e^{\prime}\alpha)_{\ell}}\delta_{i^{\prime}_{\ell+k}(\epsilon^{\prime}{y})_{\ell}}\delta_{j_{\ell+k+1}^{\prime}(\epsilon^{\prime}\lambda)_{\ell}}$
$\displaystyle=\langle{i}^{\prime}_{k}|K_{\tilde{\imath}\tilde{\jmath}}|i_{k}\rangle\\!\langle{j}_{0}|\rho|j^{\prime}_{k+1}\rangle\\!\langle{i}^{\prime}_{2k+1}|K^{\dagger}_{\tilde{\imath}\tilde{\jmath}}|i_{2k+1}\rangle|j_{0}^{\prime}\rangle\\!\langle{j}_{k+1}|\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\prod_{n=0}^{2k+1}\delta_{i_{n}i^{\prime}_{\sigma(n)}}\delta_{j_{n}j^{\prime}_{\tau(n)}}$
$\displaystyle\hskip
108.405pt\prod_{\ell=1}^{k}\delta_{i_{\ell-1}(ex)_{\ell}}\delta_{j_{\ell}(e\beta)_{\ell}}\delta_{i_{\ell+k}(\epsilon{y})_{\ell}}\delta_{j_{\ell+k+1}(\epsilon\theta)_{\ell}}\delta_{i^{\prime}_{\ell-1}(e^{\prime}{x})_{\ell}}\delta_{j_{\ell}^{\prime}(e^{\prime}\alpha)_{\ell}}\delta_{i^{\prime}_{\ell+k}(\epsilon^{\prime}{y})_{\ell}}\delta_{j_{\ell+k+1}^{\prime}(\epsilon^{\prime}\lambda)_{\ell}}$
$\displaystyle=\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\langle{i}^{\prime}_{k}|K_{\tilde{\imath}\tilde{\jmath}}|i^{\prime}_{\sigma(k)}\rangle\\!\langle{j}^{\prime}_{\tau(0)}|\rho|j^{\prime}_{k+1}\rangle\\!\langle{i}^{\prime}_{2k+1}|K^{\dagger}_{\tilde{\imath}\tilde{\jmath}}|i^{\prime}_{\sigma(2k+1)}\rangle|j_{0}^{\prime}\rangle\\!\langle{j}^{\prime}_{\tau(k+1)}|$
$\displaystyle\hskip
72.26999pt\prod_{\ell=1}^{k}\delta_{i^{\prime}_{\sigma(\ell-1)}(ex)_{\ell}}\delta_{j^{\prime}_{\tau(\ell)}(e\beta)_{\ell}}\delta_{i^{\prime}_{\sigma(\ell+k)}(\epsilon{y})_{\ell}}\delta_{j^{\prime}_{\tau(\ell+k+1)}(\epsilon\theta)_{\ell}}\delta_{i^{\prime}_{\ell-1}(e^{\prime}{x})_{\ell}}\delta_{j_{\ell}^{\prime}(e^{\prime}\alpha)_{\ell}}\delta_{i^{\prime}_{\ell+k}(\epsilon^{\prime}{y})_{\ell}}\delta_{j_{\ell+k+1}^{\prime}(\epsilon^{\prime}\lambda)_{\ell}}.$
(79)

We now replace $i^{\prime}\to\varepsilon\varsigma$ and
$j^{\prime}\to\varepsilon^{\prime}\varsigma^{\prime}$ in order to split the
$E$ and $S$ parts explicitly, leaving

$\displaystyle\omega_{k:0}^{(\rho)}$
$\displaystyle=\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\langle\varepsilon_{k}\varsigma_{k}|K_{\tilde{\imath}\tilde{\jmath}}|\varepsilon_{\sigma(k)}\varsigma_{\sigma(k)}\rangle\\!\langle\varepsilon^{\prime}_{\tau(0)}\varsigma^{\prime}_{\tau(0)}|\rho|\varepsilon^{\prime}_{k+1}\varsigma^{\prime}_{k+1}\rangle\\!\langle\varepsilon_{2k+1}\varsigma_{2k+1}|K^{\dagger}_{\tilde{\imath}\tilde{\jmath}}|\varepsilon_{\sigma(2k+1)}\varsigma_{\sigma(2k+1)}\rangle|\varepsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle\\!\langle\varepsilon^{\prime}_{\tau(k+1)}\varsigma^{\prime}_{\tau(k+1)}|$
$\displaystyle\hskip
36.135pt\prod_{\ell=1}^{k}\delta_{(\varepsilon\varsigma)_{\sigma(\ell-1)}(ex)_{\ell}}\delta_{(\varepsilon^{\prime}\varsigma^{\prime})_{\tau(\ell)}(e\beta)_{\ell}}\delta_{(\varepsilon\varsigma)_{\sigma(\ell+k)}(\epsilon{y})_{\ell}}\delta_{(\varepsilon^{\prime}\varsigma^{\prime})_{\tau(\ell+k+1)}(\epsilon\theta)_{\ell}}\delta_{(\varepsilon\varsigma)_{\ell-1}(e^{\prime}{x})_{\ell}}\delta_{(\varepsilon^{\prime}\varsigma^{\prime})_{\ell}(e^{\prime}\alpha)_{\ell}}\delta_{(\varepsilon\varsigma)_{\ell+k}(\epsilon^{\prime}{y})_{\ell}}\delta_{(\varepsilon^{\prime}\varsigma^{\prime})_{\ell+k+1}(\epsilon^{\prime}\lambda)_{\ell}}$
$\displaystyle=\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\langle\varepsilon^{\prime}_{\tau(0)}\varsigma^{\prime}_{\tau(0)}|\rho|\varepsilon^{\prime}_{k+1}\varsigma^{\prime}_{k+1}\rangle|\varepsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle\\!\langle\varepsilon^{\prime}_{\tau(k+1)}\varsigma^{\prime}_{\tau(k+1)}|\delta_{\varepsilon_{k}\varepsilon_{\sigma(2k+1)}}\delta_{\varepsilon_{\sigma(k)}\varepsilon_{2k+1}}$
$\displaystyle\hskip
36.135pt\prod_{\ell=1}^{k}\delta_{\varepsilon_{\sigma(\ell-1)}\varepsilon^{\prime}_{\tau(\ell)}}\delta_{\varepsilon_{\sigma(\ell+k)}\varepsilon^{\prime}_{\tau(\ell+k+1)}}\delta_{\varepsilon_{\ell-1}\varepsilon^{\prime}_{\ell}}\delta_{\varepsilon_{\ell+k}\varepsilon^{\prime}_{\ell+k+1}}\delta_{\varsigma^{\prime}_{\ell}\alpha_{\ell}}\delta_{\varsigma^{\prime}_{\tau(\ell)}\beta_{\ell}}\delta_{\varsigma^{\prime}_{\tau(\ell+k+1)}\theta_{\ell}}\delta_{\varsigma^{\prime}_{\ell+k+1}\lambda_{\ell}}\prod_{n=0}^{2k+1}\delta_{\varsigma_{\sigma(n)}\varsigma_{n}}.$
(80)

The ancillary part is simply
$\Omega_{k:0}^{(AB)}=d_{S}^{-k}|\alpha_{1}\beta_{1}\ldots\alpha_{k}\beta_{k}\rangle\\!\langle\theta_{1}\lambda_{1}\ldots\theta_{k}\lambda_{k}|$,
so we get for the analogue of Eq. (63),

$\displaystyle\Omega_{k:0}=\omega_{k:0}^{(\rho)}\otimes\Omega_{k:0}^{(AB)}$
$\displaystyle=\frac{1}{d_{S}^{k}}\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\delta_{\varepsilon_{k}\varepsilon_{\sigma(2k+1)}}\delta_{\varepsilon_{\sigma(k)}\varepsilon_{2k+1}}|\varepsilon^{\prime}_{0}\varsigma^{\prime}_{0}\varsigma^{\prime}_{1}\varsigma^{\prime}_{\tau(1)}\cdots\varsigma^{\prime}_{k}\varsigma^{\prime}_{\tau(k)}\rangle\langle\varepsilon^{\prime}_{\tau(k+1)}\varsigma^{\prime}_{\tau(k+1)}\varsigma^{\prime}_{\tau(k+2)}\varsigma^{\prime}_{k+2}\cdots\varsigma^{\prime}_{\tau(2k+1)}\varsigma^{\prime}_{2k+1}|$
$\displaystyle\hskip
36.135pt\langle\varepsilon^{\prime}_{\tau(0)}\varsigma^{\prime}_{\tau(0)}|\rho|\varepsilon^{\prime}_{k+1}\varsigma^{\prime}_{k+1}\rangle\prod_{n=0}^{2k+1}\delta_{\varsigma_{\sigma(n)}\varsigma_{n}}\prod_{\ell=1}^{k}\delta_{\varepsilon_{\sigma(\ell-1)}\varepsilon^{\prime}_{\tau(\ell)}}\delta_{\varepsilon_{\sigma(\ell+k)}\varepsilon^{\prime}_{\tau(\ell+k+1)}}\delta_{\varepsilon_{\ell-1}\varepsilon^{\prime}_{\ell}}\delta_{\varepsilon_{\ell+k}\varepsilon^{\prime}_{\ell+k+1}},$
(81)

and thus

$\displaystyle\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon_{k:0}^{2})]$
$\displaystyle=\operatorname{tr}[(\rho\otimes\Psi^{\otimes{k}})\,\Omega_{k:0}]$
$\displaystyle=\frac{1}{d_{S}^{2k}}\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\mathrm{Wg}(\tau\sigma^{-1})\langle\varepsilon^{\prime}_{\tau(0)}\varsigma^{\prime}_{\tau(0)}|\rho|\varepsilon^{\prime}_{k+1}\varsigma^{\prime}_{k+1}\rangle\\!\langle\varepsilon^{\prime}_{\tau(k+1)}\varsigma^{\prime}_{\tau(k+1)}|\rho|\varepsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle\delta_{\varepsilon_{k}\varepsilon_{\sigma(2k+1)}}\delta_{\varepsilon_{2k+1}\varepsilon_{\sigma(k)}}$
$\displaystyle\hskip 36.135pt\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{\varepsilon_{\ell-1}\varepsilon^{\prime}_{\ell}}\delta_{\varepsilon_{\sigma(\ell-1)}\varepsilon^{\prime}_{\tau(\ell)}}\delta_{\varsigma^{\prime}_{\ell}\varsigma^{\prime}_{\tau(\ell)}}\prod_{n=0}^{2k+1}\delta_{\varsigma_{\sigma(n)}\varsigma_{n}},$
(82)

as stated in Eq. (6), where we make the definitions

$\tilde{\Delta}_{k,\sigma,\tau}^{(d_{E})}=\delta_{e_{k}e_{\sigma(2k+1)}}\delta_{e_{2k+1}e_{\sigma(k)}}\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{e_{\ell-1}e^{\prime}_{\ell}}\delta_{e_{\sigma(\ell-1)}e^{\prime}_{\tau(\ell)}},\quad\tilde{\tilde{\Delta}}_{k,\sigma,\tau}^{(d_{S})}=\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{s^{\prime}_{\ell}{s}^{\prime}_{\tau(\ell)}}\prod_{n=0}^{2k+1}\delta_{s_{n}s_{\sigma(n)}}$
(83)

####  F.2.1 Small subsystem limit

In the $d_{E}\to\infty$ limit the only term that remains is the one with

$\displaystyle\sigma=\tau=(0,k+1)(1,k+2)\cdots(k,2k+1),$ (84)

as expressed in cycle notation, which simply means $\sigma(0)=k+1=\tau(0)$,
$\sigma(k+1)=0=\tau(k+1)$, and so on. This then leads to the correct limit, as

$\displaystyle\tilde{\tilde{\Delta}}_{k,\sigma,\tau}^{(d_{S})}\to\sum_{\begin{subarray}{c}\varsigma^{(\prime)}_{i}=1\\\
\varsigma^{\prime}\neq\\{\varsigma^{\prime}_{0},\varsigma^{\prime}_{k+1}\\}\end{subarray}}^{d_{S}}\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{\varsigma^{\prime}_{\ell}\varsigma^{\prime}_{\tau(\ell)}}\prod_{n=0}^{2k+1}\delta_{\varsigma_{\sigma(n)}\varsigma_{n}}$
$\displaystyle={d_{S}}^{2k+1},$ (85)
$\displaystyle\tilde{\Delta}_{k,\sigma,\tau}^{(d_{E})}\to\sum_{\begin{subarray}{c}\varepsilon^{(\prime)}_{i}=1\\\
\varepsilon^{\prime}\neq\\{\varepsilon^{\prime}_{0},\varepsilon^{\prime}_{k+1}\\}\end{subarray}}^{d_{E}}\delta_{\varepsilon_{k}\varepsilon_{\sigma(2k+1)}}\delta_{\varepsilon_{2k+1}\varepsilon_{\sigma(k)}}\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{\varepsilon_{\ell-1}\varepsilon^{\prime}_{\ell}}\delta_{\varepsilon_{\sigma(\ell-1)}\varepsilon^{\prime}_{\tau(\ell)}}$
$\displaystyle={d_{E}}^{2k+2},$ (86)

while keeping $\sigma\tau^{-1}=\mathbf{1}^{2k+2}$, i.e. with the argument in
the $\mathrm{Wg}$ function an identity, which generates the least denominator
powers because $\\#\bm{1}^{n}=\\#[(1)(2)\cdots(n)]=n$, i.e. the identity
generates the most number of cycles. When considering other permutations for
$\sigma$ and $\tau$ we see that the $d_{E}$ powers in the numerator can only
decrease while those in the denominator can only increase (as the identity is
the permutation that generates the most cycles) so indeed in the large
subsystem limit this is the only pair of permutations that survive; the terms
in $\rho$ yield traces when summed over basis vectors and thus give factors of
one,

$\displaystyle\sum_{\varepsilon^{\prime}_{k+1},\varsigma^{\prime}_{k+1}}\langle\varepsilon^{\prime}_{k+1}\varsigma^{\prime}_{k+1}|\rho|\varepsilon^{\prime}_{k+1}\varsigma^{\prime}_{k+1}\rangle=\sum_{\varepsilon^{\prime}_{0},\varsigma^{\prime}_{0}}\langle\varepsilon^{\prime}_{0}\varsigma^{\prime}_{0}|\rho|\varepsilon^{\prime}_{0}\varsigma^{\prime}_{0}\rangle=\operatorname{tr}\rho=1.$
(87)

This leads then to the conclusion

$\displaystyle\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon_{k:0}^{2})]$
$\displaystyle\sim\frac{d_{S}^{2k+1}d_{E}^{2k+2}}{d_{S}^{2k}}\mathrm{Wg}(\mathbf{1}^{2k+2})=d_{S}d_{E}^{2k+2}\frac{1}{(d_{E}d_{S})^{2k+2}}=\frac{1}{d_{S}^{2k+1}},\hskip
7.22743pt\text{when}\,\,d_{E}\to\infty.$ (88)

####  F.2.2 Long time limit

For the case when $k\to\infty$, the $\mathrm{Wg}$ function behaves dominantly
as in the small subsystem limit and the only pair of permutations that
identify all dimension powers in the numerator without dependence of $k$ are
identities $\sigma=\tau=\mathbf{1}^{2k+2}$. In this case

$\displaystyle\sum_{\begin{subarray}{c}\varsigma^{(\prime)}_{i}=1\\\
\varsigma^{\prime}\neq\\{\varsigma^{\prime}_{0},\varsigma^{\prime}_{k+1}\\}\end{subarray}}^{d_{S}}\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{\varsigma^{\prime}_{\ell}\varsigma^{\prime}_{\ell}}\prod_{n=0}^{2k+1}\delta_{\varsigma_{n}\varsigma_{n}}$
$\displaystyle={d_{S}}^{4k+2},$ (89)
$\displaystyle\sum_{\begin{subarray}{c}\varepsilon^{(\prime)}_{i}=1\\\
\varepsilon^{\prime}\neq\\{\varepsilon^{\prime}_{0},\varepsilon^{\prime}_{k+1}\\}\end{subarray}}^{d_{E}}\delta_{\varepsilon_{k}\varepsilon_{2k+1}}\delta_{\varepsilon_{2k+1}\varepsilon_{k}}\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq{k+1}\end{subarray}}^{2k+1}\delta_{\varepsilon_{\ell-1}\varepsilon^{\prime}_{\ell}}\delta_{\varepsilon_{\ell-1}\varepsilon^{\prime}_{\ell}}$
$\displaystyle={d_{E}}^{2k+1},$ (90)

giving

$\displaystyle\mathbb{E}_{{}_{U}}[\operatorname{tr}(\Upsilon_{k:0}^{2})]$
$\displaystyle\sim\frac{d_{S}^{4k+2}d_{E}^{2k+1}}{d_{S}^{2k}}\operatorname{tr}(\rho^{2})\mathrm{Wg}(\mathbf{1}^{2k+2})=\frac{d_{S}^{2k+2}d_{E}^{2k+1}}{(d_{E}d_{S})^{2k+2}}=\frac{1}{d_{E}},\hskip
7.22743pt\text{when}\,\,k\to\infty,$ (91)

as $\operatorname{tr}(\rho^{2})=1$ by assumption.