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 group manifold. As a by-product we compute the metric tensor and its inverse as well as the Laplace-Beltrami for .
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 . 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
and in a fully entangled Bell state of the form
where 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 and
with a randomly chosen interaction Hamiltonian.
Throughout this paper we consider two different types of randomness, namely
Quenched randomness, where 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 -invariant distribution.
Temporal randomness, where the dynamical evolution is generated by a time-dependent Hamiltonian which fluctuates randomly as a function of time . On a single system the resulting temporal evolution of the state vector can be understood as a unitary random walk in .
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 . Each of point on the sphere represents a certain unitary transformation acting on the 4-dimensional Hilbert space of the two-qubit system. Starting with , which may be located e.g. at the ‘north pole’ of the sphere, the temporal evolution 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 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 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 
in the limit and , respectively.
Furtheremore, the Rényi entanglement entropy 
is of interest. Also this entropy measure generalizes the von-Neuann entropy when setting .
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 . 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 group manifold as well as the corresponding Laplace-Beltrami operator (see supplemental material), which to our knowledge have not been published before.
2 Random unitary transformations in four dimensions
2.1 Representation of transformations
In what follows we use a particular representation of the group which was originally introduced by Tilma et alin Ref. . As reviewed in the Appendix, the group elements are generated by 15 Gell-Mann matrices , allowing one to parametrize unitary transformations by
where the 15 Euler-like angles vary in certain ranges specified in (80). Applying such a unitary transformation to the non-entangled initial state one obtains the density matrix
with the components
Remarkably, for this density matrix depends only on six angles 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 , taking advantage of the low number of parameters in this particular case.
2.2 Computing averages on the -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 . A statistical ensemble of trajectories is therefore characterized by a probability density to find a unitary transformation with the Euler angles at a given time .
For a given probability density one can compute the ensemble average of any function (such as the density matrix or the entanglement )) by integrating over the complete parameter space of the manifold weighted by :
Here is the integrated group volume which serves as a normalization factor while
denotes the volume element on the group manifold. The actual integration measure is defined by the function which depends on the chosen representation. Here we use the uniform measure, also known as Haar measure , which is by itself invariant under unitary transformations.111For example, in spherical coordinates the invariant measure of the rotational group would be given by with . In the present case of with the parametrization defined above the Haar measure is defined by 
The total group volume, first computed by Marinov , is then given by
In summary, averages over the manifold weighted by the probability density can be carried out by computing the 15-dimensional integral
The transformed state 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.
2.3 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 this is a -matrix with the elements
In general the reduced density matrix is no longer pure and its von-Neumann entropy quantifies the entanglement between the two qubits. In order to compute the entropy we determine the eigenvalues of which are given by
Having determined these eigenvalues, the entanglement of is given by
Finally, the entanglement has to be averaged over all trajectories (see Eq. (51)), 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.
2.4 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  defined by
where are the decreasingly sorted square roots of the eigenvalues of the matrix
In this expression 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
3 Quenched random interactions
In the case (a) of quenched randomness each element of the ensemble is associated with a time-independent random Hamiltonian . 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 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 and the eigenvectors factorize and thus can be treated independently. More specifically, the eigenvalues are known to be distributed as
where 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 . 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
3.1 Entanglement of the averaged density matrix
Let us now compute the average density matrix
First we compute the average over the energies
where is the normalization factor. This leads us to the result
where we defined the scaled time
and the function
Thus Eq. (32) reduces to
What remains is to determine the operators
Obviously, the first sum in Eq. (37) is given by
As for the second sum in Eq. (37), we note that the distribution of eigenvectors is invariant under a permutation of the basis vectors , hence the four operators coincide. Moreover, under a unitary transformation they transform as
where we have used the invariance of the GUE-eigenvectors under the replacement . This means that is invariant under if and only if commutes with the initial state . For a pure initial state this implies that has to be a linear combination of the identity and the initial state itself, i.e. . The linear coefficients and can be determined as follows. On the one hand, the identity
implies that , hence . On the other hand, we note that
is invariant under unitary transformations of and independent of , hence we may choose and to obtain
giving . Therefore, we arrive at the convex combination of and
with given in Eq. (36) and , which holds for any pure initial state . As expected, the averaged state lies on the segment between the initial state and the maximally mixed state due to the symmetries of the Haar measure of .
Having computed the mixed state of the ensemble we can now compute the corresponding entanglement of formation as a function of time. For a non-entangled initial pure state we find that for all times. However, if we start from the Bell state (2) with the initial entanglement we find numerically that the entanglement first decreases and vanishes at the finite scaled time (see Fig. 2).
3.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 (19)), we did not succeed to compute the average. For this reason let us consider the GUE average of the linear entropy
where denotes the reduced density matrix of the first qubit. In the qubit basis (29) this can be rewritten as
Inserting (27) and exploiting again that the GUE average factorizes, we get
where . For the initially non-entangled state this expression reduces further to
with the scaled time . As shown in D, the averages can be computed directly by integration over the given probability distributions in GUE, leading us to the final result
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 at , 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 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 .
4 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 manifold. In this case the quantity of interest is the probability distribution to find the time evolution operator with the Euler angles at time . This probability distribution allows one to compute the ensemble average of any function (such as the density matrix or the entanglement ) by integration over the complete volume weighted by , i.e. we have to compute the integral
4.1 Expected average entanglement of a uniform distribution
Before studying the temporal evolution in detail, let us consider the limit , 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  conjectured a closed expression for the expected average entanglement of an arbitrary bipartite quantum system with Hilbert space dimensions and , 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 and :
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.
4.2 Heat conduction equation on the manifold
In order to compute the probability distribution one has to solve the heat conduction equation
on the curved group manifold. Here denotes the diffusion constant while is the so-called Laplace-Beltrami operator which generalizes the ordinary Laplacian on a curved space. As stated by  this Laplace-Beltrami operator is the Markov generator of the unitary brownian motion (see F). On a Riemannian manifold with the metric tensor the Laplace-Beltrami-Operator is given by
where are the components of the inverse metric tensor and .
To our best knowledge the explicit expressions for , 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 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.
4.3 Early-time expansion
The solution of the heat conduction equation (54) and as well the averaged function can be expanded as a Taylor series around :
Using (51) we can compare the coefficients of the both Taylor series. giving
The r.h.s. is an integral over derivatives of the probability density evaluated at . If all trajectories start at it is easy to see that this probability density at is given by
Inserting this expression into (59) the integral can be evaluated by partial integration, giving
4.4 Average of the density matrix
Using (10) we calculate the derivatives . We find that
Therefore we can express all higher derivatives of in terms of the first derivative222Although this remarkable property calls for a deeper reason, we have no convincing explanation so far.
Hence, the solution of the averaged density matrix can be written as
Using the non-entangled initial state (which corresponds to taking ) this result specializes to
As can be seen, this density matrix relaxes exponentially and becomes fully mixed in the limite .
4.5 Average of the linear entropy
The same calculation can be carried out for the linear entropy defined in (21). Here we find that
For this reason, the calculation is completely analogous, giving
The result for an non-entangled initial state is therefore
For a fully entangled initial state we get