Temperature effects on a network of dissipative quantum harmonic oscillators: collective damping and diffusion processes

# Temperature effects on a network of dissipative quantum harmonic oscillators: collective damping and diffusion processes

M. A. de Ponte, S. S. Mizrahi, and M. H. Y. Moussa. Departamento de Física, Universidade Federal de São Carlos, Caixa Postal 676, São Carlos, 13565-905, São Paulo, Brazil Instituto de Física de São Carlos, Universidade de São Paulo, Caixa Postal 369, 13560-590 São Carlos, SP, Brazil
###### Abstract

In this article we extend the results presented in Ref. [Phys. Rev. A 76, 032101 (2007)] to treat quantitatively the effects of reservoirs at finite temperature in a bosonic dissipative network: a chain of coupled harmonic oscillators whichever its topology, i.e., whichever the way the oscillators are coupled together, the strength of their couplings and their natural frequencies. Starting with the case where distinct reservoirs are considered, each one coupled to a corresponding oscillator, we also analyze the case where a common reservoir is assigned to the whole network. Master equations are derived for both situations and both regimes of weak and strong coupling strengths between the network oscillators. Solutions of these master equations are presented through the normal ordered characteristic function. We also present a technique to estimate the decoherence time of network states by computing separately the effects of diffusion and the attenuation of the interference terms of the Wigner function. A detailed analysis of the diffusion mechanism is also presented through the evolution of the Wigner function. The interesting collective diffusion effects are discussed and applied to the analysis of decoherence of a class of network states. Finally, the entropy and the entanglement of a pure bipartite system are discussed.

###### pacs:
PACS numbers: 03.65.Yz; 05.10.Gg; 05.40.-a

## I Introduction

The subject of networks of interacting quantum systems has acquired an important role in the emerging field of quantum information theory. Since a realistic quantum logic processor must ultimately be composed of a large number of interacting quantum systems, it becomes mandatory to understand processes as like as perfect state transfer from one to another system of the network, and even to compute the fidelity of such a state transfer when the action of the environment is taken into account. A significant amount of result has recently been derived on the subject of perfect state transfer in optical lattices Feder () and networks of spin Spin () and harmonic oscillators Plenio (). Perfect state transfer has also been considered in networks of arbitrary topology and coupling configuration Topology () and even under random fluctuations in the couplings of a quantum chains Bose ().

Apart from state transfer, the process of decoherence of a network state has also attracted attention and interesting properties of collective damping effects, as the nonadditivity of decoherence rates, have been discussed in different contexts as in superconducting qubits Brito (), two-atom systems Ficek (), and chains of dissipative harmonic oscillators Mickel1 (); Mickel2 (); MickelDFS (); MickelGeral (). Still regarding the process of collective decoherence, the emergence of decoherence-free subspaces (DFSs) has also instigated several interesting results when considering the particular case of a composite system interacting with a common reservoir ZR (), or the more realistic situation where each system interacts with its own reservoir MickelDFS (). We call the attention to the fact that all Refs. Mickel1 (); Mickel2 (); MickelDFS () envisage such realistic cases of networks where each oscillator interacts with its own reservoir, also addressing the particular case where a common reservoir is considered.

To better understand the results in Refs. Mickel1 (); Mickel2 (); MickelDFS (), which are crucial to introduce the subject of the present work, we remember that, apart from the distinct reservoirs, the network of dissipative harmonic oscillators could present direct an indirect dissipative channels. Through the direct channels each oscillator loses excitation to its own reservoir, whereas through the indirect channels it loses excitation to all the other reservoirs but to its own. When we consider distinct reservoirs for each network oscillator, the indirect dissipative channels — intrinsically associated with the nonadditivity of decoherence rates and the emergence of DFSs MickelDFS () — are significant only in the strong coupling regime where , i.e., the number of network oscillators multiplied by their coupling strengths are about their natural frequencies . Therefore, the strong coupling regime, which brings together the collective damping effects, depends on the number of network oscillators as much as on their coupling strengths. For Markovian white noise reservoirs, however, where the spectral densities of the reservoirs are invariant over translation in frequency space, the indirect channels becomes null, except for the case Mickel1 ().

