Experimental optical simulator of reconfigurable and complex quantum environment

TL;DR

This paper presents an optical simulator capable of mimicking complex, reconfigurable quantum environments, enabling experimental exploration of open quantum system dynamics and their applications in quantum information and thermodynamics.

Contribution

The authors develop and experimentally demonstrate an optical platform that simulates arbitrary, reconfigurable quantum environments using continuous variables, advancing the study of open quantum systems.

Findings

Successfully simulated spectral density and non-Markovianity effects

Reconfigurable environment allows testing various open quantum system scenarios

Platform enables experimental investigation of quantum information processes

Abstract

No quantum system can be considered totally isolated from its environment. In most cases the interaction between the system of interest and the external degrees of freedom deeply changes its dynamics, as described by open quantum system theory. Nevertheless, engineered environment can be turned into beneficial effects for some quantum information tasks. Here we demonstrate an optical simulator of a quantum system coupled to an arbitrary and reconfigurable environment built as a complex network of quantum interacting systems. We experimentally retrieve typical features of open quantum system dynamics like the spectral density and quantum non-Markovianity, by exploiting squeezing and entanglement correlation of a continuous variables optical platform. This opens the way to the experimental tests of open quantum systems in reconfigurable environments that are relevant in, among others, quantum information, quantum thermodynamics, quantum transport and quantum synchronization.

Figures

Tables

Table 1. Table 1: Non Markovian witnesses for network 1 and 4.

${\omega }_S$	$𝒩$ of network 1	$𝒩$ of WS network
0.4	-	0.013
0.58	0.041	-
0.7	0.012	-
0.75	-	0.06
0.9	-	0.012

Table 2. Table 2: Mapping of experimental implementation for the quantum network for open quantum system.

Network component	Quantum network	Experimental implementation
Node	Quantum harmonic oscillator	Optical mode
Link	Coupling strength	Entanglement/basis change
Addressing a node	Local measurement	Pulse shaping and projective measurement

Equations

$$H_E={\mathbf{p}}^T{\mathbf{\Delta }}_{\mathbf{\omega }}\mathbf{p}+{\mathbf{q}}^T\sqrt{{{\mathbf{\Delta }}_{\mathbf{\omega }}}^{-1}}\mathbf{A}\sqrt{{{\mathbf{\Delta }}_{\mathbf{\omega }}}^{-1}}\mathbf{q};H_S={\omega }_S\left(\frac{p_S^2}{2}+\frac{q_S^2}{2}\right);H_I=k{q}_S{q}_l$$

$$\mathbf{x}(t)=\mathcal{S}(t)\mathbf{x}(0).$$

$$\begin{array}{rl}\mathbf{x}(t)& ={\mathcal{S}}_{eff}(t){\mathbf{x}}_v(0)\\ & ={\mathbf{R}}_1(t)\mathbf{\Delta }(t){\mathbf{R}}_2(t){\mathbf{x}}_v(0)={\mathbf{R}}_1(t)\mathbf{\Delta }(t){\mathbf{x}}_v(0).\end{array}$$

$$\ddot{q_S}+{\tilde{\omega }}_S^2{q}_S+{\int }_0^{t}d\tau \gamma (t-\tau )\dot{q_S}=\xi (t),$$

$$J({\omega }_S)={\omega }_S{\int }_0^{t_{max}}\mathrm{d}t\gamma (t)\cos ({\omega }_{S}t).$$

$$\mathcal{N}=\max_{\rho_{1},\rho_{2}}\int_{\frac{\partial F}{\partial t}<0}\mathrm{d}t\Big{(}-\frac{\partial F}{\partial t}\Big{)},$$

$$J(\omega_{s})=\dfrac{\omega_{s}}{t_{max}}\ln\Big{(}\dfrac{N(\omega_{s})-\langle n_{S}(0)\rangle}{N(\omega_{s})-\langle n_{S}(t_{max})\rangle}\Big{)},$$

$$F=\frac{2}{\sqrt{2(1+\cosh 2r_1\cosh 2r_2-\cos {\phi }_0\sinh 2r_1\sinh 2r_2)}}.$$

