Dynamics and manipulation of entanglement in coupled harmonic systems with many degrees of freedom

We study the entanglement dynamics of a system consisting of a large number of coupled harmonic oscillators in various configurations and for different types of nearest-neighbour interactions. For a one-dimensional chain, we provide compact analytical solutions and approximations to the dynamical evolution of the entanglement between spatially separated oscillators. Key properties such as the speed of entanglement propagation, the maximum amount of transferred entanglement and the efficiency for the entanglement transfer are computed. For harmonic oscillators coupled by springs, corresponding to a phonon model, we observe a non-monotonic transfer efficiency in the initially prepared amount of entanglement, i.e. an intermediate amount of initial entanglement is transferred with the highest efficiency. In contrast, within the framework of the rotating-wave approximation (as appropriate, e.g. in quantum optical settings) one finds a monotonic behaviour. We also study geometrical configurations that are analogous to quantum optical devices (such as beamsplitters and interferometers) and observe characteristic differences when initially thermal or squeezed states are entering these devices. We show that these devices may be switched on and off by changing the properties of an individual oscillator. They may therefore be used as building blocks of large fixed and pre-fabricated but programmable structures in which quantum information is manipulated through propagation. We discuss briefly possible experimental realizations of systems of interacting harmonic oscillators in which these effects may be confirmed experimentally.


Introduction
Quantum information processing requires as a basic ingredient the ability to transfer quantum information between spatially separated quantum bits, either to implement a joint unitary transformation or, as a special case, to swap quantum information between the qubits. For a transfer over larger distances it is usually imagined that some stationary qubits, for example in the form of trapped ions inside an optical resonator, are coupled with a quantized mode of the electromagnetic field that propagates between the spatially separated cavities [1]- [4]. Other specific realizations are possible, but the basic principle always relies on the use of some continuous degree of freedom between the qubits and their manipulation by external fields. While this appears to be the most realistic mode of transport over long distances, one may conceive other modes over shorter distances. Instead of using a quantum field one may study the possibilities offered by a discrete set of interacting quantum systems. This might involve spin degrees of freedom [5]- [9] or infinite-dimensional systems such as harmonic oscillators [10,11]. In the present paper, we explore the dynamics of entanglement in a chain of coupled harmonic oscillators [12]- [14]. Apart from its obvious relevance to quantum optical systems including photonic crystals, such a model also describes phonons in a crystal and we therefore hope that the results presented here will also have applications in condensed matter systems as well. It should be noted that the entanglement properties of a system of harmonic oscillators in the static regime have been studied in some detail [12]. This paper is organized as follows. In section 2, we present the basic physical models and their Hamiltonians. We use analytical tools from the theory of Gaussian states in continuous variable systems where some rapid development has been achieved recently (see e.g. [10] for a tutorial overview). We briefly reiterate those results that will be employed in the present investigation. Section 3 will then present the basic equations of motion in compact form for two types of interactions, namely (i) harmonic oscillators coupled by springs and, resulting from this, (ii) a model which corresponds to a rotating-wave approximation (RWA) as is appropriate in a quantum optical setting. Section 4 employs these equations for the propagation of entanglement along a chain of harmonic oscillators which might be realized by coupled nanomechanical oscillators [15,16] or optical cavities. The time evolution of the entanglement between a pair of oscillators is given analytically in a compact form. Properties such as the speed of propagation, the amount of entanglement and the transfer efficiency are then obtained from these expressions. In section 5, we present results utilizing the equations from sections 3 and 4. First, we study a method for the creation of entanglement in such a system that does not require detailed control of the interaction strength between individual oscillators but only the ability for changing the interaction strength globally [11]. The influence of imperfections such as finite temperatures or randomly varying coupling constants on such a scheme are studied. We also consider the propagation of some initially prepared entangled state along the chain. Surprisingly, for harmonic oscillators coupled by springs, we observe a non-monotonic transfer efficiency in the initially prepared amount of entanglement, i.e. an intermediate amount of entanglement is transferred with the highest efficiency. Conversely, in the RWA, the transfer efficiency is monotonic. Whereas most of these results assume a position-independent and stationary coupling, we also show that with carefully chosen position-dependent coupling the transfer efficiency in this system may be increased to unity. Finally, we study geometrical configurations that are analogous to quantum optical devices such as beamsplitters and interferometers and observe characteristic differences when initially thermal or squeezed states are entering these devices. We show that these devices may be switched on and off by changing the properties of an individual oscillator and may therefore be building blocks of large fixed but programmable structures. In section 6, we summarize the results of this paper and suggest possible experimental realizations of systems of harmonic oscillators in which these effects may be confirmed.

Models and methods
In this section, we present the systems under consideration, namely coupled harmonic oscillators, together with the Hamiltonians that describe the various models for their interaction. We will restrict our attention to Hamiltonians that are quadratic in position and momentum operators. This will be crucial for the following analysis as it permits us to draw on the results and techniques from the theory of Gaussian continuous variable entanglement. The most important results from this theory will be reviewed here briefly.

