The study of open quantum systems [1–3] has grown increasingly important as quantum technology continues to advance. The quality of engineered quantum devices hinges on our ability to understand and control those devices in a noisy environment [4]. An open environment is typically understood to cause undesirable effects to a given system of interest. For example, random noise due to an open environment can cause quantum decoherence which undermines quantum processing and reduces the reliability of quantum devices [5,6]. Thus, researchers must develop novel methods of active quantum control and manipulation [7–11].
Befittingly, an open environment can also facilitate desirable effects, e.g., a non-Markovian environment with finite memory time can regenerate quantum coherence and entanglement in discrete variable systems such as qubits [12] or generate entanglement between an optical cavity field and a movable mirror [13]. It can also be utilized for quantum memory [14] by using an exotic system such as an oscillating bound state that is in a one-dimensional photonic lattice coupled to a ”giant atom” [15] or a means to autonomously distribute entanglement [16], which is highly desirable in quantum computing [17], networks [18], metrology, and optics. Furthermore, many studies involving continuous variables systems, e.g., harmonic oscillators [19–28], have also treated environmental noise not as a hindrance but as an asset.
Despite that many harmonic oscillator chain models have been studied, many were limited in the number of oscillators, e.g., N = 2 or N = 3, and some were analyzed only in either a Markovian or non-Markovian environment and did not compare the environments. For those that were of size N and analyzed in a non-Markovian environment, their focus did not include how entanglement dynamics were governed by system and environmental parameters and did not suggest applications. In our research, we analyze the dynamics of a coupled harmonic oscillator ring system that is embedded in an environment with varying memory time. This allows us to cross over both Markovian and non-Markovian domains and compare their entanglement pathways.
It is important to understand these pathways to improve quantum information distribution in the presence of environmental noise, in order to enable applications such as entanglement networks or enhanced quantum state transfer [29]. We accomplish this by deriving two sets of differential equations that govern both Markovian and non-Markovian environments from first principles for N oscillators in continuous variables. Using these equations, we can study entanglement dynamics by systematically varying system and environmental parameters, which are the nearest neighbor and environmental coupling, the environmental memory and the central frequency.
In the results that follow, we have found many benefits of a non-Markovian environment over a Markovian one. Firstly, it can improve quantum communications by boosting the magnitude of bipartite entanglement of oscillator pairs. Secondly, by tuning system oscillator frequencies to resonate with the central frequency of the environment, long-lived bipartite entanglement can persist within a ring of oscillators connected by nearest neighbor coupling without active quantum control. Thirdly, after removing all nearest neighbor coupling, a non-Markovian environment can passively generate bipartite entanglement within initially unentangled oscillator pairs while a Markovian one cannot. We find that this can be controlled by choosing suitable parameters for the environmental memory time, the central frequency, and the system oscillator frequencies.
Lastly, having extracted the governing parameters, this allows us to engineer the entanglement dynamics and propose a way to create a “hidden” quantum network in a proof of principle example consisting of N = 6 oscillators. This is a network of pairwise entangled oscillators controlled by multiple parties, where each sub-party can control a(n) oscillator(s) that is a part of an entangled pair without any party knowing who is entangled with whom, provided that they conceal the frequency their oscillator(s).
This work is organized as follows. Section 1 is our introduction, where we establish the relevance of our work. In Section 2, we present our model and the equations that govern our system in a Markovian and a non-Markovian environment. In Section 3, we present our results that compare the entanglement dynamics of our system when simulated in a closed system, a Markovian environment and a non-Markovian one while varying the nearest neighbor coupling and environmental coupling independently for a simple network. Afterwards, we remove the nearest-neighbor coupling to isolate non-Markovian effects. Lastly, we showcase how a system of independent harmonic oscillators can establish a hidden quantum network by using the environmental memory effect and resonance with the central frequency of a non-Markovian environment. In Section 4, we summarize and discuss our findings.
The total Hamiltonian is defined in Eq.(1a), where Eq.(1b), Eq.(1c), and Eq.(1d) is the Hamiltonian of the system of interest, the environment or common bath, and the interaction between the system and environment, respectively.
The system of interest consists of N coupled harmonic oscillators described with position and momentum operators,