$$F(t)\approx \frac{2}{\sqrt{2(1+\cosh 2r_1(t)\cosh 2r_2(t)-\cos {\phi }_0\sinh 2r_1(t)\sinh 2r_2(t))}}.$$

$$F(t)\approx \frac{2}{\sqrt{{\cosh }^2(r_1(t)-r_2(t))}},F(t)\approx \frac{2}{\sqrt{{\cosh }^2(r_1(t)+r_2(t))}}.$$

$$F(t)\approx \frac{2}{\sqrt{\cosh (4r(t))(1-\cos {\phi }_0)+\cos {\phi }_0+3}},F(t)\approx \mathrm{sech}2r(t).$$

Full text

1 Introduction

Quantum information technologies are nowadays getting to the regime where noisy intermediate scale quantum systems (NISQ) [1] show quantum advantage when compared with classical equivalents [2, 3, 4, 5].

Nevertheless, whatever the platform considered, both decoherence and losses are still obstacles that can be mitigated but not completely avoided [6, 7], so any platform should be considered an open system. The theory of open quantum systems largely explored the role of the environment and showed that, when opportunely engineered, it can be promoted to be an ally of the open system [8, 9, 10, 11, 12, 13, 14].

Beyond quantum information platforms, the study of open quantum systems and a structured environment is essential for understanding biological systems [15, 16, 17], boosting quantum thermal machines [18, 19] and machine learning protocol [20], achieving collective phenomena such as dissipative phase transitions [21] or synchronization in the quantum realm [22, 23] and explaining the emergence of the classical world from quantum constituents [24, 25, 26]. It is thus crucial to have experimental platforms to put to the test the variety of open quantum systems for different purposes. In this work we demonstrate an optical simulator of arbitrary quantum environments interacting with an open system.

The fine-grained structure of the environment can be described in many cases as a network of quantum harmonic oscillators. We experimentally reproduce the dynamics of the open system coupled to networks whose interaction structure can take an arbitrary shape. More generally in this work we consider quantum complex networks, mimicking real-world ones [27], as examples of complex environments [28, 29]. The role of the open system and of the network (environment) are played in the experiment by optical spectral modes interacting in a multimode non-linear process pumped by a femtosecond laser.

We show that the quantum optical platform, by using squeezing and entanglement correlations as resources along with continuous variables (CV) measurement [30, 31], is able to reproduce two crucial features of the system-environment dynamics. The first one is the energy exchange/dissipation, characterized by the spectral density of the environmental coupling [32]. The second property is the quantum non-Markovianity (QNM) [33, 34, 35, 36, 37]. Specifically we test the QNM as introduced by Breuer et al [35], consiting of a back-flow of information from the environment to the system during the dynamics [38].

While there exists number of experiments that have shown the ability to emulate open quantum systems and in particular to control Markovian - non-Markovian transition or dephasing effects [39, 40, 41, 42, 43, 44], these have been mainly implemented via single or few qubits. Only preliminary experimental studies [45] of continuous variable (CV) systems have been reported and without any control of the open system dynamics. Here we take a significant step beyond these regimes by not only studying and controlling the dynamics of open system coupled to a multiparty CV environment but also shaping the fine-grained complex structure of the latter and finally experimentally testing the probing schemes for the environmental properties. Reservoir control and engineering are essential in quantum information tasks [46, 47, 48] and the probed quantities, spectral density and quantum non-Markovianity, are the key features to look at [49, 50, 51, 52, 53, 54, 55].

2 Results

2.1 Mapping of the open quantum system into CV optical systems

The open quantum system we emulate (S, also named the probe) is coupled to an environment (E), as shown in Fig. 1 A. S is a harmonic oscillator of frequency $\omega_{s}$. The environment is modelled by an ensemble of $N$ other harmonic oscillators coupled with each other via spring-like interactions. The coupling strength between oscillators $i$ and $j$ is denoted as $g_{ij}$. Without loss of generality, we consider the case where the frequencies of the oscillators in E are the same, $\omega_{0}$. We assume that S is coupled to only one node of the environmental network labelled $l$ (see Fig. 1). The Hamiltonians of the network ($H_{E}$), of the system ($H_{S}$) and of their interaction ($H_{I}$) are [56, 27]:

[TABLE]

where $\bm{A}$ is a symmetric and real matrix of size $(N)$ called the adjacency matrix of the environmental network, such that $A_{ij}=\delta_{ij}\omega_{i}^{2}/2-(1-\delta_{ij}g_{ij}/2)$, $\bm{q}=(q_{1},...q_{N})^{T}$, $\bm{p}=(p_{1},...p_{N})^{T}$ and $\left(p_{S},q_{S}\right)$ stands for renormalized quadrature operators111With $\left(\bm{q^{\prime}},\bm{p^{\prime}},q_{S}^{\prime},p_{S}^{\prime}\right)$ the usual quadrature operators of N+1 harmonic oscillators, they are defined as $\bm{q}^{T}=\bm{q^{\prime}}^{T}\sqrt{\bm{\Delta_{\omega}}}$; $\bm{p}^{T}=\bm{p^{\prime}}^{T}\sqrt{\bm{\Delta_{\omega}}^{-1}}$; $p_{S}=p_{S}^{\prime}\sqrt{\omega_{S}^{-1}}$; $q_{S}=q_{S}^{\prime}\sqrt{\omega_{S}}$..

Hamiltonians in Eq. (1) are given with $\hbar=1$, and $m=1$. Also, in the rest of the article numerical values of couplings and frequencies are specified relative to a fixed (arbitrary) frequency unit. This is because on the one hand any possible value and unit can be chosen in the simulations implemented via the optical system, and on the other hand the properties of the open system in presence of different environments are driven by the ratio between involved frequency and coupling terms $\{\omega_{S},\omega_{0},g,k\}$ rather than their absolute value.

We can easily derive the temporal dynamic of the full system, described by the evolution matrix $\mathscr{S}(t)$ such as:

[TABLE]

where $\bm{x}(t)=(q_{s},\bm{q},p_{s},\bm{p})$. More details about the derivation of $\mathscr{S}$ can be found in [27]. $H_{E}$,$H_{S}$ and $H_{I}$ being quadratic Hamiltonians leading to Gaussian processes, at each time $t$ the evolution matrix $\mathscr{S}(t)$ is symplectic. To study the energy transfer between the network and the system, some of the involved oscillators should be initialized in a non-zero energy state; such preparation can be included in the symplectic operation, $\mathscr{S}_{eff}(t)=\mathscr{S}(t)\mathscr{S}_{in}$, so that the global process can be written as:

[TABLE]

The second line is the Bloch-Messiah decomposition of the symplectic transformation where $\bm{R}_{1}$ and $\bm{R}_{2}$ are orthogonal matrices, corresponding to linear optics operations (or more generally to basis change), while $\bm{\Delta}$ is a diagonal matrix corresponding to squeezing operations [59]. As the initial state preparation is included in $\mathscr{S}_{eff}$, $\bm{x}_{v}(0)$ are quadratures of a collection of oscillators in the vacuum states and $\bm{R}_{2}$ can be discarded.

We can obtain the same transformation of Eq. (3) on quadratures of optical modes, using multimode squeezing ($\bm{\Delta}$) and mode basis change ($\bm{R}_{1}$). Experimentally, these can be implemented in optical parametric process pumped by an optical frequency comb and measured via ultrafast shaped homodyne detection [31, 30, 27]. We can thus emulate the evolution of the open quantum system at a time $t$ by implementing the transformation $\bm{R}_{1}(t)\bm{\Delta}(t)$ on optical modes of different optical spectrum that play the role of the harmonic oscillators composing the network and the system. Starting from a bunch of initial modes to which the squeezing operation $\bm{\Delta}(t)$ is applied, the linear optics operation $\bm{R}_{1}(t)$ corresponds to a basis change that can be realised by measuring the quadratures $\bm{x}(t)$ via the appropriate local oscillator shape in homodyne detection [30, 27]. This platform can simulate the dynamic of an open quantum system in environments with tunable spectral features, induced by environment. The environment can have any complex structure of correlations, as long as the number of optical modes that can be detected in the mode-selective homodyne - i.e. the number of addressable harmonic oscillators- is large enough. In the following we will show how our experimental system is able to simulate environments of different shape and size by recovering their specific features in the interaction with the open system. In particular we will recover the spectral density $J(\omega)$ [32], that reveals the energy flow between the system and the environment, and the quantum non-Markovian behaviours [35], i.e. an information back-flow from the environment to the system.