The physical models
The general set-up consists of a chain of M coupled harmonic oscillators, where the coupling is assumed to be such that the corresponding Hamiltonian is at most quadratic in position and momentum. We will number the harmonic oscillators from 1 to M with periodic boundary conditions such that the (M + 1)th oscillator is identified with the first. The choice of periodic boundary conditions yields exact and compact analytical solutions since we can employ normal co-ordinates straightforwardly. A similar approach is less successful in the non-periodic boundary case. In addition, we allow for the existence of a distinguished decoupled oscillator with index 0 which will be a convenient notation for some of the later studies. Arranging the position and momentum operators in the form of a vector R = (q 0 ,q 1 , . . . ,q M ,p 0 ,p 1 , . . . ,p M ), (1) we can then write the general Hamiltonian in the form where V is the potential matrix and T the kinetic matrix. We will consider three basic settings for which we now provide the matrices T and V explicitly. In all these cases, we assume that the oscillators in the chain are all identical with a mass m = 1 and eigenfrequency ω = 1.
(a) Uncoupled oscillators. If the oscillators are not coupled with each other, then the potential energy of the kth oscillator is simply given byq 2 k /2 while its kinetic energy isp 2 k /2. As a consequence, both the potential and the kinetic matrices are diagonal and identical, namely Oscillators coupled by springs. If neighbouring oscillators (except for the 0th oscillator) are coupled via springs that obey Hooke's law, the Hamiltonian is given by where c denotes the coupling strength and we have used periodic boundary conditions, i.e. q M+1 = q 1 . Keeping in mind that we wish to leave the oscillator with index 0 uncoupled, the potential matrix becomes whereas the kinetic matrix is given by the identity matrix T = 1 M+1 .
(c) Oscillators in the RWA. An interaction that provides simpler dynamics is obtained via the RWA in quantum optical systems. Indeed, if we define annihilation and creation operatorŝ then we observe that the interaction in case (b) includes terms of the formâ † kâ † k+1 , i.e. interaction terms for which both harmonic oscillators are being excited simultaneously. Such terms are not energy-conserving, and in quantum optics they are usually be neglected in the framework of the RWA. Following this practice amounts to considering the Hamiltonian In terms of position and momentum operators this can be written as so that in this case both the potential and the kinetic matrices are given by Note that the matrix V can be conceived as the adjacency matrix of a weighted graph G(v, e) encoding the interaction pattern between the systems in the canonical co-ordinates corresponding to position. Vertices of the graph correspond to physical systems, i.e. the individual harmonic oscillators, whereas the weight associated with each of the edges quantifies the coupling strength [17]. The main diagonal corresponds to loops of the weighted graph. This intuition is in immediate analogy to graph states for spin systems with an Ising interaction between the constituents [18]- [20] and can be useful in the study of more complex geometries. Before we study the dynamical properties of these systems, we provide in the following subsection a brief overview over the main technical tools that we will employ.

Analytical tools
Analysing the entanglement properties of infinite-dimensional systems such as harmonic oscillators is generally technically involved unless one restricts attention to specific types of states. Indeed, in recent years, a detailed theory of so-called Gaussian entangled states has been developed. As we will employ some of its basic results in the subsequent analysis, we are providing a brief review of some useful results. A more detailed tutorial review can be found, e.g. in [10], and more technical details concerning Gaussian states can be found in [21].
The central variables in the analysis of harmonic oscillators are the canonical operators for position and momentum. Let us assume a system with n harmonic oscillators. As stated above it is convenient to arrange these in the form of a vector R T = (q 1 , . . . ,q n ,p 1 , . . . ,p n ).
The characteristic feature distinguishing the quantum harmonic oscillator from its classical counterpart is the canonical commutation relation (CCR) between position and momentum. Employing the vector R these can be written in the particularly convenient form [ R j , R k ] = iσ j,k , where the real skew-symmetric block diagonal (2n × 2n) matrix σ, the symplectic matrix, is given by assuming units where h = 1 and k = 1, a choice that will be adopted for the rest of this paper. Instead of referring to states, i.e. density operators, one may equivalently refer to functions that are defined on phase space. Whereas there are many equivalent choices for phase-space distributions, for the purposes of this work it is most convenient to introduce the (Wigner-) characteristic function. Using the Weyl operator W ξ = e iξ T σR for ξ ∈ R 2n , we define the characteristic function as The state and its characteristic function are related to each other according to a Fourier-Weyl relation Gaussian states are exactly those states for which the characteristic function χ ρ is a Gaussian function in phase space [21], i.e. if the characteristic function is of the form As is well known, Gaussians are completely specified by their first and second moments, d and γ, respectively. As the first moments can be always made zero utilizing appropriate local displacements in phase space, they are not relevant in the context of questions related to squeezing and entanglement and will be ignored in the following. The second moments can be collected in the real symmetric 2n × 2n covariance matrix γ defined as With this convention, the covariance matrix of the nth-mode vacuum is γ = 1 2n . As the covariance matrix encodes the complete information about the entanglement properties of the system, we will use it to quantify the amount of entanglement between two groups of oscillators. There is no unique way to quantify entanglement for mixed states, and several different measures grasp entanglement in terms of different operational interpretations. For the purposes of this work, we settle for the logarithmic negativity which is comparatively easy to compute and possesses an interpretation as a cost function [22]- [26]. Given two parties, A and B, the logarithmic negativity is defined as where ρ T B is the state that is obtained from ρ via a partial transposition with respect to system B and · 1 denotes the trace norm. As we focus attention on Gaussian states which we characterize via the covariance matrix γ rather than the density matrix ρ, we need to provide a prescription for the evaluation of the logarithmic negativity directly in terms of the covariance matrix. To this end, it is important to note that on the level of covariance matrices the transposition is reflected by time reversal which is a transformation that leaves the positions invariant but reverses all momenta,q →q,p → −p. The partial transposition is then the application of this time-reversal transformation to a subsystem, i.e. one party. Let us now consider a system made up of m + n oscillators, where m oscillators are held by party A and n oscillators by party B. Applying time reversal to the n oscillators held by party B, the covariance matrix will be transformed to a real symmetric matrix γ T B given by where The (2n × 2n) matrix γ T B is the matrix collecting the second moments of the partial transpose ρ T B of ρ. The logarithmic negativity is then given by where the γ j are the symplectic eigenvalues of γ T B . For a general covariance matrix, γ j arises in the diagonalization of γ using symplectic transformations, i.e. transformations S that preserve the CCR so that SσS T = σ. The resulting diagonal matrix is the Williamson normal form of a covariance matrix whose diagonal elements are the symplectic eigenvalues. It is sometimes useful to know that the symplectic eigenvalues can be obtained directly without explicit diagonalization of the matrix as the positive square roots of the usual eigenvalues of −σγσγ [24]. For all Hamiltonians that are quadratic in the canonical position and momentum operators the ground state is an important example of a Gaussian state. For a Hamilton operator of the form we find that the covariance matrix of the ground state is given by which reduces to when T = 1 n . If, on the other hand, as for the interaction in case (c), we have T = V , then the ground state is given by γ = 1 n ⊕1 n , which is the same as the ground state of n non-interacting harmonic oscillators. The primary aim of this work is the investigation of the dynamical properties of the system of harmonic oscillators and the evolution of entanglement properties under such dynamics. The dynamics of the covariance matrix under a Hamiltonian quadratic in position and momentum operators can be obtained straightforwardly from the Heisenberg equation For our time-independent Hamiltonian (18), this leads to the covariance matrix at time t as Equipped with these tools we can now proceed to the analysis of the entanglement dynamics of the harmonic chain.