A model of N harmonic oscillators connected with nearest neighbor coupling, ηk, arranged as a 1-dimensional linear ring, where k ∈ {1, …, N}.
We model the Hamiltonian of our environment as the sum of a set of independent harmonic oscillators described with position and momentum operators,
In general,
We use the approximate non-Markovian master equation to evolve our system in a weak non-Markovian environment [30–33], which is given in Eq.(2), (setting ħ = 1). The operators
For this model, we define Ō(t) in Eq.(3), where Ô(t, s) takes the form
After we satisfy Eq.(A1), we can now reformulate Eq.(3) as
The Ornstein-Uhlenbeck correlation function for the environment is α (t, s) [34–36]. The inverse of the correlation time, τc = γ−1, or the environmental memory, is the parameter γ that defines the time window for quantum information leaked from the system into the environment to return back to the system. The central frequency, Ω, is the dominate frequency value within a spectrum of frequencies that environmental oscillators could possess, characteristic of colored noise. Together, γ and Ω are the key parameters that determine the non-Markovian effects.
If we take the limit as γ → ∞, then τc → 0, and
Therefore, we can recover all system dynamics in the Markovian regime and make comparisons with our non-Markovian regime.
In all our simulations, we initialize the system harmonic oscillators as Gaussian states [37]. These states can be entangled [38–40] by violating the Peres-Horodecki separability criterion [41] or satisfying the inseparability criteria [42,43] for continuous variables. We closely follow the procedure used by Plenio [44] to calculate the logarithmic negativity within our oscillators as the primary method to characterize the entanglement within them. While there exist many entanglement measures, we choose logarithmic negativity because our Hamiltonian uses continuous variables and is linear, which preserves Gaussian states. Furthermore, our linear Hamiltonian allows for the computation of the covariance matrix relatively easily, which is the input in the logarithmic negativity function, thereby making it also simple to compute. For further technical details regarding entanglement with continuous variables or Gaussian states we recommend the references [45,46].
To obtain the full entanglement dynamics of our system oscillators, we calculate the logarithmic negativity pairwise for every possible oscillator pair combination. The procedure for the calculation of the logarithmic negativity between oscillators i and j, where i ≠ j ∈ {1, …, N}, starts by constructing their covariance matrix, C, with elements Ck,l = 2Re[〈RkRl〉 – 〈RkRl〉], where the list
Afterwards, we take the partial transpose of
Finally, given that the individual or joint expectation values of the position and momentum operators are time dependent, 𝒩i,j is also time dependent. For dimensionless units, we let τ = ω0t, where we let ω0 = 1 to be an arbitrary reference frequency. Hence, our results will show 𝒩i,j as a function of τ instead of t.
For simplicity, we let our system consist of N = 4 harmonic oscillators. Since we aim to study entanglement propagation within our ring system, we initialize the total state as |Ψ〉 = |0,0〉 ⊗ S2(ζ)|0,0〉, where S2(ζ) is the two-mode squeezing operator, ζ = reiϕ, ϕ is the angle in optical phase space, and r is the squeezing parameter, which we set to r = 1. To express |Ψ〉 in the coordinate basis, Ψ(x1, x2, x3, x4, r) = ψ0(x1) ⊗ ψ0(x2) ⊗ ψs(x3,x4, r). The first two oscillators, 1 and 2, are initialized as vacuum states given in Eq.(8a) while the remaining oscillators, 3 and 4, are initialized together as a two-mode squeezed state [47] given in Eq.(8b).
For all results that follow, we use Ψ(x1, x2, x3, x4, r) unless otherwise specified. Lastly, we further simplify by setting mi = 1 and letting all nearest neighbor coupling, environmental coupling, and oscillator frequencies to vary identically, i.e. ηi = η, βi = β, ωi = ω, where i = {1, …,4}, for all results except for the results in section 3.4, it varies ωi differently.
As we analyze our ring of oscillators in a closed system (non-existent noise), a Markovian environment (white noise), and a non-Markovian environment (colored noise), the entanglement within any oscillator pair are colored with links as shown in Figure 2. This color scheme applies for all results unless otherwise stated. Note, for N = 4, oscillator pairs 1,3 and 2,4 are colored Yellow because we find that their entanglement dynamics are identical, i.e. 𝒩1,3 (τ) = 𝒩2,4 (τ). Similarly, this is true for oscillator pairs 2,3 and 1,4, i.e. 𝒩2,3 (τ) = 𝒩1,4(τ), and we color them Green. If their entanglement dynamics are unique, we use different colored lines. i.e. 𝒩1,2 (τ) is Blue and 𝒩3,4 (τ) is Red.

