text:
\section{Introduction}

The quest towards understanding how thermodynamics emerges from quantum theory has seen a great deal of recent progress~\cite{Gogolin, Goold}. Most prominently, the conundrum of how to recover irreversible phenomena that obey the second law of thermodynamics, starting from reversible and recurrent Schr\"odinger dynamics, has been rectified by considering an analogous dynamical equilibrium. This is known generically as \emph{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~\footnote{ 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 \emph{thermalisation} and the equilibrium state with a \emph{thermal state}. For detail see e.g.~\cite{Gogolin}}.

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.~\cite{GarciaPintos}). In a typical setup, the closed dynamics of a system and environment, as illustrated in Fig.~\ref{Fig interventions}(a), is considered. Remarkably, it can be shown that the expectation value for \textit{any} observable $\mathcal{O}$ on the system alone lies in the neighbourhood of its time-averaged value for almost all times [Fig.~\ref{Fig interventions}(b)]. Yet, such equilibration results, in full generality, remain unsolved~\cite{Gogolin}; in particular, it is unknown whether they extend to correlations between observables at different times.

\begin{figure}[t]
      \begin{tikzpicture}[yscale=0.2, xscale=0.15, every node/.style={scale=0.925}]
      \shade[outer color=white, inner color=red!20!white, draw=black] (0,0) [rounded corners] rectangle (8,10);
      \shade[inner color=red!90!white, outer color=red!25!white, draw=black] (2,3) [rounded corners] rectangle (6,7);
      \node[scale=1.5] at (11,4) {$\longrightarrow$};
      \node[above] at (11,4.5) {$e^{-iH\tau}$};
      \shade[xshift=5.5in,outer color=white, inner color=red!10!white, draw=black] (0,0) [rounded corners] rectangle (8,10);
      \shade[xshift=5.5in,inner color=white, outer color=red!5!white, draw=black] (2,3) [rounded corners] rectangle (6,7);
      \node[below] at (4.5,0) {$t=0$};
      \node[below] at (18.5,0) {$t=\tau$};
      \node at (1,8) {$E$};
      \node at (3,6) {$S$};
      \node at (15,8) {$E$};
      \node at (17,6) {$S$};
      \node at (11,-4) {(a)};
      \end{tikzpicture}
~
    \begin{tikzpicture}[yscale=0.25, xscale=0.18, every node/.style={scale=0.925}]
    \draw[<->] (0,-4) -- (0,6);
    \draw[-] (0,0) -- (8,0);
    \draw[-, dashed] (8,0) -- (12,0);
    \draw[->] (12,0) -- (21,0);
    \draw [thick, blue] plot [smooth, tension=1] coordinates { (0,2) (1,5) (2,0) (3,0.3) (4,0) (5,0.1) (6,0) (7,-0.2) (8,0) (9,0)};
    \draw [thick, blue] plot [smooth, tension=1] coordinates { (12,0) (13,-0.1) (14,0.1) (15,0) (16,0.2) (17,0) (18,-3) (19,0) (20,0) (21,0)};
    \node[label={[label distance=0.5cm,text depth=-1ex,rotate=90]right:$\langle\mathcal{O}\rangle_{\rho(t)}-\langle\mathcal{O}\rangle_{\omega}$}] at (-2.5,-5) {};
    \node[right] at (21,0) {$t$};
    \node[below] at (1,0) {$0$};
    \draw[-, dashed] (1,0) -- (1,5);
    \node at (11,-6) {(b)};
    \end{tikzpicture}

\vspace{0.15 in}