In the weak coupling regime where , the indirect channels always disappears. However, these indirect channels, coming from the strong coupling regime, remains in the case where all network oscillators interacts with a common reservoir MickelDFS (), even assuming a common Markovian white noise reservoir. This is due to the fact that a common reservoir induces an additional correlation between the network oscillators, restoring the indirect decay channels.

Recently, a generalization of Refs. Mickel1 (); Mickel2 () has been presented through a comprehensive treatment of networks of dissipative quantum harmonic oscillators, whichever its topology, i.e., whichever the way the oscillators are coupled together, the strength of their couplings and their natural frequencies MickelGeral (). Focusing on the general more realistic scenario where each oscillator is coupled to its own reservoir, the case where all the network oscillators are coupled to a common reservoir was also addressed. However, after deducing the master equation for the case where all the reservoirs are at finite temperatures, all further analysis of the dynamics of the network states was restricted to the case where the reservoirs are at K. Whereas a quantitative analysis of the decoherence and the evolution of the linear entropy of representative states of the network were given at K, only a brief qualitative analysis of the equilibrium states of the network was presented at finite temperatures. In the present manuscript we extend the treatment in Ref. MickelGeral () given a detailed analysis of the temperature effects on networks of dissipative quantum harmonic oscillators.

The present extension of Ref. MickelGeral () that accounts for the temperature effects coming from thermal reservoirs is not only interesting due to its more realistic approach but also from the mathematical development here achieved. In fact, we present an alternative approach to previous results in the literature Gardiner () regarding the obtainment of the solution of the master equation and the estimation of decoherence times through the Wigner distribution function. To circumvent noise effects, many of the nowadays experiments demonstrating quantum logic operations through atom-field interactions occur in cryogenic environments where temperature effects are negligible. In cavity quantum electrodynamics, the setup is cooled to around 0.5 K by a He-He refrigerator to avoid blackbody radiation in the High-Q superconducting cavity. Under such a specific condition, the temperature effects on the decoherence process are almost negligible. However, when the setup is scaled from one single cavity to a network of High-Q cavities, major questions arise due to temperature effects. First of all, would the DFSs survive despite the temperature effects? Apart from the special class of states composing the DFSs, how the temperature affects other states of the network, as for example initial entangled states? Evidently, these questions present no obvious answers, even under the assumption that all network cavities are cooled at low temperatures. On this regard, we expect that the collective damping effects coming from the indirect dissipative channels to play a major role for the answers to the above questions.

Apart from providing the mathematical treatment of the temperature effects on a network of dissipative harmonic oscillators, in the present manuscript we also analyze the role played by temperature on the evolution of particular states of the network other than those composing DFSs. We reserve the analyses of the emergence of DFSs under temperature effects to an specific work Mickel4 () where the mechanism of construction of such privileged states is also discussed along with decoherence.

This paper is organized as following: in Section II we revisit our model of a bosonic dissipative network MickelGeral () and present the derivation of the master equation governing the dynamics of the associated density operator. In Section III we present the solution of the normal ordered characteristic equation obtained from the master equation for the density operator of the network. In Section IV we analyze the evolution of two general classes of initial states of the network, given by mixtures of coherent and number states, through the normal ordered characteristic equation, the Glauber-Sudarshan -function, and the Wigner distribution. A detailed analyses of the diffusion processes is presented in Section V and the collective decoherence rates of a family of states of the network is analyzed in Section VI. In Section VII we discuss the entropy and the entanglement degree of a pure bipartite system and, finally, in Section VIII we present our concluding remarks.

## Ii The master equation of a bosonic dissipative network