The equations of motion
In two separate subsections, we provide the explicit solutions to the equations of motion for the two coupling models (b) and (c) that will be investigated both analytically and numerically in the remainder of the paper.

Harmonic oscillators coupled by springs
This model is characterized by a Hamiltonian of the form Note that for the moment we neglect the decoupled additional 0th oscillator. In the following, we will provide an explicit form for the equations of motion for the canonical positions and momenta in the Heisenberg picture. To this end, we can diagonalize this Hamiltonian by introducing the normal co-ordinateŝ This leads to where Here we have used the fact that Q u = Q † −u , P u = P † −u . We introduce the annihilation and creation operatorŝ so that the Hamiltonian takes the form In the Heisenberg picture the annihilation and creation operators evolve according toâ s (t) = e −iω s tâ s (0) andâ † s (t) = e +iω s tâ † s (0). Separating the real and imaginary parts, we get Substituting these into equation (24), we obtain the time evolution for the original position and momentum operatorŝ q n (t) = M r=1 (q r (0)f r−n (t) +p r (0)g r−n (t)), where we have defined the useful functions Throughout the paper, dots will denote time derivatives. Defining the covariance matrix elements to be (once we ignore the displacements q i ρ ) we find their values at time t as where the γ on the right-hand side are the initial values of the covariance matrix elements.

This model is specified by the Hamiltonian
In the following, we can carry out the analysis along the same lines as in the previous subsection. Again, employing the normal co-ordinates, equations (24), we obtain where we now have Introducing the annihilation and creation operatorŝ the Hamiltonian assumes a form In the Heisenberg picture, the annihilation and creation operators then evolve in time asâ s (t) = e −i 2 s tâ s (0) andâ † s (t) = e +i 2 s tâ † s (0). Separating again real and imaginary parts we obtain Transforming back to the original position and momentum operators we find where we have defined the functions F k and G k as Note that these functions are slightly simpler than the corresponding ones in equations (31) as they lack the frequency s in the denominator. The covariance matrix elements vary in time as

Propagation of entanglement along the chain
We investigate the capacity of the harmonic chain for transmission of quantum information. The clearest signature for the ability to transmit quantum information and coherence is the proof of the ability to transmit one constituent part of an entangled pair of oscillators through the chain.
To analyse the situation, we require the equations of motion for the covariance matrix. Let us assume the existence of a distinguished oscillator 0 which is entirely decoupled from the others. We imagine that at time t = 0 this oscillator is prepared in a two-mode squeezed state with the first oscillator of the chain, whereas all other oscillators are assumed to be in their respective ground state (assuming no interaction). The initial conditions then read The 0th oscillator for the interaction via springs will obey a free time evolution which is given by (using equations (29) and noting thatq 0 =Q 0 andp 0 =P 0 ) Similarly, for the RWA interaction, Note, however, that they correspond to local unitary rotations on the 0th oscillator which do not affect the entanglement between this oscillator and the remaining ones. To simplify the expressions, we will therefore omit this time evolution in the following. Again we will treat the two types of interactions described by H Spring and H RWA separately. The next two subsections establish the analytical expressions for the time evolution of the entanglement between the 0th oscillator and the nth oscillator.

Harmonic oscillators coupled by springs
In this case, employing the special form of the initial conditions for the system (43), the elements of the covariance matrix describing the 0th and the nth oscillator at time t are given by and In the limit of a chain of infinite length, i.e. when M → ∞, we find and Here we have employed the definitions ζ = c/(1 + 2c) and = √ 1 + 2c, and introduced the functions and a nn (t) = While the above set of equations determines the time evolution exactly, they do not provide very much direct insight into the dynamics of the system. In the following, we will show, however, that we can obtain very good and compact approximations to the above exact solution in terms of elementary functions. Although the following derivation is not rigorous in the sense that it does not provide error bounds, a numerical comparison between the approximations and the exact results shows the impressive quality of the approximate solution. As a first step, we expand the functions f k and g k to first order in ζ, and also drop a factor of 1/ √ 1 − 2ζ cos φ in g k . In the following, we will employ Bessel functions which satisfy the relations On using trigonometrical addition theorems one finds that in first order A further crucial approximation replaces the time-dependent quantities a nn (t), b nn (t), c nn (t), d nn (t), e nn (t) by their time averages, i.e.
With all these approximations we finally obtain Numerical comparisons in later sections will show that these approximate solutions are very good approximations when r is not too small. We have so far collected all the elements of the covariance matrix involving the 0th and the nth oscillators. Since we will investigate the entanglement properties between the two oscillators, we can trace out the rest of the oscillators which leaves us with the covariance matrix of the reduced system comprising only two oscillators. Employing the ordering (q 0 ,p 0 ,q n ,p n ) we find where where we have used the abbreviations s (t) = t − π(n − 1)/2 and J s (t) = J n−1 (ζ t). From this explicit form for the covariance matrix of the 0th and nth oscillators we can now determine the symplectic eigenvalues as solutions of the polynomial [24] The solution is given by where we have Note that we have dropped the time argument in J s (t) to make the expressions appear more compact. The logarithmic negativity, finally, is in this approximation given by Note that the other symplectic eigenvalue is greater than unity.