\begin{tikzpicture}[xscale=0.28,yscale=0.3, every node/.style={xscale=0.8,yscale=0.9}]
    \draw[-, dashed, thick, red!50!blue, rounded corners, fill=red!50!blue!4!white] (-2.5,0) -- (-2.5,2.5) -- (25,2.5) -- (25,-0.25) -- (23,-0.25) -- (23,-2.5) -- (20,-2.5) -- (20,0.5) -- (16.5,0.5) -- (16.5,-2.5) -- (10.5,-2.5) -- (10.5,0.5) -- (7.5,0.5) -- (7.5,-2.5) -- (4.5,-2.5) -- (4.5,0.5) -- (1.5,0.5) -- (1.5,-2.5) -- (-2.5,-2.5) -- (-2.5,0);
    \draw[-, line width=0.5mm, white] (16.01,-2.5) -- (14,-2.5);
    \node at (15.1,-2.55) {$\cdots$};
    \shade[inner color=white, outer color=black!10!white, draw=black] (-1,0) ellipse (1cm and 2cm);
    \node at (-1,0) {$\rho$};
    \draw[-] (-0.15,1) -- (5,1);
    \draw[-] (-0.15,-1) -- (2,-1);
    \node[above] at (0.6,0.8) {\scriptsize${E}$};
    \node[above] at (0.6,-1.2) {\scriptsize${S}$};
    \shade[top color=blue!50!white, bottom color=white, draw=black] (2,-2.5) rectangle (4,0);
    \node at (3,-1.25) {$\mathcal{M}_{0}$};
    \draw[-] (4,-1) -- (5,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (5,-1.5) rectangle (7,1.5);
    \node at (6,0) {$U_{1}$};
    \draw[-] (7,1) -- (11,1);
    \draw[-] (7,-1) -- (8,-1);
    \shade[top color=blue!50!white, bottom color=white, draw=black] (8,-2.5) rectangle (10,0);
    \node at (9,-1.25) {$\mathcal{M}_{1}$};
    \draw[-] (10,-1) -- (11,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (11,-1.5) rectangle (13,1.5);
    \node at (12,0) {$U_2$};
    \node at (15,0) {$\cdots$};
    \draw[-] (13,1) -- (14,1);
    \draw[-] (13,-1) -- (14,-1);
    \draw[-] (16,1) -- (20.5,1);
    \draw[-] (16,-1) -- (17,-1);
    \shade[top color=blue!50!white, bottom color=white, draw=black] (17,-2.5) rectangle (19.5,0);
    \node[rotate=45] at (18.25,-1.25) {$\mathcal{M}_{_{k-1}}$};
    \draw[-] (19.5,-1) -- (20.5,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (20.5,-1.5) rectangle (22.5,1.5);
    \node at (21.5,0) {$U_k$};
    \draw[-] (22.5,1) -- (24,1);
    \node at (24,1) {\Large\textbackslash};
    \draw[-] (22.5,-1) -- (24,-1);
    \node at (25,-1.2) {$\rho^{(S)}_k$};
    \node at (12,-5) {(c)};
    \end{tikzpicture}
\caption{\textbf{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\"odinger'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.}\label{Fig interventions}
\end{figure}

To see this, consider the case depicted in Fig.~\ref{Fig interventions}(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 \emph{memory} built up within the dynamics itself, independently of the external interventions, which are the hallmark of so-called \emph{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~\cite{Srednicki1999}. 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~\cite{processtensor2, PhysRevA.92.022102}, 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~\cite{Schlosshauer2007, VegaAlonso}; this is mathematically equivalent (through an analogue of Stinespring's dilation theorem) to continually refreshing (discarding and replacing) the environment's state~\cite{CollisionModels}, i.e., artificially throwing away information from the environment (see Fig.~\ref{Markov approx circuit}). In reality, this is not usually how a closed system evolves, leading to the question: \textit{In the general case, how does dissipative Markovian dynamics emerge from closed Schr\"odinger 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~\cite{spohn}.

\begin{figure}[t]
\begin{tikzpicture}[xscale=0.3,yscale=0.34, every node/.style={xscale=0.8,yscale=0.9}]
    \draw[-, line width=0.5mm, white] (16,-3) -- (14,-3);
    \shade[inner color=white, outer color=black!10!white, draw=black] (-1,1.1) ellipse (0.9cm and 1cm);
    \shade[inner color=white, outer color=black!10!white, draw=black] (-1,-1.1) ellipse (0.9cm and 1cm);
    \node at (-1,1.1) {$\rho_{_E}^{(1)}$};
    \node at (-1,-1.1) {$\rho_{_S}$};
    \draw[-] (-0.15,1) -- (5,1);
    \draw[-] (-0.15,-1) -- (2,-1);
    \node[above] at (0.6,0.8) {\scriptsize${E}$};
    \node[above] at (0.6,-1.2) {\scriptsize${S}$};
    \shade[top color=blue!50!white, bottom color=white, draw=black] (2,-2.5) rectangle (4,0);
    \node at (3,-1.25) {$\mathcal{M}_{0}$};
    \draw[-] (4,-1) -- (5,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (5,-1.5) rectangle (7,1.5);
    \node at (6,0) {$U_{1}$};
    \draw[-, rounded corners] (7,1) -- (8.25,1);
    \node at (8.25,1) {\Large\textbackslash};
    \draw[-, rounded corners] (9.75,1) -- (11,1);
    \shade[inner color=white, outer color=black!10!white, draw=black] (9.75,1.25) ellipse (0.75cm and 0.8cm);
    \node at (9.75,1.25) {$\scriptstyle{\rho_{_E}^{(2)}}$};
    \draw[-] (7,-1) -- (8,-1);
    \shade[top color=blue!50!white, bottom color=white, draw=black] (8,-2.5) rectangle (10,0);
    \node at (9,-1.25) {$\mathcal{M}_{1}$};
    \draw[-] (10,-1) -- (11,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (11,-1.5) rectangle (13,1.5);
    \node at (12,0) {$U_2$};
    \node at (15,0) {$\cdots$};
    \draw[-, rounded corners] (13,1) -- (14.25,1);
    \node at (14.25,1) {\Large\textbackslash};
    \draw[-] (13,-1) -- (14,-1);
    \draw[-] (19.25,1) -- (20.5,1);
    \shade[inner color=white, outer color=black!10!white, draw=black] (19.25,1.25) ellipse (0.75cm and 0.8cm);
    \node at (19.275,1.25) {$\scriptstyle{\rho_{_E}^{(k)}}$};
    \draw[-] (16,-1) -- (17,-1);
    \shade[top color=blue!50!white, bottom color=white, draw=black] (17,-2.5) rectangle (19.5,0);
    \node[rotate=45] at (18.25,-1.25) {$\mathcal{M}_{_{k-1}}$};
    \draw[-] (19.5,-1) -- (20.5,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (20.5,-1.5) rectangle (22.5,1.5);
    \node at (21.5,0) {$U_k$};
    \draw[-] (22.5,1) -- (24,1);
    \node at (24,1) {\Large\textbackslash};
    \draw[-] (22.5,-1) -- (24,-1);
    \node at (25,-1.2) {$\rho^{(S)}_k$};
    \end{tikzpicture}
    \caption{\textbf{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. \ref{Fig interventions}c)}\label{Markov approx circuit}
\end{figure}

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~\cite{Ledoux, milman, boucheron2013concentration, PopescuWinter, Li} to determine when a subsystem dynamics might typically be close to Markovian. Our approach is broadly similar to that of Refs.~\cite{PopescuWinter, Masanes} 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 \emph{process tensor} framework~\cite{PhysRevLett.101.060401, PhysRevA.80.022339, Costa}, 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~\cite{Modi, processtensor}.

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~\cite{Breuer2015}. It allows for a precise quantification of memory effects (including across multiple time points)~\cite{processtensor2}, 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.

\section{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.~\ref{Fig interventions}(c). Formally, a $k$-step process tensor is a completely positive causal~\footnote{ 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\l{}kowski isomorphism~\cite{Choi1975, Jamiol1972,processtensor}. 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_iB_i}\equiv|\psi^+\rangle\!\langle\psi^+|\equiv\sum_{j,\ell=1}^{d_S}|jj\rangle\!\langle{\ell\ell}|/d_S$ at each time $i$~\cite{Simon2017}. 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$~\footnote{ 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
\begin{gather}
\begin{split}
&\Upsilon_{k:0}=\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\mbb1)\mathcal{S}_k\cdots(U_1\otimes\mbb1)\mathcal{S}_1(U_0\otimes\mbb1),
\end{split}
\label{choi evolution}
\end{gather}
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.~\ref{PT circuit diagram}.

\begin{figure}[ht]
\centering
    \begin{tikzpicture}[thick, scale=0.3, every node/.style={scale=0.95}]
    \shade[inner color=white, outer color=black!10!white, draw=black] (-1,0) ellipse (1cm and 2cm);
    \node at (-1,0) {$\rho$};
    \draw[-] (-0.15,1) -- (2,1);
    \draw[-] (-0.15,-1) -- (2,-1);
    \node[above] at (1,0.8) {\scriptsize${E}$};
    \node[above] at (1,-1.2) {\scriptsize${S}$};
    \shade[bottom color=red!50!white, top color=white, draw=black] (2,-1.5) rectangle (4,1.5);
    \node at (3,0) {$U_{0}$};
    \draw[-] (4,1) -- (6,1);
    \draw[-] (4,-1) -- (6,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (6,-1.5) rectangle (8,1.5);
    \node at (7,0) {$U_{1}$};
    \draw[-] (8,1) -- (10,1);
    \draw[-] (8,-1) -- (10,-1);
    \node at (11,0) {$\cdots$};
    \draw[-] (12,1) -- (14,1);
    \draw[-] (12,-1) -- (14,-1);
    \shade[bottom color=red!50!white, top color=white, draw=black] (14,-1.5) rectangle (16,1.5);
    \node at (15,0) {$U_{k}$};
    \draw[-] (16,1) -- (17,1);
    \node at (17,1) {\Large\textbackslash};
    \draw[-] (16,-1) -- (18,-1);
    \draw[-] (3,-3) -- (10,-3);
    \draw[-] (3,-5) -- (10,-5);
    \node at (11,-4) {$\cdots$};
    \draw[-] (12,-3) -- (18,-3);
    \draw[-] (12,-5) -- (18,-5);
    \draw (3,-3) arc (90:270:1);
    \node[left] at (2.25,-4) {$\Psi_{A_1B_1}$};
    \draw[-] (3,-7) -- (10,-7);
    \draw[-] (3,-9) -- (10,-9);
    \node at (11,-8) {$\cdots$};
    \draw[-] (12,-7) -- (18,-7);
    \draw[-] (12,-9) -- (18,-9);
    \draw (3,-7) arc (90:270:1);
    \node[left] at (2.25,-8) {$\Psi_{A_2B_2}$};
    \node[left] at (2.25,-10.5) {$\vdots$};
    \draw[-] (3,-12) -- (10,-12);
    \draw[-] (3,-14) -- (10,-14);
    \node at (11,-13) {$\cdots$};
    \node[left] at (11.5,-10.5) {$\vdots$};
    \draw[-] (12,-12) -- (18,-12);
    \draw[-] (12,-14) -- (18,-14);
    \draw (3,-12) arc (90:270:1);
    \node[left] at (2.25,-13) {$\Psi_{A_kB_k}$};
    \shade[outer color=red!50!blue!4!white, inner color=white, draw=black] (18,0) [rounded corners] rectangle (20.5,-15);
    \node[rotate=55] at (19.25,-7) {$\Upsilon_{k:0}$};
    \draw[-, blue!60!white] (5,-1) -- (5,-3);
    \node[blue!60!white] at (5,-1) {$\times$};
    \node[right] at (4.8,-2.1) {\scriptsize$\mathcal{S}_1$};
    \node[blue!60!white] at (5,-3) {$\times$};
    \draw[blue!60!white,-] (9,-1) -- (9,-7);
    \node[blue!60!white] at (9,-1) {$\times$};
    \node[right] at (8.8,-6.1) {\scriptsize$\mathcal{S}_2$};
    \node[blue!60!white] at (9,-7) {$\times$};
    \draw[blue!60!white,-] (13,-1) -- (13,-12);
    \node[blue!60!white] at (13,-1) {$\times$};
    \node[right] at (12.8,-11.1) {\scriptsize$\mathcal{S}_k$};
    \node[blue!60!white] at (13,-12) {$\times$};
    \node[above] at (3.25,-3.2) {\scriptsize$A_1$};
    \node[above] at (3.25,-5.2) {\scriptsize$B_1$};
    \node[above] at (3.25,-7.2) {\scriptsize$A_2$};
    \node[above] at (3.25,-9.2) {\scriptsize$B_2$};
    \node[above] at (3.25,-12.2) {\scriptsize$A_k$};
    \node[above] at (3.25,-14.2) {\scriptsize$B_k$};
    \end{tikzpicture}
\caption{\textbf{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.}\label{PT circuit diagram}
\end{figure}

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

\section{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 \emph{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}$~\cite{Simon2017}.

In~\cite{processtensor2}, 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\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
\begin{gather}
\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
\label{def nonMarkov trDistance}
\end{gather}
where $\|X\|_1\equiv\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})}})\geq2\,\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~\cite{PopescuWinter} and equilibration~\cite{Gogolin, ShortFarrelly, Masanes} (or lack thereof, e.g.~\cite{PhysRevLett.106.040401}). It has an experimental importance for comparing quantum processes~\cite{PhysRevA.71.062310} 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~\cite{Breuer2015}. The measure we use here captures \textit{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.

\section{Typicality of Markovian processes}\label{sec: Main Typicality}
In order to quantify non-Markovianity, as defined in Eq.~\eqref{def nonMarkov trDistance}, 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~\cite{ZhangMInt, Li, GuMoments, Collins, Puchala} and leads to the relatively straightforward application of concentration of measure results~\cite{Ledoux, milman, boucheron2013concentration}. Specifically, we use the Haar measure to sample two distinct types of $ES$ evolution:
\begin{enumerate}[leftmargin=*]
\item Random interaction: All $U_i$  independently chosen,
\item Constant interaction: $U_i=U_j,\quad \forall \leq{i}\neq{j}\leq{k}$,
\end{enumerate}
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.~\ref{Fig interventions}(c); the entire set enters into the process tensor through Eq.~\eqref{choi evolution}. 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.

\subsection{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_Ed_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$,
\begin{gather}
    \mathbb{P}[\mathcal{N}\geq\mathcal{B}_k(d_E,d_S)+\epsilon]\leq\mathrm{e}^{-
\mathcal{C}(d_E,d_S)\epsilon^2},\label{eq: main result}
\end{gather}
where $\mathcal{C}(d_E,d_S)=c \, d_Ed_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.

\begin{figure}[t]
\centering
\begin{tikzpicture}[xscale=1.8, yscale=0.8, every node/.style={scale=1}]
    \shade[thick, scale=0.65, inner color=gray!5!white, outer color=white, draw=black]  (-2.225,0.4) ellipse (2.75cm and 2.5cm);
    \draw[-, qmark, very thick, dotted, decorate,decoration={snake, amplitude=.4mm, segment length=0.25in}, rotate around={150:(-0.75,0)}] (-0.7725,0) -- (0.7725,0);
    \node[qmark, right, above] at (-1,0.15) {$\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}$};
    \node at (-1.425,-1.025) {$\Upsilon$};
    \draw[qorange,|<->|] (-1.32,-0.875) -- (-1.26,-0.65);
    \node[right, below, qorange] at (-1.2,-0.625) {$\epsilon$};
    \node[green!50!black] at (-2.33,0) {$\mathcal{B}_k$};
    \node[red!80!black] at (-1.9,0.2) {$\mathbb{E}[\mathcal{N}]$};
    \draw[-, red!80!black, very thick, dotted, rotate around={150:(-0.75,0)}] plot[smooth cycle, tension=.7] coordinates {(-1,0) (-0.8,0.3) (-0.6,0.4) (-0.3,0.6) (0,0.5) (0.4,0.8) (0.7,0.3) (1,0) (0.7,-0.3) (0.4,-0.6) (0,-0.5) (-0.3,-0.6) (-0.6,-0.4) (-0.8,-0.3)};
    \draw[-, green!50!black, very thick, dotted, rotate around={150:(-0.75,0)},scale=1.4] plot[smooth cycle, tension=.7] coordinates {(-1,0) (-0.8,0.3) (-0.6,0.4) (-0.3,0.6) (0,0.5) (0.4,0.8) (0.7,0.3) (1,0) (0.7,-0.3) (0.4,-0.6) (0,-0.5) (-0.3,-0.6) (-0.6,-0.4) (-0.8,-0.3)};
    \draw[blue!80!black, semithick, |<->|] (-1.415,-0.83) -- (-1.15,0.16);
    \node[right] at (-1.5, -0.04) {$\color{blue!80!black}{\mathcal{N}}$};
    \node[scale=1] at (-1.5,-2) {(a)};
    \end{tikzpicture}
    ~
    \begin{tikzpicture}[xscale=1.8, yscale=0.8, every node/.style={scale=1}]
    \shade[thick, scale=0.65, inner color=white, outer color=white, draw=black]  (-2.225,0.4) ellipse (2.75cm and 2.5cm);
    \foreach\i in {0,0.02,...,0.7} {
        \fill[opacity=\i*0.1, blue, rotate around={60:(-1,-1)}] (0,0) ellipse ({0.7-\i} and {1.8-1.8*\i});
      }
    \draw[qmark, -, very thick, dotted, decorate,decoration={snake, amplitude=.4mm, segment length=0.25in}, rotate around={150:(-0.75,0)}] (-0.7725,0) -- (0.7725,0);
    \node[qmark, right, above] at (-0.9,0.35) {$\Upsilon^{\scriptscriptstyle{(\mathrm{M})}}$};
    \node at (-1.5,-1) {$\Upsilon$};
    \node at (-1.5,-2) {(b)};
    \end{tikzpicture}
    \caption{\textbf{Concentration of Markovian processes in large environments.} (a) A geometric \emph{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})}}$.}\label{main cartoon}
\end{figure}

The proof is presented in full in Appendix~\ref{Appendix: Main Result}, 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~\cite{Ledoux,boucheron2013concentration,PopescuWinter} 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.~\eqref{eq: main result}. 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 \emph{small subsystem} or \emph{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^2d_E\gg1\gg\epsilon$; this is fulfilled for $\epsilon=d_E^{-1/3}$~\footnote{Similar to~\cite{PopescuWinter}, here we look for an $x>0$ such that $\epsilon=d_E^{-x}$ and $\epsilon^2d_E=d_E^x$}. A geometrical cartoon to illustrate the result is presented in Fig.~\ref{main cartoon}.

This is the case because the upper bound $\mathcal{B}_k \geq \mathbb{E}(\mathcal{N})$ is given by
\begin{gather}
\mathcal{B}_k(d_E,d_S)\equiv\begin{cases}
\frac{\sqrt{d_E\,\mathbb{E}[\tr(\Upsilon^2)]-x}+y}{2}
&\text{if}\quad d_E<d_S^{2k+1}
\\[0.1in]
\frac{\sqrt{d_S^{2k+1}\mathbb{E}[\tr(\Upsilon^2)]-1}}{2}&\text{if}
\quad d_E\geq{d}_S^{2k+1},\end{cases}\label{avg upper bound}\end{gather}
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}[\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~\cite{chaosYoshida, Brandao2016}, as we will discuss below) in Appendix~\ref{sec:Purity Averages}.

The purity $\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~\ref{appendix: avg processes}, 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 \emph{noisiness} of the quantum process itself and entanglement between system and environment. The expressions take the form
\begin{gather}
    \mathbb{E}_{_{U_i}}[\tr(\Upsilon^{2})]=\frac{d_E^2-1}{d_E(d_Ed_S+1)}\left(\frac{d_E^2-1}{d_E^2d_S^2-1}\right)^k+\frac{1}{d_E},\label{average purity ergodic}
\end{gather}
for the random interaction picture and
\begin{gather}
    \mathbb{E}_{_U}[\tr(\Upsilon^2)]=d_S^{-2k}\hspace*{-0.1in}\sum_{\sigma,\tau\in\mathcal{S}_{2k+2}}\hspace*{-0.1in}\mathrm{Wg}(\tau\sigma^{-1})\,\varrho\,_\tau\,\Delta_{k,\sigma,\tau}^{(d_E,d_S)},\label{average time independent purity}
\end{gather}
in the constant interaction case, where $\mathcal{S}_n$ is the symmetric group over $n$ elements, $\mathrm{Wg}$ is known as the Weingarten function~\cite{GuMoments}, $\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$.

\subsection{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.~\eqref{average time independent purity} 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~\cite{GuMoments,Collins}, 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,
\begin{gather}
    \lim_{d_E\to\infty}\mathds{E}[\tr(\Upsilon^{2})]=\frac{1}{d_S^{2k+1}},
\end{gather}
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:
\begin{gather}\lim_{d_E\to\infty}\mathds{E}[\mathcal{N}]=\lim_{d_E\to\infty}\mathcal{B}_k(d_E,d_S)=0.\end{gather}
Furthermore, from Eq.~\eqref{average purity ergodic} 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
\begin{gather}
    \lim_{k\to\infty}\mathds{E}[\tr(\Upsilon^{2})]=\frac{1}{d_E},
\end{gather}
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
\begin{gather}\lim_{k\to\infty}\mathds{E}[\mathcal{N}]\leq\lim_{k\to\infty}\mathcal{B}_k(d_E,d_S)=1,\end{gather}
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.~\eqref{eq: main result}, and the bound on non-Markovianity, Eq.~\eqref{avg upper bound}, 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_Ed_S/4$ (see e.g.~\cite{AlmostAllPureMME}). As we also detail in Appendix~\ref{sec:Purity Averages}, when we take the Haar measure, in either case, we recover the average purity
\begin{equation}
    \mathds{E}[\tr(\rho_S^{2})]=\frac{d_E+d_S}{d_Ed_S+1},
\end{equation}
where $\rho_S\equiv\Upsilon_{0:0}=\tr_E(U\rho\,U^\dagger)$, as found in~\cite{Lubkin, PageEntropy, LloydPagels, ScottCaves, Li, Giraud, Pasquale}. This then leads to
\begin{gather}
    \mathds{E}\left[\mathcal{D}\left(\rho_S,\frac{\mbb1}{d_S}\right)\right]\leq\frac{1}{2}\sqrt{\frac{d_S^2-1}{d_Ed_S+1}},
\end{gather}
which goes as $\sqrt{d_S/4d_E}$ when $d_E\gg{d}_S$, and is also a standard result \cite{AlmostAllPureMME}.

\subsection{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.~\ref{numPlot}. 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.

\begin{figure}[t]
     \includegraphics[width=0.475\textwidth]{numeric_avgN_DIF}
     \caption{\textbf{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 $\lfloor40/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~\cite{Mezzadri}.}\label{numPlot}
\end{figure}

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~\cite{chaosYoshida}. The mechanism for equilibration is precisely that of dephasing~\cite{Gogolin}, 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~\cite{Cotler2017}, that reproduce the Haar distribution up to the $n$-th moment. As mentioned above, our expression for the upper bound in Eq.~\eqref{avg upper bound} 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.~\eqref{eq: main result} would also hold for $n$-designs with large enough $n$~\cite{Low}. This type of distribution has recently been shown to be approximated by a wide class of physical evolutions~\cite{PhysRevX.7.021006}. In particular, the unitary circuits involved in complex quantum computations typically look Haar random~\cite{Brandao2016}; 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\geq0}$ for time-independent Hamiltonians $H$ cannot strictly become Haar random~\cite{chaosYoshida}, certain complex Hamiltonian systems can be treated as random up to a certain Haar-moment via unitary designs~\cite{PhysRevX.7.021006}.

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.

\section{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~\cite{milz_kolmogorov_2017}.

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.~\eqref{def nonMarkov trDistance}. 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
\begin{gather}
\mathcal{N}_\text{coarse}\leq\mathcal{N},
\label{eq:coarsegraining}
\end{gather}
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 \textit{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} |\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 \textit{eigenstate thermalization hypothesis}~\cite{PhysRevLett.120.150603}, 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~\cite{arXiv:1704.06204}. 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.

\section{Discussions and Conclusions}\label{sec: discussion}
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~\cite{Schlosshauer2007, VegaAlonso}, 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~\cite{tamascelli2018}; it can always be bounded by a function of time scales in the system-environment Hamiltonian~\cite{Luchnikov2018}, 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.~\eqref{eq:coarsegraining} 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~\cite{Gogolin}, 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~\cite{fixedpoint}. Even if we do not characterise the equilibrium state here, e.g. as done for quantum states under randomized local Hamiltonians~\cite{Cramer}, 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.~\ref{Fig interventions}(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.~\cite{PopescuWinter} 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.

\begin{acknowledgments}
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.
\end{acknowledgments}