A pictorial representation of N = 4 oscillators arranged in a ring configuration embedded in a shared-Closed (non-existent noise), Markovian (white noise), and non-Markovian (colored noise) environment, respectively. We use different colored lines connecting oscillator i and j to represent unique entanglement dynamics or logarithmic negativity 𝒩i,j (τ). 𝒩1,2 (τ) is Blue, and 𝒩3,4 (τ) is Red. If the entanglement dynamics are identical within an oscillator pair, we use the same colored lines. 𝒩1,3 (τ) = 𝒩2,4 (τ) are colored Yellow and 𝒩2,3 (τ) = 𝒩1,4 (τ) are colored Green.
It is evident in the left plots in Figure 3a,b, and c that the nearest neighbor coupling η is directly proportional to the rate in which entanglement within each oscillator pair reaches successive relative maximums regardless of its environment. For instance, the spacing in τ between relative maximums or peaks are smaller for the Blue and Red curves for larger values of η. Moreover, the peaks for 𝒩1,3 (τ) < 𝒩1,4(τ) or 𝒩2,4 (τ) 𝒩2,3 (τ) which is due to the nearest neighbor coupling between oscillator pair 1,4 or pair 2,3 having a larger effect than the indirect or effective next-nearest neighbor coupling between oscillator pair 1,3 or pair 2,4. Additionally, the magnitude of the peaks of all oscillator pairs changes with τ, which can be seen most clearly with the Blue and Red curves.

Entanglement dynamics in different environments η is varied. For all plots ω = 1. 𝒩1,2 (τ) is Blue, and 𝒩3,4 (τ) is Red. 𝒩1,3 (τ) = 𝒩2,4 (τ) are colored Yellow and 𝒩2,3 (τ) = 𝒩1,4 (τ) are colored Green. Left side shows η dependence for 0.01 ≤ η ≤ 0.07. Right side is only for η = 0.07 and the dashed Gray line approximates an enveloping function following the relative maximums of 𝒩1,2 (τ) and 𝒩3,4 (τ).
For the Blue and Red curves in the right plot of Figure 3a, the alternating peaks approximately follow a sinusoidal enveloping function (dashed Gray line). Also, the smaller Green curve peaks before and after the larger Blue and Red curves and the Yellow curve subtly peaks between the curves just mentioned. The order in which these peaks occur suggests that the overall entanglement ”flows” within the ring system in the following manner. Entanglement initially exists between oscillator pair 3,4, then it drastically drops towards zero while simultaneously developing between oscillator pair 1,3 and pair 2,4 causing them to both briefly peak and then go to zero. Following this, oscillator pair 1,4 and pair 2,3 peaks and goes to zero and then very briefly oscillator pair 1,3 and pair 2,4 peaks again and drops to zero before finally oscillator 1,2 peaks then drops to zero. This order then happens in reverse and the cycle repeats, which indicates that entanglement flows roughly in both clockwise and counterclockwise directions simultaneously around the ring system, especially for small η. For larger values, this circular flow of entanglement may become more convoluted.
In Figure 3b, the magnitude of entanglement within any oscillator pair substantially decays as a function of τ in a Markovian environment, which starkly contrasts the results in Figure 3a. However, a shred of similarity is that the Yellow and Green curves are peaking in a similar manner but only initially. Entanglement ”tries” to flow circularly, but the drastic decay overwhelms the tiny Yellow and Green peaks, making them approach zero much faster than the Red and Blue peaks. Hence, entanglement only exists within oscillator pair 1,2 and pair 3,4, where they take turns reaching their peaks. As we follow their enveloping function (dashed Gray line) it drops towards zero as expected, which reinforces the fact that a Markovian environment poorly facilitates entanglement.
Interestingly, while one expects the entanglement in oscillator pairs to exponentially decay in a Markovian environment it does not exactly appear to do so given the non-monotonic dashed Gray line. The reason is ultimately due to nearest neighbor coupling with the remaining system oscillators. When we calculate the logarithmic negativity for a pair of oscillators, we trace everything else out besides those oscillators, so the remaining system oscillators act as an “effective non-Markovian environment” and is responsible for the non-monotonic decay. If we remove all nearest neighbor-coupling in the forthcoming results, then this non-monotonic exponential decay should disappear, which is exactly what we see in section 3.2.2.
In Figure 3c, the entanglement dynamics flows circularly for small η, but it may become more complicated with increased η, e.g., η ≥ 0.04 for certain durations. Also, we find that the entanglement within oscillator pairs that are next-nearest neighbors diminishes as η increases, i.e., the tiny Yellow peaks. As for the small Green peaks, they remain for all τ and η. Furthermore, the enveloping function (dashed Gray line) that follows the Blue and Red curves is roughly periodic, but it decays overall. Within the parameter ranges we considered, a non-Markovian environment typically allows for entanglement to persist in our ring system much longer than a Markovian one.
To further compare how our ring system behaves in different environments, we vary the environmental coupling parameter β and analyze the resulting entanglement dynamics as shown in Figure 4. We choose to show when τ = {400, …,500} because this range shows the entanglement peaks of oscillator pairs alternating regularly and clearly as opposed to showing the full graph, which would be too visually cluttered.