Hamiltonian in the RWA
For the Hamiltonian H RWA one proceeds along very similar lines, i.e. taking the limit M → ∞ to find We find the covariance matrix elements to be Note that now the terms corresponding to a nn (t), b nn (t), c nn (t), d nn (t) and e nn (t) do not need to be approximated as all time dependences cancel each other out conveniently in the expressions as opposed to the spring case. Again we can write the covariance matrix of the reduced system comprising only the 0th and the nth oscillator employing the ordering (q 0 ,p 0 ,q n ,p n ), and find where we have used the abbreviations RWA (t) = 2 RWA t − π(n − 1)/2 and J RWA (t) = J n−1 (ct). Note that these expressions are very similar to those for the first interaction type. The symplectic eigenvalues are given by the solutions of This gives rise to with Again we have dropped the explicit time dependence in J RWA for brevity of notation. The logarithmic negativity is then given by Having prepared all the analytical work we need, we can now proceed to investigate the entanglement dynamics of the harmonic chain.

Study of the entanglement dynamics of the harmonic chain
In this section, we study numerically and analytically, various aspects of the entanglement dynamics of the harmonic chain with various initial states and geometrical arrangements for both types of interactions. We begin by briefly revisiting the effect of the spontaneous creation of entanglement which is obtained when the interaction strength between the oscillators is globally changed suddenly [11]. This effect can only be observed in a model in which the oscillators are coupled via springs. In the RWA, this effect does not occur as both the coupled and uncoupled chains have an identical ground state. Then we move on to consider the propagation of entanglement through a harmonic chain. Surprisingly, for harmonic oscillators coupled with springs corresponding to a phonon model, we observe a non-monotonic transfer efficiency in the initially prepared amount of entanglement, i.e. an intermediate amount of entanglement is transferred with the highest efficiency. In the RWA, the transfer efficiency is monotonic although equally surprising. We provide analytical expressions for the propagation speed of the entanglement through the chain together with approximate analytical expressions for the transfer efficiency. We also study the influence of imperfections such as finite temperatures and different coupling constants. Finally, we study different geometrical configurations that are analogous to quantum optical devices such as beamsplitters and interferometers and observe characteristic differences when initially thermal or squeezed states are entering these devices. We propose ways in which these devices may be in and out of action, thereby allowing for the creation of pre-fabricated quantum networks that can be programmed by external switches.

Spontaneous creation of entanglement
We consider all the harmonic oscillators to be in the ground state and initially uncoupled, i.e. γ q n q m = δ nm , γ q n p m = 0, γ p n p m = δ nm .
We suddenly switch on the interaction at time t = 0 and observe the dynamical evolution of entanglement. We do not observe any entanglement with the H RWA interaction (because we have not included a mechanism to produce a spontaneous excitation, i.e. we do not have terms such as a † k a † k+1 ), so we will exclusively deal with H Spring . We limit our scope to numerical results because the analytical expressions, while they can be provided, are too complicated to yield any insight. Note also that the approximations used below in the squeezed state case cannot be applied here. A typical example of the time evolution of entanglement is shown in figure 1.
In an open chain of length 30, we study the time evolution of the entanglement between the first and the last oscillator. We observe no entanglement for a finite time until at a time t 0 one first encounters a build-up of entanglement between the two oscillators. This time t 0 is approximately given by i.e. the time t 0 is approximately linear in n, the separation of the oscillators. It should be noted that t 0 is half as large as the time that is required for a perturbation at the first oscillator to travel to the nth oscillator. This suggests that the origin of the entanglement between the first and the nth oscillator arises from the interaction of those oscillators exactly half-way in between. Their entanglement is generated by the initial sudden switching of the interaction and then propagates through the chain. This idea will be further corroborated in the later subsections when the propagation of pre-prepared entanglement is considered. Furthermore one finds that the dependence of the maximal degree of entanglement as measured by the logarithmic negativity is approximately proportional to n −1/3 for large n until it reaches values of about 10 −2 when it begins to drop quite rapidly to vanish entirely. It should be noted that for any parameters of the model, there exists a finite n such that the state of the first and the nth oscillator is separable at all times. For the coupling value of c = 0.1 this will happen for n ∼ = 20 000. Therefore, the value of largest separation is very large indeed for reasonable coupling strengths. To appreciate that this is a somewhat surprising behaviour, it should be contrasted with the entanglement structure in the ground state of the chain. Then it can be shown that for any chosen coupling strength two distinguished oscillators are never in an entangled state, unless they are immediately neighbouring [12].
Here we have studied the special case of the entanglement between the endpoint of an open chain. It should be noted that this is a particularly favourable configuration. For a given distance between oscillators one always obtains the largest amount of entanglement when one places them at the opposite ends of an open chain. Two oscillators in a very long chain with the same distance and with positions well away from the ends of the chain will lead to considerably smaller amounts of entanglement. Indeed, the amount of entanglement will differ by approximately a factor of 4. This discrepancy is due to the fact that at the ends of the chain the oscillators possess fewer neighbours with which they can become entangled. As we discard all oscillators other than two any entanglement with other oscillators will deteriorate the entanglement between the distinguished oscillators. While this does not explain the factor of 4 quantitatively, it gives an intuitive picture for the decrease of entanglement that will be discussed in more detail later on.

Thermal state case
In the previous section, we have studied the entanglement dynamics in an environment that is at zero temperature. This is reflected by the fact that the initial state of our harmonic oscillators is assumed to be the vacuum. We will now move one step further towards a more realistic description by setting the initial state as a thermal equilibrium state. As with the ground-state case, we do not establish entanglement for the H RWA interaction as a thermal state can be represented as a mixture of displaced vacuum states which do not lead to any spontaneous entanglement in the RWA. Therefore, we shall again only consider the H Spring interaction. The thermal equilibrium state is given by The equations of motion are then γ q n q m (t) = 1 + 2 e ω/T − 1 (a nm (t) + d nm (t)),