In order to show the reconfigurability of our system we implement the networks shown in Fig. 1. We set 3 different linear networks of 16 nodes with different periodic coupling strength [32]. Then we set networks of 50 nodes derived from complex-network models like the Watts-Strogatz model, characterized by short average path lengths and high clustering where we took rewiring probability $p_{WS}=0.1$ [60], and the Barabasi-Albert model characterized by a power law distribution of the degree with connection parameter set to $\kappa=2$ [61].

2.2 Spectral density for a given network

The reduced dynamic of the system in presence of the environment can be derived by tracing out the environmental degrees of freedom from the evolution given by the total Hamiltonian $H_{E}+H_{S}+H_{I}$. In general the evolution of the open system takes the form of a non-unitary master equation for the state, or an equivalent generalized Langevin equation for the position observable, [62, 63, 64, 32, 65] such as:

[TABLE]

where $\tilde{\omega}_{S}$ is a renormalized system frequency and $\xi(t)$ is Langevin forcing of the system. The dissipation and memory effect of the system are featured by the damping kernel $\gamma(t)$. This can be equivalently characterized via the spectral density defined as:

[TABLE]

The value of the spectral density at a given frequency $\omega_{S}$ shows the strength of the energy flow between the system and the environment, e.g. its damping rate. Specific network structures of the environment are characterized by different shapes of $J(\omega_{S})$ that can be easily recovered by calculating the evolution of the system plus the environment via the total Hamiltonian in Eq. (1) and getting $\gamma(t)$ from the reduced dynamics of the system. It should be noted that $J(\omega_{S})$ in Eq. (5) is normally defined with $t_{max}=\infty$ but as here we are dealing with finite environments we set a finite $t_{max}$, which can be considered as the time the open system interacts with all the elements of the networks before the revival dynamics arises due to the finite size effects of the environment (see Methods).

For unknown network structures, the shape of $J(\omega_{S})$ can be recovered by probing the excitation number of the system $S$ that interacts with the environment [32, 56, 27] at different frequencies $\omega_{S}$ (see Eq. (7) in Methods). This is the approach we follow in the experimental simulation: given all the network parameters, the experimental setup composed by different optical modes can implement the quadrature evolution $\bm{x}(t)$ of the network plus the system in Eq. (3). We get the value of $J(\omega_{S})$ by monitoring the excitation number $\langle\hat{n_{S}}\rangle$, that can be recovered from homodyne measurements of $\langle\hat{q_{S}}^{2}\rangle$ and $\langle\hat{p_{S}}^{2}\rangle$ and we compare the results with the expected theoretical shape. The protocol is shown in Fig. 2. The experimental data in Fig. 3 are obtained from homodyne measurements of the mode corresponding to (S) having interacted with (E) until $t_{max}$. Each dot is the value of $J(\omega_{S})$ recovered when the environment interacts with the system at a specific frequency $\omega_{S}$. This corresponds to a given measurement setting, i.e. to a given basis change $\bm{R}_{1}(t_{max},\omega_{S})$ in Eq. (3), and in particular to a given Local Oscillator (LO) spectrum in the homodyne measurement, set by the SLM (Spatial Light Modulator) mask. One dot on each curve is the average of 20 measurements and the error bars are obtained from the standard deviation. The theoretical curves are calculated from Eq. (5) given the parameters of the networks shown in Fig. 1 (see Methods). Except for some noise due to the instabilities of the experimental system, the experimental data shown in Fig. 3 match the shapes of the theoretical curves for all complex environments.