We present here a brief review of the steps for the derivation of the master equation of a bosonic network, as developed in Ref. MickelGeral (). We start from the general case of a network of interacting oscillators, where each one interacts with each other, from which all other topologies can be recovered. As depicted in Fig. 1, we also consider the case where each oscillator interacts with its own reservoir due to this more realistic approach for most of the physical systems. However, as pointed in Ref. MickelDFS (), despite the realistic scenario of the case of distinct reservoirs, the case of a common reservoir is more general from the technical point of view. In fact, as discussed at the end of this section, the master equation for the case of distinct reservoirs can be deduced from the case of a common reservoir.

We start by considering a general Hamiltonian for a bosonic network, , composed by a set of coupled oscillators

 HS=ℏN∑m=1⎡⎣ωma†mam+12N∑n(≠m)=1λmn(a†man+ama†n)⎤⎦, (1)

distinct reservoirs, modeled by a set of modes,

 HR=ℏN∑m=1∑kωmkb†mkbmk, (2)

and the coupling between the network oscillators and their respective reservoirs

 HI=ℏN∑m=1∑kVmk(b†mkam+bmka†m), (3)

where () is the creation (annihilation) operator for the th bath mode coupled to the th network oscillator whose creation (annihilation) operator reads (). The coupling strengths between the oscillators are given by the set , while those between the oscillators and their reservoirs by . We assume, from here on, that and run from to .

Before addressing the dissipative process through Hamiltonian (3), we focus first on Hamiltonian to show how to derive different topologies of a nondissipative network of coupled harmonic oscillators. Rewriting in a matrix form

 HS=ℏ(a†1⋯a†N)⎛⎜ ⎜ ⎜⎝H11⋯H1N⋮⋱⋮HN1⋯HNN⎞⎟ ⎟ ⎟⎠⎛⎜ ⎜⎝a1⋮aN⎞⎟ ⎟⎠, (4)

