PAPER • OPEN ACCESS Recent citations Two-particle quantum correlations in - Quantum transport in non-Markovian stochastically-coupled networks dynamically disordered photonic lattices Ricardo Román-Ancheyta et al To cite this article: Roberto de J León-Montiel et al 2019 New J. Phys. 21 053041 View the article online for updates and enhancements. This content was downloaded from IP address 158.109.49.69 on 23/07/2021 at 16:00 New J. Phys. 21 (2019) 053041 https://doi.org/10.1088/1367-2630/ab1c79 PAPER Two-particle quantum correlations in stochastically-coupled OPEN ACCESS networks RECEIVED 9 January 2019 Roberto de J León-Montiel1,6 , Vicenç Méndez2,6 , Mario A Quiroz-Juárez1 , Adrian Ortega3, REVISED 25 March 2019 Luis Benet3 , Armando Perez-Leija4,5 and Kurt Busch4,5 1 Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Apartado Postal 70-543, 04510 Cd. Mx., México ACCEPTED FOR PUBLICATION 2 25 April 2019 Grup de Física Estadística, Departament de Física, Universitat Autònoma de Barcelona, E-08193 Barcelona, Spain 3 Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Apartado Postal 48-3, 62251 Cuernavaca, Morelos, México PUBLISHED 4 29 May 2019 Max-Born-Institut, Max-Born-Straße 2A, D-12489 Berlin, Germany 5 Humboldt-Universität zu Berlin, Institut für Physik, AG Theoretische Optik Photonik, Newtonstraße 15, D-12489 Berlin, Germany 6 These authors contributed equally to this work. Original content from this work may be used under E-mail: roberto.leon@nucleares.unam.mx the terms of the Creative Commons Attribution 3.0 Keywords: many-particle quantum correlations, quantum networks, off-diagonal dynamical disorder licence. Any further distribution of this work must maintain attribution to the Abstract author(s) and the title of the work, journal citation Quantum walks in dynamically-disordered networks have become an invaluable tool for under- and DOI. standing the physics of open quantum systems. Although much work has been carried out considering networks affected by diagonal disorder, it is of fundamental importance to study the effects of fluctuating couplings. This is particularly relevant in materials science models, where the interaction forces may change depending on the species of the atoms being linked. In this work, we make use of stochastic calculus to derive a master equation for the dynamics of one and two non-interacting correlated particles in tight-binding networks affected by off-diagonal dynamical disorder. We show that the presence of noise in the couplings of a quantum network creates a pure-dephasing-like process that destroys all coherences in the single-particle Hilbert subspace. Moreover, we show that when two or more correlated particles propagate in the network, coherences accounting for particle indistinguishability are robust against the impact of off-diagonal noise, thus showing that it is possible, in principle, to find specific conditions for which many indistinguishable particles can traverse stochastically-coupled networks without losing their ability to interfere. 1. Introduction The study of quantum random walks in noisy environments have played a fundamental role in understanding non-trivial quantum phenomena observed in an interdisciplinary framework of studies ranging from biology [1, 2], chemistry [3], materials science [4] and electronics [5], to photonics [6–9] and ultracold matter [10, 11]. For many years, most of the research efforts had been focused on the propagation of single particles [12]; however, a great interest in describing the dynamics of correlated particles in noisy systems has recently arisen [13–16], mainly because it has been recognized that many-particle quantum correlations can be preserved in noisy networks by properly controlling the initial state of the particles, their statistics, indistinguishability or their type of interaction [17, 18]. In general, the interesting features in the dynamics of quantum correlated particles traversing noisy networks are due to the tunneling amplitudes in the associated Hamiltonians. Therefore, including noise into the off-diagonal elements of the Hamiltonian allows one to assess the effects of decoherence and noise. On many occasions, when describing the evolution of correlated particles in network systems affected by non-dissipative noise, a physically accurate result can be obtained after averaging over many realizations of the noisy walks. In other words, in most cases, one does not have a master equation to analytically describe the phenomenon under study. Indeed, this represents a serious problem, specially in cases where the number of particles or network sites is extremely large. In such scenarios, computing the evolution of the system quickly becomes a computationally © 2019 The Author(s). Published by IOP Publishing Ltd on behalf of the Institute of Physics and Deutsche Physikalische Gesellschaft New J. Phys. 21 (2019) 053041 R de J León-Montiel et al demanding task, which can only be tackled by developing sophisticated computer algorithms [19]. Consequently, most of the work is generally focused on optimizing numerical approaches, and the physical interpretation of the noise effects are sometimes overlooked. In the present work we introduce an exact analytical approach to study quantum walks in noisy systems. We use stochastic calculus to derive a master equation for the propagation of two correlated particles in a quantum network affected by off-diagonal dynamical disorder. By using our results, we show that off-diagonal noise produces an effective pure-dephasing-like process that destroys all coherences in a single-particle quantum walk. Remarkably, we find that when two or more indistinguishable particles propagate in a noisy system, coherences accounting for particle indistinguishability are robust against the dephasing-like process. These results elucidate the role of particle indistinguishability in the preservation of quantum coherence in systems that interact with a noisy environment. 2. Single-particle dynamics We start by describing the dynamics of a single particle in a quantum network affected by random fluctuations in the coupling between sites. In this situation, the time evolution of the single-particle wavefunction at the nth site, ψn, is given by the stochastic Schrödinger equation (with ÿ=1) dyn = - iwn yn - i å knm (t ) ym, (1) dt m¹n where ωn stands for the energy of the nth site, and the coupling between them is given by knm (t ) = knm + fnm (t ), with fnm (t ) = fmn (t ) describing a white-noise process with zero average, that is, á fnm (t )ñ = 0, and á fnm (t ) fjl (t ¢)ñ = gnm dnm, jld (t - t ¢). Here dnm, jl = dnj dml + dnl dmj , with δnm being the Kronecker delta. γnm denotes the noise intensity, that is, how strong the stochastic fluctuations are, and áñ denotes averaging over the noise realizations. Following a treatment equivalent to the one used in [20, 21], where fluctuations are introduced in the site- energies rather than the couplings, we can obtain a master equation for a stochastically-coupled network by taking the time derivative of rnm (t ) = á yn y*mñ. Thus, by using equation (1), we can write drnm dy* dy = yn m + y*m n , dt dt dt = - i (wn - wm ) rnm + i å kmj rnj - i å knj rjm j j - i å gmj á yn y*j hmj (t )ñ + i å gnj á yj y*mhnj (t )ñ , (2) j j where we have defined a new stochastic variable hnm (t ) = -fnm (t ) gnm , which satisfies the conditions á hnm (t )ñ = 0, and á hnm (t ) hjl (t ¢)ñ = dnm, jld (t - t ¢). Notice that equation (2) is not yet complete, as it remains to compute the correlation functions of the last two terms. To do so, we employ the Novikov’s theorem [22], which for the fourth term on the right hand side of equation (2) takes the form d [yn (t ) y*j (t )] á yn y*j hmj (t )ñ = å pq ò dt ¢á hmj (t ) hpq (t ¢)ñ dhpq (t ¢) , 1 d [yn (t ) y*j (t )] = 2 å dmj,pq dhpq (t ) , (3) pq where the operator d dhpq (t ) stands for the functional derivative with respect to the stochastic process. At this point, it is worth remarking that the Novikov’s theorem can be used for solving different types of Gaussian stochastic processes, such as non-Markovian colored noise [23]. Now, to solve the functional derivative in equation (3), we make use of the expression t ⎡ yn (t ) y*m (t ) = ò0 dt ¢ ⎢ f (yn y*m,¼) - i å gmr yn y*r hmr (t ) ⎣ r ⎤ + i å gnr yr y*mhnr (t ) ⎥ , (4) r ⎦ which corresponds to the formal integration of equation (2) without the noise average. Note that, in equation (4), the function f (yn y*m,¼) contains all terms that do not depend on stochastic variables. 2 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al Thus, by using equation (4) we obtain d [yn (t ) y*j (t )] = - i å gjr yn y*r djr , pq + i å gnr yr y*j dnr , pq , (5) dhpq (t ) r r where we have used of the relation dhjr dhpq = djr , pq . We can now substitute equation (5) into (3) to find i i á yn y*j hmj (t )ñ = - 2 å gjr rnr djr , mj + 2 å gnr rrj dnr , mj . (6) r r Similarly, the fifth term on the right hand side of equation (2) is found to be i i á yj y*mhnj (t )ñ = - 2 å gmr rjr dmr , nj + 2 å gjr rrm djr , nj . (7) r r Finally, by substituting equations (6) and (7) into (2), we obtain drnm ⎡ ⎤ 1 = - ⎢i ( w n - w m ) + å (gnj + gmj ) ⎥⎥ rnm dt ⎢⎣ 2 j ⎦ + i å (kmj rnj - knj rjm) + gnm rnm + dnm å gnj gmj rjj , (8) j j which corresponds to a master equation for the time evolution of a single particle in a stochastically-coupled quantum network. Note that equation (8) can be cast into a Lindblad equation [12] by noticing that the deterministic part of equation (1) corresponds to the dynamics of a quantum network described by a tight- binding Hamiltonian of the form Hd = å wn∣nñá n ∣ + å knm∣ nñá m ∣ , (9) n m¹n with ∣nñ denoting the particle (or excitation) being at the nth site of the network. Then, by making use of equation (9), we can write equation (8) as drnm = - i [Hd, r ]nm + deph [r ]nm , (10) dt where [,] stands for the commutator. Notice that the effect of the fluctuating coupling is captured by the pure- dephasing-like operator ⎡ ⎤ 1 deph [r ]nm = - ⎢ å (gnj + gmj ) - gnm ⎥⎥ rnm + dnm å gnj gmj rjj , (11) ⎢⎣ 2 j ⎦ j with γnm describing the dephasing rate introduced in the coupling between the nth and the mth site. To elucidate the effect of the stochastic coupling between sites, we now compute the dynamics of a single excitation in a fully connected network composed by three sites. We have chosen this configuration, because it constitutes the simplest quantum network that one can investigate both theoretically and experimentally [5, 24]. The energies of the sites are arbitrarily chosen to be ω1=ω2=ω3=5 ps−1, whereas the coupling between them are set to κ12=1 ps−1, and κ13=κ23=0.5 ps−1. Figure 1 shows some examples of platforms where single-excitation stochastic networks have been successfully implemented, namely optical tweezers [25], waveguide arrays [24], superconducting circuits [26–28], and electrical-circuit arrays [29]. Quite recently, it has been shown that stochastic quantum networks can also be implemented in fascinating platforms, such as nuclear magnetic resonance systems [30], and ion traps [31, 32]. The time evolution of the diagonal (populations) and off-diagonal (coherences) elements of the system’s density matrix, solved by means of equation (8), is shown in figure 2. In all figures, the dephasing rate is set to γ12=γ13=γ23=0.38 ps−1. For the sake of comparison, we have included the numerical solution (dashed lines) of equation (1), which corresponds to the average of 10 000 random realizations, where the dephasing coefficient is defined by means of the relation [33, 34]: gnm = s nm 2 Dt , with s nm 2 being the variance of the Gaussian distribution containing the values of the stochastic variable fnm (t ), and Δt the correlation time. Notice that the effect of the fluctuating coupling is a pure-dephasing-like process, described by the operator deph [r ], that destroys the coherence between the network’s sites. This phenomenon is equivalent to the quantum Zeno effect [35], a process in which the system is driven by its environment into a steady state where the regular hopping of the wavefunction is no longer sustained, i.e. a state in which the initial excitation becomes incoherently delocalized [36, 37]. 3 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al Figure 1. Schematic representation of photonic and electronic platforms where single-excitation stochastic networks have been investigated: (a) optical tweezers, (b) waveguides, (c) superconducting circuits, and (d) electrical-circuit arrays. Figure 2. Dynamics of a single excitation injected into site 1 of a stochastically-coupled three-site quantum network. (a) Time evolution of the population in each of the sites; (b) and (c) show the real and imaginary parts of the coherence (off-diagonal) terms, respectively. Note that the effect of the fluctuating coupling is a non-dissipative pure-dephasing-like process that destroys coherence between the network’s sites. The solid line corresponds to the solution using our derived master equation (equation (8)); whereas the dashed line shows the numerical solution of equation (1) obtained by averaging 10 000 realizations. In both cases, we have set the dephasing rates to γ12=γ13=γ23=0.38 ps−1. 4 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al 3. Two-particle wavefunction dynamics We now turn our attention to the description of two-particle correlation dynamics. To this end, we use the concept of two-particle probability amplitude [24, 38], and derive the corresponding equations of motion for finite tight-binding networks comprising N sites. We start by noting that the probability amplitudes for a quantum particle, initialized at a site n, are governed dUp, n by the equations [24, 38]: dt = -iwn Up, n - iå rN= 1 kpr (t ) Ur , n, where Up,n stands for the impulse response of the system, that is, the unitary probability amplitude for a single particle traveling from site n to site p. As in the previous section, the coupling kpr (t ) represents a Gaussian Markov process with zero average. We can then write, in terms of single-particle probability amplitudes, the two-particle probability amplitudes at sites p and q as: yp, q (t ) = åm = 1, n = 1xm, n [Up, n (t ) Uq, m (t )  Up, m (t ) Uq, n (t )], where xm, n is the initial probability amplitude profile that fulfills the conditions åm = 1, n = 1∣ xm, n ∣2 = 1. Notice that the sign±determines whether the particles are bosons (+) or fermions (−), respectively. Then, by taking the time derivative of the two-particle wavefunction, we obtain the equation dy p , q = - i (wp + wq ) yp, q - i å [kpr (t ) yr , q + kqr (t ) yp, r ] , (12) dt r which describes the dynamics of two-particle quantum correlations. Notice that two-particle quantum states evolve in a Hilbert space composed by a discrete set of N2-mode states occupied by the two particles. One important fact to highlight regarding equation (12) is the presence of the term (wp + wq ) yp, q , which implies that during evolution the wavefunction ψp,q acquires a phase that a single particle acquires when it traverses the same network twice [39]. Indeed, such effects can be expected since we are dealing with two correlated particles [40]. Finally, we remark that the modulus squared of the two-particle wavefunction gives the probability of finding one particle at site p and the other at q [41–46]. We can now follow the same procedure as in the previous section to obtain a master equation for the two- particle wavefunction dynamics by taking the time derivative of rpq, p ¢ q ¢ = á ypq y*p ¢ q ¢ñ. Thus, by using equation (12), we obtain (see appendix A for details) drpq, p ¢ q ¢ ⎡ = ⎢ - i (wp + wq - wp ¢ - wq ¢) - gpq - gp ¢ q ¢ dt ⎣ 1 ⎤ - å (glp + glq + glp ¢ + glq ¢ ) ⎥ rpq, p ¢ q ¢ 2 l ⎦ - i å (klq rpl, p ¢ q ¢ + klp rlq, p ¢ q ¢) l + i å (klq ¢ rpq, p ¢ l + klp ¢ rpq, lq ¢) l - å (dpq glq glp rll, p ¢ q ¢ + dp ¢ q ¢ glp ¢ glq ¢ rpq, ll ) l + å (dqq ¢ glq glq ¢ rpl, p ¢ l + dqp ¢ glq glp ¢ rpl, lq ¢) l + å (dpq ¢ glp glq ¢ rlq, p ¢ l + dpp ¢ glp glp ¢ rlq, lq ¢) l + gqq ¢ rpq ¢ , p ¢ q + gqp ¢ rpp ¢ , qq ¢ + gpp ¢ rp ¢ q, pq ¢ + gpq ¢ rq ¢ q, p ¢ p, (13) which is the master equation that describes the time evolution of two correlated particles in a stochastically- coupled quantum network. Note that we can write equation (13) in its equivalent Lindblad form by following the procedure described above and noticing that the Hamiltonian for the two-particle system is given by  = Hd ⨂ + ⨂Hd , where ⨂ stands for the tensor product and  is the identity matrix, whose dimension is the same as the single-particle Hamiltonian. Before considering particular examples, it is worth noting that in the following we will use the compact notation ∣1n , 1m ñ to represent the states where one particle is populating the site n and another the site m, i.e. ∣ 1n ñ Ä ∣ 1m ñ, whereas states µ(∣ 1n , 1m ñ + ∣ 1m , 1n ñ) are symmetrized wavefunctions. We now examine the evolution of two-particle correlations in a three-site fully-connected network. As initial states we consider three different bosonic cases: (i) two indistinguishable particles in the separable state ∣ y (0)ñ = (∣ 11 , 12 ñ + ∣ 12 , 11 ñ) 2 , (ii) an incoherent two-distinguishable-particle state represented by r (0) = (∣ 11 , 12 ñá 11 , 12 ∣ + ∣ 12 , 11ñá 12 , 11 ∣) 2, and (iii) two particles in an entangled state 5 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al Figure 3. Density matrices (absolute value) for (a), (d) separable, (b), (e) incoherent and (c), (f) entangled states at t=0 ps and t=1 ps, respectively. Note that the stochastic coupling affects the system in such a way that when indistinguishable particles are injected into the system (figures 3(a) and (c)) the probability of finding both particles in the same site is the largest (bunching effect); whereas for distinguishable photons (figure 3(b)), the probability of finding them in different sites becomes larger (anti-bunching effect). The parameters used for the quantum network—namely site-energies, couplings and dephasing rates—are set to be ω1=ω2=ω3=5 ps−1; κ12=2 ps−1, κ13=κ23=1 ps−1, and γ12=γ13=γ23=0.38 ps−1. ∣ y (0)ñ = (∣ 11 , 11ñ + ∣ 12 , 12 ñ) 2 . Figure 3 shows the evolution of the initial states at t=1 ps. Notice that the stochastic fluctuations affect the system in such a way that, when indistinguishable particles (figures 3(a), (d) and (c), (f)) are injected in the system, the probability of finding both particles in the same site is the largest, that is, the photons bunch in all sites with the same probability. This effect could be thought of as a generalized Hong- Ou-Mandel effect produced by the pure-dephasing-like process. In striking contrast, when distinguishable photons are injected in the system (figures 3(b), (e)), the probability of finding them in different sites becomes larger, thus leading to an anti-bunching effect. An important aspect to point out regarding the results shown in figure 3 is the computation time required for obtaining them. Remarkably, solving the master equation in equation (13) takes ∼0.521 s; whereas the pure numerical solution of equation (12) takes ∼2.4 h. This implies that our derived master equation improves the computation time by at least four orders of magnitude, while providing the maximum accuracy possible (see appendix B for details). Recently, it has been shown that coherences arising from particle indistinguishability are robust against noise [24, 38]. By making use of our model, we have verified that in the steady-state, coherences accounting for particle indistinguishability do survive the impact of stochastic fluctuations in the coupling between sites (see appendix C for details). These results imply that it is possible, in principle, to find specific conditions for which many indistinguishable particles can traverse stochastically-coupled networks without losing their ability to interfere. Finally, notice that the generalization of our results to N correlated particles is straightforward following similar steps as above by introducing the N-particle probability amplitude N Yp, q, r , ¼(t ) = å ¼ ja, b, c , ¼[cap,,bq,,cr,,¼ + caper , b, c , ¼ + ...] , (14) a, b, c , ¼ ¼ with cap,,bq,,cr,,¼ = Up, a (t ) Uq, b (t ) Ur , c (t )..., where Um,n represents the probability amplitude for each particle at site n when it is injected into channel m. The superscript ‘per’ stands for the cyclic permutations of the subscripts p, q, r, K in the corresponding transition amplitudes. 6 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al 4. Conclusions In this work, we have derived a master equation for the propagation of correlated particles in quantum networks affected by off-diagonal dynamical disorder. Unlike commonly-used computational methods, where many stochastic trajectories are needed, our equation allows one to find the average trajectory of correlated particles in a single calculation. By using our results, we showed that the effect of introducing noise in the couplings of a quantum network leads to a dephasing-like process that destroy all coherences in the single-particle Hilbert subspace. Interestingly, we found that when two or more correlated particles propagate in a stochastically- coupled network, coherences accounting for the indistinguishability of the particles endure the impact of noise. These results may help elucidating the role of particle indistinguishability to preserve quantum coherence and entanglement propagating through complex dynamically-disordered systems. Acknowledgments This work was supported by DGAPA-UNAM under the project UNAM-PAPIIT IA100718, and by CONACYT under the project CB-2016-01/284372. Armando Perez-Leija and Kurt Busch acknowledge financial support by the Deutsche Forschungsgemeinschaft within the priority program SPP 1839 ‘Tailored Disorder’ (PE2602/2-2 and BU 1107/12-2). VM acknowledges financial support by Ministerio de Economia y Competitividad through Grant No. CGL2016-78156-C2-2-R. Appendix A. Derivation of the two-particle master equation We start by writing the expression for the probability amplitude dynamics of a quantum particle initiated at site n dUq, n = - iwq Uq, n - i å krq (t ) Ur , n, (A.1) dt r where ωn stands for the energy of the nth site, and the coupling between the rth and qth sites is given by krq (t ) = krq + frq (t ), with frq (t ) = fqr (t ) describing a Gaussian Markov process with zero average, that is á frq (t )ñ = 0, (A.2) á frq (t ) fjl (t ¢)ñ = grq drq, jld (t - t ¢). (A.3) Here drq, jl = drj dql + drl dqj , with δrq being the Kronecker delta. γrq denotes the noise intensity, that is, how strong the stochastic fluctuations are, and áñ denotes stochastic averaging. By defining the stochastic variable frq (t ) = - grq xrq (t ), we can write dUq, n = - iwq Uq, n - i å krq Ur , n + i å grq xrq (t ) Ur , n, (A.4) dt r r with the properties of the stochastic variable ξrq given by á xrq (t )ñ = 0, (A.5) á xrq (t ) x jl (t ¢)ñ = drq, jld (t - t ¢). (A.6) Notice that because noise (dynamic disorder) is introduced in the couplings, we must keep in mind that r ¹ q and, consequently, j ¹ l . Now, to compute the evolution of the two-particle density matrix rpq, p ¢ q ¢ = á ypq y*p ¢ q ¢ñ, with yp, q (t ) = åm = 1, n = 1xm, n [Up, n (t ) Uq, m (t )  Up, m (t ) Uq, n (t )], we first write d (ypq y*p ¢ q ¢) = - i [wp + wq - wp ¢ - wq ¢ ] ypq y*p ¢ q ¢ dt - i å klq ypl y*p ¢ q ¢ - i å klp ylq y*p ¢ q ¢ l l + i å klq ¢ ypq y*p ¢ l + i å klp ¢ ypq y*lq ¢ l l - iå glq ypl y*p ¢ q ¢xlq (t ) - i å glp ylq y*p ¢ q ¢xlp (t ) l l + iå glq ¢ ypq y*p ¢ lxlq ¢ (t ) + i å glp ¢ ypq y*lq ¢xlp ¢ (t ). (A.7) l l 7 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al We can formally integrate equation (A.7), and obtain t ⎧ ypq y*p ¢ q ¢ = ò0 dt ¢ ⎨ f (ypq y*p ¢ q ¢,¼) ⎩ - i å glq ypl (t ¢) y*p ¢ q ¢ (t ¢) xlq (t ¢) l - i å glp ylq (t ¢) y*p ¢ q ¢ (t ¢) xlp (t ¢) l + i å glq ¢ ypq (t ¢) y*p ¢ l (t ) xlq ¢ (t ¢) l ⎫ + i å glp ¢ ypq (t ¢) y*lq ¢ (t ¢) xlp ¢ (t ¢) ⎬ , (A.8) l ⎭ where f () is a function that contains all terms that do not depend on the stochastic variables. Concurrently, we can write the average of equation (A.7) as d á ypq y*p ¢ q ¢ñ = - i [wp + wq - wp ¢ - wq ¢] á ypq y*p ¢ q ¢ñ dt - i å klq á ypl y*p ¢ q ¢ñ - i å klp á ylq y*p ¢ q ¢ñ l l + i å klq ¢ á ypq y*p ¢ lñ + i å klp ¢ á ypq y*lq ¢ñ l l - iå glq á ypl y*p ¢ q ¢xlq (t )ñ l - i å glp á ylq y*p ¢ q ¢xlp (t )ñ l + i å glq ¢ á ypq y*p ¢ lxlq ¢ (t )ñ l + i å glp ¢ á ypq y*lq ¢xlp ¢ (t )ñ. (A.9) l It is clear that in order to obtain the master equation for rpq, p ¢ q ¢ (t ), we must evaluate the correlation functions in the last four terms of equation (A.9). To do so, we invoke the Novikov’s theorem [22, 47], which for the first correlation function in equation (A.9) takes the form d [ypl (t ) y*p ¢ q ¢ (t )] á ypl y*p ¢ q ¢xlq (t )ñ = å rs ò dt ¢áxlq (t ) xrs (t ¢)ñ dxrs (t ¢) , d [ypl (t ) y*p ¢ q ¢ (t )] =å rs ò dt ¢dlq,rsd (t - t ¢) dxrs (t ¢) , 1 d [ypl (t ) y*p ¢ q ¢ (t )] = 2 å dlq,rs dxrs (t ) . (A.10) rs Here, we have taken into account the fact that, in the Stratonovich interpretation [48], ò d (t ) = 1/2. We can then use equation (A.8) to write the functional derivative as d [ypl (t ) y*p ¢ q ¢ (t )] = - i å gsl yps (t ) y*p ¢ q ¢ (t ) dsl, rs dxrs (t ) s - i å gsp ysl (t ) y*p ¢ q ¢ (t ) dsp, rs s + i å gsq ¢ ypl (t ) y*p ¢ s (t ) dsq ¢ , rs s + i å gsp ¢ ypl (t ) y*sq ¢ (t ) dsp ¢ , rs, (A.11) s 8 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al where we used the relation dxsl dxrs = dsl, rs . By substituting this result into equation (A.10), we can write i á ypl y*p ¢ q ¢xlq (t )ñ = - å dsl, lq gsl rps, p ¢ q ¢ 2 s i - å dsp, lq gsp rsl, p ¢ q ¢ 2 s i + å dsq ¢ , lq gsq ¢ rpl, p ¢ s 2 s i + å dsp ¢ , lq gsp ¢ rpl, sq ¢. (A.12) 2 s Similarly, the remaining correlation functions are given by i á ylq y*p ¢ q ¢xlp (t )ñ = - å dsq, lp gsq rls, p ¢ q ¢ 2 s i - å dsl, lp gsl rsq, p ¢ q ¢ 2 s i + å dsq ¢ , lp gsq ¢ rlq, p ¢ s 2 s i + å dsp ¢ , lp gsp ¢ rlq, sq ¢, (A.13) 2 s i á ypq y*p ¢ lxlq ¢ (t )ñ = - å dsq, lq ¢ gsq rps, p ¢ l 2 s i - å dsp, lq ¢ gsp rsq, p ¢ l 2 s i + å dsl, lq ¢ gsl rpq, p ¢ s 2 s i + å dsp ¢ , lq ¢ gsp ¢ rpq, sl , (A.14) 2 s i á ypq y*lq ¢xlp ¢ (t )ñ = - å dsq, lp ¢ gsq rps, lq ¢ 2 s i - å dsp, lp ¢ gsp rsq, lq ¢ 2 s i + å dsq ¢ , lp ¢ gsq ¢ rpq, ls 2 s i + å dsl, lp ¢ gsl rpq, sq ¢. (A.15) 2 s Finally, by substituting equations (A.12)–(A.15) into equation (A.9) we obtain drpq, p ¢ q ¢ = - i (wp + wq - wp ¢ - wq ¢ ) rpq, p ¢ q ¢ dt 1 - å [(g lp + g lq + g lp ¢ + g lq ¢ ) - gpq - gp ¢ q ¢ ] rpq, p ¢ q ¢ 2 l - i å (klq rpl, p ¢ q ¢ + klp rlq, p ¢ q ¢ - klq ¢ rpq, p ¢ l - klp ¢ rpq, lq ¢ ) l - å (dpq g lq g lp rll, p ¢ q ¢ + dp ¢ q ¢ g lp ¢ g lq ¢ rpq, ll ) l + å (dqq ¢ g lq g lq ¢ rpl, p ¢ l + dqp ¢ g lq g lp ¢ rpl, lq ¢ ) l + å (dpq ¢ g lp g lq ¢ rlq, p ¢ l + dpp ¢ g lp g lp ¢ rlq, lq ¢ ) l + gqq ¢ rpq ¢ , p ¢ q + gqp ¢ rpp ¢ , qq ¢ + gpp ¢ rp ¢ q, pq ¢ + gpq ¢ rq ¢ q, p ¢ p , (A.16) which is the result shown in equation (13) of the main text. Appendix B. Comparison between master equation and the direct stochastic numerical simulation We now provide a quantitative comparison between the time evolution of a two-particle state obtained by means of our derived master equation and by directly implementing the stochastic equations. Figure B1 shows the 9 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al Figure B1. Density matrices (absolute value) for a separable state, ∣ y (0)ñ = (∣ 11 , 12 ñ + ∣ 12 , 11 ñ) 2 , at t=1 ps, t=3 ps, and t=5 ps, obtained by means of the derived master equations (a)–(c), and by the direct numerical evaluation of the stochastic equations (d)– (f). Figures 1(g)–(i) show the absolute difference between both solutions, Dr = ∣ r(pq , p ¢ q ¢ - r pq, p ¢ q ¢ master) (numerical) ∣, at the corresponding evolution times. evolution of a separable state, ∣ y (0)ñ = (∣ 11 , 12 ñ + ∣ 12 , 11 ñ) 2 , propagating in a dynamically-disordered three-site network. The parameters used for the quantum networks—namely site-energies, couplings and dephasing rates—are the same as those used for obtaining figure 3 of the main text. Figures B1(a)–(c) show the results obtained by using the derived master equation (equation (13) of the main text), whereas figures B1(d)–(f) show the results obtained by numerically solving equation (12) of the main text using the Taylor Integration package [49]. The latter were obtained by averaging over 10 000 different realizations of the two-particle random walk. It is important to highlight that the computation time required for each case was Tc(master) = 0.521 s, and Tc(numerical) = 2.4 h for the master equation and direct stochastic evaluation, respectively. Clearly, our derived equation improves the computation time by at least four orders of magnitude, while providing the maximum accuracy possible. For the sake of completeness, in figures B1(g)–(i), we have included the absolute difference between the absolute value of the density matrix elements obtained from the master equation and the numerical solution, i.e. Dr = ∣∣ r(pq , p ¢ q ¢ ∣ - ∣ r pq, p ¢ q ¢ master) (numerical) ∣∣. Finally, we would like to remark that while the derived master equation provides the exact solution, the accuracy of the stochastic-computation solution strongly depends on the number of realizations being used for the average, which implies that many realizations (and therefore longer computation times) are required in order to obtain reliable numerical results. This is the reason why, when possible, one should use master equations instead of direct stochastic numerical simulations. Appendix C. Quantitative analysis of two-particle coherence preservation To quantify the amount of surviving coherence in the steady states, we use two different coherence measures, the physically intuitive norm of coherence [50]: CN (r ) = å i ¹ j∣ ri, j∣, and the relative entropy of coherence [51] CE (r ) = S (rdiag ) - S (r ), with S representing the von Neumann entropy and ρdiag the matrix obtained from the density matrix ρ after removing all off-diagonal elements. Note that, in both measures, a totally mixed (or incoherent) state is signaled by a vanishing coherence measure. Figure C1 shows the evolution of the norm of coherence and the entropy of coherence for the initial separable (solid line), path-entangled (dashed line), and 10 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al Figure C1. Time evolution of (a) the norm of coherence CN (r ), and (b) the relative entropy of coherence CE (r ), considering three different initial states, namely separable (solid line), path-entangled (dashed line), and incoherent (dotted line) states. incoherent (dotted line) states. Interestingly, when distinguishable particles (incoherent states) are injected into the system (CN=CE=0), coherence due to propagation-induced indistinguishability rapidly emerges and the system evolves into a steady state where CN=0.250 and CE=0.061. Remarkably, the same results as those shown in figure C1 can be observed for any value of the dephasing rates. Therefore, we can convincingly say that under the influence of fluctuating couplings, identical particles always evolve into a steady state in which coherences due to indistinguishability perpetually prevail. ORCID iDs Roberto de J León-Montiel https://orcid.org/0000-0002-3045-3604 Vicenç Méndez https://orcid.org/0000-0002-2605-6606 Mario A Quiroz-Juárez https://orcid.org/0000-0002-5995-9510 Luis Benet https://orcid.org/0000-0002-8470-9054 Armando Perez-Leija https://orcid.org/0000-0002-6218-295X References [1] Rebentrost P, Mohseni M, Kassal I, Lloyd S and Aspuru-Guzik A 2009 New J. Phys. 11 033003 [2] Chin A, Datta A, Caruso F, Huelga S and Plenio M 2010 New J. Phys. 12 065002 [3] Park H et al 2016 Nat. Mater. 15 211 [4] Ziman J M 1979 Models of Disorder: The Theoretical Physics of Homogeneously Disordered Systems (Cambridge: Cambridge University Press) [5] León-Montiel R J, Quiroz-Juárez M A, Quintero-Torres R, Domínguez-Juárez J L, Moya-Cessa H M, Torres J P and Aragón J L 2015 Sci. Rep. 5 17339 [6] Saikin S K, Eisfeld A, Valleau S and Aspuru-Guzik A 2013 Nanophotonics 2 21 [7] Viciani S, Lima M, Bellini M and Caruso F 2015 Phys. Rev. Lett. 115 083601 [8] Biggerstaff D N, Heilmann R, Zecevik A A, Grafe M, Broome M A, Fedrizzi A, Nolte S, Szameit A, White A G and Kassal I 2016 Nat. Commun. 7 11282 [9] Caruso F, Crespi A, Ciriolo A G, Sciarrino F and Osellame R 2016 Nat. Commun. 7 11682 [10] Schonleber D W, Eisfeld A, Genkin M, Whitlock S and Wuster S 2015 Phys. Rev. Lett. 114 123005 [11] Trautmann N and Hauke P 2018 Phys. Rev. A 97 023606 [12] Breuer H P and Petruccione F 2002 The Theory of Open Quantum Systems (Oxford: Oxford University Press) [13] Rigovacca L and Franco C D 2016 Sci. Rep. 6 22052 [14] Siloi I, Benedetti C, Piccinini E, Piilo J, Maniscalco S, Paris M G A and Bordone P 2017 Phys. Rev. A 95 022106 [15] Rossi M A C, Benedetti C, Borrelli M, Maniscalco S and Paris M G A 2017 Phys. Rev. A 96 040301 [16] de Vega I and Alonso D 2017 Rev. Mod. Phys. 89 015001 [17] Bromley T R, Cianciaruso M and Adesso G 2015 Phys. Rev. Lett. 114 210401 [18] Beggi A, Buscemi F and Bordone P 2016 Quantum Inf. Process. 15 3711 [19] Piccinini E, Benedetti C, Siloi I and Paris M G A 2017 Comput. Phys. Commun. 215 235 [20] Eisfeld A and Briggs J S 2012 Phys. Rev. E 85 046118 [21] León-Montiel R J and Torres J P 2013 Phys. Rev. Lett. 110 218101 [22] Novikov E A 1965 Sov. Phys.—JETP 20 1290 [23] Strunz W T and Yu T 2004 Phys. Rev. A 69 052115 [24] Perez-Leija A, Guzman-Silva D, León-Montiel R J, Graefe M, Heinrich M, Moya-Cessa H, Busch K and Szameit A 2018 njp Quantum Inf. 4 024002 [25] León-Montiel R J and Quinto-Su P A 2017 Sci. Rep. 7 44287 [26] Mostame S, Rebentrost P, Eisfeld A, Kerman A J, Tsomokos D I and Aspuru-Guzik A 2012 New J. Phys. 14 105013 [27] Mostame S, Huh J, Kreisbeck C, Kerman A J, Fujita T, Eisfeld A and Aspuru-Guzik A 2017 Quantum Inf. Process. 16 44 [28] Potocnik A et al 2018 Nat. Commun. 9 904 11 New J. Phys. 21 (2019) 053041 R de J León-Montiel et al [29] Quiroz-Juárez M A, Aragón J L, León-Montiel R J, Vázquez-Medina R, Domínguez-Juárez J L and Quintero-Torres R 2016 Europhys. Lett. 116 50004 [30] Wang B X, Tao M J, Ai Q, Xin T, Lambert N, Ruan D, Cheng Y C, Nori F, Deng F G and Long G L 2018 njp Quantum Inf. 4 52 [31] Gorman D J, Hemmerling B, Megidish E, Moeller S A, Schindler P, Sarovar M and Haffner H 2018 Phys. Rev. X 8 011038 [32] Maier C, Brydges T, Jurcevic P, Trautmann N, Hempel C, Lanyon B P, Hauke P, Blatt R and Roos C F 2019 Phys. Rev. Lett. 122 050501 [33] Laing C and Lord G J 2008 Stochastic Methods in Neuroscience (Oxford: Clarendon) [34] León-Montiel R J, Svozilík J and Torres J P 2014 Phys. Rev. E 90 012108 [35] Misra B and Sudarshan E C G 1977 J. Math. Phys. 18 756 [36] Levi F and Mintert F 2014 New J. Phys. 16 033007 [37] León-Montiel R J, Vallés A, Moya-Cessa H M and Torres J P 2015 Laser Phys. Lett. 12 085204 [38] Perez-Leija A, León-Montiel R J, Sperling J, Moya-Cessa H, Szameit A and Busch K 2017 J. Phys. B: At. Mol. Opt. Phys. 51 024002 [39] Giuseppe G D, Martin L, Perez-Leija A, Keil R, Dreisow F, Nolte S, Szameit A, Abouraddy A F, Christodoulides D N and Saleh B E A 2013 Phys. Rev. Lett. 110 150503 [40] Klyshko D N 1982 Sov. Phys.—JETP 56 753 [41] Abouraddy A F, Saleh B E A, Sergienko A V and Teich M C 2001 Phys. Rev. Lett. 87 123602 [42] Saleh B E A, Teich M C and Sergienko A V 2005 Phys. Rev. Lett. 94 223601 [43] Abouraddy A F, Saleh B E A, Sergienko A V and Teich M C 2002 J. Opt. Soc. Am. B 19 1174 [44] Bromberg Y, Lahini Y, Morandotti R and Silberberg Y 2009 Phys. Rev. Lett. 102 253904 [45] Lebugle M, Gräfe M, Heilmann R, Perez-Leija A, Nolte S and Szameit A 2015 Nat. Commun. 6 8273 [46] Weimann S, Perez-Leija A, Lebugle M, Keil R, Tichy M, Gräfe M, Heilmann R and Nolte S 2016 Nat. Commun. 7 11027 [47] Méndez V, Horsthemke W, Mestres P and Campos D 2011 Phys. Rev. E 84 041137 [48] van Kampen N G 1981 J. Stat. Phys. 24 175 [49] Perez J A and Benet L 2017 Taylorintegration.jl https://github.com/PerezHz/TaylorIntegration.jl [50] Streltsov A, Adesso G and Plenio M B 2017 Rev. Mod. Phys. 89 041003 [51] Baumgratz T, Cramer M and Plenio M B 2014 Phys. Rev. Lett. 113 140401 12