Entanglement formation under random interactions
Abstract
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 timeindependent 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 byproduct we compute the metric tensor and its inverse as well as the LaplaceBeltrami for .
1 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 twoqubit system is initially prepared in a welldefined pure state. As examples we consider two different initial states, namely, a nonentangled pure state
(1) 
and in a fully entangled Bell state of the form
(2) 
where denotes the canonical qubit configuration basis. Starting with the given initial state the system then evolves unitarily as
(3) 
where the time evolution operator is determined by and
(4) 
with a randomly chosen interaction Hamiltonian.
Throughout this paper we consider two different types of randomness, namely

Quenched randomness, where is timeindependent. In this case one considers an ensemble of twoqubit 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 timedependent Hamiltonian 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 .
The difference between the two cases is visualized in Fig. 1. In this figure the big red sphere stands symbolically for the 15dimensional group manifold of . Each of point on the sphere represents a certain unitary transformation acting on the 4dimensional Hilbert space of the twoqubit 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 vonNeumann entropy of the reduced states
(5) 
where denotes the timedependent reduced density matrix of the respective qubit. In cases where the logarithm is too difficult to evaluate we resort to the socalled linear entropy
(6) 
as an alternative entanglement measure. Note that both measures can be obtained from the more general Tsallis entanglement entropy [3]
(7) 
in the limit and , respectively.
Furtheremore, the Rényi entanglement entropy [4]
(8) 
is of interest. Also this entropy measure generalizes the vonNeuann 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 2qubit 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 byproduct, we compute the metric tensor and its inverse on the group manifold as well as the corresponding LaplaceBeltrami 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. [9]. As reviewed in the Appendix, the group elements are generated by 15 GellMann matrices , allowing one to parametrize unitary transformations by
(9)  
where the 15 Eulerlike angles vary in certain ranges specified in (80). Applying such a unitary transformation to the nonentangled initial state one obtains the density matrix
(10) 
with the components
(11) 
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 timedependent 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 :
(12) 
Here is the integrated group volume which serves as a normalization factor while
(13) 
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 [10], which is by itself invariant under unitary transformations.^{1}^{1}1For 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 [9]
The total group volume, first computed by Marinov [11], is then given by
(14) 
In summary, averages over the manifold weighted by the probability density can be carried out by computing the 15dimensional integral
(15) 
with the measure (2.2) and the integration ranges specified in (80). The uniform Haar measure corresponds to taking .
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
(16) 
For the initial state this is a matrix with the elements
(17) 
In general the reduced density matrix is no longer pure and its vonNeumann entropy quantifies the entanglement between the two qubits. In order to compute the entropy we determine the eigenvalues of which are given by
(18) 
Having determined these eigenvalues, the entanglement of is given by
(19) 
Finally, the entanglement has to be averaged over all trajectories (see Eq. (51)), i.e.
(20) 
However, if the average of the vonNeumann entropy is too difficult to compute, we will also use the linear entropy
(21) 
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 [12] defined by
(22) 
where are the decreasingly sorted square roots of the eigenvalues of the matrix
(23) 
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
(24) 
where .
3 Quenched random interactions
In the case (a) of quenched randomness each element of the ensemble is associated with a timeindependent random Hamiltonian . Since the spectral decomposition
(25) 
of a randomly chosen Hamiltonian is always nondegenerate, the time evolution operator can be written as
(26) 
Hence the state of the system evolves as
(27) 
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
(28) 
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 fourdimensional Hilbert space according to Haar measure, independent of the eigenvalues. If one defines the qubit basis
(29) 
this average can be carried out by setting
(30) 
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
(31) 
3.1 Entanglement of the averaged density matrix
Let us now compute the average density matrix
(32) 
First we compute the average over the energies
(33) 
where is the normalization factor. This leads us to the result
(34) 
where we defined the scaled time
(35) 
and the function
(36) 
Thus Eq. (32) reduces to
(37) 
What remains is to determine the operators
(38) 
Obviously, the first sum in Eq. (37) is given by
(39) 
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 GUEeigenvectors 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
(41) 
implies that , hence . On the other hand, we note that
(42) 
is invariant under unitary transformations of and independent of , hence we may choose and to obtain
(43) 
giving . Therefore, we arrive at the convex combination of and
(44) 
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 nonentangled 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 vonNeumann entanglement entropy of the individual pure states would be straightforward 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
(45) 
where denotes the reduced density matrix of the first qubit. In the qubit basis (29) this can be rewritten as
(46) 
Inserting (27) and exploiting again that the GUE average factorizes, we get
where . For the initially nonentangled 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
(49)  
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
(50) 
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 longtime 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 Timedependent 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
(51) 
over the ranges specified in (80), where denotes the volume element according to the Haar measure defined in Eq. (13).
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 [5] 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 :
(52) 
Applying this formula to a random twoqubit system, one obtains
(53) 
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
(54) 
on the curved group manifold. Here denotes the diffusion constant while is the socalled LaplaceBeltrami operator which generalizes the ordinary Laplacian on a curved space. As stated by [17] this LaplaceBeltrami operator is the Markov generator of the unitary brownian motion (see F). On a Riemannian manifold with the metric tensor the LaplaceBeltramiOperator is given by
(55) 
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 Earlytime expansion
The solution of the heat conduction equation (54) and as well the averaged function can be expanded as a Taylor series around :
(56) 
(57) 
Using (51) we can compare the coefficients of the both Taylor series. giving
(58) 
where denotes the volume element defined in Eq. (13). Using the heat equation (54) we can replace the partial derivative, obtaining
(59) 
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
(60) 
Inserting this expression into (59) the integral can be evaluated by partial integration, giving
(61) 
4.4 Average of the density matrix
Using (10) we calculate the derivatives . We find that
(62) 
Therefore we can express all higher derivatives of in terms of the first derivative^{2}^{2}2Although this remarkable property calls for a deeper reason, we have no convincing explanation so far.
(63) 
Hence, the solution of the averaged density matrix can be written as
(64) 
Using the nonentangled initial state (which corresponds to taking ) this result specializes to
(65) 
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
(66) 
For this reason, the calculation is completely analogous, giving
(67) 
The result for an nonentangled initial state is therefore
(68) 
For a fully entangled initial state we get