2.3 Quantum non-Markovianity

In addition we simulate quantum Markovian or non-Markovian behaviour of the studied networks for some parameters range. In this work we use the definition of quantum non Marknovianity (QNM) introduced by Breuer et al in [35], where the memory effect in the system is associated to a back flow of information from E to S. A Markovian process continuously tends to reduce the distinguishability between any two quantum states of a given system, in a non-Markovian behavior it tends to increase, so that the flow of information about such distinguishability is reversed from the environment to the system. The original definition of QNM based on trace distance between the states can be expressed, in the case of Gaussian states, via the fidelity as proposed in [35, 66]. Thereforewe use the following QNM witness:

[TABLE]

where $(\rho_{1},\rho_{2})$ is a pair of states of S and $F$ their fidelity. Experimentally only a finite set of quantum states can be accessed, we have then chosen the two experimentally accessible states $\rho_{1},\rho_{2}$ that minimize the fidelity at $t=0$, in order to have the maximal sensitivity in the QNM witness $\mathcal{N}$. Such states are two vacuum squeezed states that are squeezed along two orthogonal directions (see Methods).

The quantum non-Markovianity can be associated to specific structures of the spectral density [32], in particular it has been shown that maximal values of $\mathcal{N}$ are reached at the edges of a band-gap in spectral density, where the band-gap is a region where $J(\omega_{S})$ is close to zero. In this protocol we focus on some given values of the probe frequency where the spectral density of specific networks (linear and WS) have particular features: large values, or being within or at the edge of the gap. For each time $t$ the protocol is performed for two different input states $\rho_{1sq}$ and $\rho_{2sq\perp}$, in order to measure their fidelity. The first state $\rho_{1sq}$ has squeezing of $-1.8$ dB along the $q$ quadrature and antisqueezing of $+2.9$ dB, and the second $\rho_{2sq\perp}$ has squeezing of $-1.3$dB along the $p$ quadrature and antisqueezing of $+2.4$ dB. The fidelity between the two states is calculated via the covariance matrix of the two recovered via the measurement of $\langle\hat{q_{S}}^{2}\rangle$ and $\langle\hat{p_{S}}^{2}\rangle$.

Fig. 4 shows the fidelity measurement as a function of time $t$ for the linear periodic network called network 1 and the Watts-Strogatz network at given values of frequency $\omega_{S}$. For both environments we can observe a back flow of the information in the system for some of the monitored frequencies of the probe. When the environment takes the form of network 1, the system is non-Markovian for $\omega_{S}=0.58$, at the edge of the gap [32], while no information exchange is perceived between the system and the environment for $\omega_{S}=0.70$. At such frequency, the dynamics of the system can be interpreted as unitary, as shown by a value of $J(\omega_{S})$ close to [math] in Fig. 3 A. In the case of the Watts-Strogatz network, non-Markovianity is observed for $\omega_{S}=0.4$ and $\omega_{S}=0.75$ with a larger information back-flow for the latter, while no information exchange is noticeable at $\omega_{S}=0.9$. In order not to overestimate the witness value because of high frequency fluctuations in the experimental data, the derivative of Eq. 6 is evaluated on averaged curves ( solid lines in Fig. 4 ). The obtained values of $\mathcal{N}$ are gathered in the Table 1.

3 Discussion