Entanglement dynamics in different environments as β varies while η = 0.07, ω = 1. Solid lines colored Blue, Red, Yellow and Green correspond to 𝒩1,2(τ), 𝒩3,4(τ), 𝒩1,3(τ) = 𝒩2,4(τ), and 𝒩2,3(τ) = 𝒩1,4(τ), respectively in a Markovian environment. Dashed lines colored Cyan, Magenta, Orange and Lime Green correspond to 𝒩1,2(τ), 𝒩3,4(τ), 𝒩1,3(τ) = 𝒩2,4 (τ), and 𝒩1,4(τ) = 𝒩2,3(τ), respectively in a non-Markovian environment. Larger to smaller peaks correspond to βM = {0.01, …, 0.04} with increments of 0.01 or βNM = {0.03, …,0.12} with increments of 0.03 and Ω = 0, γ = 0.5.
The entanglement within oscillator pairs in a Markovian environment are depicted as solid colored lines. As we increase βM = {0.01, …, 0.04} with increments of 0.01, the overall magnitude of entanglement within all oscillator pairs decreases, which corresponds to the larger peaks to smaller peaks viewed top to bottom. For a non-Markovian environment, the entanglement within oscillator pairs are dashed lines and are colored Cyan, Magenta, Lime Green, and Orange for 𝒩1,2 (τ), 𝒩3,4 (τ), 𝒩1,4(τ) = 𝒩2,3(τ), and 𝒩1,3(τ) = 𝒩2,4(τ), respectively. Proportionally, as β increases the overall magnitude of entanglement within all oscillator pairs also decreases. Larger peaks to smaller peaks correspond to βNM = {0.03, …, 0.12} with increments of 0.03. As evident from the plot, each 𝒩i,j(τ)|βM ≈ 𝒩i,j(τ) |βNM. Thus, for the parameter ranges we considered, we found a case where a non-Markovian environment may permit entanglement dynamics that are near identical to the dynamics in a Markovian one but with βNM = 3βM.
In Figure 5, we set η = 0 to isolate the effects of the environment. This decouples all oscillators in our ring system into independent oscillators. The colors Dark Purple, Red, and Magenta correspond to β = {0.04,0.08,0.12} in the plot and as β increases, we find that the entanglement initially set within oscillators 3 and 4 monotonically decays faster. This decay also persists regardless of environment, but a Markovian environment loses entanglement substantially more in terms of magnitude and duration compared to a non-Markovian one, which can be seen with the solid or dashed lines respectively.

Comparison of entanglement preservation within oscillator pair 3,4 in different environments. 𝒩3,4 (τ)|Markovian are solid lines while 𝒩3,4 (τ) |non–Markovian are dashed lines. For dashed lines, Ω = 0 and γ = 0.5. For either dashed or solid lines, η = 0, ω = 1 and the colors Dark Purple, Red, and Magenta corresponds to β = {0.04,0.08,0.12}.
A Markovian and non-Markovian environment can have near identical entanglement decay, e.g., the solid Dark Purple line overlaps the dashed Magenta line. Given the parameter range we considered, this can occur when
Interestingly, when system oscillator frequencies are resonant with the central frequency of the non-Markovian environment, e.g., ω = Ω = 1, the peaks of the Blue and Red curves are larger and more regular compared to the non-resonant case. This contrast is apparent when comparing the left plot of Figure 6 to the right plot of Figure 3c. If we extend the simulation time, we can see that the enveloping function that follows the Blue and Red peaks for the resonant case (solid Black line) outlasts the non-resonant case (dashed Gray line) in the right plot of Figure 6. Hence, by maintaining resonance we can further improve entanglement preservation in our ring system in a manner where it flows back and forth between oscillator pairs 1,2 and 3,4 indefinitely.

Entanglement dynamics in a non-Markovian environment involving resonance. For both plots, η = 0.07, β = 0.06, ω = Ω = 1, γ = 0.5. Left plot shows the resonance case where 𝒩1,2 (τ) is Blue, and 𝒩3,4 (τ) is Red. 𝒩1,3 (τ) = 𝒩2,4(τ) are colored Yellow and 𝒩2,3 (τ) = 𝒩1,4 (τ) are colored Green. Right plot shows an extended simulation of the enveloping functions for the resonance/non-resonance cases which are shown as a solid Black line and a dashed Gray line, respectively. Note, the dashed Gray line originates/continues from Figure 3c.
We can understand the resonance effect in more detail by varying Ω/ω = {0.1,…,2} with increments of 0.1 while we set ω = 1. We again decouple our system oscillators by letting η = 0 and show the entanglement dynamics of our system at resonance and as it passes through it in the left and right plots of Figure 7, respectively.

