On the qubits dynamics in random matrix environment

We consider the entanglement evolution of two qubits embedded into disordered multiconnected environment. We model the environment and its interaction with qubits by large random matrices allowing for a possibility to describe environments of meso- and even nanosize. We obtain general formulas for the time dependent reduced density matrix of the qubits corresponding to several cases of the qubit-environment interaction and initial condition. We then workout an analog of Born–Markov approximation to find the evolution of the widely used entanglement quantifiers: the concurrence, the negativity and the quantum discord. We show that even in this approximation the time evolution of the reduced density matrix can be non-Markovian, thereby describing certain memory effects due to the backaction of the environment on qubits. In particular, we find the vanishing of the entanglement at finite moments and its subsequent revivals (entanglement sudden death and entanglement sudden birth). Our results, partly known and partly new, can be viewed as a manifestation of the universality of certain properties of decoherent qubit evolution which have been found previously in various versions of bosonic macroscopic environment.


Introduction
The time evolution of a system interacting with a 'large' in certain sense environment (reservoir, bath) is a basic theme of statistical mechanics (macroscopic systems), quantum optics (small systems) and adjacent branches of physics, commonly known as the theory of open systems [1][2][3]. A new surge of interest to the theme has been motivated by quantum information theory, where the robustness of quantum phenomena, entanglement in particular, with respect to the decoherence effects of the environment is of primary importance [4,5].
One usually assumes that a state of the composite system consisting of a given system  and its environment  evolves unitarily via the Liouville-Neumann equation. The reduced density matrix r ( ) t S of  obtained by tracing out the environment degrees of freedom is usually dissipative and dephasing. The 'partial tracing' procedure makes the analysis of  r ( ) t quite non-trivial since it is rather difficult if possible to obtain in general an exact and efficient 'governing' equation for  r ( ) t . Among various approximations of the field one has to mention first the so-called Markovian approximation arising on the time intervals proportional to the inverse square of the system-reservoir interaction [1-3, 6, 7]. In this case one obtains a well defined local in time and memoryless dynamics in term of certain differential equation for the reduced density matrix. There is, however, a number of important problems where one needs to study the non-Markovian dynamics in which the environment acts as a memory allowing the earlier states of system to affect its later dynamics, see [8][9][10][11] for recent reviews.
One of these problems arises in quantum information theory and is a subject of considerable activity over the last two decades. It concerns the time evolution of quantum correlations, in particular, entanglement, a counterintuitive key property of quantum systems whose state cannot be written as the product of states of their subsystems, thereby being a quantum analog of the independency property of probability theory. The time evolution of quantum systems embedded into an environment degrades, in general, the quantum correlations.
Quantum information theory deals often with bipartite systems: 1 . 2

A B
A basic entity of quantum information theory is a two-state quantum system, known as the qubit. Its counterpart in the theory of open systems is a two-level system arising in a number of physical situations [1][2][3]6]. Thus, an important problem of the quantum information theory is the analysis of time evolution of (1.1), (1.2) where a are the qubits. The problem has been actively studied in the last decades in the frameworks of various models of the bosonic environment described by the so-called spin-boson type Hamiltonian and its approximations see [8,10,11,12,14,15] for reviews.
In this paper we study the entanglement dynamics of two qubits by using large random matrices that have been widely used to describe complex quantum systems of large but not necessary macroscopic size, in particular, quantum coherence effects, see e.g. [16][17][18][19][20] for results and references. We consider Hamiltonians which are similar to the two-qubit versions of spin-boson type Hamiltonians but have large random matrices instead of their bosonic counterparts, i.e., the environment Hamiltonian and the qubit-environment interaction Hamiltonian. The quantum dynamics and their analysis in our work generalize those introduced in [21,22], where one qubit embedded into a 'random matrix' environment has been considered Note that the random matrix models of the dynamics of two qubits and of more general quantum systems embedded into an environment have been already considered, see e.g. [23][24][25][26] and references therein. The main difference between our work and the above works is that they are confined to the Heisenberg time scale pertinent for the quantum chaos studies, while in our work we deal with a shorter time scale commonly used in nonequilibrium statistical mechanics and the theory of open systems. Besides, the above works focus mostly on the entanglement characteristics known as the purity and the fidelity, while we carry out the simultaneous comparative study of three entanglement characteristics: the concurrence, the negativity and the quantum discord, widely used in quantum information theory and allowing for the study of quite interesting quantum coherence phenomena. This makes essentially complementary the results of [23][24][25][26] and our results.
The paper is organized as follows. In section 2 we first briefly outline certain known facts on the qubit dynamics for one and two qubits and on the three widely used entanglement characteristic (quantifiers): the negativity, the concurrence and the quantum discord. We then describe our random matrix model of the twoqubit dynamics which can be both Markovian and non-Markovian depending on the value of parameters of the model and we comment on related earlier results in the literature. In section 3 we present our results both analytic and numeric on the time evolution of the above entanglement quantifiers paying special attention to the entanglement vanishing and revival (known as early stage disentanglement or entanglement sudden death (ESD) and entanglement sudden birth (ESB)) and their dependence on the Markovianity of the model in question.
We find a variety of cases of the entanglement evolution, partly known but found in different models of entanglement and partly new. This shows that the proposed random matrix models describe sufficiently well interesting properties of entanglement, thus can be used for further studies of quantum coherence (see, e.g. [27]). This also shows that the properties possess a certain amount of universality (or robustness) to the models used for their description.