In summary we have experimentally demonstrated the simulation of an open system coupled to complex network environments of different shape. Any network shape can be engineered and probed in the actual platform. In particular we have shown probing techniques for the spectral density of the environmental coupling and of the quantum non-Markovianity. Our platform is the first experimental setup where continuous variable open systems with engineered environment are tested. It goes beyond the few-qubits implementation by controlling a multipartite system with up to 50 components. The environment and system size can be increased in future experiments by considering both spectral and time multiplexing [67], moreover non-Gaussian interaction can be added [68, 69]. Applications are relevant in the context of quantum information technologies. Dissipation phenomena in energy transfer and in particular vibronic dynamic can be mapped via the demonstrated experimental apparatus [15, 16, 17] opening the way to the test of artificial light-harvesting architecture. Moreover we can test engineered environment to enhance quantum thermal machines [18, 19]. Finally, we can explore different probing schemes like multiple harmonic oscillators (measured modes) coupled with different partitioning of the environment, via weak or strong interaction. We can then test paradigmatic collective phenomena, like quantum phase transitions and quantum synchronization [21, 22, 23] and the emergence of classical world from the quantum one as the effect of the interaction with a structured environment [24, 25, 26].

4 Methods

4.1 The experimental setup

A general scheme of the experiment is shown in Fig. (5). The main goal is to produce highly multimode non classical light. Frequency doubled pulses from a Ti:Sapphire laser are propagated in an OPO cavity with round trip time matched to the pump pulse train cycle time. The non linear crystal inside the cavity is a BiBO crystal of length $0.2$ mm. The repetition rate of the pulse train is $75$MHz and pulses have a duration of about $100$ fs. Correlations appear among many spectral modes of the frequency comb of the down converted light, giving rise to a squeezed vacuum state at $795$ nm central wavelength with a highly multimode structure [31, 30]. This multimode squeezing structure can be measured via a spectrally resolved homodyne detection. A pulse shaper, located upstream, enables to set the spectral mode of the local oscillator in the homodyne detection, which is then mode probed in the detection. We find that the squeezed modes spectra are described by a set of Hermite-Gaussian functions. The first one of the series is a Gaussian function with FWHM of around 6.5 nm.

The quadratures evolution of the involved optical modes in the parametric process is then described by an equation of the same form of the Eq. (3). The harmonic oscillators initially in a vacuum state and whose quadratures are squeezed as $\bm{\Delta}(t)\bm{x}(0)$ are, in the actual setup, the squeezed modes.

The quadratures $\bm{x}(t)=\bm{R}_{1}(t)\bm{\Delta}(t)\bm{x}(0)$ are the ones accessed via homodyne detection. The transformation $\bm{R}_{1}(t)$ describes the basis change from the modes with Hermite-Gauss spectral shapes to the measured modes [31, 30, 70]. In order to exactly match the evolution of the multimode optical system with the one of the network plus the probe, we have to set both the squeezing values in $\bm{\Delta}(t)$ and the detected pulse shapes given by $\bm{R}_{1}(t)$. In the present experiment, numerical analysis showed that the probing of $J(\omega_{S})$ and of $\mathcal{N}$ are not very sensitive to changes in the squeezing values in $\bm{\Delta}(t)$. In these protocols the numbers and values of non-zero diagonal terms in $\bm{\Delta}$ set the number of oscillators that are initially in a not-vacuum state and the value of their excitation numbers. Interactions between the oscillators and their evolution can be established via $\bm{R}_{1}(t)$. Thus only $\bm{R}_{1}(t)$ is modified accordingly to the dynamics we have to simulate. In both protocols (probing of spectral density and quantum non-Markovianity) we follow the spread of energy/information from few initially populated harmonic oscillators to a large number of harmonic oscillators in the networks.

So if the number of the harmonic oscillators that can be initially set in a not-vacuum state is limited by the number of produced squeezed modes, the number of harmonic oscillators that can be reached by some excitation (i.e. the number of the total oscillators in the networks) is limited by the number of spectral modes we can measure, as this is what limit the size of the matrix $\bm{R}_{1}(t)$. A larger number of spectral modes to measure means a larger spectrum to be shaped for the local oscillator field. So in the end the number of the simulated harmonic oscillators depends on the capability of the pulse shaper which is here limited by the optical complexity [71] and the spectral width of the field that is used as local oscillator. In order to not be limited by the latter, the local oscillator, that is derived from the main laser source is broadened with a 2 cm-long all-normal dispersion photonic crystal fiber (PCF), before entering the pulse shaper stage. In such a fiber, the broadening mechanism only relies on self-phase modulation, which is known to be low-noise [72]. If needed in future experiments and protocols, $\bm{\Delta}(t)$ can be opportunely controlled via the shaping of the pump in the parametric process [73].