c nm (t) + a nm (t)).
Appropriate values for the temperatures have to be obtained from experiment in the particular context under consideration, but it appears possible nowadays to achieve ratios T/ω 1 in different physical systems (note that we have taken h = 1 and k = 1) such as nanomechanical oscillators. We consider again an open chain instead of a closed ring. This renders the analytical treatment more difficult but makes no difference for the numerics. The entanglement depends little on the temperature as long as the mean thermal photon number is well below unity. Figure 2 shows the temperature dependence of entanglement evolution. We observe that down to temperatures corresponding to x = 10, the entanglement evolution is almost exactly the same as in the groundstate case. Only when x < 10 we see an effect of a finite temperature. Even for x = 6, we still have a significant portion of entanglement albeit a small delay in the arrival of entanglement. A more realistic scenario is dealt with in [11]. Decoherence mechanisms may be included without leaving the harmonic setting. Often, the high-temperature limit of Ohmic quantum Brownian motion is appropriate with an independent heat bath for each oscillator in the limit of negligible friction and under the assumption of product initial conditions. Such a decoherence mechanism can be accounted for by adding terms of the form [27,28] −ξ[q n , [q n , ρ]] to the idealized unperturbed generator of the dynamical map for each of the oscillators, where the real number ξ specifies the decoherence time scale. However, in cases where product initial conditions are inappropriate or unrealistic, decoherence may still be modelled using, for example, time-convolutionless projection operator techniques [29]. In small systems, nonproduct conditions may be incorporated by explicitly appending heat baths to each of the oscillators with a linear coupling [30], according to Hamiltonians of the form with real numbers ξ j , where theq (i) j denote the canonical co-ordinates corresponding to position of the ith oscillator of the jth heat bath consisting of m oscillators. Assuming a particular form of the spectral density, the coupling strength to the finite heat baths may be chosen in a way that is consistent with empirically known values for energy dissipation. Often, Q-factors are approximately known for resonators, which quantify the number of radians of oscillations necessary for the energy to decrease by a factor of 1/e. Hence, on the basis of these Q-factors, the appropriate coupling may be evaluated. Figure 3 shows the influence of decoherence in case of an open chain with the same parameters as in figure 1 for an Ohmic heat bath, i.e. ξ j = j /m, where > 0 is the cut-off frequency of the environment modes. One finds that the created entanglement by suddenly switching on the interaction is surprisingly robust against decoherence within this model.

Entanglement transport through the harmonic chain
In the previous section, we considered the case where no entanglement was present in the initial state of the system. Entanglement emerged as a consequence of a sudden change in the coupling constant between neighbouring harmonic oscillators. In this subsection, we will investigate a different situation and consider the transmission of entanglement through a onedimensional chain. To this end, we initialize two harmonic oscillators in a two-mode squeezed state. We assume that one of these oscillators is decoupled from the rest and assign it the index 0, while the other oscillator, with index 1, forms part of a chain of harmonic oscillators with nearest-neighbour interaction. The remainder of the chain starts out in the ground state corresponding to zero temperature (assuming no interactions). By evolving the initial state, we expect the entanglement to travel along the chain so that with increasing time more and more distant oscillators will be entangled with the 0th harmonic oscillator. There are a number of free parameters that can be varied: the coupling strength c, the amount of initial entanglement quantified by the two-mode squeezing parameter r, the time t and the position of the oscillator to be entangled with the 0th oscillator. To simplify the analytical work, we shall be dealing with the limit M → ∞.  While in the previous section, the interaction in the RWA does not lead to the spontaneous creation of entanglement, here it allows for the propagation of entanglement. We give an example for the time evolution of the logarithmic negativity between the 0th and the 30th oscillator for both interactions in figure 4. For both interactions we obtain qualitatively the same behaviour but we observe that under the RWA interaction the entanglement propagates somewhat faster but as expected this difference decreases with decreasing coupling constant c. Another difference is the fact that the entanglement under the RWA interaction does not exhibit the small-amplitude oscillations that the interaction due to harmonic oscillators coupled by springs exhibits due to the existence of counter-rotating terms of the formâ kâk+1 . The propagation of the quantum entanglement can be seen even more clearly in figure 5, where for a ring composed of 40 oscillators and a coupling constant of c = 0.1 the time evolution of the logarithmic entanglement between an uncoupled oscillator and the nth oscillator is shown when initially the uncoupled oscillator and the first oscillator are coupled. One observes that with increasing time more and more distant oscillators are becoming entangled. Entanglement propagates both clockwise and anti-clockwise around the ring. After a sufficiently long time it becomes important that the ring has a finite size and the two counter-rotating 'entanglement waves' meet at the opposite end of the ring and we observe some entanglement enhancement. Both figures 4 and 5 suggest that entanglement can be distributed to distant oscillators. It will therefore be interesting to study the efficiency for this transfer when we vary the amount of entanglement provided initially by varying r. In particular, we will be interested in the first local maximum in the amount of entanglement N f as quantified by the logarithmic negativity. We separate the study of the entanglement transfer efficiency for the two interactions as they exhibit distinctly different behaviours. We begin with the interaction describing oscillators interacting via springs. Figure 6  shows the amount of entanglement at the first local maximum as squeezing is varied. We observe the remarkable fact that for large initial entanglement, the value of N f saturates. We can obtain an analytic expression for the saturation value as a function of c and n by taking the r → ∞ limit. Thus, we approximate η 2 Spring ∼ = −w/(8y 1 ) of equation (60) and, discarding all the terms except cosh 2 (r) ∼ = e 2r /4 and cosh(2r) ∼ = e 2r /2, we obtain where J max = 0.6748851(n − 1) −1/3 [31] is the value of the first maximum of the nth Bessel function of the first kind. This substitution provides a very good approximation as the first local maximum of the logarithmic negativity coincides with that of the nth Bessel function. Therefore, we find N sat = −log 2 (|η Spring |) and is shown in figure 7. We observe that the saturation value decreases for both increased coupling strength and increased distance. While the latter is intuitive the former might be surprising as one could have thought that it is advantageous to increase the coupling strength to facilitate the transfer of entanglement but this intuition is clearly contradicted by figure 7. Even more striking is that the entanglement vanishes entirely when the interaction strength becomes too high. We believe that this is due to the fact that the initial entanglement disperses across several oscillators, and we will discuss this at the end of this subsection.  If we translate these findings into an entanglement transfer efficiency defined by then we observe that this efficiency exhibits a non-monotonic behaviour. Indeed, we observe a maximum in the efficiency as shown in figure 8 for the same parameters as in figure 6. We have not yet found a convincing and intuitive explanation for the occurrence of a maximum in the transfer efficiency. 3 In fact, we will shortly see that the phenomenon of a non-monotonic transfer efficiency is absent in the RWA interaction. The surprising implication of this non-monotonic behaviour of the transmission efficiency is that it is advantageous to transmit entanglement in intermediate size portions rather than in one very large packet.
Let us now consider the same problem of the entanglement transfer efficiency in the RWA interaction. Indeed, we can show that there is still saturation in the amount of entanglement that can be transmitted and we find the value of the saturation to be (after taking r → ∞) where again J max = 0.6748851(n − 1) −1/3 [31] is the value of the first maximum of the nth Bessel function of the first kind. Note that this expression is independent of the coupling strength c. Unlike the case of oscillators interacting with springs there is no maximum in the efficiency for the RWA interaction as can be seen clearly in figure 9. Indeed, for large r, the efficiency is tending to zero while for r approaching 0 the efficiency tends to T eff = J max . Since the efficiency is not equal to unity for both interactions, the question arises as to where the rest of the entanglement is located. The most obvious place is to search in the neighbourhood of the nth oscillator. Since we always determine the entanglement between individual oscillators, we ignore many others that have interacted with it and have thereby become entangled. Any entanglement between these two oscillators will therefore deteriorate as they are being entangled with other oscillators that we choose to ignore. This viewpoint is corroborated by determining the entanglement between the 0th oscillator and a whole group of neighbouring oscillators instead of a single one. The result of this can be seen in figure 10. This shows the change in the amount of entanglement as we increase the number of oscillators in the second group. The graph supports the idea that the missing entanglement between the 0th and the nth oscillator is due to the creation of entanglement between the nth oscillator and its neighbours because we start to recover entanglement as we compute the entanglement between the 0th oscillator and the neighbourhood of oscillators surrounding the nth oscillator. This spread of entanglement is not dissimilar to a dispersion of the energy of a wavepacket as it experiences different group velocities. However, the effect on the entanglement can be considerably stronger as the energy will only decrease linearly with the width of the wave-packet, whereas the entanglement can drop much more rapidly and become zero at finite spreading.