Entanglement dynamics in a non-Markovian environment for oscillator pair 3,4 and pair 1,2 at resonance and non-resonance. Entanglement within all other pairwise combinations are effectively zero. For all plots shown, η = 0, β = 0.30, ω = 1, and γ = 0.5. Left plot has Ω/ω = 1 and shows 𝒩1,2 (τ) as Blue, and 𝒩3,4(τ) as Red‥ Solid Cyan and Magenta lines are the average values of 𝒩1,2(τ) and 𝒩3,4(τ) for τ = {200,…, 1000}, respectively. Zoomed plots within the left plot for τ = {950,…, 1000} reveal that both 𝒩1,2 (τ) and 𝒩3,4 (τ) behave sinusoidally in the long time limit. Right plots shows the average values at saturation of 𝒩1,2 (τ) and 𝒩3,4 (τ) as Blue and Red dots. The Blue and Red lines that join the dots approximate functions with a maximum occurring when Ω/ω = 1 and are highlighted in Cyan and Magenta.
Contrary to Figure 5, entanglement within oscillator pair 3,4 is not just preserved longer than the Ω = 0 case, it can be preserved indefinitely with larger magnitude and with larger environmental coupling, i.e. β = 0.30. Also, the entanglement dynamics resembles the latter half of an inverted logistic function. This is evidenced by the left plot where the Red curve decreases initially and then behaves sinusoidally with constant amplitude, which can be seen in the zoomed plots within the left plot. Moreover, the Blue curve increases simultaneously and then also behaves sinusoidally with constant amplitude, resembling the later half of a logistic function. While these curves may suggest that entanglement is transferred from oscillator pair 3,4 to pair 1,2, we find that entanglement is actually generated simultaneously within oscillator pair 1,2 and 3,4 and they reach saturation at the same time because all oscillators have identical β, which will be explained in the following sections. Also, we find that entanglement within all other pairwise combinations are near or exactly zero, but if we change their initial states, the dynamics for all pairwise combinations are present and we will elaborate on this immediately after this section.
For each increment of Ω where Ω/ω ≠ 1, 𝒩1,2 (τ) and 𝒩3,4 (τ) increase and decrease respectively just like the left plot in Figure 7 but with extra oscillations initially before saturating to a different average value. For clarity, we do not show their plots for every Ω increment, but we do show their average values at saturation and depict them as dots with their corresponding color in the right plot. At resonance, the magnitude of the average values of 𝒩1, 2 (τ) and 𝒩3,4 (τ) at saturation maximizes and are shown as Cyan and Magenta lines/dots and we note that the resonance curves are asymmetric, which has been seen in similar work [48]. Crucially, we find that the mere presence of a non-zero central frequency regardless of resonance is responsible for the entanglement within oscillator pairs to be preserved in the long time limit. However, the non-resonance ratio Ω/ω ≠ 1 cannot be too far from unity, otherwise both 𝒩1,2 (τ) and 𝒩3,4 (τ) ≈ 0. In sum, the central frequency’s role enables indefinite preservation and the generation of entanglement within oscillator pairs despite not satisfying resonance exactly, but if it does then the amount of entanglement maximizes simultaneously.
Entanglement generation within oscillator pairs in a non-Markovian environment is not only dependent upon the number of oscillators, but also their initial states. We examine N = 4 identical and independent oscillators initialized as vacuum states, Ψ(x1, x2, x3, x4)separable = ψ0 (x1) ⊗ ψ0(x2) ⊗ ψ0(x3) ⊗ ψ0(x4) and vary Ω = {1, …, 2} with increments of 0.05 as we maintain resonance, Ω = ω. We show the results in Figure 8 and find that all possible pairwise oscillator combinations have entanglement generated within them identically and have the same dynamics, 𝒩i,j (τ), where i ≠ j ∈ {1,…,4 }.