4.2 Protocol

4.2.1 Spectral density measurement

The protocol is based on the bosonic resonator network mapping established by J. Nokkala et al [27]. In table 2 are gathered the equivalent items involved in the mapping to emulate network of bosonic oscillator.

As shown in Fig. 2, the protocol is carried out in two stages: first, basis change and accordingly the mask characteristics setting the probe measurement are computed by a Mathematica code and second, the variances $\langle\hat{q_{S}}^{2}\rangle$ and $\langle\hat{p_{S}}^{2}\rangle$ are measured. The process is as follows:

• •

Interaction time $t_{max}$, environment structure, coupling strength $g$ and $k$ and the frequency $\omega_{0}$ are set and entered in the code. At first, the damping kernel $\gamma(t)$ is numerically calculated and then we selected as $t_{max}$ a value where the resulting solution is flat and close to 0. The $\gamma$ computation for the considered environments are in the supplementary material.

For sufficiently short times the system cannot resolve the different frequencies of the network, making the spectral density a continuous function of frequency in this regime, as seen in Fig. 3. For some networks the spectral density additionally assumes a constant form for a transient where the shape is not sensitive to small differences in interaction time. The set interaction time for each network are : $t_{max,1}=t_{max,2}=t_{max,3}=150$, $t_{max,4}=90$ and $t_{max,5}=250$.

• •

A set of matrices $\mathscr{S}_{eff}(t_{max},\omega_{s})$ is evaluated for 120 values of $\omega_{s}$ in the range $\{0.2,0.7\}$ for network 1, 2 and 3. Then 100 matrices are also computed in the ranges $\{0.1,1.1\}$ and $\{0.5,0.8\}$ for respectively the network 4 and 5. The matrices $\bm{R}_{1}(t_{max},\omega_{s})$ and $\bm{\Delta}(t_{max},\omega_{s})$ are obtained from the Bloch Messiah decomposition.

• •

From $\bm{R}_{1}(t_{max},\omega_{s})$ we can derive the spectral mode corresponding to the system/probe of frequency $\omega_{s}$ having interacted with the network for the time $t_{max}$. The corresponding optical spectrum of the Local Oscillator is shaped via the SLM masks. The average values $\langle\hat{q_{S}}^{2}\rangle$ and $\langle\hat{p_{S}}^{2}\rangle$ are obtained via homodyne detection.

• •

The average photon number $\langle n_{S}(t_{max})\rangle$ of the system is derived so that we can get $J(\omega_{s})$ from the following equation [27]

[TABLE]

where $N(\omega_{s})=\left(\mathrm{e}^{\omega_{S}/T}-1\right)^{-1}$ is the thermal average boson number with $T$ being the temperature of the environment.

4.2.2 Quantum non-Markovianity

Although the emulated total system dynamic remains unchanged, the way to highlight QNM is slightly different than the way to recover the spectral density function.

• •

Environment structure, coupling strength $g$ and $k$, the frequency $\omega_{0}$ and the probe frequency $\omega_{s}$ are set and entered in the code.

• •

A set of matrices $\mathscr{S}_{eff}(t,\omega_{s})$ is evaluated for 251 values of $t$ in the range $\{0,500\}$. The same set is applied to two different input states for the probe/system oscillator, $\rho_{1sq}$ and $\rho_{2sq\perp}$ consisting in two vacuum states squeezed along two orthogonal directions. The two are naturally encoded in the first two modes $HG_{0}$,$HG_{1}$ of the Hermite-Gauss series that diagonalize the parametric down conversion Hamiltonian, $\rho_{HG_{0}}(0)=\rho_{1sq}$; $\rho_{HG_{1}}(0)=\rho_{2sq\perp}$.

• •

The set of SLM masks corresponding to the temporal evolution of the two initially squeezed oscillators are evaluated.