Speed of entanglement propagation
We have found that the propagation of half the two-mode squeezed state through the chain takes a finite time as the entanglement between the 0th and the nth oscillator is exactly zero for a finite time interval (see e.g. figure 4). After a certain amount of time, the two oscillators in question become entangled and the logarithmic negativity reaches a temporary maximum. We are able to determine this time analytically for both types of interactions that we are considering. To make the analysis tractable, we consider an infinitely long chain. We find that the first maximum of the logarithmic negativity coincides with the first maximum of a Bessel function J m (x). We know the position of this maximum occurs at x = m + 0.8086165m 1/3 [31]. Noting that m = n − 1, x = ζ t = ct/ √ 1 + 2c for Hooke's law interaction and x = ct in the RWA interaction, we obtain We observe that the time that is required for entanglement to be established between the 0th and the nth oscillator is a function of the coupling strength and the position n. For large separations n, it quickly becomes linear in n and, as expected, the larger the coupling c the faster entanglement is established. We also see that the RWA interaction produces faster entanglement since t Spring /t RWA = √ 1 + 2c. As n − 1 is the separation of the first and the nth oscillator, we can define the speed of propagation to be Clearly, for large n, the speeds approach a constant dependent on c. For the RWA interaction, the propagation velocity increases linearly with c. This is an attractive feature, because, unlike the case of interaction via springs, the efficiency under the RWA does not decrease as we increase c because its efficiency is independent of c.

Optimization of entanglement transfer and generation
In the previous section, we have studied the entanglement transfer along a chain of identical harmonic oscillators as this will be the situation that is most easily implemented experimentally. However, we observe that the transmission efficiency decreases with distance. One might expect that one can improve this efficiency by tuning the couplings and the eigenfrequencies of the harmonic oscillators suitably. Indeed, in this section, we will show what can be achieved in this more general setting. For simplicity, we consider the task of transmitting one-half of a two-mode squeezed state from one end of an open chain to the other. We assume, as usual, one decoupled harmonic oscillator with index 0 and a chain of length M through which the other half of the two-mode squeezed state is transmitted. Perfect transmission from one end of the chain to the other is possible in the RWA of nearest neighbours if we choose the interaction strength and V n,n = 1, with the real number c being sufficiently small for V to be positive. The choice of the diagonal elements being all equal to 1 is equivalent to the requirement that we choose the eigenfrequencies ω n of the uncoupled oscillators as That the transmission is perfect can be shown by first realizing that in an interaction picture with respect to H 0 = i (q 2 i +p 2 i )/2 in which the diagonal elements of V vanish (this interaction picture will leave all entanglement properties unaffected as it is of direct sum form) we can replace where Q and P are column vectors, by the complex notation Q − iP so that Now we need to realize that V I is a quantum mechanical representation of a rotation. This allows the evaluation of the matrix elements of e iV I t . In particular, we find that so that one can generate an interchange between the first and the Mth co-ordinate by waiting for a time t = π/c (see [8] for an analogous argument in spin chains). Without the assumption of the RWA, i.e. choosing the Hamiltonian H Spring to describe the time evolution the above simple argument fails and indeed it is not possible to tune the nearestneighbour couplings alone to generate perfect transfer of entanglement for M > 2. However, if one chooses the couplings as above and decreases the value of the constant c, then, for a fixed distance, one can obtain arbitrarily good transfer efficiency at the expense of an increased delay time. This should not come as a surprise, as in the case of c → 0 the RWA becomes exact as the terms that are neglected in the RWA are of order O(c 2 ). Therefore, for entanglement distribution over a fixed distance they will play a decreasing role as the time of arrival for the entanglement is of the order O(c −1 ).
The case M = 2 is an exception, where one may realize an exact swap of the state of the 1st to the 2nd oscillator. To show this, note that specific covariance matrix elements of the 0th and the 2nd oscillator are given by Considering the functions f 1 and g 1 as specified in equations (31), we find that there exist a real number c and a time t such that simultaneously can be satisfied. This is the case when we chose the real c such that there exist natural numbers k and l such that and t = lπ. Then, it follows that, as the 0th oscillator is invariant and the state of the 0th and the 2nd oscillator necessarily corresponds to a pure Gaussian state, the 1st oscillator is necessarily decoupled from the other two. In this sense, the state can be swapped from one oscillator to the other, while retaining the entanglement with the 0th oscillator. Hence, one has a perfect channel for appropriate times.