Generalities
We will use the setting that has been worked out in a number of recent papers on the dynamics of a pair of qubits embedded in a sufficiently 'large' environment, see [8,12,14,15]  . In other words, every qubit has its own environment and its own interaction with the environment, hence, the qubits are dynamically independent. Here the entanglement between the qubits for > t 0 arises only because they are initially entangled (see the initial conditions (2.12)-(2.14). The Hamiltonian (2.3) can describe two initially entangled and excited two-level atoms spontaneously emitting into two different cavities, two sufficiently well separated impurity spins, say, nitrogen vacancy centers in a diamond microcrystal, etc.
(ii) Hamiltonian H F . Here  H is given by (2.2) and i.e., the first qubit is free (  ), but the second qubit is as in (2.3). Here also the qubits do not interact and their entanglement for > t 0 is possible only if they are entangled initially (see initial conditions (2.12)-(2.14) below). The free qubit is known as the ancilla or spectator in certain contexts of quantum information theory, see, e.g. [4,11,25,28].
(iii) Hamiltonian H C . Here again  H is given by i.e., the qubits have a common environment. Here the entanglement between the qubits for > t 0 is due the initial conditions (see the next subsection) but also due to the Hamiltonian providing an effective interaction between the qubits.
The Hamiltonian H F can be viewed as a particular case of H I in which the coupling of one of the spin with its environment is switched off. Likewise, the Hamiltonians H I and H C can be viewed as two limiting cases of a more general Hamiltonian describing two qubits which interact with an environment via a short range interaction. Indeed, denoting  r the radius of the interaction and d the distance between the qubits, we can say that the cases   r d and   r d of the general Hamiltonian are given by H I and H C (assuming additionally that the radius  r of qubit-qubit interaction is much smaller than  r ), see, e.g. [29] for more. The the density matrix of the composite (two qubits plus environment) is where  r ( ) 0 is a 4×4 positive definite and trace one matrix, the initial density matrix of the qubits, and  P is the projection on a pure state of the environment, often an eigenstate of  H . Note that we consider here and below the zero-temperature case, hence pure states for the environment appropriate for the quantum information theory setting, although mixed (thermal) states, the Gibbs distribution  in particular, can also be considered.
Define the reduced density matrix of  (qubits) as where  Tr denotes the partial trace with respect to the degrees of freedom of  . We have where the linear superoperator F( ) t is known as the influence (Feynman-Vernon) functional in the theory of open systems and as the quantum channels in quantum information theory.
An important property of the above three general two-qubit Hamiltonians (2.1)-(2.5) is that they are essentially determined by the one-qubit Hamiltonians is the free superoperator for the qubit A, see (2.4), and according to (2.7), (2.8) and (2.10) the superoperators of reduced evolution (quantum channels) are: . 2.11 This simplifies considerably the analysis of the evolution governed by H I and H F which is the case in several physical situations, see [1,8,12,15,24].
We will discuss the concrete forms of the above Hamiltonians in section 2.4. Here we just note that from the point of view of statistical mechanics and condensed matter theory the Hamiltonians H I and H F of (2.3) and (2.4) seem less interesting than the Hamiltonian H C of (2.5 ), since H I and H F describe non-interacting quantum systems. They are, however, of considerable interest for quantum information theory, since the dynamics determined by H I and H F allow for the study of the emergence of quantum correlations in a 'pure' form, i.e., without dynamical correlations due to the interaction between the qubits. In particular, it seems that the Hamiltonian H I could be a simple model appropriate for quantum computing, where qubits are independent in typical solid state devices. Besides, the dynamical independence of qubits can describe the absence of non-local operations in the quantum information protocols.
Dynamical correlations are the case only in the time evolution determined by Hamiltonian H C . The corresponding quantum channel F ( ) t C cannot be written directly via F A and F B , see (2.11). However, it follows from (2.1)-(2.5) that it is constructed from the same operators as in F A and F B . This facilitate considerably the study of the corresponding evolution (see [9,12,27] for this channel with bosonic environment and related quantum correlation).

Initial conditions
Here we describe three families of initial conditions  r ( ) 0 for the reduced density matrix (2.7) of two qubits. The conditions have been used in a variety of recent papers (see, e.g. reviews [8,12] and references therein). We write below ñ = |a a a , 1 2 1,2 for the vectors ñ Ä ñ | | a a 1 2 of the standard product basis of the state space of two qubits where ñ =  |a a , are the basis vectors of the state space of one qubit. We will also omit the subindex  in the reduced density matrices below.  . For a r = 1 W 3 k reduces to r Y k . According to (2.12)-(2.14) the corresponding states are not entangled if a = = ( ) n 0, 1, 1, 2, 3.

2.15
n It is easy to find that in the basis ñ = ++ñ ñ = +-ñ ñ = -+ñ all the above initial condition have the so-called X-form [8,12,13] * * r r r r r r r r which arises in a number of physical situations and is maintained during widely used dynamics (see [8,12,15,30] for reviews). The form is also maintained during the quantum dynamics governed by our random matrix Hamiltonians, see (2.47 )-(2.49) below.
In what follows we will call the model of the two-qubit evolution the pair consisting of one of Hamiltonians (2.1)-(2.5) and one of initial conditions (2.12)-(2.14). Thus, a particular model is denoted , , , 1, 2, 3 , 1, 2, 2.18 and for n=3 the value of k = 1, 2 from (2.14) has to be indicated.

One-qubit random matrix Hamiltonian
Since the system Hamiltonian  H in the total Hamiltonian (2.1) is always the Hamiltonian (2.2) of two qubits, we have to indicate the environment and the interaction Hamiltonians in general formulas (2.3)-(2.5). They are essentially determined by the one qubit Hamiltonians (2.9).
A well known and fairly realistic example of such a Hamiltonian is the spin-boson Hamiltonian describing a spin (two-level system) coupled to a 'large' (usually macroscopic) collection of quantum oscillators (bosons), e.g. the cavity radiation field in quantum optics, see [1,3,6]for this and other physical realizations. The Hamiltonian is is assumed to be smooth enough. The spin-boson model (2.19) proves to be quite non-trivial in general and admits a number of approximations, among them the Markov approximations is one of often used in various physical situations [1,3,6]. As in classical dynamical systems the approximation is the case if the environment characteristic time is much smaller than that of the system [31]. The rigorous analysis Markov approximation is given in [2,7].
In this approximation the superoperator t F( ) of (2.8) has the exponential form where τ is the 'slow' time (see, e.g. (2.29) below) and the superoperator can be viewed as a quantum analog of the Fokker-Planck operator of classical Markov processes [31]. There is, however, a wide variety of cases, including those of quantum information theory, where the Markov approximation is not adequate enough (see, e.g. recent reviews [9,11,12,15]). One can use in these cases the spin-boson Hamiltonian (2.19) and its two-spin (two-qubit) versions which are, however, rather complicated in general. Several interesting cases are based on exact particular solutions of the spin-boson model corresponding to one or two excitation in the composite system and on the Lorentzian form of the spectral density (2.20), see [12] for a review. Being rather special, this model spectral density provides an interpolation between the Markovian (l - V 1 1 , flat spectrum and weak coupling) and non-Markovian , narrow spectrum and strong coupling) evolutions in the frameworks of the same model. In this paper we will use a different asymptotically exactly solvable model for the qubit-environment evolution which is a two-qubit extension (in the sense of (2.1)-(2.5)) of the one-qubit model proposed and studied in [21,22]. In these works H bos and H dip in (2.19) are replaced by random matrices of large size. The resulting formulas for the asymptotically exact solution are also rather complicated even in the one-qubit case. Thus we will confine ourselves to a certain asymptotic regime, where the time is inverse proportional to the square of the coupling constant, see (2.29). This regime is an essential ingredient to obtain the Markov approximation for the spin-boson and similar models [1, 2, 7] but it yields in our case both the Markovian and the non-Markovian evolution depending on the form (flat or narrow) of the density of states of the environment (see formulas (2.23) below) which plays the role of the spectral function (2.20) in our model.
We will briefly describe the model following [21,22] and adding its properties not given in these works but important in the present work. Let M N be a N×N Hermitian matrix (random or not), be its density of states where n 0 assumed to be smooth enough (see (2.20)).
Furthermore, let W N be a random N×N Hermitian matrix distributed according to the matrix Gaussian law - In other words, the entries of Gaussian random variables such that where the symbol á ñ ... denotes the expectation. This is known as the Gaussian unitary ensemble (see, e.g. [16,20]).
We define the Hamiltonian of our composite system (a qubit plus an environment) as the following randoḿ where 1 l is the unit l×l matrix. Thus, this is a particular case of (2.1) in which (see (2.19 . It should be noted that the model yields the same results for any real symmetric or Hermitian random matrix in (2.25) whose entries    W j k N , 1 jk are independent and identically distributed random variables satisfying (2.24). However, in this case the techniques are more involved requiring an extension of those of random matrix theory for so-called Wigner matrices (see [20]).
In addition, our results are also valid in the case, where M N is random and independent of W N for all N. In this case we have to assume that the sequence { } M N of random matrices is defined for all  ¥ N on the same probability space (as well as the sequence { }) W N , is independent of { } W N and satisfies (2.23) for all typical realizations (with probability 1).
Having M N and W N random, we can view the second term in (2.25) as the Hamiltonian of a 'typical' N-level reservoir and the third term as a 'typical' interaction between the qubit and the reservoir. Note that the randomness (frozen disorder) is a basic ingredient of the theory of disordered extended systems. Its successful and efficient use is justified by establishing the selfaveraging property of the corresponding results, i.e., their validity for the overwhelming majority of realizations of randomness for macroscopically large systems, see, e.g. [32,33]. This guarantees the representativity of expectations of extensive observables allowing one to obtain an adequate description of the system in question by using only expectations. For quantum information aspects of the selfaveraging property see recent review [34].
In our case the selfaveraging property follows from the bound [22,35] for the variance of the entries of 2×2 reduced density matrix r r = ab a b= where C is a numerical constant. According to the bound the selfaveraging property of r ( ) ( ) t N holds for large N and any fixed t and even for large t if  t N . The above suggests that in our model the limit  ¥ N plays the role of the macroscopic limits in statistical mechanics. Note however that in statistical mechanics the density of states of a macroscopic many-body system (macroscopic environment in particular) is exponential in its volume and/or in its number of degrees of freedom since the exponent is the entropy of macroscopic system, the bosonic environment in our case. On the other hand, the density of states of our random matrix environment is linear in N, see (2.23), a common situation in one-body problems of condensed matter theory. Thus, interpreting N as the volume of reservoir, we have to admit that our environment, being much larger than our system (just one qubit), is much 'less' macroscopic than that of statistical mechanics origin. Note that large random matrices have been widely used to model a variety of complex quantum systems, including heavy nuclei and atoms, meso-and nanoscopic particles, quantum networks, etc, i.e., not necessarily macroscopic quantum systems. Thus, our model can describe a qubit in a mesoscopic (or even nanoscopic) cluster, a quantum dot, a quantum network, etc [16,17,19,20,24,25]. In other words, we propose the model in which the bosonic (quantum Gaussian) noise that models macroscopically large and many-body quantum environments is replaced by the random matrix noise that models smaller quantum environments of mesoscopic or even nanoscopic size.
Returning to (2.26), we conclude that it allows us to confine ourselves to the mean value r Denoting the limiting mean value r ( ) t , we have for the one-qubit version of (2.8) for the model (2.23)-(2.25) [21,22]: , where ab ( ) A t and ab ( ) B t are given by contour integrals of analytic functions solving a system of two selfconsistent equations, see formulas (4.1)-(4.7) of [21]. This is reminiscent to the mean field theories of statistical mechanics where the selfconsistent equations are quite common (recall the Curie-Weiss and van der Waals equations). In fact, it is widely believed that random matrices of large size provide a kind of mean field models for one body disordered quantum systems.
The coefficients ab ( ) A t and ab ( ) B t depend on numerical parameters s v , , see (2.25), the parameter E, the limiting value of the eigenvalue of M N corresponding to its eigenstate in the one qubit analog of  P of (2.6), and the functional parameter n 0 (see (2.23)). This allows for a variety of asymptotic regimes of statics and dynamics of the model (see [21,22,35]).
One of the properties of the random matrix model (2.23)-(2.25), which can be attributed to its non-Markovianity, is that the infinite time limit r ¥ ( )of (2.27) is not an equilibrium measure, moreover, r ¥ ( ) depends on the initial conditions r ( ) 0 , see formula (4.11) of [21]. This can be viewed as a manifestation of the role of the 'size' of the environment in the 'forgetting' the initial conditions during the reduced matrix evolution. Indeed, according to the above, random matrices can model extended systems of intermediate size (between microscopic and macroscopic) which are not necessarily large enough to provide the complete loss of memory on the initial conditions. Note that this property is unusual (and undesirable) for equilibrium statistical mechanics, since it implies that the long time evolution does not equilibrates the system. On the other hand, the property could be of interest for quantum information theory, since systems, which remember 'forever' information about their initial conditions, can be a resource for storing the quantum information.
It follows from (2.27) that the evolution of diagonal entries of the limiting reduced density matrix (controlling the amplitude damping and dissipation) and that of its off-diagonal entry (controlling the dephasing and decoherence) are dynamically independent, i.e., r r r a r r r a

a a a a a a a a a
-- In view of the complex form of ab ( ) A t and ab ( ) B t (see [21]) the asymptotic (Bogolyubov-van Hove) regime was considered in [21,22] yielding, upon passing to the interaction representation in which r a a , are the same but r a a -, get the factor a e ts 2i : and n 0 is the limiting density of states (2.23) of the reservoir. Note that the dynamics of the diagonal and off-diagonal entries of the reduced density matrix are again independent, (2.28), however, it cannot be purely dephasing, i.e., to have zero rate of decay of diagonal entries and non-zero rate of off-diagonal entries (see, e.g., [36] for pure dephasing evolution in a Markovian model). Note also that one obtain a simple model with pure dephasing dynamics by replacing s x by s z in (2.25) or in (2.19) (see [37] for a rather general version of analogous static model).
It follows from (2.28) and (2. Thus the dependence of r ¥ ( )on r ( ) 0 is generic in the asymptotic regime (2.29). It is also of interest that formulas (2.28) and (2.30)-(2.32), as well as more general formulas (2.27), contain the both the energy E of the reservoir and the energy (or, rather, level spacing) s 2 of the qubit. This is one more manifestation of the above mentioned fact, according to which our random matrix environment can be viewed as an one-body system, similar to those of condensed matter theory, where the energies of various constituents can be of the same order of magnitude. If, however, we assume that  E s or that n 0 of (2.23) is 'flat' enough in the interval - 2 , hence n n a  0 , then we obtain from (2.28) and (2.30)-(2.32) r t r r pn a , 0 The first term on the right is the infinite temperature Gibbs (hence equilibrium) distribution of a two-level system (see (2.34)). Moreover, the whole rhs solves the equation In the studies of open quantum systems the memory effects are often ignored, and this leads to dynamical semigroups (2.21) and the memoryless Markovian dynamics [2,7,9] governed by the Born-Markov-Lindblad master equation. However, actual dynamics of quantum systems can deviate from the Markovian dynamics. This motivated a considerable activity in the study of various aspects of non-Markovian dynamics in the last decades [9][10][11][12].
Following the activity, we will discuss briefly the non-Markovianity of the dynamics defined by (    The idea behind the definition (as well as several other, see, e.g [11]) is that if the dynamics is Markovian, then the distance (2.42) between the density matrices r ( ) t 1 and r ( ) t 2 , resulting from a pair of distinct initial conditions via the dynamics in question, is monotone decreasing in time ( t ( ) d is negative) signalling the one way flow of information from the system to the environment, hence the memoryless dynamics. Thus, the positivity of (2.41) can be viewed as a manifestation of the non-Markovianity of the dynamics in question, hence memory effects resulting from the backflow of the information from the environment to the system. It can be shown that for the dynamics (2.29)-(2.32)  n n n n n n n n n n n F = - We will also mention an interesting class of unital superoperators (quantum channels) which are 'quantum' analogs of doubly stochastic matrices and are defined by the property F = ( ) 1 1 [5]. It is easy to check that the channel defined by

Two-qubit random matrix Hamiltonians
We will present now the versions of general two-qubit Hamiltonians (2.1)-(2.5) in which the environment Hamiltonians and those of the qubits-environment interaction are the random matrices (2.23), (2.24) that are used in the one-qubit Hamiltonian (2.25). We will use the same sub-symbols I, F and C as in (2.1)-(2.5) for the corresponding Hamiltonians. The Hamiltonian of two qubits will always be given by ( are two independent random matrices with the same probability distribution as W N in (2.24). We obtain the Hamiltonian of two qubits, each of them interacts via independent random matrices with its own random matrix environment: One can also use distinct matrices  = M a A B , , N a provided that their limiting density of states of (2.23) are the same (otherwise the final formulas are too complicated).
(ii) Hamiltonian H F . Here, in full correspondence with general formulas (2.4) and the case (i) above we set: is the Hamiltonian of two identical qubits where only one of them interacts via a random matrix with a random matrix environment.
(iii) Hamiltonian H C . Here, in view of (2.5) and the case (i) above, we choose M N as the the common environment Hamiltonian and  to obtain the Hamiltonian of two identical qubits interacting via a random matrix with a common random matrix environment. Note that random matrices, in particular, Hamiltonians similar to our Hamiltonians (2.3)-(2.5), have been used in the description of the qubits dynamics, although in a more generic form which does not possesses the Pauli matrix structure, see [24][25][26] and references therein. The Pauli matrix structure of the one-qubit random matrix Hamiltonian (2.25) has been considered in [21] and later in [22,23,35].
First of all, [23][24][25][26] deal with a regime where the evolution time varies on the so-called Heisenberg time scale t H of the environment. The scale is of the order of magnitude of the inverse level spacing of the environment and is a quantum analog of Poincaré time scale where the classical statistical mechanics description fails. In the case of random matrix models t~N H , where N is the size of the corresponding matrix. This time scale is widely used in the studies of quantum chaos and related delicate phenomena (Loschmidt echo, correlation hole, etc), see, e.g. [18] for a review. On the other hand, [21,22,35] exploit a 'complementary' regime of much shorter than the Heisenberg (Poincaré) time scale common in non-equilibrium statistical mechanics (known there as kinetics scale) and in the theory of open systems. We believe that this time scale is of considerable importance in quantum information theory. Technically the difference is that [23][24][25] , , while we use the limit  ¥ N assuming that the evolution time is N-independent, or, at most, that t  N , see (2.29).
The two different time scales mentioned above and the corresponding energy scales~-E N 1 and E is N-independent are well known in the random matrix theory as the global (or macroscopic) scale (or regime) and the local (or microscopic) scales (or regime), see e.g. [16,20]. The local regime requires a quite involved and sophisticated techniques and less feasible by now. This is why [23][24][25] confine themselves to the first order approximation in the qubit-environment coupling for the reduced density matrix. The work [26] extends the result of [23][24][25] on an essentially larger interval of time and the system-environment interaction by using the Grassmann integration techniques but still remaining in the frameworks of the local regime, i.e., within the Heisenberg time scale t~N H . On the other hand, the global regime is fairly well understood and allow us to find closed expression for the reduced density matrix, applicable within the kinetics time scale. We believe that the fact that our models reproduce several recent results on the evolution of entanglement, e.g. its disappearances and revivals (see figures 1-3) found before for bosonic environments, see e.g. reviews [12,15], provides a serious justification of our approach. Moreover, it yield several new facts on interrelations between the negativity and concurrence (see e.g. figures [3][4][5]. Another difference between our paper and [23][24][25][26] is that the latter deal mainly with the quantum coherence (in particular, entanglement) characteristics known as the purity and the fidelity and only partly with the concurrence, while we carry our rather detailed comparative and asymptotically exact in N analysis of three entanglement quantifiers: the concurrence, the negativity and the quantum discord (see section 2.5 for their definition).

Entanglement quantifiers
Entanglement, having a short but highly non-trivial mathematical definition (a state of a bipartite system (1.2) is entangled if it is not a tensor product of the states of its parties  A and  B ), is a quite delicate and complex quantum property admitting a wide variety of physical manifestations and potential applications. This motivated the introduction and the active study of a number of its quantitative characteristics (entanglement quantifiers or entanglement witnesses) which are functionals of the corresponding state and determine its 'amount' of entanglement, see reviews [4,8,14,28,30,40]. We will consider in this paper three widely used entanglement quantifiers of bipartite states: the negativity, the concurrence and the quantum discord. Since    there is a number of reviews and a considerable amount of original works treating these characteristics, we will give here only their brief description.
(i) Negativity r [ ] N , see reviews [4,8,28]. If a two-qubit state ρ acted upon with the partial transposition with respect to one of the two subsystems has negative eigenvalues, then it is entangled and vice versa. In other words, the reservoir is not able to create entanglement of two qubits if and only if the partial transposition preserves its positivity.
The corresponding numerical characteristic is known as negativity and defined as the sum of absolute values of negative eigenvalues of the partially transposed state r PT . If ρ has the X-form (2.17), then r PT has the same form and and we have in this case r r r r r r r r r r r = -- (ii) Concurrence r [ ] C , see reviews [4,8,12,28]. It is quite popular entanglement quantifier of two-qubit states, closely related to another entanglement quantifier, known as the entanglement of formation and applicable in general to multiqubit systems. We have r l  The concurrence varies from 0 for separable states to 1 for maximally entangled states and its positivity is a necessary and sufficient condition for the state to be entangled.
(iii) Quantum discord r [ ] D , see reviews [14,30]. Quantum discord has a rather non-trivial definition. The idea behind it is to view 'genuine' quantum correlations as those that are not classical. Accordingly, the quantum discord is defined as the difference between the quantum generalizations of two classically equivalent definitions of the mutual information of the classical information theory. Quantum discord is non-negative in general and is positive for the entangled states. However, there exist not entangled states having a positive discord, hence not classical. In other words, the quantum discord 'feels' a subtle difference between product states and classical states and can be viewed as a measure of total non-classical (quantum) correlations of bipartite systems and, more generally, pair correlations in quantum many-body systems including those that are not captured by the concurrence and the negativity (2 qubits) and the entanglement of formation (many qubits). In general quantum discord is more robust than entanglement against decoherence, and it does not present sudden death phenomena (see below).
Unfortunately, we are not aware of a compact formula for the quantum discord of an arbitrary X-state (2.17) similar to (2.50) and (2.51) for the negativity and concurrence. However, for the states arising in our models we found an ansatz that simplifies considerably the numerical analysis:

Time evolution of entanglement quantifiers
We will present now our results on the time evolution of the negativity, the concurrence and the quantum discord for the random matrix models = = Nn N I F C n , , , , 1, 2, 3 given by initial conditions (2.12)-(2.14) and Hamiltonians (2.47) for two identical qubits each interacting with its own environment, (2.48) for two identical qubits but only one of them interacts with an environment and (2.49) for two identical qubits both interacting with the same environment.
Let us mention first that for all these models we have an analog of bound (2.26) providing the selfaveraging property (typicality) of the reduced density matrices in question for large N. The corresponding proofs are analogous to that of (2.26) for the one-qubit case (see [22]), but are rather tedious and we do not present them here. The property guarantees the validity of the description of the model in question by using only the mean reduced density matrix (the representativity of means). Recall that the representativity of means is a basic general property of statistical mechanics (where the expectations are with respect to the Gibbs measure in the equilibrium case and with respect to the corresponding non-equilibrium distribution in general) and of the theory of disordered systems where one has one more expectation with respect to the frozen macroscopic disorder.
Let us also mention that one can consider models with classical environment. A seemingly simple model could be that with a Langevin-type stochastic dynamics of the density matrix, i.e., say, the N=1 version of the one-qubit Hamiltonian (2.25) (and the corresponding two-qubit Hamiltonians (2.47)-(2.49)) whose parameters are classical noises (Ornstein-Uhlenbeck random process, random telegraph process, etc) describing the environment, see, e.g. [12,13,41] and references therein. An analog of the total density matrix (2.6) in this approach is the random 4×4 positive definite matrix solving a Langevin-type stochastic differential equation which play the role of the quantum (Liouville-Neumann) equation. The role of the reduced density matrix (2.7) plays the expectation of this random solution for which one often obtains a governing (master) equation. However, in this approach the variance of the solution is of the same order of magnitude as its expectation (recall the standard Langevin equation for the Ornstein-Uhlenbeck random process). As a result, in this case there is no a counterpart of bound (2.26) guaranteeing the representativity of means, thus the complete description of the solution (the random density matrix) is given by its whole probability distribution.
By using (2.11), the reduced density matrices for models In and Fn, = n 1, 2, 3 can be easily obtained for all three families of initial conditions (2.12)-(2.14) from the one-qubit density matrix given by (2.28) and (2.30)-(2.32). For example, we have for the model I1 The formulas for all other models = In Fn n , , 123 are similar. As for the models = Cn n , 1, 2, 3 of two identical qubits interacting with a common environment, their reduced density matrix cannot be obtained as the tensor product those for the one-qubit case. It turns out, however, that in this case the reduced density matrix can be obtained by generalizing the techniques of [21,22] yielding (2.26), (2.27). Corresponding calculations and results are rather involved and diverse and will be published elsewhere [27]. Here we will only give several interesting, we believe, examples of the corresponding dynamics for the model C2 in the asymptotic regime (2.29) and for the maximally entangled case a = 1 2 2 (the Bell case) of the initial conditions (2.13) (see figures 4 and 5).
We will present now the result of numerical analysis of our random matrix models. It will be convenient to use the energy units where the qubit amplitude s of (2.46) is set to 1.
We will start with the measure (2.41), (2.42) of the non-Markovianity of our reduced dynamics of two qubits. Unfortunately, we were not able to find an analog of the one-qubit formula (2.44) for two-qubit, i.e., to carry out explicitly the maximization procedure over all pairs of initial conditions (4 × 4 positive definite matrices of trace 1). Since, however, we deal only with three families (2.12)-(2.14) of initial conditions, it seems reasonable to perform the maximization only over the initial states n while considering the model Nn for a given = n 1, 2, 3. This gives a lower bound for the maximum (2.41) over all the states, thus allows us to conclude that the reduced two-qubit dynamics for the model Nn is not Markovian if the bound is positive. In particular, it can be shown that for the model F1, where the initial conditions (2.12) are parameterized by the parameter a Î [ ] 0, 1 1 , the corresponding restricted maximum is attained for the pairs of the parameter, i.e., for the unentangled cases of (2.12). Figure 1(a) gives the plot of the distance (2.42) for the model F1 (one qubit is free and one interacts with its own environment) and the maximizing pairs (3.2) of initial conditions. According to the figure, the corresponding integral in (2.41) is 0.092, thus the reduced two-qubit dynamics with parameters indicated in the figure is non-Markovian by the criterion (2.43). Note that we consider the case where n n n = < + -0.8 1 0 2 , i.e., the one-qubit dynamics is Markovian (the distance (2.42) is monotone decreasing and the rhs of (2.44) is zero). On the other hand, in the two-qubit model F1 the distance may be both monotone and not monotone in the both cases:  n n n + -0 2 and n n n > + -0 2 for different a a ¢  ( ) , 1 1 . This can be viewed as a manifestation of the 'increase' of the amount of non-Markovianity of the two-qubit reduced dynamics comparing with that of the one-qubit reduced dynamics with the same one-qubit Hamiltonian.
Next, we will discuss the negativity and the concurrence. We will consider first the models In and correspond to the Markovian and the non-Markovian one-qubit dynamics according to (2.44 )).
It is important that all the plots become zero at a finite moment t ESD and remain zero afterwards, i.e., an initially entangled state becomes disentangled starting from t t = . ESD This is the ESD phenomenon [42,43], see also [12,15] for reviews. It is of interest that according to figure 1(b) in the both models F1 and I1 the ESD moments for the concurrence and negativity coincide. Moreover, it can be shown analytically that In In while we have for any ρ of X form Hence, in spite of rather different definitions (2.50) and (2.51) of these entanglement quantifiers, they can be closely related in certain cases, see, e.g. [28] for a similar fact. This is why we give the plot of N 2 I1 and N 2 F1 on figure 1(b), where the above equality is not valid. Figure 2 shows that t ESD depends rather non-trivially on the choice of the family of the initial conditions and on of the parameter of the family. It is also of interest that the plot of concurrence for a = 0.7 2 (solid line) of figure 2(b) is quite similar to a plot obtained experimentally, see [8].
We have also proved analytically the relation Fn n n n 'separating' the dependence of the concurrence on the initial conditions ( a ( ) C n n is the concurrence at t = 0) and on the Hamiltonian ( f is a smooth function such that = ( ) f a b 0, , 1 and a and b are given in (2.44) and (2.31)). The relation implies that t ESD is independent of inial conditions (parameter a = n , 1,2 n in (2.12) and (2.13)).
Passing to the models = Cn n , 1, 2, 3 determined by Hamiltonian H C of (2.49) and describing the time evolution of two qubits embedded in a common environment, we find that their reduced dynamics are more diverse than those of models In and Fn. Figure 3 depicts the concurrence and the negativity in the model C2 with maximally entangled initial condition (2.13) and with the Lorentzian density of states (2.23) of the environment Note that a more general function (2.22) is widely used as a model spectral density (2.20) for special versions of the spin-boson dynamics of two qubits, see, e.g. [10,12] and references therein. It follows from figure 3 that for (3.6) and the parameters indicated in the figure caption the plots of the concurrence and the negativity have two disjoint branches.
The branch of figure 3(a), from t = 0 to t t = ESD , is quite similar to those of figures 1(b) and 2 and typical for the Markovian evolution (monotone decreasing curves). However, in the model C2 there is the second branch of these quantifiers, although with the amplitude which is at most 18% of the maximum amplitude of the first branch. The 'revival' of the concurrence and the negativity, hence the entanglement, after the ESD is known as the ESB and was found in special versions of the spin-boson dynamics, where the ESD and ESB can even appear several times, although often with decaying amplitudes, see, e.g. [12]. On the other hand, in our case the second branches of the concurrence and the negativity tend to a non-zero value at infinity, thereby exhibiting the entanglement trapping phenomenon [12]. It is also of interest that for the first branches we have an analog of of (3.3) while for the second branches the strict analog of (3.4) is valid except the endpoints of the branch, where the both plots equal zero, see figure 1(b). We have also found a case of the model C2, where the second branch becomes zero at certain finite t t ¢ > ESD ESB and is zero afterwars, i.e., one more ESD. Note that since even partial entanglement of large number of copies of a bipartite quantum system may be distilled to a smaller number of maximally entangled systems [4,5,8], many quantum information processing and communication tasks may be carried out as long as there is still some entanglement left. This is why the result of figure 3(b) seems to be of interest.
The time evolution of the quantum discord is presented on figure 4. (non-Markovian one-qubit regime according to (2.44)) the plot (solid line) attains its minimum which is about ´-1.5 10 2 of its initial value, and then grows slowly to its infinite time value which is about´-2.7 10 2 of its initial value.
And only for n n n = + _ 1 0 2 the plot (dotted-dashed line) decays monotonically to zero at infinity. Recall that this is the 'most genuine' Markovian regime of our model, since, in this case, we have (2.43) (but not yet (2.35)-(2.37)) and the infinite time reduced density matrix (2.30)-(2.32) is independent of initial conditions (see (2.33)) although is not necessarily of the Gibbs form (2.34).
The time evolution of quantum discord in the model C2 corresponding to the maximally entangled Bell state Y a = | 2 1 2 2 of (2.13), is presented on figure 4(b). The plot is more structured. In particular, according to the solid lines on figures 4(a) and (b) the quantum discord grows monotonically for large times. This seems somewhat counterintuitive since one would expect that quantum correlations decrease with time. It is also of interest that the plot on figure 4(b) has a plateau which could be interpreted as the 'freezing' of the discord for a finite time interval found in several non-Markovian models with bosonic environment and in certain experiments [12,30].
We conclude that the quantum discord can be preserved and even amplified on certain time intervals of its evolution in the dissipative environment, see [12,14,30] for discussions and references on the properties of quantum discord in other environments.
Note also that the vanishing of the quantum discord, i.e., an analog of the ESD phenomenon of figures 1(a)-3 would have in general a number of important implications, in particular, providing an efficient indicator of the classicality of state in question, see, e.g. [4], especially sections IV.C and VIII.A1. We have not found an analog of the ESD for quantum discord in our models but only either its asymptotic vanishing at infinity or identical vanishing for non-entangled initial conditions (see (2.15) and the solid lines on figure 5), thereby confirming the 'larger' robustness of the quantum discord with respect to the decoherent influence of the environment. This is well seen from figure 5 depicting the value of quantum discord at infinity for models F2 and ( ) F3 2 as functions of the parameters of the models.

Conclusion
In this paper we have considered several schemes of the time evolution of two qubits embedded into an environment. The schemes are determined by certain entangled initial conditions and one of three Hamiltonians describing: (a) each qubit interacts with its own environment, (b) one of qubits interacts with its own environment and one is free and (c) both qubits interact with a common environment.
A concrete realization of these quite natural schemes is determined by the choice of the Hamiltonians of the environment and of its interaction with qubits. An often used choice corresponds to the so-called spin-boson Hamiltonian describing a macroscopic and many body environment, see e.g. [12] for a review. Another possibility is to use large random matrices to model a complex disordered and not necessary macroscopic environment. One has to choose then the time (the energy) scale of random matrix theory. We choose the socalled global regime which provides the description analogous to that of non-equilibrium statistical mechanics and the theory of open systems (see, however, works [23][24][25][26] for another choice known as the local regime, pertinent to the quantum chaos studies and leading to complementary results).
Obtained models are asymptotically exactly solvable in the limit of large matrix size and allow us to analyze a variety of the qubit dynamics interpolating between the Markovian (memoryless) and non-Markovian (including the environment backaction) dynamics.
We have examined three widely used numerical characteristics (quantifiers) of quantum correlations, entanglement in particular: the negativity, the concurrence and the quantum discord and we found a number of patterns of time dependence of these quantifiers. A part of patterns is qualitatively similar to those obtained previously in the frameworks of bosonic models (the vanishing and the revival of entanglement, its freezing and locking) while another part consists of new patterns (connections between the negativity and the concurrence, two seemingly different entanglement quantifiers, the monotone growth of the discord for large times, relevant parameters determining a pattern). We believe that these findings manifest the universality (the independence of the model) of certain fundamental properties of entanglement as well as indicate new directions of its studies.