Entanglement generation within all pairwise oscillator combinations 𝒩i,j (τ), where i ≠ j ∈ {1, …,4} with η = 0, β = 0.30, γ = 0.5, and Ω = ω. 𝒩i,j(τ) is shown as pastel colored lines of Pink, Yellow, Green, and Blue and correspond to Ωsub = {1.25,1.5,1.75,2}. Pastel colored dots correspond to the average value of 𝒩i,j (τ) for the duration τ = {50, …,200 } with different values of Ω = {1, …,2 } incremented at 0.05.
The entanglement dynamics, 𝒩i,j(τ), for a subset of Ω values, Ωsub = {1.25,1.5,1.75,2 }, are depicted as pastel colored lines of Pink, Yellow, Green, and Blue, respectively. They are read with the vertical and bottom horizontal axis. Pastel colored dots correspond to the average value of 𝒩i,j(τ)|Ω for the duration τ = {50, …,200} are read with the vertical and top horizontal axis. We find that 𝒩i,j (τ)|Ω resembles the latter half of a logistic function but has small oscillations with constant amplitude for a given Ω. As Ω increases, the average values for 𝒩i,j (τ)|Ω at saturation increases but the amplitudes decrease.
While we have not determined exactly how entanglement generation depends on the number of oscillators and each of their initial states, we suspect that these preliminary results point to a underlying physical mechanism that involves a complicated interdependence among oscillators and warrants future work.
In Figure 9, we explore the parameter space of β by Ω, where β = {0.30,…,1.05} incremented by 0.05 and Ω = ω = {1, …,8} incremented by 1. Within this parameter space, we find similar entanglement dynamics shown in Figures 7 and 8, i.e., they resemble the latter half of a logistic function or its inverted counterpart. As β increases the duration required to reach saturation, τsat, decreases until it reaches a minimum, τsat → min(τsat). If β is allowed to continue increasing even after reaching min(τsat), then the amplitude of oscillations within the overall entanglement dynamics begins growing and fluctuating wildly and the saturated value obtained prior begins reducing in value until it reaches zero for oscillator pairs 1,2 and 3,4. Therefore, within the parameter space there exists an optimum environmental coupling, βopt|Ω, where we can simultaneously have min(τsat) and the maximum average value of entanglement within oscillator pair 1,2 and 3,4.

Entanglement dynamics for 𝒩1,2(τ)|Ω and 𝒩3,4( τ)|Ω using βopt|Ω, where Ω = ω = {1,…, 8} with increments of 1, which correspond to a Dark Blue to a faded Light Blue and a Dark Red to a faded Light Red, respectively. 𝒩1,2 (τ)|Ω=2 and 𝒩3,4 (τ)|Ω=2 using β = 0.30 are shown as overlaid Blue and Red dashed curves overlaid. Remaining parameters used are η = 0 and γ = 0.5.
We show 𝒩1,2 (τ)|Ω and 𝒩3,4 (τ)|Ω with solid lines following a color range. They correspond to a Dark Blue to a faded Light Blue and a Dark Red to faded Light Red and both use βopt|Ω. If we do not use βopt|Ω, e.g., β = 0.30 with Ω = 2, then the entanglement dynamics for oscillator pair 1,2 and 3,4 would show saturation occurring much later as evidenced by the overlaid Blue and Red dashed curves, respectively. Given resonance and η = 0, β governs τsat and there exists βopt|Ω. Therefore, the simultaneous increase and decrease of entanglement within oscillator pair 1,2 and 3,4 is due to using identical environmental coupling, which also occurred in Figure 7. The takeaway is that if we desire maximum entanglement in the shortest t possible, we must use βopt|Ω and maximize Ω while maintaining resonance.
In Figure 10, we examine how generated entanglement within oscillator pair 1,2 and 3,4 depends on the environmental memory. We use the same parameters associated with 𝒩1,2 (τ)|Ω=2 and 𝒩3,4(τ)|Ω=2 in Figure 9, but extend the simulation to τ = {0, …, 1000} and vary γ = {0.5,…, 15} with increments of 0.5. Once maximum entanglement saturation is reached, their values are plotted as Blue and Red dots for each value of γ.

Entanglement saturation value within oscillator pair 1,2 and 3,4 during τ = {0, …1000} as a function of γ = {0.5, …, 15} incremented with 0.5. Blue dots are for oscillator pair 1,2 and Red dots are for oscillator pair 3,4 and their linearly fitted functions are shown as Blue and Red lines. Resonance is maintained, Ω = ω = 2, η = 0, and β = 0.30.
We find a negative linear correlation between γ and the entanglement saturation value. This follows our physical intuition in that as γ → ∞, the correlation time τc → 0 and the environment becomes Markovian which does not facilitate entanglement generation. The difference between the Blue and Red lines are approximately constant when γ < 4 and their slopes are also approximately equal which indicates that the non-Markovian environment acts on the oscillator pairs impartially. Lastly, we expect that any underlying analytical function that governs entanglement generation would have γ dependence.
In this section, we demonstrate how to utilize environmental memory and resonance with the central frequency of a non-Markovian environment to engineer a “hidden” quantum network. The network is hidden because entanglement is generated through the environment itself, which is not explicitly observable. Additionally, it is also hidden in another sense. Parties participating within this network have no knowledge of who is entangled with whom by concealing the frequency of oscillators held in their possession, thus their entanglement connection is also “hidden”.
We showcase this with six possible network configurations each consisting of N = 6 independent system oscillators, in a shared non-Markovian environment with γ = 0.5 and Ω = 4. We can refer to each network configuration and its entanglement dynamics as the row and column of plots, Config(row, column), shown in Figure 11, e.g., the top left plot corresponds to Config(1,1).