• •

Homodyne measurements are used for the evaluation of the states $\rho_{HG_{0}}(t)$ and $\rho_{HG_{1}}(t)$ at time t and their fidelity.

Acknowledgements

Funding Statement

This work was supported by the European Research Council under the Consolidator Grant COQCOoN (Grant No. 820079). S.M. acknowledge financial support from the Academy of Finland via the Centre of Excellence program (Project No. 336810 and Project No. 336814). J.N. acknowledges financial support from the Turku Collegium for Science, Medicine and Technology as well as the Academy of Finland under project no. 348854. R.Z. acknowleges funding from the Spanish State Research Agency, through the María de Maeztu project CEX2021-001164-M and the QUARESC project PID2019-109094GB-C21 (AEI /10.13039/501100011033), and CAIB QUAREC project (PRD2018/47).

Competing interests

The authors declare that they have no competing interests.

Supplementary Files

4.3 Witnessing non-Markovianity with pure squeezed states

The used witness $\mathcal{N}$ of quantum non-Markovianity (QNM), specified in Sec. $(2.3)$ of the main manuscript, involves a maximization over pairs of distinct initial states. Experimental quantification will naturally limit the possible pairs to experimentally accessible states, such as pure squeezed states characterized by the magnitude $r_{i}$ and phase $\varphi_{i}$ of squeezing where the index $i\in\{1,2\}$ indicates the two states.

Here plausibility is provided for the following claims: non-Markovianity is witnessed independently of the initial phase difference $\varphi_{0}:=\varphi_{1}-\varphi_{2}$; squeezing opposite quadratures, as we have done, can be expected to give the highest value for $\mathcal{N}$ in this set of states. In particular, if the open system dynamics is found to be (non-)Markovian for one choice of $\varphi_{0}$ it will be found (non-)Markovian for all values $\varphi_{0}\neq 0$; only the numerical value of $\mathcal{N}$ will change.

The key argument is that the choice of $\varphi_{0}$ maximizing distinguishability at $t=0$ maximizes it also for any time $t>0$ whereas the qualitative behavior is independent of $\varphi_{0}$. As will be seen, this is a consequence of the following points, suggested both by intuition and numerical simulations:

the phase difference remains nearly constant: $\varphi(t)\approx\varphi_{0}$; 2. 2.

$r_{i}(t)$ is nearly independent of $\varphi_{i}$; 3. 3.

qualitative behavior of $r_{i}(t)$ is not sensitive to $r_{i}(0)$.

We start from the fidelity $F$ between the two states. It reads

[TABLE]

As anticipated only the relative phase matters. We turn our attention to $F(t)$. Because of point $1.$ we may substitute $\varphi(t)$ with $\varphi_{0}$. Then the squeezing parameters are the only source of non-monotonicity.

[TABLE]

We consider two limiting cases: identical phase $\varphi_{0}=0$ and opposite phase $\varphi_{0}=\pi$. Because of point $2.$ the resulting expressions can be directly compared. We get for the two, respectively,

[TABLE]

As expected, these are the maximal and minimal values of fidelity, respectively, in the interval $\varphi_{0}\in[0,\pi]$ because there $F(t)$ is monotonically decreasing. This follows from the non-negativity of the hyperbolic functions and continuity in $\cos{\varphi_{0}}$ which is monotonous in $[0,\pi]$. Given point $3.$, full contribution of non-monotonicity from both squeezings is achieved only in the latter case and therefore $\mathcal{N}$ is maximized when $\varphi_{0}=\pi$.

In the special case $r_{1}=r_{2}:=r$ point $3.$ becomes unnecessary. Here fidelity in general and in the case $\varphi_{0}=\pi$, respectively, simplifies to

[TABLE]

There are caveats. First of all this is clearly not a proof but rather justifies why our choice for the initial states and the interpretation of the results are reasonable. Second, the equations above implicitly assume pure states. Since we have used a weak coupling and pure state for the network purity can be expected to remain high for the considered times and fidelity should behave as above. This has been checked with numerical simulations.