we identify the elements of the matrix as

 Hmn={ωmform=nλmnform≠n, (5)

whose values characterize whichever the network topology, i.e., whichever the way the oscillators are coupled together, the set of coupling strengths , and their natural frequencies .

To obtain the master equation of the network we first diagonalize the Hamiltonian (within the physical regime where the normal modes are assumed to be positive) through the canonical transformation

 Am=∑nCmnan, (6)

where the coefficients of the th line of matrix define the eigenvectors associated to the eigenvalues of matrix . With being an orthogonal matrix, its transposed turns out to be exactly its inverse , resulting in the commutation relations and , which enable the Hamiltonian to be rewritten as a sum , where

 H0 =ℏ∑mϖmA†mAm+ℏ∑m∑kωmkb†mkbmk, (7a) V =ℏ∑m,n∑kC−1mnVmk(b†mkAn+bmkA†n). (7b) With the diagonalized Hamiltonian H0 we are ready to introduce the interaction picture, defined by the transformation U0(t)=exp(−iH0t/ℏ), where
 V(t)=ℏ∑m,n[Omn(t)A†n+O†mn(t)An], (8)

and . Next, we assume that the interactions between the resonators and the reservoirs are weak enough to allow a second-order perturbation approximation. We also assume a Markovian reservoir such that the density operator of the global system can be factorized as . Under these assumptions the reduced density operator of the network of dissipative coupled resonators satisfy the differential equation

 dρS(t)dt=−1ℏ2∫t0dτTrR[V(t),[V(τ),ρS(t)⊗ρR(0)]]. (9)

Since for a thermal reservoir , we have to solve the integrals appearing in Eq. (9), related to correlation functions of the form

 ∫t0dτ⟨Omn(t)O†m′ℓ(τ)⟩ =C−1mnCℓm′∫t0dτ∑k,k′VmkVm′k′⟨bmkb†m′k′⟩ ×exp{−i[(ωmk−ϖn)t−(ωm′k′−ϖℓ)τ]}. (10)

Considering that the reservoir frequencies are very closely spaced to allow a continuum summation, we obtain

 ∫t0dτ⟨Omn(t)O†m′ℓ(τ)⟩=Nδmm′C−1mnCℓmγm(ϖℓ)+¯nm(ϖℓ)~γm(ϖℓ)2ei(ϖℓ−ϖn)t, (11)

where we have defined the average excitation of the reservoir associated to the th oscillator as through the relation , apart from the damping rates

 γm(ω) =∫t0dτ∫∞0dνNπ[Vm(ν)σm(ν)]2e−i(ν−ω)(τ−t), (12a) ~γm(ω) =∫t0dτ∫∞0dνNπ[Vm(ν)σm(ν)]2¯nm(ν)¯nm(ω)e−i(ν−ω)(τ−t), (12b) with σm(ν) being the density of states of the mth reservoir. In the context of Markov approximation, where Vm(ϖn),σm(ϖn) and ¯nm(ϖn) are slowly varying functions around the normal modes ϖn we can simplify the expressions (12) to their usual forms
 γm(ω)=~γm(ω)=1N[Vm(ω)σm(ω)]2. (13)

Back to the Schrödinger picture and to the original field operators , we finally obtain from the steps outlined above, the master equation

 dρS(t)dt=iℏ[ρS(t),HS]+∑m,n[Γmn+Υmn2LmnρS(t)+Υmn2LmnρS(t)], (14)

where we have defined the damping and the diffusion matrix elements and in the forms

 Γmn =N∑ℓCℓnγm(ϖℓ)C−1mℓ, (15a) Υmn =N∑ℓCℓn~γm(ϖℓ)¯nm(ϖℓ)C−1mℓ, (15b) whereas the Liouville operators accounting for the direct (m=n) and indirect (m≠n) dissipative channels, are given by
 LmnρS(t) ≡[anρS(t),a†m]+[am,ρS(t)a†n], (16a) LmnρS(t) ≡[a†nρS(t),am]+[a†m,ρS(t)an]. (16b) As mentioned in the Introduction and discussed in Ref. MickelGeral (), the oscillators lose excitation to their own reservoirs through the direct dissipative channels, whereas through the indirect channels they lose excitation to all the other reservoirs but not to their own. Although in Ref. MickelGeral () we have obtained the master equation for the general case of reservoirs at finite temperatures, all further analysis was carried out for reservoirs at 0 K where the diffusion matrix elements Υmn are null. Next, considering the case of reservoirs at finite temperatures, we must discuss the master equation (14) under the weak and strong coupling regimes between the network oscillators.
 Γmn =Nγm(ωm)δmn, (17a) Υmn =N~γm(ωm)¯nm(ωm)δmn, (17b) where, evidently, we have also approximated the normal modes by their original natural frequencies. Under the above considerations, the master equation (14) becomes
 dρS(t)dt =iℏ[ρS(t),HS]+N∑m[γm(ωm)+~γm(ωm)¯nm(ωm)2LmmρS(t) +~γm(ωm)¯nm(ωm)2LmmρS(t)], (18)

where, essentially, the indirect dissipative channels disappear, establishing the additivity of the decoherence rates.

### ii.2 Strong coupling regime

The strong coupling regime means that , i.e., at least one of the coupling between two network oscillators must be of the order of any natural frequency . In this case, the indirect dissipative channels become effective, inducing collective damping and diffusion effects that we must investigate. As pondered in the Introduction, would the collective effects of nonadditivity of the decay rates and the emergence of DFSs still survive despite the temperature effects?

It must be mentioned that Markovian white noise reservoirs washes out the collective damping effects introduced by the strong coupling regime since the spectral densities are invariant over translation in frequency space, i.e., , rendering the same matrix elements as in Eq. (17a). However, Markovian white noise reservoirs do not washes the collective diffusion effects. In fact, only under the additional assumption that for whatever normal mode , we recover Eq. (17b) for the diffusion matrix , erasing the collective effects completely. Next we discuss the case where the whole network is under the action of a common reservoir.

### ii.3 A common reservoir

When all the network oscillators are coupled to a single reservoir, the master equation, derived in Ref. MickelDFS (), is similar to that in Eq. (14), replacing the damping and the diffusion matrix elements by

 Γmn =N∑ℓ,n′Cℓnγmn′(ϖℓ)C−1n′ℓ, (19a) Υmn =N∑ℓ,n′Cℓn~γmn′(ϖℓ)C−1n′ℓ, (19b) where the damping rates γmn(ω) and ~γmn(ω) for the case of a single common reservoir are given by MickelDFS ()
 γmn(ω) =∫t0dτ∫∞0dνNπVm(ν)Vn(ν)σ2(ν)e−i(ν−ω)τ, (20a) ~γmn(ω) =∫t0dτ∫∞0dνNπVm(ν)Vn(ν)σ2(ν)¯n(ν)¯n(ω)e−i(ν−ω)τ. (20b) As mentioned above and discussed in Ref. MickelDFS (), the master equation for the case of distinct reservoirs can be deduced from the case of a single common reservoir. The above deduction of the master equation (14), where we started from the case of distinct reservoirs, was entirely due to its broad application in many physical systems. To demonstrate how to derive the case of distinct reservoir from that of a common one, we remember that Vm(ν) gives the distribution function of the reservoir modes coupled to the mth oscillator. Therefore, in the absence of overlap between the distribution functions, i.e., ∫dνVm(ν)Vn(ν)=0 for m≠n, Eqs. (20) reduce to those in Eqs. (12). In this case, the occurrence of the indirect-decay channels follows entirely from the strong coupling between the oscillators, as discussed in the subsections presented above. When there is a significant overlap between the distribution functions, i.e., ∫dνVm(ν)Vn(ν)≠ 0 for at least one m≠n, we get the indirect-decay channels even when the network oscillators do not interacts at all. The strength of the damping and the diffusion matrix elements being defined by the amount of the overlap, i.e., when the overlap between the distributions Vm(ν) and Vn(ν) is maximum, the strengths Γmn and Υmn equals Γmm and Υmm.
 ddtχ({ηm},t)=−∑m,n[ηmΥmn2η∗n+ηm(HDmn)∗∂∂ηn+c.c.]χ({ηm},t), (21)

where we defined the matrix elements

 HDmn=Γmn/2+iHmn. (22)

As noted in Ref. MickelGeral (), the matrix is an extension of the free evolution in Eq. (5), which takes into account the dissipative mechanisms of the network.

Starting with the assumption that Eq. (21) admits a solution of the form , we obtain two differential equations, one accounting for the dynamic process, given by

 ddtϕ({ηm},t)=−∑m,n[ηm(HDmn)∗∂∂ηn+c.c.]ϕ({ηm},t), (23)

and the other accounting for the stationary solution of the characteristic function, given by

 (24)

If we perform the substitution in the first differential equation (23), it turns out to be exactly that appearing in Ref. MickelGeral () for the derivation of the solution of the Glauber-Sudarshan -function. Therefore, following the steps outlined in Ref. MickelGeral (), the solution of Eq. (23) can be written as

 ηm(t)=∑ℓ,nηn(0)D∗nℓexp(−Ω∗ℓt)(D−1ℓm)∗, (25)

where we employed the diagonal form of following from the transformation . Note that writing the solution (25) in a matrix form, it becomes

 η(t) =η(0)∙D∗∙exp(−Ω∗t)∙(D−1)∗ =η(0)∙exp[−(D∙Ω∙D−1)∗t] =η(0)∙exp[−(HD)∗t] (26)

such that

 dη(t)dt=−η(t)∙(HD)∗, (27)

or, equivalently,

 dηn(t)dt=−∑mηm(t)(HDmn)∗, (28)

representing a system of coupled differential equations which follows from Eq. (23) under the assumption that , with MickelGeral ().

The second differential equation (24) can be solved assuming a general Gaussian form

 φ({ηm})=exp(−12∑m,nηmΠmnη∗n), (29)

where the elements of matrix are the coefficients to be determined. Substituting (29) into Eq.(24) and changing conveniently the labels and of the involved matrices, we verify that the differential equation (24) reduces to a matrix equation of the form

 (HD)∗∙Π+Π∙(HD)⊤=Υ+Υ⊤, (30)

which is explicitly written as

 ∑ℓ(HDmℓ)∗Πℓn+∑ℓΠmℓHDnℓ=Υmn+Υnm. (31)

with the superscript in Eq. (30) standing for transposed. It is worth noting that for identical reservoirs, where and so , we obtain a symmetric dissipative matrix , i.e., , making Eq. (30) the well-known Lyapunov equation. The solution of Eq. (30), namely the determination of , can be obtained by converting the matrix equation into a system of algebraic equations, i.e., into a new matrix equation of the simplified form , with the elements of matrix being the unknown variables. To this end, it is useful to define the column vector

 vec(Π)≡(Π11Π21⋯ΠN1Π12⋯ΠN2⋯Π1N⋯ΠNN)⊤, (32)

where the first elements of correspond to the first column of matrix , whereas the next elements correspond to the second column of and so on. As so, the matrix equation (30) can be rewritten into the form livroSalomon ()

 [I⊗(HD)∗+HD⊗I]∙vec(Π)=vec(Υ+Υ⊤)%, (33)

where is an identity matrix. From the mathematical properties presented in Appendix A for the matrix , we verify that the elements of matrix can be written as

 Πℓℓ′=∑m,n,m′,n′Υm′n′+Υn′m′Ωm+Ω∗nDℓ′mD−1mm′(DℓnD−1nn′)∗, (34)

finally leading to the solution of Eq. (24) through Eq.(29). In fact, substituting the expression (34) into the left hand side of Eq. (31), and using the relation , we obtain

 ∑ℓ[(HDmℓ)∗Πℓn+ΠmℓHDnℓ]=∑m′,n′(Υm′n′+Υn′m′)δmn′δnm′=Υmn+Υnm, (35)

which is exactly the right hand side of Eq. (31).

We thus verify that the solution of the characteristic equation is of the form

 χ({ηm},t)=φ({ηm})[ϕ({ηm},t=0)|{ηm}⇒{ηm(t)}. (36)

Since for we get , such that , we end up with the solution of the characteristic function

 χ({ηm},t)=φ({ηm})φ({ηm(t)})[χ({ηm},t=0)|{ηm}⇒{ηm(t)}, (37)

given in terms of its initial state. An interesting point to be noted is that the dynamics of the problem, given by , takes into account only the dissipative rates together with the free evolution Hamiltonian in Eq. (5), leaving aside the diffusive process associated with . Such a diffusive process, appearing in the ratio is however, modified, by the dissipative mechanisms.

## Iv Dynamics of the network states: characteristic function, Glauber-Sudarshan P-function, Wigner distribution and density operator

Starting from two general classes of initial network states, given by mixed superpositions of coherent and number states, we next analyze the evolution of such states through the characteristic function, the Glauber-Sudarshan -function, and the Wigner distribution. We also compute the network density operator for the case of mixed superposition of Fock states.

### iv.1 A mixed superposition of coherent states

Considering that the initial network state comprehends a mixture of superpositions of coherent states like , the initial density operator becomes

 (38)

where is the probability associated with the state . The parameters () represent a set of variables defining the probability density function , and stands for a product of coherent states, where represents the state associated with the th network oscillator. In the particular case where , the pure state becomes the discrete superposition

 (39)

where we have defined . Through the definition of the time-dependent vector elements and matrix elements

 Θmn(t) =∑ℓDmℓexp(−Ωℓt)D−1ℓn, (40a) Jmn(t) =Πmn−∑m′,n′Πm′n′Θ∗mm′(t)Θnn′(t), (40b) we verify, after a rather lengthy calculation, that the evolution of the initial network state (38) can be described either through the characteristic function
 χ({ηm},t) ×exp{∑m[ηmK∗m(rȷ;t)−η∗mKm(sȷ;t)]−12∑m,nηmJmn(t)η∗n}, (41)

either by the Glauber-Sudarshan -function

 P({ξm},t) =(2/π)NdetJ∑ȷpȷ∫dsȷΛ(sȷ)∫drȷΛ∗(rȷ)⟨{βm(rȷ)}∣∣{βm(sȷ)}⟩ ×exp{−2∑m,nJ−1mn(t)[ξm−Km(sȷ;t)][ξn−Kn(rȷ;t)]∗}, (42)

or even by the Wigner distribution function

 W({ξm},t) =(2/π)Ndet~J∑ȷpȷ∫dsȷΛ(sȷ)∫drȷΛ∗(rȷ)⟨{βm(rȷ)}∣∣{βm(sȷ)}⟩ ×exp{−2∑m,n~J−1mn(t)[ξm−Km(sȷ;t)][ξn−Kn(rȷ;t)]∗}. (43)

Note that the difference between the Glauber-Sudarshan -function and the Wigner distribution comes from the time-dependent function associated with the width of the their Gaussian function. Consequently, the Wigner function can be obtained from the Glauber-Sudarshan -function through the substitution . Whereas the Glauber-Sudarshan -function diverges when there is no diffusion process such that (with all the reservoirs at K), the width of the Wigner function presents an additional term inhibiting any singularity.

For the case of K reservoirs MickelGeral (), the density operator of the network, to be used below, is given by

 ρS(t) =∑ȷpȷ∫dsȷΛ(sȷ)∫drȷΛ∗(rȷ) (44)

### iv.2 A mixed superposition of Fock states

We now assume the initial network state to be a mixture of superposition of Fock states , where the parameter indicates the number of photons in the th oscillator while the coefficient represents the probability amplitude associated with each state composing the whole superposition. The initial density operator is thus given by

 (45)

where is the probability associated with the state . Since the Fock state , of the th oscillator, can be expanded as a superposition of coherent states of the form

 |xm⟩=Nm∫2π0dθme−ixmθm∣∣βmeiθm⟩, (46)

it is easy to note that the initial state (45) can be obtained by Eq. (38), identifying

 Λ(rȷ) →Λ(θȷ)=∑{xm}C(ȷ){xm}∏mNme−ixmθm, (47a) ∫drȷ →∫2π0dθm, (47b) βm(rȷ) →βmeiθm%, (47c) such that we can use the results of the previous subsection to obtain the characteristic function, the Glauber-Sudarshan P-function and Wigner distribution for a mixed superposition of pure Fock states. Alternatively, such functions may be directly computed from the initial state (45). Their expressions are presented in Appendix B, where the density operator for a mixed superposition of pure Fock states is also presented.

where the unitary operation satisfies . From this matrix relation, we obtain the evolved diffusion coefficients

 Dm(t)=∑n,n′U†mn(t)~Jnn′(t)Un′m(t) (49)

as the elements of diagonal matrix . In this framework, the rotated Wigner distribution, written as

 W({~ξm},t)=∑ȷpȷ∫dsȷ∫drȷW({~ξm};rȷ,sȷ,t), (50)

are composed by diagonal () and off-diagonal () elements defined by

 W({~ξm};rȷ,sȷ,t) =(2/π)Ndet~JΛ∗(rȷ)Λ(sȷ)⟨{βm(rȷ)}∣∣{βm(sȷ)}⟩ ×exp{−∑m2Dm(t)[~ξm−~Km(sȷ;t)][~ξm−~Km(rȷ;t)]∗}, (51)

where and . The vector gives the excitation intensity of the th normal-mode oscillator through . We stress that the larger or smaller values of depend on the network topology (contained within the matrix elements ), apart from the regime of coupling strengths between the oscillators (contained within the matrix elements ). As a particular example of this dependence, we consider a degenerate symmetric network, i.e., a degenerate network of oscillators, all of them interacting with each other, where , and . In this case, in the strong coupling regime, we obtain the expression

 ~Jmn(t)=δmn+2¯nN(1−e−Nγt), (52)

and the diffusion coefficients

 Dm(t)=⎧⎨⎩~Jmm(t)−~Jmn(t)=1form≠N,~Jmm(t)+(N−1)~Jmn(t)=1+2¯n(1−e−Nγt)form=N, (53)

showing that can reduce or enhance the strength of the diffusion coefficients associated with the normal-mode oscillators.

### v.1 Directional and mean diffusion times

From the above TD diffusion coefficients (53) we define the directional diffusion time

 1τ(m)diff=ddtDm(t)∣∣∣t=0, (54)

displaying a tendency to a significant spread of the peak — common to all elements (the diagonal and off-diagonal) of the Wigner function — associated with the th normal-mode oscillator. Since each normal-mode oscillator defines a direction in the coordinate frame , we are naturally led to define the mean diffusion time, associated with all the dimensions of the space, as the average value

 1τdiff=1N∑m1τ(m)diff=1NddtTrD(t)∣∣∣t=0. (55)

The average diffusion time becomes useful to compute the decoherence time of any network state when complemented with the estimated time for a significant decay of the peaks associated with the interference terms of the Wigner function (), to be defined below as ,

As an illustrative example of the above theory, below we analyze the diffusion coefficients for the weak and strong coupling regimes considering the case of a degenerate symmetric network.

#### v.1.1 The weak coupling regime

In the weak coupling regime, the matrix , already in a diagonal form, is defined by the elements , such that . In this regime, all the diffusion coefficients equal to

 Dm(t)=D(t)=1+2¯n(1−e−γt). (56)

The average diffusion time becomes

 τdiff=12¯nγ, (57)

showing, as expected, that the larger the temperature, the smaller the time required for a significant diffusion rate. In this case, the coefficients are mode independent and assume a common value, such that the spreads of the peaks associated with the diagonal terms of the Wigner function occurs homogeneously in all directions.

#### v.1.2 The strong coupling regime

In the strong coupling regime, the elements of matrix are given by Eq. (52) and the diffusion coefficients by Eq. (53), showing that only the th normal-mode oscillator undergoes the diffusion process. For all the normal-mode oscillators but the th, the diffusion coefficients are counterbalanced by the diffusion rates and coming from the direct- and indirect-decay channels, respectively. The diffusion coefficients in this regime lead to the same mean diffusion time as that in Eq. (57), showing that the average diffusion effect comes entirely from the temperatures of the reservoirs. As to be demonstrated in the next section, this interesting result is not limited to the degenerate symmetric topology.

### v.2 Diffusion and topology

Starting from Eq. (55) and noting that (), with the elements of matrix given by Eq. (40b), we obtain the general expression

 τ−1diff=2NTrΥ, (58)

applicable to whatever the network topology and the strength coupling regime between the oscillators, where

 TrΥ=N∑m,n~γm(ϖn)¯nm(ϖn)CnmC−1mn. (59)

We note that the information regarding the topology of the network is contained only in the product which acts as a normalized distribution function () when computing the average value of the diffusion rate given by Eq. (58).

We identify two general situations where, as in the case of a degenerate symmetric network, the diffusion mechanism becomes independent of the topology of the network. The first situation occurs when identical reservoirs are assumed, such that and, consequently, , making the mean diffusion

 τ−1diff=2∑m~γ(ϖm)¯n(ϖm)%, (60)

independent of the network topology. The second situation arises from the assumptions of Markovian white noise reservoirs and low-temperature regime, where the normal-mode frequencies satisfy the relation