Six different entanglement network configurations each involving N = 6 oscillators are depicted as colored circles. The colors Red, Yellow, Green, and Blue are unrelated to all previous Figures. They instead correspond to the system oscillator frequencies ω = {2,3,4,5}. A solid line linking a pair of oscillators means entanglement was initially set within them at τ = 0, i.e., solid Green or Yellow line between oscillators 3 and 4. Dotted lines linking a pair of oscillators means entanglement was not present within them initially but generates within them later, τ > 0. A stand alone oscillator means that entanglement is neither generated or initialized involving that oscillator with any other one for τ ≥ 0, i.e., oscillator 6. The corresponding entanglement dynamics for each network configuration are shown below with solid or dashed lines, which means that entanglement was either initialized or generated within them. The parameters used in all plots are η = 0, β = 0.30, Ω = 4, and γ = 0.5.
If oscillator pairs are initialized with entanglement, they are shown with solid lines in the configurations and plots. If oscillator pairs are not initialized with entanglement but has it generated within them later, they are shown with dotted lines in the configurations and dashed lines in the plots. We set all oscillators to have identical mass and environmental coupling, and remove nearest neighbor coupling, i.e., mi = 1, βi = 0.30, ηi = 0, where i = {1, …, 6}. Also, we allow oscillators take on any of the frequency values in list, ω = {2,3,4,5 } and are color coded Red, Yellow, Green, or Blue, respectively. Note, these colors used here have no relation to all previous Figures. The initial state of our system for each network is Φ(x1, x2, x3, x4, r, x5, x6) = ψ0(x 1) ⊗ ψ0(x2) ⊗ ψs (x3, x4, r) ⊗ ψ0(x5) ⊗ ψ0(x6).
Starting with Config(1,1), we see that the entanglement dynamics for 𝒩1,2(τ)|ω=2, 𝒩3,4(τ)|ω=3, and 𝒩5,6 (τ)|ω=5 all resemble damped sinusoidal functions and undulate independently from each other which coincides with oscillator pairs possessing different frequencies. The dampening for oscillator pair 3,4 or 5,6 happens quicker than for oscillator pair 1,2 because their frequencies are closer in resonance to the central frequency, Ω = 4. The entanglement within each of these oscillator pairs do eventually reach saturation due to Ω > 0 and not being too far off resonance, but simply has more undulations and takes a longer duration to do so, which echos the results found in Figure 7.
If the frequency of an oscillator pair is resonant with the central frequency and it had no initial entanglement, then its entanglement dynamics resembles the later half a logistic function. This is shown in Config(l, 2), where 𝒩1,2 (τ)|ω=4 increases in value and then saturates. Meanwhile, 𝒩3,4 (τ)|ω=3 and 𝒩5,6 (τ)|ω=5 are identical to Config(1, 1). Additionally, we again find that possessing a higher frequency can result in more entanglement generation than lower ones. This is evidenced by Config(1,1) and Config(1,3), where the Blue curve is larger than the Red curve, and by Config(1,2), where the Blue curve is larger than the Green curve after saturation. It is also shown when comparing 𝒩3,4(τ)|ω=4 in Config(1,3) to 𝒩3,4(τ)|ω=3 in Config(1,1). These results are consistent with our findings in Figure 8 and Figure 9.
In Config(1, 3), we find that if an oscillator pair had initial entanglement and is resonant, i.e. 𝒩3,4(τ)|ω=4 now shown as a Green solid curve, then its dynamics resembles the later half of a inverted logistic function. Also regardless of resonance, possessing initial entanglement does not allow for entanglement generation involving it with any other oscillator(s) in general, i.e. 𝒩3,j(τ)|ω = 𝒩4,j(τ)|ω = 0 where j ∈ {1,2,5,6}. Thus, these curves are not shown in all the plots nor do we see connections involving oscillator 3 and 4 with any other oscillator in all configurations. These results seem to support entanglement monogamy and warrants future work. However, in order to have a strong proof of genuine entanglement monogamy, we need to consider multi-partite entanglement [49,50].
In Config(2,1), 𝒩3,4(τ)|ω=3 behaves just like in Config(1,1) and Config(1,2). However, we now have identical entanglement dynamics for 𝒩1,2(τ)|ω=2, 𝒩1,5(τ)|ω=2, and 𝒩2,5(τ)|ω=2 due to them all having the same frequency. They are all shown as a single Red dashed curve and behave similar to 𝒩1,2(τ)|ω=2 in Config(1,1) but much less pronounced. The lower overall magnitude of entanglement generated within these three pairs suggests that it depends on number of oscillators that share the same frequency or the number of pairwise combinations. This reflects the findings found in Figure 8 that suggested an interdependency among the oscillators. Crucially, if an oscillator has an off-resonant frequency with the others, then no entanglement can be generated involving it. This is shown with oscillator 6 colored in Blue and standing alone, unable to link with any other oscillators. In its corresponding plot, all possible pairwise combinations involving oscillator 6 are shown as a Blue dashed curve and they all have zero entanglement generated within them, i.e. 𝒩6,i|ω=5 = 0 where i ∈ {1, …,5}.
Likewise in Config(2,2), we have identical entanglement dynamics to Config(2,1), but 𝒩1,2(τ)|ω=4, 𝒩1,5(τ)|ω=4, and𝒩2,5(τ)|ω=4 are now resonant with the central frequency and are shown as a single Green dashed curve. They behave similar to 𝒩1,2(τ)|ω=4 in Config(1,2), but with a lower saturation value likely because of the same reason when we had compared the Red dashed curve in Config(2,1) to the one in Config(1,1).
In the last configuration, Config(2,3), all oscillator pairs are resonant with the central frequency. This results in identical entanglement generation within six oscillator pairs, i.e. 𝒩1,2(τ)|ω=4, 𝒩1,5(τ)|ω=4, 𝒩1,6(τ)|ω=4, 𝒩2,5(τ)|ω=4, 𝒩2,6(τ)|ω=4, and 𝒩5,6(τ)|ω=4 and are all shown as a single Green dashed curve. Their saturation value is slightly less than𝒩1,2(τ)|ω=4, 𝒩1,5(τ)|ω=4, and𝒩2,5 (τ)|ω=4 in Config(2,2). Meanwhile, 𝒩3,4(τ)|ω=4 is shown as a Green solid curve and we see that its saturation value is slightly less than 𝒩3,4(τ)|ω=4 in Config(1,3) despite having the same frequency. It also is approximately equal to 𝒩3,4(τ)|ω=3 in Config(2,2) despite being a higher frequency. These results again show an interdependency among oscillators depending upon their initial states.
In summary, if the frequencies of oscillator pairs are resonant with the central frequency, regardless of initial entanglement, then those oscillator pairs will reach entanglement saturation without any extra oscillations resembling the results found in Figure 7, Figure 8, and Figure 9. Furthermore, if we refer to the resonant cases for pairs of oscillators that did not have initial entanglement, which are the Green dashed curves in Config(1,2), Config(2,2), and Config (2,3), specifically in that order, we see that the magnitude of entanglement saturation decreases when there are more pairwise combinations. This further supports how entanglement dynamics depends on the initial states of oscillators and the number of them as discussed in Figure 8. Furthermore, if two or more oscillators share the same frequency and is approximately resonant with the central frequency, then entanglement will be generated within all possible pairwise combinations of them, although they will have extra undulations before saturating and be lesser in magnitude compared to if they were exactly resonant.
Lastly, an intuitive reason as to why changing the frequency of the oscillators affects the entanglement dynamics is because it changes the effective coupling between the system and environment. This is due to each oscillator interacting predominantly with the bath mode they are resonant with in a non-Markovian environment, i.e., a bath with a non-flat spectrum. An analytical reason is warranted, but we leave this for future work.
Harmonic oscillator systems are among the most extensively studied models in physics. We have presented a detailed analysis of a ring of oscillators coupled to a noisy environment by numerically simulating the entanglement dynamics in a Markovian and a non-Markovian environment. This included how to manipulate entanglement evolution by varying physical parameters such as coupling, memory times, and resonant conditions. Through this analysis, we have shown that passive entanglement storage may be enhanced with the memory effect and resonance.
We also showed how to passively distribute entanglement within a small group of independent oscillators to establish a “hidden quantum network”. Specifically, we showed that initially unentangled oscillator pairs may develop entanglement without cross-entangling facilities provided they all share a common non-Markovian environment. To not overstate our results about scalability, we emphasize that our presented results only apply for system sizes where N ≤ 6. However, given that our analytical set of differential equations is solved for 𝒩 oscillators, future work involving system sizes with N > 6 can be readily carried out. Additionally, we show that resonance may maximize and help to stabilize the shared entanglement. These results suggest a scheme for hidden quantum networks, where participating parties may distribute information via a common environment.
In conclusion, non-Markovian environments and resonance may be useful for enhancing long-lived entanglement, facilitating hidden entanglement networks, and potentially enabling secure quantum information processing.