Sensitivity to random variations in the coupling
In the preceding subsections, we have discussed an ideal model in which all experimental parameters can be determined perfectly. Any real experimental set-up however will suffer small variations in parameters such as the coupling strength between neighbouring oscillators. To confirm that the effects that have been found in this work can be observed in real experiments, we consider in the following the impact of random-position-dependent variations in the coupling strength between neighbouring oscillators. As an example, we consider an open chain with potential matrix where c i,j = c + c i,j is the position-dependent coupling between the ith and jth oscillators, where c i,j is a realization of a random variable distributed according to a normal N(0, c) distribution. For a chain of length 10, an average coupling constant of c = 0.1 and an initial two-mode squeezed state with squeezing parameter r = 0.8. Figure 11 shows the ratio of the first maximum for the case of slightly perturbed couplings over the idealized case versus the perturbation size c. We observe that for c/c 0.25, the achieved entanglement is greater than 95% of the degree of entanglement in the unperturbed case. Similar results apply for the RWA interaction. These results indicate that sending quantum information along the chain is stable under perturbations. Similar considerations can be made for the spontaneous creation of entanglement which show that the results are much more sensitive to perturbations. Indeed, with the same specifications as above, we find the entanglement at the first maxima to be between a small fraction and twice the amount for the non-perturbed case when c/c = 0.5. This suggests that the experimental demands for the verification of the spontaneous creation of entanglement are considerably higher than for the distribution of entanglement.

Other geometrical arrangements: beamsplitters and interferometers
So far we have studied only a linear chain of harmonic oscillators through which quantum entanglement can be propagated. However, it might be interesting to consider more complicated structures which may be used as building blocks for more complicated networks, in principle, any arrangement corresponding to an arbitrary weighted graph. In this subsection, we will study briefly two possible extensions of the linear chain, namely a Y-shaped configuration which can be used for the generation of entanglement and a configuration resembling an interferometer. We show furthermore how such configurations may be switched on and off thereby controlling the transport of quantum information in such a structure. A more detailed discussion of such structures and their optimization will be presented elsewhere. The material in this subsection should merely serve as examples for possible alternative ways of creating and manipulating entanglement through propagation in pre-fabricated structures. We begin by considering a chain ofY-shape which is shown in figure 12. One arm consisting of M in oscillators is connected to two further arms each consisting of M out oscillators.As usual, we consider nearest-neighbour interactions only and, for simplicity and the clearest demonstration of the effects, we restrict attention to the RWA interaction. We assume that the structure is initially in the ground state, i.e. at temperature T = 0. At time t = 0, we perturb the first harmonic oscillator exciting it either to a thermal state characterized by covariance matrix elements for some z, or a pure squeezed state characterized by covariance matrix elements As an example we choose a coupling constant of c = 0.2, and let the arms of the Y-shape contain M out = 30 oscillators each while the base contains M in = 10 oscillators. In figure 13, we present the results for the choice of a squeezed state with γ q 1 q 1 = 10 = 1/γ p 1 p 1 and a thermal state with γ q 1 q 1 = 10 = γ p 1 p 1 .  Squeezed state Figure 13. In a chain with Y-shape and nearest-neighbour interaction of RWA type the first oscillator at the foot of theY-shape is either excited to a squeezed state with γ qq = 10 = 1/γ pp or a thermal state with γ qq = 10 = γ pp . The remaining oscillators are in the ground state. The perturbation propagates along the chain into both arms of the Y-shape. For an initial thermal state excitation no entanglement is ever found between the ends of the two arms of theY-shape while entanglement is generated when the initial state is a squeezed state. The coupling constant is chosen as c = 0.2, the arms of the Y-shape contain 30 oscillators each while the base contains 10 oscillators. We observe that for an initial thermal state excitation no entanglement is ever found between the ends of the two arms of the Y-shape. This can be understood because a thermal state is a mixture of coherent states, i.e. displaced vacuum states. If the system is initialized in the vacuum state it will evidently not lead to any entanglement in the RWA and therefore an initialization in a thermal state cannot yield entanglement either. On the other hand, considerable entanglement is generated when the initial state is a squeezed state. It is possible to optimize the generation of entanglement by adjusting the strength of the nearest-neighbour couplings but this will be pursued elsewhere. These two observations resemble closely optical beamsplitters which do not create entanglement from thermal state input but can generate entanglement from squeezed inputs (see [34] for a comprehensive treatment of the entangling capacity of linear optical devices).
Another interesting set-up is shown in figure 14. We will henceforth call this the interferometric set-up. Let the number of oscillators on the left (including the junction), upper arm, lower arm and on the right be M L , M U , M D and M R , respectively. If we prepare a two-mode squeezed state between a decoupled oscillator and the leftmost oscillator of the interferometric set-up, then we are interested in how much entanglement propagates through the set-up depending on the properties of the two arms. One may vary different parameters such as the length of one of the arms, the coupling strength or eigenfrequency of the oscillators in one arm. We will focus on how the change in eigenfrequency ω of the harmonic oscillators in one of the arms affects the propagation of entanglement through the interferometric device. We change the eigenfrequencies of the oscillators smoothly across one arm following ω i = 1 + (ω − 1) × min(i, M U + 1 − i)/(M U /2) so that the oscillator half-way through the arm has frequency ω. Figure 15 shows the logarithmic negativity between the decoupled oscillator and the last one in the interferometric configuration at the time t = 250 plotted against ω. The other parameters are M U = M D = 30, M L = M R − 1 = 9 and c = 0.2. One clearly observes interference fringes in the frequency ω that are related to the effective path-length difference between the upper and the lower arms. The interference fringes do not have full amplitude and their amplitude is reduced for increasing ω. More sophisticated choices for the coupling parameters in the interferometric structure can improve on these imperfections. This demonstrates that in interferometric structure the transmission through the device will be strongly influenced by changes of the properties of one arm of the structure.
This shows that more complicated structures such as the Y-shape or the interferometric may be used to create entanglement from an initially unentangled system and transport it. There is a distinct analogy here to quantum optical networks which might be used for information processing either employing the polarization degree of freedom or, as we did here, the excitation number degree of freedom. This suggests that one could construct similar 'hardwired' networks on the level of interacting quantum systems that could then perform certain quantum information processing or communication tasks. This might involve structures such as the Y-shape presented here but may also implement structures such as interferometer structures shown in figure 14 or multi-input devices. If one were to consider hardwired structures, then it would be necessary to devise methods by which these structures could be switched on and off. Here we explore two possibilities. First, one might change the coupling strength c Junction of the oscillator at the three-way junction in the Y-shape. Apart from the obvious fact that they remain disentangled for c Junction = 0 (i.e. uncoupled), we find the first maximum for the entanglement decreases to roughly half the value for large coupling strength c Junction = 0.8. This is shown in figure 16. A further increase of the coupling strength to c Junction = 5 does not lead to further significant change. A different approach would be to change the eigenfrequency or the mass of the junction oscillator while keeping the coupling strength the same as the other oscillators. Indeed, if we increase the eigenfrequency ω Junction , the entanglement can be reduced to an arbitrarily small amount for both RWA and spring interactions. A decrease of the eigenfrequency is less efficient but would also allow a significant reduction of the amount of entanglement generated in the device. Figure 17 demonstrates these effects achieved by changing the eigenfrequency of the junction oscillator in theY-shape. It should be noted that the dependence of the logarithmic negativity with the eigenfrequency ω Junction is almost perfectly fitted by a Lorentzian line shape.
The above examples suggest that it is possible to switch on and off pre-fabricated devices such as the Y-shape shown above by adjusting either the coupling strength (decreasing it, i.e. approach decoupling) or the eigenfrequency (increasing it). Such a manipulation of the junction  By increasing and decreasing ω Junction respectively in the Y-shape we observe a noticeable change in the amount of entanglement that is generated in the device. The other parameters are chosen as in figure 13.
oscillator dictates the quantum information flowing through the junction. An example for a possible implementation of such a switch from an optical setting would be coupled optical cavities which are filled with atoms. Laser irradiation of these atoms would then lead to a shift in the resonance frequency of the cavity which corresponds to a change in the eigenfrequency of a harmonic oscillator in the above examples. In this way, individual cavities might be decoupled. A detailed study of such a scheme will be presented elsewhere.
Finally, we would like to briefly mention an analogy with the monotonic and non-monotonic behaviour of the efficiency. We find that by increasing the mass of the junction oscillator in the beamsplitter configuration, one can obtain monotonically decreasing entanglement between the two ends for the RWA interaction, whereas the Hooke law interaction produced non-monotonic behaviour.

