{'fraction_non_alphanumeric': 0.040213965039640844, 'fraction_numerical': 0.01103257235648104, 'mean_word_length': 4.199528184752918, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Markovian and non-Markovian dynamics induced by a generic environment

We study the open dynamics of a quantum two-level system coupled to an environment modeled by random matrices. Using the quantum channel formalism, we investigate different quantum Markovianity measures and criteria. A thorough analysis of the whole parameter space, reveals a wide range of different regimes, ranging from strongly non-Markovian to Markovian dynamics. In contrast to analytical models, all non-Markovianity measures and criteria have to be applied to data with fluctuations and statistical uncertainties. We discuss the practical usefulness of the different approaches.


INTRODUCTION
Open quantum systems have been of interest for a long time [1]. The interest stems from the natural separation of a quantum system into a central system of interest, and an uninteresting or uncontrollable part, usually denoted by environment. In 1967, Lindblad [2], Gorini, Kossakowski and Sudarshan [3] arrived at the so-called Lindblad master equation to describe the evolution of a central system weakly interacting with a memoryless environment. This equation has been of paramount importance in the field as one can, both, analytically solve several instances of the equation [4], and describe accurately a wide range of experimental situations [5]. The dynamics produced by such an equation is called quantum Markovian dynamics. Recently there has been an effort to classify and understand systematically open quantum systems which lie outside this description.
Definitions and measures for quantum non-Markovianity (NM) have received considerable interest in the last 10 years or so [6,7]. They are meant to characterize quantum processes (viz. the dynamics of open quantum systems) which cannot be described by a master equation with constant Lindblad operators. For such systems, there might be hierarchy of stronger/weaker non-Markovianity in the sense that it is easy/not-easy to describe the quantum channel by some effective evolution equation within the systems state space alone. A second quite different idea is that NM is a feature which might be taken advantage of in order to perform certain tasks. Both are valid points of view, but it is not clear to what extend some of the popular definitions and measures of NM provide relevant information with respect to these questions. Several reviews of the field can be found in Refs. [8][9][10][11].
Non-Markovianity has been studied extensively, for quantum processes with analytical solution, implying no fluctuations or uncertainties, applications to more realistic processes are relatively rare [12,13].
In the present paper, we study a quantum channel derived from the coupling of a two-level system to a "generic" quantum environment. We use random matrix theory (RMT) to describe that environment. Choosing the Hamiltonian in the environment from the Gaussian unitary ensemble (GUE), we find a variety of different behaviors in the relevant parameter space, ranging from Markovian (Lindblad-Dynamics) to strongly NM behavior. This model is ideal to discuss the questions raised above. Since the model has no known analytical solution, all criteria and measures must be calculated numerically, with unavoidable statistical errors, which resembles an experimental situation in which finite statistics come into play.
In Sec. II we present the model and find the structure of the quantum channel describing the dynamics of the system. In Sec. III, we introduce the NM measures, which we will use in Sec. IV for analysis. In Sec. IV we compare and interpret the different measures for the whole available region in the parameter space. We finish with some closing remarks in Sec. V.

II. MODEL
The model to be used has been introduced in Ref. [14], which focussed on the derivation of an analytical description in the linear response regime. We find it suitable for the purpose of the present work, since it is a generic model, which, nontheless shows a broad range of different behaviors with respect to quantum NM. In this section, we describe the Hamiltonian of our system, and the quantum channel formalism used for a complete description of its dynamics.

A. The Hamiltonian
We consider a two level system (qubit) coupled to an environment. The Hilbert space of the qubit will be labeled by the subindex c and that of the environment by arXiv:1906.10750v1 [quant-ph] 25 Jun 2019 the subindex e . We assume that the dynamics in the whole Hilbert space is unitary, with the evolution governed by the Hamiltonian All subindices of operators in the right hand side indicate the subspace in which they act, except for σ z , which is a Pauli matrix acting on the qubit. ∆ is the level splitting in the qubit and the parameter λ controls the strength of the coupling between the central system and the environment. The first two terms in Eq. (1) represent the free evolution of both central system and environment, while the third term provides the coupling which is assumed to be separable. The Hamiltonian of the environment H e shall be chosen from the Gaussian unitary ensemble to provide generality to the results discussed here [15]. We measure time in units of the Heisenberg time of H e and energy in units of the means level spacing in the center of the spectrum of H e ; in such units, = 1. The density matrix of the central system for a time t is given by where we have chosen a product state as the initial state of central system and environment. Just as H e , V e is also chosen from the Gaussian unitary ensemble, and the angular brackets denote an ensemble average over both random matrices. The magnitude of the matrix elements The most general form of v c has a parallel and perpendicular component, with respect to the internal Hamiltonian σ z . If there is only a parallel component (v c ∝ σ z ) the qubit dynamics becomes dephasing. The channel acting on the qubit can be obtained in terms of the fidelity amplitude for the Hamiltonians H e ± λV e . This is already a very rich case; however, it has been considered before [16,17]. The other limiting case is when the coupling is perpendicular to the internal Hamiltonian, which is the case we are studying here; see also [14]. Thus we simply set, without loosing any generality, where σ x is one of the three Pauli matrices.

B. Quantum channel formalism
We describe the reduced dynamics of the qubit with the quantum channel formalism, which means that the evolution of the system state is described in terms of a linear time dependent map Λ t acting on the space of density matrices S(C 2 ) of the central system. The map Λ t takes an arbitrary initial state, and returns the state evolved according to Eq. (2) for a time t, Since the image of this map is a density matrix, Λ t has two properties: (i) it preserves the trace of the argument and (ii) is completely positive. Maps with these characteristics will be referred to as quantum channels. We will also be interested in more general linear operators that preserve the trace but, even though map hermitian operators to hermitian operators, are not necessarily completely positive. In fact, we shall consider maps that are generally non positive. We will call such maps quantum maps. In this language, a quantum channel is a quantum map, but not necessarily the other way around. A quantum map K, and any linear map, is determined by its action on a basis; consider the computational basis {|a b| } a,b∈{0,1} and arrange the resulting elements in the matrix This is the so called Choi-matrix representation [18,19].
Since K maps hermitian matrices to hermitian matrices, C K must also be hermitian. It can be seen that the Choi matrix can be obtained by applying the extended map id ⊗ K to a Bell state in the Hilbert space of two qubits: where |Bell = (|00 + |11 / √ 2. Thus, if K is a quantum channel, 1/2 C K is a two qubit density matrix, which in turn implies that C K ≥ 0, i.e. all of its eigenvalues must be positive or equal to zero. We will now outline the procedure to construct the Choi matrix, based on the dynamics. The procedure will rely on two properties of our particular channel. First, it is unital, which means that the identity is mapped onto the identity. Second, the evolution of a diagonal operator remains diagonal, i.e. if c is diagonal, so is c (t) in Eq. (4). The proof of both properties is found in the Appendix A. Let ρ x,y,z c be the density matrices associated with the +1 eigenvalues of the corresponding Pauli matrices. From the second property and trace conservation, we find that and from the first property, plus linearity of quantum maps, we get for some time dependent functions r, z x and z y . From these equations, one can see that the dynamics in the z axis is decoupled from the ones in the xy plane. If we define z 1,2 = (z x ∓ iz y )/2, we can write directly Notice that r| t=0 = z 1 | t=0 = 1, and z 2 | t=0 = 0.

III. NON-MARKOVIANITY MEASURES
We consider NM as a property of a quantum process Λ t , which is a one-parameter family of quantum channels with t ∈ R + 0 and Λ 0 = id. The NM criteria and measures used here are based on two different concepts, (i) divisibility and (ii) contractivity. Both of them require knowledge of the intermediate quantum map In Appendix B 1, we calculate the Choi representation of this intermediate quantum map with the following result: . The parameters r , z 1 and z 2 are the same as r, z 1 and z 2 but calculated at a time t + ε. When D = 0 or 2r − 1 = 0, Λ t is not invertible, and therefore Λ t+ε,t may not exist.

A. Divisibility
A quantum process Λ t is divisible if and only if for any t, ε > 0 it holds that Λ t+ε can be written as the composition Λ t+ε = Λ x • Λ t , with Λ x being a valid quantum channel. Here, Λ x can be identified with the intermediate quantum map Λ t+ε,t given in Eq. (10). Hence the divisibility of a quantum process is equivalent to all intermediate quantum maps being valid quantum channels. Formally speaking, the quantum process Λ t is divisible if and only if (a) Λ t is invertible for almost all t ∈ R + 0 , and (b) ∀t, ε > 0 : Λ t+ε,t = Λ t+ε • Λ −1 t is a valid quantum channel if it exists.
In condition (a), we allow Λ t to be non-invertible at a finite (countable) number of points in time. In condition (b), we check the complete positivity of the intermediate map only for those t, where Λ t is invertible.
a. RHP-Markovianity One of the definitions of NM is given in terms of the divisibility of the quantum process under consideration. Following Rivas et al. [20], we call a quantum process Λ t , RHP-Markovian if and only if the two conditions above are fulfilled.
To check for the complete positivity of the intermediate quantum map, Rivas et al. consider the trace norm of the associated Choi matrix defined as C t+ε,t 1 = j |λ j |, where λ j are the eigenvalues of C t+ε,t . Since the sum of the eigenvalues always is equal to the dimension of the Hilbert space of the physical system, due to trace preservation, any negative eigenvalue will necessarily lead to C t+ε,t 1 being larger than two. In addition, since the composition of two quantum channels is again a quantum channel, it is sufficient to check complete positivity for infinitesimal ε, only. Hence, Rivas et al. define the function which is zero if the intermediate map is completely positive, and greater than zero otherwise. Moreover, to define a measure of the degree of non-Markovianity of the process, the authors integrate this function over time. We shall label this quantity as Notice that N RPL (t) = 0 if and only if the process is divisible for all times up to t. b. Application to our model In order to check the RHP-Markovianity (divisibility) of Λ t as defined in Sec. II, we use the Choi representation, Eq. (11) of the intermediate quantum map Λ t+ε,t . The eigenvalues of C t+ε,t are Hence, the eigenvalues are non-negative if and only if (i) |Z 1 | ≤ q and (ii) |Z 2 | ≤ 1 − q.
Since we can limit to infinitesimal ε, we expand the different functions in Eq. (11) around t and obtain, to first order, where we have used the fact that r = r + εṙ, z 1 = z 1 + εż 1 , and z 2 = z 2 + εż 2 . The two conditions (i) and (ii) can now be written as where we introduced and Additionally we used the fact that for any complex number c, |1 + εc| = 1 + εc + O(ε 2 ) on the expression for |Z 1 |. Finally, we combine the two inequalities into We can now relate our inequalities to the criterium of Rivas et al. as follows. In our case, the trace norm of the Choi matrix, can be written as This yields We showed that non-negativity of the eigenvalues is equivalent to the double inequality δ 2 ≤ δ q ≤ δ 1 . Then we saw that it is also equivalent to g(t) = 0. This means that the double inequality holds if and only if g(t) = 0.

B. Contractivity
Markovianity of classical stochastic processes imply that probabilities distributions decrease their Kolmogorov distance with time [21]. This is interpreted as a loss of information of the initial conditions. Carrying this ideas to a quantum level results in a definition of Markovianity [22]. Let ρ 1,2 (t) denote the evolution of two states ρ 1,2 . We define where T [ 1 (t), 2 (t)] = tr(| 1 (t) − 2 (t)|)/2 is the trace distance, which is directly related with the probability of distinguish the state 1 (t) from the state 2 (t), i. e., it is their distinguishability [23], and |A| = √ AA † . In other words, σ is the derivative of the distance between the evolved states. We say that a process is contractive if for all ρ 1,2 and all t ≥ 0, we have that σ(ρ 1 , ρ 2 , t) ≤ 0. A process is said to be non-Markovian if it is not contractive. Breuer et al. then define the following quantity as a measure for the degree of non-Markovianity: The calculation of this measure is greatly simplified when the process acts on a qubit. To perform this maximization, one should consider only pure, orthogonal initial states [24]. Indeed, we found that T [ 1 (t), 2 (t)] only depends on the vector difference between the representations of initial states in the Bloch ball (see Appendix B 2). Moreover, the distance between the two points representing the initial states enters as a homogeneous scale factor. It therefore possible, restricting the maximum search to such cases, where ρ 1 is a pure state, and ρ 2 the uniform mixture. If the pure state ρ 1 is parametrized in spherical coordinates by the angles θ and φ, we obtain [cf. with In what follows, we derive a criteria for NM based on the contractivity, which can be compared to the criteria Eq.
for all times. Condition (27) becomes −ṙ(2r − 1) = δ q (2r − 1) 2 ≥ 0 which in turn is equivalent to δ q ≥ 0 (at least as long as 2r − 1 = 0, i.e. away from the points where Λ t is not invertible). To consider condition (28), where From this, it is clear that the largest time derivative of To summarize, for the process to be Markovian (in the sense of contractivity), it is required that both, δ q ≥ 0 andṀ max ≤ 0. This is equivalent to where This double inequality is the analog of Eq. (21), which has been derived as criterium for Markovianity in the sense of divisibility. For unital maps like the one considered here, contractivity of the trace distance can be identified with positivity. That means that our quantum process Λ t is contractive if and only if all intermediate maps are positive. Note that divisibility is defined as all intermediate maps being completely positive. Thus, there may be processes which are contractive but not divisible [25]. Therefore, we may find regions in the parameter space of our system, where the dynamics is contractive (i.e. BLP-Markovian) but not divisible (i.e. RHP-Markovian). A recent, comprehensive discussion on different criteria for divisible and contractive processes can be found in Ref. [26].

C. Maximal recovery
This quantifier of non-Markovianity can be based on any capacity-like property of the channel. In fact, we use distinguishability, maximized over states, just as the BLP-measure. However, instead of summing up all the small increments, we search for the maximum distinguishability recovery over the whole quantum process [27]: where T ( 1 , 2 ) is defined in the text below Eq. (24). The BLP-measure has been related to the backflow of information, which can be quantified indeed in terms of the recovery of distinguishability. However, in order to quantify the amount of information recovered, it is much more sensible to use N MDR , than integrating over all backflow in a process, where information is fluctuating back and forth between system and environment. Of course there is a price to pay. It is quite more expensive to compute N MDR , than it is to compute N BLP , because of the additional degrees of freedom, t 1 and t 2 .

IV. NUMERICAL SIMULATIONS
In this section, we apply different methods to characterize the (non-)Markovian dynamics of a generic open quantum system. The quantum process to be studied is obtained by numerical simulations of the Hamiltonian in Eq. (1) with the initial state of the environment taken as the maximally mixed state. This includes a Monte-Carlo sampling of a random matrix ensemble. Therefore, the numerical data are contaminated by residual statistical fluctuations due to the finite size of the sample. We stress that residual fluctuations would be present, in experimental situations, also. The dimension of the Hilbert space of the environment is set to N = 200 and the sample size is fixed to N sam = 2400 unless otherwise stated. In our model system, defined in Eq. (1), we can identify three different energy scales. The average level spacing in H e (which is set equal to one), the spacing ∆ between the two levels of the qubit, and the coupling strength λ between qubit and environment. In an effort to explore the properties of our model as thoroughly as possible, we consider a rectangular region in the parameter space (∆, λ). For ∆, the interesting range reaches from values much smaller than the mean level spacing to values much larger, thus we choose the limits 0.016 ≤ ∆ ≤ 16, covering a range of three orders of magnitude. For λ, we choose the lower limit where λ is much smaller than the average level spacing in H e , such that perturbation theory would be applicable for the Hamiltonian H e + λV . The upper limit, by contrast, is dictated by the requirement that our model reproduces the behavior of the random matrix model in the limit N → ∞, in such a way that finite size effects are still negligible. For that to hold, the spreading of eigenstates of H e + λV expanded in the basis of H e must be small compared to the spectral range of H e . Both considerations lead us to the limits 1/32 ≤ λ ≤ 1/2. A similar parameter range has been explored in Ref. [14].
In order to obtain reliable values for the two measures, we establish a finite ending time for the processes to be studied. At that time, the system will have relaxed so much that the remaining dynamics is unusable for any practical purposes. Our approach for finding a sensible definition for the ending time t End is described in Sec. IV A. In Sec. IV B we present and discuss the results for three non-Markovianity measures. Moreover we select several representative parameter sets that are analyzed in detail in Sec. IV C. In particular, we study the dependence of the measures on the number of samples considered. The motivation is to understand the statistical significance of the results presented in Sec. IV B. At the end, this could be useful for identifying a NM measure which is accurate, robust and significant. We finish this section by analyzing the local (in time) criteria for non-Markovianity in Sec. IV D. We consider divisibility and contractivity, via the corresponding inequalities Eq. (21) and Eq. (32). Both conditions are tested for the infinitesimal intermediate quantum map Λ t+ε,t as defined in Eq. (10), and examine carefully the usefulness of such kind of expression under statistical fluctuations.

A. Process ending time
For the NM measures to be considered below (Sec. IV B) it is essential to define an ending time t End for the quantum process in question. However, two conflicting requirements arise. On the one hand the ending time should be sufficiently large, such that the dynamics of the process is completely contained, but on the other hand it should also be sufficiently short, such that the contaminating contribution from residual statistical fluctuations remain small.
The quantum process studied here, has the convenient property that if one chooses as initial state an eigenstate of σ y , the system converges to the uniform mixture in the limit of long times. Along this process, the purity P (t) = tr[ c (t) 2 ] decays from P (0) = 1 to P (∞) = 1/2. We choose the ending time for the process at that time, where the purity of Λ t [ρ y c ] is equal to 0.51, which means that the purity has decayed to 2% above its minimum value. Of course other values of the same order are equally possible, but they do not change the general results of our study.
In Fig. 1 we show the process ending time as a function of the parameters ∆ and λ, color coded over the parameter space. While we use a linear scale for λ, ∆ is varied on a log-scale. The resulting ending time varies over several orders of magnitude, so we also use a log-scale for the color mapping. While λ increases, the t End becomes smaller exponentially fast, since for large λ, the Fermi golden rule approximation applies [14]. It is however quite remarkable that for small λ, the largest t End can be found near ∆ = 0.16, which is approximately equal to t −1 H , where t H = 2π is the Heisenberg time in the random matrix environment. In other words, at ∆ = 0.16 the period of the Rabi oscillation is equal to the Heisenberg time.
For all non Markovianity measures we choose t = t End unless otherwise stated, and we shall thus drop the time dependence.

B. Three measures for non-Markovianity
Analyzing Fig. 2, we find that both measures reach smallest values in the upper left corner of the parameter space, shown. Indeed, due to the following argument, we expect the quantum process to be at least close to Markovian in this region. The standard prescription for deriving a quantum master equation via the Born-Markov approximation consists in the following steps [28]: (i) Couple the central system weakly to each of the many degrees of freedom in the environment (Born approximation), (ii) let the number of degrees of freedom in the environment go to infinity, and (iii) assume the environment correlation functions to decay almost instantaneously on the time scale of the reduced dynamics (Markov condition). In terms of level density and average local level spacing, condition (i) and (ii) lead to a wide range in energy with an exponentially high level density, which means that the perturbation strength will be large as compared to the level spacing. This regime is known as the Fermigolden-rule regime [29]. Note that in this parameter region, ∆ 1, which results in a slow system dynamics so that condition (iii) is fulfilled.
By contrast, for sufficiently small coupling λ < 7/32 and not too small ∆, the dynamics is clearly NM. It is clear that in this region at least some of the conditions mentioned above are not fulfilled. The region of strongest NM behavior is around ∆ ≈ t −1 H and small values of λ. An interesting area is in the upper right region of the NM maps in Fig. 2. There, the BLP measure of Eq. (25) tends to very small values, while the RHP measure, Eq. (13), remains constant. As explained at the end of Sec. III B, the RHP criterion is more restrictive than the BLP criterion, as it requires complete positivity and not just positivity for the intermediate maps. It is thus possible that a given quantum process is BLP-Markovian but not RHP-Markovian. Due to residual statistical fluc- tuations, a definite judgement is difficult. Fig. 3 shows N MDR (t) as defined in Eq. (34) up to the ending time t = t End as a function of ∆ and λ. Its behavior is more similar to N BLP than to N RHP , which may be related to its common origin. Note though that the boundary between the regions of Markovian and strongly non-Markovian behavior is sharper. The region of strong NM is also somewhat larger, located in an area parallel to the M-NM boundary, at values for λ below 1/4, of about 12 blocks in size. Finally, note the region in the upper right corner. There, the MDR measure leads to (relatively) lower values than the BLP measure, in distinction to the RHP measure which reaches much larger values.

C. Robustness and convergence of the NM measures
In the previous section, we presented the results of three different measures for NM. For our model, we found a wide range of different behavior, depending on the choice of the parameters, ∆ and λ. Here, we study the robustness and accuracy of the measures in more detail. For that purpose we select three points in parameter space, where the behavior of the quantum process is quite different: • Point P1 (∆ = 10 −1.4 , λ = 3/8), where the dynamics is Markovian, or at least very close to it.
• Point P3 (∆ = 10 0.8 , λ = 3/8), where the dynamics looks like being more NM according to RHP than according to BLP or MDR -compare the two heat maps in Fig. 2.
Let us now discuss the numerical results. In Fig. 4, the NM-measures are shown as a function of time. In other words, we compute the NM-measures as if the quantum process would end at time t, instead of t End . It is clear from the definition of all three measures, that they must be monotonously increasing: For unambiguously whether the process is Markovian or non-Markovian. For P2 (strongly NM; middle panel), the RHP-measure increases to very large values of the order of 10 2 . The BLP measure also increases along the full time range up to values of the order 10 0 , while the MDR measure quickly saturates at a value of the order of 10 −1 (note that the MDR measure is by definition limited to values below one). The different behavior between BLP and MDR will be discussed below, where we consider the criteria for contractivity. In any case, all three measures clearly show the NM of the process.
For P3 the RHP-measure increases continuously as in the previous cases, reaching values of the order of 10 0 . By contrast, the other two measures, BLP and MDR remain below a value of the order of 10 −4 . This may hint towards the possibility that here, the quantum process is P-divisible but not CP-divisible.
In Fig. 5, we plot the NM-measures versus the sample size, N sam , where we expect that the NM-measures approach a limit value for N sam → ∞, the true ensemble average. For point P1 (top panel), as the ensemble size increases, all measures tend algebraically to zero. For point P2 (middle panel) all measures converge to a finite values, whereas for point P3, the measures based on distinguishability seem to drop to zero while the one based on divisibility attains a finite value. In other words, this result suggests that the dynamics is P-divisible but not CP-divisible at that point [25].

Environment size
We also experimented with different environment sizes, but the results did not change significantly. In fact, it can be shown that this may have an effect at short times only; in general, one would expect finite size effects of N e at times of the order of 1/N e in units of the Heisenberg time. These are not our concern, since we define our model in the limit N e → ∞.

Non-equivalence of NM measures
A second interesting question is that of "quantitative equivalence". Two measures M 1 , M 2 for a physical property may be called (quantitatively) equivalent, if and only if for any two states A, B of some system. For example, if one thermometer (calibrated according to the empirical temperature M 1 ) finds that a body A is colder than a body B, any other thermometer (calibrated according to some different empirical temperature M 2 ) should find the same relation. Analyzing carefully our results in Fig. 2, we can indeed find pairs of points in parameter space, where the two NM-measures N RHP and N BLP violate this condition. For instance, N RHP is clearly larger at P4 than at P2, while in the case of N BLP it is just the other way round.

D. Time evolution of Markovianity criteria
Here, we study the behavior of the two time-local criteria for non-Markovianity which are based on divisibility and contractivity, described in Sec. III.
For the divisibility, we consider the requirements given in Eq. (21) on the one hand, and the condition g(t) = 0 with g(t) given in Eq. (12). While formally, both criteria are equivalent [see Eq. (23)], we will see that in the presence of experimental/statistical uncertainties, one might be easier to verify than the other. In Fig. 6 we present the aforementioned expressions for points P1, P2 and P3 for times ranging from zero to t End . For point P1 (top two panels), we see that the inequalities are saturated, in the sense that δ 1,2,q are apparently all equal. Despite that, the function g(t) is not identically equal to zero, however it is small. The inequalities have allowed to correctly identify the point as Markovian, in agreement with the analysis of Fig. 5. In fact, for larger ensemble sizes the value of g(t) diminishes.
The behavior of δ 1,2,q and g(t) for point P2 can be seen in the middle panels of Fig. 6. Here, the different curves corresponding to δ 1,2,q cross each other in a systematic and ordered fashion, indicating non-Markovianity beyond statistical fluctuations. This is also seen in the behaviour of the corresponding g(t), which oscillates regularly around values of the order of 0.1. Notice that we identified this point as displaying non-Markovian behavior, again with the aid of Fig. 5.
Finally, point P3 is studied in the lower panels of Fig. 6. We can see that the curve corresponding to δ q is not between the ones corresponding to δ 1,2 . The conclusions are confirmed by the behavior of g(t), where one finds notorious fluctuations on top of a smooth curve which increases systematically as a function of time. Thus, according to the divisibility criterion, the system is non-Markovian, as concluded from the lower panel of Fig. 5.
In the case of the Markovianity measure based on the contractivity of the process, the condition σ ≤ 0 translates to Eq. (32) for the channels here considered, see Eq. (9). For a given time t, a certain initial pair of states ρ (max) 1,2 will yield the maximum N BLP (t) as defined in Eq. (25). For these states, one can calculate and see how the final value of the measure is built with time. Notice that σ max is different from the derivative of N BLP (t), recall Eq. (25), as for each ending time t, the states that maximize the quantifier are different, whereas in Eq. (35) we fix the ending time; however, for t equal to the ending time, they coincide. On the top left panel we can see that Eq. (32) is apparently fulfilled during the whole process, however, since the two curves for δ C 1 and δ C 2 lie on top of each other, the inequality δ C 1 ≥ δ C 2 might be violated on a smaller scale. On the top right panel, we see that indeed the measure is close to zero; for t < 0.7 it is numerically zero, and afterwards, close to 10 −8 . From this evidence, we arrive at the conclusion that the point is Markovian, with respect to contractivity, in agreement with the same case studied in Fig. 6.
The point P2 is analyzed in the middle panels. The function δ q oscillates around zero with decreasing amplitude, while δ C 2 provides an upper bound for δ C 1 . Indeed, whenever this bound is saturated, δ q has an node, and σ max has a minimum, making the system "very Markovian". On the other hand, when the difference between δ C 1 and δ C 2 is largest, σ max has a maximum, and the system becomes very non-Markovian. Thus, as can be seen on the right, σ max oscillates regularly and with a relatively high amplitude around zero. The NM measure, N BLP adds up the areas below the positive parts of σ(t), which we expect to have a similar behavior as σ max (t), shown here. Therefore, the point P2 shows a genuine non-contractive behavior and the dynamics is NM in this case.
The point P3 is analyzed in the bottom panels. All three functions δ q and δ C 1,2 display a similar behaviour as for point P1, but for δ q we observe stronger statistical fluctuations. In this case, however, σ max oscillates, and has a maximum in 0. The measure picks up small statistical fluctuations which diminish as we increase the sample size and may therefore be regarded as spurious. We can conclude that the process for point P3 is contractive, in agreement with previous conclusions.
The criteria provided in Eqs. (21) and (32), indeed provide a usefull tool to understand if a certain non-zero value for one of the NM measures should be regarded as statistically significant or not. In some cases it is helpful to analyze the behavior of the measures under variation of the sample size in order to arrive at the correct decision.

V. CONCLUSIONS
In this article we studied a qubit coupled to a generic environment modeled by random matrices. The model Hamiltonian contains a factorizable interaction between qubit and environment, and provides the qubit with an internal dynamics perpendicular (in the Bloch representation) to the one induced by the interaction. This induces a channel structure for which we were able to derive analytical conditions for several criteria of Markovianity. In spite of its simplicity, the model displays rich dynamics in the qubit, beyond pure depolarization or dephasing.
We then applied the criteria to determine for which parameters the model yield Markovian dynamics in the qubit. We found several difficulties with verifying Markovianity criteria for the numerical data, which we expect to appear also in real experimental situations. Fluctuations due to noise, and/or due to finite sample sizes may contribute notably to the finite value of a non-Markovianity measure and thereby suggest non-Markovian behavior, whereas the clean system really is Markovian. An analysis like the present one, where the ensemble size is increased such that residual fluctuations diminish, eventually reveals the true behavior.
Acknowledgements-Support by projects CONACyT 285754, UNAM-PAPIIT IG100518, Benemérita Univer- In this section, we prove that the quantum channel, describing the dynamics of the qubit under coupling to the RMT environment, has the form of a X-state [30] as postulated in Eq. (9). For doing so, we use the entire expansion of the evolution of system and environment in a Born series. We consider the more general case of an arbitrary mixed initial state e in the environment, and only at the end specialize to the case e = 1 1/N e .
In the interaction picture, a solution to the Hamiltonian may be written as Here, we stick to the convention chosen in the main part of this paper, where the time variable t measures time in units of the Heisenberg time, which results in = 1. The echo-operator M (t) describes the evolution of the state χ(t), such that χ(t) = M (t) χ(0). In the original Schrödinger picture, it thus holds: As a formal solution of the evolution equation in the interaction picture, the echo operator fulfills the following integral equation Here, the second line represents the afore mentioned Born series, where the interaction has been written in blockmatrix notation:σ Here, we will calculate the average of˜ (t) = M (t) c ⊗ e M (t) † , with respect to the random matrix ensemble for V e . This defines the quantum map Λ t as and thereby its Choi-matrix representation, in Eq. (9). Below, we will also consider this quantum map in the interaction picture, defined asΛ t : c (0) →˜ c (t) = tr e ˜ (t) . (A3) a. Simplified notation For our purpose, it is convenient to introduce the following more compact notation for the multi-dimensional time integrals: Note that we enumerate the different time variables from smallest to largest starting at the time closest to zero. Next, we will introduce an independent notation for the integrands. Let us consider the caso of an even number of product terms, first: sinceṼ e (τ ) is Hermitian. Note that the product terms on the LHS of the first equation, must be ordered according to decreassing time arguments, just as in Eq. (A1). From an explicit computation we find Note that the product terms must be ordered such that time increases from right to left. For an odd number of terms: With this, we may write for the echo operator b. Ensemble averaged quantum channel Averaging over the random matrix V e implies that only such terms survive, which contain an even power of matrices A(τ k ) and A(σ k ). This means that the indices of summation must either be both even or both odd. Therefore, For the quantum channel in the interaction picture, Eq. (A3) it is now easily verified that This and Eq. (A2) then show that Λ t [ |0 0| ] is indeed of the form postulated in Eq. (9), and it yields the following expressions for r(t): The second equation results from the fact that the reduced evolution of the qubit conserves the trace. Naturally, it is difficult to prove this directly, from the expressions derived here.
Again, the resulting qubit state is diagonal, however unless we specialize to the case e = 1 1/N , the functionr(t) is different from r(t) corresponding to the previous case.
We continue with the off-diagonal blocks of the Choi-matrix representation of the quantum channel: This result confirms again the general X-state structure of the Choi-representation of our quantum channel. In terms of the parametrization in Eq. (9), we find: where the phases e ±i∆t arise from returing to the Schrödinger picture, according to Eq. (A2). Finally, For any Hermiticity conserving linear map, the Choi-representation itself must be Hermitian, in particular also for our quantum channel, as can be seen from its definition in terms of the reduced dynamics in Eq. (2). This implies that In the case of z 1 (t) the equivalence of the Eqs. (A15) and (A17) is rather obvious. One simply has to exchange the variable names k and k along with τ and σ. In the case of z 2 (t), we do not see any simple way of proving the equivalence in terms of the present approach.
c. Collecting results for the case e = 1 1/N From the considerations in the previous paragraph, we found that the Choi-matrix representation of the quantum channel defined by the Eqs. (2) and (A2) is given by The general X-state structure (i.e. all the zeros in this matrix) follows directly from the considerations in this section. For simplicity we omitted the time argument in this representation. The functions r(t), z 1 (t), z 2 (t), andr(t) are as defined above. Some of the dependencies between these matrix elements could not be proven within our derivation, but are valid due to elementary properties of the evolution equation (2). This is the case for trace conservation (diagonal blocks) and Hermiticity (off-diagonal blocks). As a last point, we show that for e = 1 1/N e , it holds that r(t) =r(t). In this case, and for e (0) = 1 1 2 , we find from Eq. (2): This means that no matte if we perform an ensemble average or not, the resulting quantum channel is unital (it maps the identity onto itself). This in turn implies: Appendix B: Two representations of a quantum channel We will use two different representations. (i) The "super operator" representation which is nothing else than the standard matrix representation of a linear map on a vector space. This representation is easy to read directly from the evolution of a standard set of states; in fact, we shall construct it in that way. (ii) The Choi-matrix representation in which verifying inherent quantum channel properties is easier.
For the super operator representation, we represent the density matrices as column vectors, in the so-called "antilexicographical" ordering [31]. For a single qubit, we have This fixes the matrix representation L t of the map Λ t , since the columns of L t must be the images of the canonical basis vectors.
The super operator representation has the convenient property, that the representation L of the composition of two quantum maps, Λ = Λ 2 • Λ 1 (where Λ 2 is applied to the result of Λ 1 ) is simply given by the matrix product of the representations L 1 and L 2 of the individual maps: L = L 2 L 1 .

Quantum map for intermediate time steps
We know that the superopertor has the following form where r and z 1,2 are functions on time.
For given quantum maps Λ t and Λ t+ε , we compute the quantum map which takes states (t) form time t to t + ε. For the moment we assume ε > 0 to be finite. The central question is whether Λ t+ε,t = Λ t+ε • Λ −1 t is a CPmap, or not. While the super-operator representation L t is appropriate to compute the composition of Λ t+ε and Λ −1 t , the Choi representation is needed for verifying the complete positivity. From L t given on (B3) of the section II we have its inverse matrix where d = r 2 − (1 − r) 2 = 2r − 1 and D = |z 1 | 2 − |z 2 | 2 . If we denote with primes the functions evaluated on t + ε, that is, r = r(t + ε) similarly for z 1 and z 2 , then We found that the super-operator representation corre-sponding to Λ t+ε,t is given by where q = (r − r − 1)/d, Z 1 = (z 1 z * 1 − z 2 z * 2 )/D and Z 2 = (z 2 z 1 − z 1 z 2 )/D. The corresponding Choi matrix turns out of reshuffle the above matrix: on subsection III A the divisibility of Λ t+ε,t is explored via the positivity (non negative eigenvalues) of Choi matrix C t+ε,t .

Trace distance and contractivity
We have introduced the trace distance and its definition through of the eq. (24). Furthermore, turns out that the trace distance of an Hermitian matrix is equal to one half of the sum of absolute values of its eigenvalues.