{'fraction_non_alphanumeric': 0.03719023945960305, 'fraction_numerical': 0.011991954764525066, 'mean_word_length': 4.471760797342193, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entanglement formation under random interactions

The temporal evolution of the entanglement between two qubits evolving by random interactions is studied analytically and numerically. Two different types of randomness are investigated. Firstly we analyze an ensemble of systems with randomly chosen but time-independent interaction Hamiltonians. Secondly we consider the case of a temporally fluctuating Hamiltonian, where the unitary evolution can be understood as a random walk on the SU (4) group manifold. As a by-product we compute the metric tensor and its inverse as well as the Laplace-Beltrami for SU (4).


Introduction
If two initially separable quantum systems are exposed to random interactions they are expected to become entangled, exhibiting random quantum correlations. How do these quantum correlations arise as a function of time? To address this question we study the entanglement between two qubits subjected to random interactions as a function of time. The study of entanglement dynamics under random environments has attracted much interest recently, for instance, the emerging entanglement between coupled quantum systems through a bosonic heat bath [1]. Although our system is oversimplified in comparison with these dissipative systems, we believe that our study may give the upper bound for the entanglement under the strong random interactions.
In what follows we assume that the two-qubit system is initially prepared in a well-defined pure state. As examples we consider two different initial states, namely, a non-entangled pure state (1) and in a fully entangled Bell state of the form where {|00 , |01 , |10 , |11 } denotes the canonical qubit configuration basis. Starting with the given initial state the system then evolves unitarily as where the time evolution operator is determined by U (0) = 1 and with a randomly chosen interaction Hamiltonian.
Throughout this paper we consider two different types of randomness, namely (a) Quenched randomness, where H(t) = H is time-independent. In this case one considers an ensemble of two-qubit systems starting from the same initial state, where each member evolves by a different but temporally constant Hamiltonian drawn from an SU (4)-invariant distribution.
(b) Temporal randomness, where the dynamical evolution is generated by a timedependent Hamiltonian H(t) which fluctuates randomly as a function of time [2]. On a single system the resulting temporal evolution of the state vector can be understood as a unitary random walk in C 4 .
The difference between the two cases is visualized in Fig. 1. In this figure the big red sphere stands symbolically for the 15-dimensional group manifold of SU (4). Each of point on the sphere represents a certain unitary transformation U (t) acting on the 4dimensional Hilbert space of the two-qubit system. Starting with U (0) = 1, which may be located e.g. at the 'north pole' of the sphere, the temporal evolution U (t) can be represented by a certain trajectory (blue line) on the group manifold.
Let us now think of an ensemble of such systems, represented by a set of statistically independent trajectories. In the quenched case (a), where a random Hamiltonian is chosen at t = 0 and kept constant during the temporal evolution, these trajectories are straight, advancing at different pace and pointing in different directions, while in case (b) they may be thought of as random walks on the group manifold. At a given final observation time the trajectories of the ensemble terminate in different points marked by small blue bullets in the figure, each of them representing a unitary transformation. Applying this transformation to a pure initial state one obtains a final pure state with a certain individual entanglement. In the sequel we are interested in the statistical distribution of these final states and their entanglement.
To quantify the entanglement we use two different entanglement measures. For a pure state the entanglement is defined as the von-Neumann entropy of the reduced states where ρ 1,2 (t) = Tr 2,1 ρ(t) denotes the time-dependent reduced density matrix of the respective qubit. In cases where the logarithm is too difficult to evaluate we resort to the so-called linear entropy as an alternative entanglement measure. Note that both measures can be obtained from the more general Tsallis entanglement entropy [3] E q (t) = 1 − Tr[ρ q 1 (t)] q − 1 (7) in the limit q → 1 + and q → 2, respectively.
Furtheremore, the Rényi entanglement entropy [4] H q (t) = log Tr[ρ q is of interest. Also this entropy measure generalizes the von-Neuann entropy when setting q → 1.
Our main results are the following. In the first case (a), where a temporally constant Hamiltonian is randomly chosen, the mean entanglement is expected saturate at a certain value for t → ∞. Our findings confirm this expectation, but surprisingly we observe that the average entanglement overshoots, i.e., it first increases, then reaches a local maximum, then slightly decreases again before it finally saturates at some stationary value.
In the second case (b), where the Hamiltonian changes randomly as a function of time, the average entanglement saturates as well, although generally at a different level. This saturation level, which is the average entanglement of a unitarily invariant distribution of 2-qubit states, has been computed previously in Refs. [5][6][7][8]. Here we investigate the actual temporal behavior of the entanglement before it reaches this plateau. As a by-product, we compute the metric tensor and its inverse on the SU (4) group manifold as well as the corresponding Laplace-Beltrami operator (see supplemental material), which to our knowledge have not been published before.
Remarkably, for ρ(0) = |11 11| this density matrix depends only on six angles α 1 , . . . , α 6 out of 15. Since the observables investigated in this paper are invariant under local unitary transformations, any pure separable initial state will give the same result. Therefore, without loss of generality we can restrict ourselves to ρ(0) = |11 11|, taking advantage of the low number of parameters in this particular case.

Computing averages on the SU (4)-manifold
In the following section we will consider an ensemble of trajectories of unitary transformations generated by random interactions. Using the above representation, each trajectory can be parametrized by a time-dependent vector of Euler angles α(t).
A statistical ensemble of trajectories is therefore characterized by a probability density p(α, t) to find a unitary transformation with the Euler angles α at a given time t.
For a given probability density p(α, t) one can compute the ensemble average of any function f (α) (such as the density matrix ρ(α) or the entanglement E(α)) by integrating over the complete parameter space of the SU (4) manifold weighted by p(α, t): Here V SU (4) is the integrated group volume which serves as a normalization factor while denotes the volume element on the SU (4) group manifold. The actual integration measure is defined by the function µ(α) which depends on the chosen representation.
The transformed state ρ(α) = U α ρ(0) U † α is still pure but generally entangled. Being interested in the average entanglement of states generated by random unitary transformations, it makes a difference whether the entanglement is computed before taking average over α or vice versa, as will be discussed in the following.

Average of the entanglement
Let us first discuss the case of computing the entanglement before taking the average over all α. In this case one has to compute the reduced density matrix of the first qubit for given ρ(α), defined as the partial trace For the initial state ρ(0) = |11 11| this is a 2 × 2-matrix with the elements In general the reduced density matrix σ(α) is no longer pure and its von-Neumann entropy quantifies the entanglement E(α) between the two qubits. In order to compute the entropy we determine the eigenvalues of σ which are given by +32 sin 2 (α 2 ) sin 4 (α 6 ) cos (4α 4 ) + 6 cos (4α 6 ) + 50 Having determined these eigenvalues, the entanglement of ρ(α) is given by Finally, the entanglement has to be averaged over all trajectories (see Eq. (52)), i.e.
However, if the average of the von-Neumann entropy is too difficult to compute, we will also use the linear entropy as an alternative entanglement measure.

Entanglement of the average
Alternatively, we may first compute the average density matrix ρ = ρ(α) α and then determine the entanglement of the resulting mixed state. To this end a suitable entanglement measure is needed. An interesting quantity in this context is Wootters concurrence [12] defined by where λ 1 , . . . , λ 4 are the decreasingly sorted square roots of the eigenvalues of the matrix In this expression σ y is the Pauli matrix while ρ * denotes the complex conjugate of σ without taking the transpose.
From the concurrence one can easily compute the entanglement of formation of the mixed state, which is given by

Quenched random interactions
In the case (a) of quenched randomness each element of the ensemble is associated with a time-independent random Hamiltonian H. Since the spectral decomposition of a randomly chosen Hamiltonian is always non-degenerate, the time evolution operator can be written as Hence the state of the system evolves as where ρ(0) denotes the initial state.
The Hamiltonian itself has to be drawn from a certain probabilistic ensemble of Hermitian random matrices [2,13]. Here the most natural choice is again the Gaussian unitary ensemble (GUE). This ensemble has the nice property that the probability distributions for the eigenvalues E j and the eigenvectors |φ j factorize and thus can be treated independently. More specifically, the eigenvalues are known to be distributed as where A = 1 2σ 2 is a constant determining the width of the energy fluctuations and therewith the time scale of the temporal evolution. In the following the corresponding average over the energies will be denoted by . . . E . On the other hand, the orthonormal set of eigenvectors is randomly oriented in the four-dimensional Hilbert space according to Haar measure, independent of the eigenvalues. If one defines the qubit basis this average can be carried out by setting and integrating over all Euler angles α according to the Haar measure (see Appendix B). This average will be denoted by . . . α . The total GUE average thus factorizes as

Entanglement of the averaged density matrix
Let us now compute the average density matrix First we compute the average over the energies is the normalization factor. This leads us to the result where we defined the scaled time and the function Thus Eq. (33) reduces to What remains is to determine the operators Obviously, the first sum in Eq. (38) is given by As for the second sum in Eq. (38), we note that the distribution of eigenvectors |φ j = U † α |j is invariant under a permutation of the basis vectors |j , hence the four operators T jj coincide. Moreover, under a unitary transformation V ∈ SU (4) they transform as where we have used the invariance of the GUE-eigenvectors under the replacement U α → U α V . This means that T jj is invariant under V if and only if V commutes with the initial state ρ(0). For a pure initial state this implies that T jj has to be a linear combination of the identity and the initial state itself, i.e. T jj = a1 + bρ(0). The linear coefficients a and b can be determined as follows. On the one hand, the identity implies that Tr 4 j=1 T jj = 1, hence 4a + b = 1/4. On the other hand, we note that is invariant under unitary transformations of ρ(0) and independent of j, hence we may choose j = 4 and ρ(0) = |4 4| to obtain giving a = b = 1 20 . Therefore, we arrive at the convex combination of 1/4 and ρ(0) with f (τ ) given in Eq. (37) and τ = t/ √ 2A, which holds for any pure initial state ρ(0). As expected, the averaged state lies on the segment between the initial state ρ(0) and the maximally mixed state 1/4 due to the symmetries of the Haar measure of SU (4) .
Having computed the mixed state of the ensembleρ(t) we can now compute the corresponding entanglement of formation as a function of time. For a non-entangled  initial pure state ρ(0) = |11 11| we find that E F (ρ(t)) = 0 for all times. However, if we start from the Bell state (2) with the initial entanglement E F (ρ(0)) = ln 2 ≈ 0.693 we find numerically that the entanglement first decreases and vanishes at the finite scaled time τ ≈ 0.4997. (see Fig. 2).

Average of the individual entanglement
Instead of computing the entanglement of the average density matrix, let us now compute the average of the individual entanglement of each trajectory, i.e. the entanglement is computed before taking the GUE average. Although the von-Neumann entanglement entropy of the individual pure states would be straight-forward to compute (see (20)), we did not succeed to compute the average. For this reason let us consider the GUE average of the linear entropy where ρ 1 (t) denotes the reduced density matrix of the first qubit. In the qubit basis (30) this can be rewritten as Inserting (28) and exploiting again that the GUE average factorizes, we get where c µβ j = φ j |µβ . For the initially non-entangled state ρ(0) = |11 11| this expression reduces further to This function is plotted in Fig. 3. As one can see, the linear entanglement entropy (black line) first increases rapidly, then reaches a local maximum L 0.199936 < 0.2 at τ max 0.817377, then decreases again and finally saturates at the value Because it would need much more effort to calculate the analytical the linear entropy for different initial states analytically, we used numerical methods. The results are compared in Fig. 3. As one can see clearly, all the lines tend to touch the limit value at the fixed time τ 0.817 and the curves do not intersect. Fig. 1 explains the meaning of this result: Each single trajectory of the ensemble on the surface is deterministic with a given initial point, direction and velocity. Since all members of the ensemble share the same initial starting point on the upper half of the sphere, the probability for finding the walkers can be slightly higher on the upper half in the long-time limit because all trajectories will periodically return to this point. This is why the limit depends on the initial state and therefore deviates from the Haar measure. The bump can be seen as a transient state in which the probability distribution seems to be almost randomly distributed before saturating. Note that in contrast to the case discussed before (see Fig. 2) the system remembers its initial state, saturating at different levels of entanglement in the limit t → ∞.

Time-dependent random interactions
Let us now consider the case (b) of a temporally varying Hamiltonian, where the state vectors of the ensemble perform a unitary random walk on the SU (4) manifold. In this case the quantity of interest is the probability distribution p(α, t) to find the time evolution operator U (t) with the Euler angles α at time t. This probability distribution allows one to compute the ensemble average of any function f (α) (such as the density matrix ρ(α) or the entanglement E α ) by integration over the complete SU (4) volume weighted by p(α, t), i.e. we have to compute the integral over the ranges specified in (A.4), where dV SU (4) denotes the volume element according to the Haar measure defined in Eq. (13).

Expected average entanglement of a uniform distribution
Before studying the temporal evolution in detail, let us consider the limit t → ∞, where we expect the state vectors to be uniformly distributed on the group manifold.
Since such an ensemble is by itself invariant under unitary transformations, the state vectors are distributed according to a Gaussian Unitary Ensemble (GUE). Starting from this observation, Page [5] conjectured a closed expression for the expected average entanglement of an arbitrary bipartite quantum system with Hilbert space dimensions m and n, which was later proven rigorously by Foong, Kanno, and Sen [6,7]. This formula describes the average entanglement of a random pure state between two subsystems with Hilbert space dimensions m and n: Applying this formula to a random two-qubit system, one obtains This is the average entanglement between two qubits in a randomly chosen pure state.

Heat conduction equation on the SU (4) manifold
In order to compute the probability distribution p(α, t) one has to solve the heat conduction equation on the curved SU (4) group manifold. Here D denotes the diffusion constant while ∆ α is the so-called Laplace-Beltrami operator which generalizes the ordinary Laplacian on a curved space. As stated by [17] this Laplace-Beltrami operator is the Markov generator of the unitary brownian motion (see Appendix F). On a Riemannian manifold with the metric tensor g ij the Laplace-Beltrami-Operator ∆ is given by where g ij are the components of the inverse metric tensor and |g| = √ det g.
To our best knowledge the explicit expressions for g ij , g ij and ∆ α have not been published before. This is perhaps due to the fact that the formulas are so complex that even powerful computer algebra systems such as Mathematica R are not able to compute the inverse metric directly. Instead one has to invert the matrix manually element by element. Our explicit results are included in the supplemental material attached to this paper.

Early-time expansion
The solution p(α, t) of the heat conduction equation (55) and as well the averaged function f (t) can be expanded as a Taylor series around t = 0: Using (52) we can compare the coefficients of the both Taylor series. giving where dV SU (4) = |g|d 15 α denotes the volume element defined in Eq. (13). Using the heat equation (55) we can replace the partial derivative, obtaining The r.h.s. is an integral over derivatives of the probability density p(α, t) evaluated at t = 0. If all trajectories start at α 0 it is easy to see that this probability density at t = 0 is given by Inserting this expression into (60) the integral can be evaluated by partial integration, giving

Average of the density matrix
Using (10) we calculate the derivatives n α ρ(α). We find that Therefore we can express all higher derivatives of n ≥ 1 in terms of the first derivative § ∂ n ∂t n ρ(t) Hence, the solution of the averaged density matrix can be written as Using the non-entangled initial state |ψ 0 = |11 (which corresponds to taking α 0 = 0) this result specializes to As can be seen, this density matrix relaxes exponentially and becomes fully mixed in the limite t → ∞.

Average of the linear entropy
The same calculation can be carried out for the linear entropy defined in (22). Here we find that For this reason, the calculation is completely analogous, giving The result for an non-entangled initial state L(α 0 ) = 0 is therefore For a fully entangled initial state L(α 0 ) = 1/2 we get In both cases the averaged linear entropy tends towards 1 5 which coincides with the value of a randomly distributed density matrix according to Haar-measure (cf. Ref. [5,14]).

Average of the Tsallis entropy
Because of computational difficulties we calculated the Tasllis and von-Neumann entropy only calculated to first order in t. Moreover, we restricted ourselves to the case of a fully entangled initial state since the Taylor expansion fails for a non-entangled initial state due to a logarithmic factor in time at t = 0.  Using the definition of the Tsallis entropy (7) and following the same lines as outlined above, one ends up with the first order approximation In the limit q → 2 this reproduces the first order terms of (70) whereas for q → 1 we obtain the first-order approximation of the von-Neumann entropy By direct approximation on can derive the first order in τ = 2Dt of the von-Neumann entropy directly for a few timesteps using an initial state |ψ 0 = (cos φ, 0, 0, sin φ) with φ ∈ (0, π/4): For a fully entangled initial state (φ → π/4) this reproduces Eq. (72). Unfortunately, the expansion does not converge for φ → 0, which is why one has approximate an unentangled initial state directly, obtaining where γ is Eulers constant with numerical value of γ ≈ 0.577.
To clarify these results Fig. 4 shows numerical calulations and the analytical firstorder approximations as dashed lines for several initial states.
It is interesting to note that the von-Neumann entropy shows a bump if one starts near the limit value of 1/3, whereas the linear entropy does not have such a behaviour. But since one can interpret the linear entropy as lowest-order expansion term of the von-Neumann entropy around a pure state, it is not surprising that the von-Neumann entropy shows a more complex behavior. Moreover, a fixed initial state with a von-Neumann entropy of 1/3 is obviously not yet randomly distributed in SU (4).

Average of the Rényi entropy
We perform similar calculations as above for the Rényi entropy as defined in (8). This leads to the following approximation in first order for small times using a fully entangled initial state: For q → 1 we obtain the von-Neumann entropy as already computed in Eq. (72).

Discussion
In this paper we have studied how the entanglement of a two-qubit system subjected to random interactions evolves in time. We have considered two different types of randomness, namely, quenched and time-dependent random interactions, starting either from a separable or from a maximally entangled initial state. The main results are the following: • Quenched random interactions: Since entanglement measures are non-linear, it makes a difference whether the average over the quenched disorder is carried out before or after evaluating the entanglement. In the first case we can compute the In general, it would be interesting to see whether there exists an inequality relating the quantities E(ρ(t)) and E( ρ(t) ). averaged density matrix explicitly (see Eq. (45)), finding that the entanglement of formation decreases monotonously, see Fig.2. In the second case we can only compute the so-called linear entropy (50), which is found to overshoot before it saturates.
• Time-dependent random interactions: As outlined in the Introduction, this problem is equivalent to solving a random walk on the SU (4) group manifold. We find exact expressions for both the averaged density matrix and the averaged linear entropy, confirming that these quantities vary exponentially with time. For the averaged von-Neumann entropy as a special case of Tsallis entropy, however, we could only obtain a first-order approximation.
It should be noted that taking the average after evaluating the entanglement leads to results which are not directly accessible in experiments. The reason is that the computation of entanglement in pure-state systems requires the knowledge of the full density matrix in one of the subsystems. This matrix can only be measured by means of repeated experiments under identical conditions, which in our case would mean to use the same realization of randomness. Having estimated the entanglement this result has then to be averaged over different realizations of randomness. In experiments, where the randomness differs upon repetition, one would instead obtain a fully mixed density matrix without any information about entanglement.
As a by-product, when solving the random walk problem on the SU (4) manifold we had to compute the corresponding Laplace-Beltrami operator which in turn required to compute the metric tensor and its inverse. The computation of the metric was extremely difficult and to our best knowledge has not been done before. The resulting expressions are lengthy (see Supplemental Material). Surprisingly, we finally arrive at very simple results (see Sect. 4.4). This indicates there might be a deeper mathematical structure behind the problem that we failed to understand so that the brute-force calculation presented here is perhaps not really necessary. It would be interesting to investigate this point in more details. Table D1. Results for the integration of R ijkl (τ )

Appendix E. Numeric calculations
Depending on the case of the randomness we use different approaches for the numerical calculation.
• In quenched randomness case it is possible to compute ρ(t) directly for a single member of the ensemble using Eq. (28), because H is random but constant.
• In the temporal case we have to choose a new random Hamiltonian at each numerical time step dt. The unitary operator is chosen as U = A −1 · A † with A = 1 + iH dt 2 to ensure a unitary transformation. To ensure scale invariance in time, one has to scale a chosen H by 1/ √ dt as expected for noise terms.
Afterwards one can carry out the average over the ensemble to compute the averaged density matrix and the entanglement of formation, or to obtain the averaged entanglement measure (linear, von-Neumann).
The ensemble of the random Hamiltonians H is the GUE. To obtain a H ∈ GUE one has to chose a Hermitian matrix containing • random numbers z with z ∈ N (0, σ 2 ) on the diagonal entries and • complex random numbers z 1 +iz 2 √ 2 with z 1 , z 2 ∈ N (0, σ 2 ) on the off-diagonal entries, where N (0, σ 2 ) is the normal distribution. These are random matrices according to the GUE whose components are random numbers out of N (0, σ 2 dt). Introducing a normalized time t → τ σ 2 the random distribution changes to N (0, dτ ), which is why dH τ generates a Wiener process.
Having a look at the differential dU (τ ) we conclude that where we expanded the exponential function in the last step for dτ → 0 using (F.5). This is a stochastic differential equation for the Ito process {U (τ )} τ ∈R + and therefore describes a unitary Brownian motion on SU (4) (see also [18,19]).
As stated by [17] the Markov generator A of this diffusion process is given by the Laplace-Beltrami-Operator on the Riemannian manifold A = 1 2 . As described above, the metric is induced by the matrix scalar product. Since the Laplace-Beltrami-Operator is independent of the choice of the basis, we use an operator α that is parametrized by the Euler angles α, getting ∂p(α, τ ) ∂τ = 1 2 α p(α, τ ). (F.7) Using the rescaled time t = τ σ 2 and the definition D ≡ σ 2 2 we end up with the diffusion equation (55).