Conclusions and discussions
We have investigated the entanglement dynamics of systems of harmonic oscillators both analytically and numerically. Particular attention has been paid to harmonic oscillators coupled by springs (Spring) and to harmonic oscillators with a linear coupling in a RWA as is appropriate in a quantum optical setting. After an introduction to the mathematical formalism and the derivation of the analytical solutions for the equations of motion for these interactions we then investigated several possible scenarios. We considered the generation of entanglement without detailed local control of individual systems. This was achieved by first switching off any interaction between the oscillators, cooling them to near the ground state and subsequently switching on the coupling suddenly. Surprisingly, entanglement will be generated over very large distances which is in stark contrast with the entanglement properties of the stationary ground state of a harmonic chain where only nearest neighbours exhibit entanglement [12]. We have also demonstrated that a linear chain of harmonic oscillators is capable of transporting quantum information and quantum entanglement for various types of nearest-neighbour coupling. For position-independent nearestneighbour coupling we observe that the transmission efficiency is a non-monotonic function in the coupling strength for Hooke's law coupling, whereas it is monotonically decreasing for the RWA coupling. In both cases, this suggests that it is advantageous to transmit entanglement in smaller portions rather than large units. However owing to the rapid decline in efficiency with the spring interaction for very small r, one should avoid sending in r too small. The propagation speed for the quantum entanglement has been provided analytically. For the above effects, we have studied their sensitivity to random variations in the coupling between the oscillators and to finite temperatures.
Finally, we have proposed more complicated geometrical structures such as Y-shapes and interferometric set-ups that allow for the generation of entanglement in pre-fabricated structures without the need for changing any coupling constants. We have also shown that these structures may be switched on and off by changing the coupling of only a single-harmonic oscillator with its neighbours. This suggests the possibility for the creation of pre-fabricated structures that may be 'programmed' by external actions. Therefore, quantum information would be manipulated through its propagation in these pre-fabricated structures somewhat analogous to modern microchips and as opposed the most presently suggested implementations of quantum information processing where stationary quantum bits are manipulated by a sequence of external interventions such as laser pulses.
All these investigations were deliberately left at a device-independent level. It should nevertheless be noted that there are many possible realizations of the above phenomena. These include nanomechanical oscillators [11], arrays of coupled atom-cavity systems, photonic crystals, and many other realizations of weakly coupled harmonic systems, potentially even vibrational modes of molecules in molecular quantum computing [35]. A forthcoming publication will discuss device-specific issues of such realization as well as improved structures (including novel topological structures as well as changes of their internal structure such as diatomic chains) that allow for better performances with fewer experimental resources. We hope that these ideas may lead to the development of novel ways for the implementation of quantum information processing in which the quantum information is manipulated by flowing through pre-fabricated circuits that can be manipulated from outside.