# Experimental optical simulator of
reconfigurable and complex quantum environment

P. Renault1,∗, J. Nokkala2, G. Roeland 1, N. Y. Joly 5,6,7,
R. Zambrini3, S. Maniscalco4, J. Piilo2 N. Treps1, V. Parigi1
1Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université PSL,
Collège de France, 75005 Paris
2Department of Physics and Astronomy, University of Turku,
FI-20014 Turun yliopisto, Finland
3Instituto de F́ısica Interdisciplinar y Sistemas Complejos IFISC (CSIC-UIB),
Campus Universitat Illes Balears,
E-07122 Palma de Mallorca, Spain
4QTF Centre of Excellence, Department of Physics, University of Helsinki, P.O.
Box 43, FI-00014 Helsinki, Finland,
5 Department of Physics, Friedrich-Alexander-Universität, Staudtstraße 2,
91058 Erlangen, Germany
6Max Planck Institute for the Science of Light, Staudtstraße 2, 91058 Erlangen
Germany
7Interdisciplinary Centre for Nanostructured Films, Cauerstr.3, 91058
Erlangen, Germany
∗ current address: Department of Physics, University of Virginia, 382
McCormick Rd, Charlottesville, 22904-4714, VA, USA

###### 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.

##  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]:

$H_{E}=\bm{p}^{T}\bm{\Delta_{\omega}}\bm{p}+\bm{q}^{T}\sqrt{\bm{\Delta_{\omega}}^{-1}}\bm{A}\sqrt{\bm{\Delta_{\omega}}^{-1}}\bm{q};\quad
H_{S}=\omega_{S}\left(\frac{p_{S}^{2}}{2}+\frac{q_{S}^{2}}{2}\right);\quad
H_{I}=kq_{S}q_{l}$ (1)

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.

Figure 1: Structure of the global system and environments. A: Scheme of a
global system where a probe (green) is attached to an environment network
(blue) via a coupled strength $k$ (orange). B-D: Environment structures
investigated in this work. B: linear structure used for the network 1, 2 and
3. These are networks of 16 nodes with periodic coupling strength such as
network 1: $g_{i,i+1}=\\{0.1;0.05;0.1;...\\}$; network 2:
$g_{i,i+1}=\\{0.1;0.1;0.05;0.1;...\\}$; network 3:
$g_{i,i+1}=\\{0.1;0.05;0.025;0.1;...\\}$. The coupling strengths between the
probe and these three environments are the same: $k_{1-2-3}=0.01$. C and D
Structure that represents the environments called respectively network 4 and
network 5 and both composed by $N=50$ nodes. They are two instances of the
Watts-Strogatz (WS) and Barabasi-Albert (BA) model as defined in complex
networks theory [57, 58]. The coupling strengths are constant within these two
networks: $g_{i,i+1}=0.08$ for network 4 and $g_{i,i+1}=0.02$ for network 5
and we have $k_{4}=0.02$ and $k_{5}=0.004$. For each network, the nodes have
the same frequency $\omega_{0}=0.25$.

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

$\bm{x}(t)=\mathscr{S}(t)\bm{x}(0).$ (2)

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:

$\begin{split}\bm{x}(t)&=\mathscr{S}_{eff}(t)\bm{x}_{v}(0)\\\
&=\bm{R}_{1}(t)\bm{\Delta}(t)\bm{R}_{2}(t)\bm{x}_{v}(0)=\bm{R}_{1}(t)\bm{\Delta}(t)\bm{x}_{v}(0).\end{split}$
(3)

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

Figure 2: Protocol for the recovery of the spectral density function for a
given environment. The parameters of a given network are used to calculate the
quadrature dynamics of the system $S$ for different frequencies $\omega_{S}$
while the decomposition of such evolution gives the measurement setting,
$\bm{R}_{1}(\omega_{S})$ is tied to a specific mask for the spatial light
modulator (SLM) in the shaping of the local oscillator (LO). Measurement of
$\langle\hat{q}_{S}^{2}\rangle$ and $\langle\hat{p}_{S}^{2}\rangle$ via the
quadrature statistics are associated to the $\langle\hat{n}_{S}\rangle$ value
and then to $J(\omega_{S})$.

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:

$\ddot{q_{S}}+\tilde{\omega}_{S}^{2}q_{S}+\int_{0}^{t}d\tau\gamma(t-\tau)\dot{q_{S}}=\xi(t),$
(4)

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:

$J(\omega_{S})=\omega_{S}\int_{0}^{t_{max}}\mathrm{d}t\gamma(t)\cos(\omega_{S}t).$
(5)

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.

Figure 3: Spectral density measurements. Spectral density as a function of the
frequency $\omega_{s}$ when the environment has a linear structure (A) or a
complex structure (B). Orange lines: theoretical calculation; blue dots:
experimental values recovered by probing the optical mode simulating the open
system. Red vertical lines show the probe frequencies of network 1 and 4 for
which QNM is investigated.

###  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:

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

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 $0$ 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.

Figure 4:  Fidelity as a function of the time parameter for several probe frequencies. (left) The environment takes the form of network 1 and pink and yellow dots are respectively for $\omega_{S}=0.7$ and $\omega_{S}=0.58$. (right) The environment takes the form of the Watts-Strogatz network, pink dots are measured for $\omega_{S}=0.9$, red dots for $\omega_{S}=0.75$ and the yellow dots for $\omega_{S}=0.4$. Solid lines are obtained via averaging over 50 points.  $\omega_{S}$ |  $\mathcal{N}$ of network 1 |  $\mathcal{N}$ of WS network
---|---|---
0.4 | - | 0.013
0.58 | 0.041 | -
0.7 | 0.012 | -
0.75 | - | 0.06
0.9 | - | 0.012
Table 1: Non Markovian witnesses for network 1 and 4.

##  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.

Figure 5: Scheme of the experimental setup. The beam from the femto-second
laser source is divided in two paths: the first one to produce multimode
squeezed light {second-harmonic generation + SPOPO cavity} and the second one
to spectrally shape the local oscillator {PCF + pulse shaper}. The two beams
are then recombined for the homodyne detection. Inset: electron micrograph of
the end face of the PCF.

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.

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
Table 2: Mapping of experimental implementation for the quantum network for
open quantum system.

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]

$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{)},$ (7)

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.

## References

  * [1] J. Preskill, “Quantum Computing in the NISQ era and beyond,” Quantum, vol. 2, p. 79, Aug. 2018.
  * [2] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, et al., “Quantum computational advantage using photons,” Science, vol. 370, no. 6523, pp. 1460–1463, 2020.
  * [3] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, et al., “Quantum supremacy using a programmable superconducting processor,” Nature, vol. 574, no. 7779, pp. 505–510, 2019.
  * [4] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, et al., “Strong quantum computational advantage using a superconducting quantum processor,” Physical review letters, vol. 127, no. 18, p. 180501, 2021.
  * [5] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, et al., “Quantum computational advantage with a programmable photonic processor,” Nature, vol. 606, no. 7912, pp. 75–81, 2022.
  * [6] C. Berdou, A. Murani, U. Reglade, W. Smith, M. Villiers, J. Palomo, M. Rosticher, A. Denis, P. Morfin, M. Delbecq, et al., “One hundred second bit-flip time in a two-photon dissipative oscillator,” arXiv preprint arXiv:2204.09128, 2022.
  * [7] R. Lescanne, M. Villiers, T. Peronnin, A. Sarlette, M. Delbecq, B. Huard, T. Kontos, M. Mirrahimi, and Z. Leghtas, “Exponential suppression of bit-flips in a qubit encoded in an oscillator,” Nature Physics, vol. 16, no. 5, pp. 509–513, 2020.
  * [8] M. J. Biercuk, H. Uys, A. P. VanDevender, N. Shiga, W. M. Itano, and J. J. Bollinger, “Optimized dynamical decoupling in a model quantum memory,” Nature, vol. 458, no. 7241, pp. 996–1000, 2009.
  * [9] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, “Quantum computation and quantum-state engineering driven by dissipation,” Nature physics, vol. 5, no. 9, pp. 633–636, 2009.
  * [10] J. T. Barreiro, P. Schindler, O. Gühne, T. Monz, M. Chwalla, C. F. Roos, M. Hennrich, and R. Blatt, “Experimental multiparticle entanglement dynamics induced by decoherence,” Nature Physics, vol. 6, no. 12, pp. 943–946, 2010\.
  * [11] V. V. Albert, B. Bradlyn, M. Fraas, and L. Jiang, “Geometry and response of lindbladians,” Physical Review X, vol. 6, no. 4, p. 041031, 2016.
  * [12] C. P. Koch, “Controlling open quantum systems: tools, achievements, and limitations,” Journal of Physics: Condensed Matter, vol. 28, p. 213001, may 2016.
  * [13] B.-H. Liu, X.-M. Hu, Y.-F. Huang, C.-F. Li, G.-C. Guo, A. Karlsson, E.-M. Laine, S. Maniscalco, C. Macchiavello, and J. Piilo, “Efficient superdense coding in the presence of non-markovian noise,” EPL (Europhysics Letters), vol. 114, no. 1, p. 10005, 2016.
  * [14] Z.-D. Liu, O. Siltanen, T. Kuusela, R.-H. Miao, C.-X. Ning, C.-F. Li, G.-C. Guo, and J. Piilo, “Efficient quantum teleportation under noise with hybrid entanglement and reverse decoherence,” arXiv preprint arXiv:2210.14935, 2022.
  * [15] A. Mattioni, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, “Design principles for long-range energy transfer at room temperature,” Phys. Rev. X, vol. 11, p. 041003, Oct 2021.
  * [16] F. Caycedo-Soler, A. Mattioni, J. Lim, T. Renger, S. Huelga, and M. Plenio, “Exact simulation of pigment-protein complexes unveils vibronic renormalization of electronic parameters in ultrafast spectroscopy,” Nature Communications, vol. 13, no. 1, pp. 1–8, 2022.
  * [17] A. Nüßeler, D. Tamascelli, A. Smirne, J. Lim, S. F. Huelga, and M. B. Plenio, “Fingerprint and universal markovian closure of structured bosonic environments,” Phys. Rev. Lett., vol. 129, p. 140604, Sep 2022.
  * [18] G. Manzano, G.-L. Giorgi, R. Fazio, and R. Zambrini, “Boosting the performance of small autonomous refrigerators via common environmental effects,” New Journal of Physics, vol. 21, p. 123026, dec 2019.
  * [19] M. Kloc, K. Meier, K. Hadjikyriakos, and G. Schaller, “Superradiant many-qubit absorption refrigerator,” Phys. Rev. Applied, vol. 16, p. 044061, Oct 2021\.
  * [20] A. Sannia, R. Martínez-Peña, M. C. Soriano, G. L. Giorgi, and R. Zambrini, “Dissipation as a resource for quantum reservoir computing,” arXiv preprint arXiv:2212.12078, 2022.
  * [21] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, “Spectral theory of liouvillians for dissipative phase transitions,” Physical Review A, vol. 98, no. 4, p. 042118, 2018.
  * [22] G. Manzano, F. Galve, G. L. Giorgi, E. Hernández-García, and R. Zambrini, “Synchronization, quantum correlations and entanglement in oscillator networks,” Scientific Reports, vol. 3, p. 1439, Mar. 2013.
  * [23] A. Cabot, G. L. Giorgi, F. Galve, and R. Zambrini, “Quantum synchronization in dimer atomic lattices,” Phys. Rev. Lett., vol. 123, p. 023604, Jul 2019\.
  * [24] F. Galve, R. Zambrini, and S. Maniscalco, “Non-Markovianity hinders Quantum Darwinism,” Scientific Reports, vol. 6, p. 19607, Jan. 2016\.
  * [25] T. P. Le and A. Olaya-Castro, “Strong quantum darwinism and strong independence are equivalent to spectrum broadcast structure,” Phys. Rev. Lett., vol. 122, p. 010403, Jan 2019.
  * [26] C. Foti, T. Heinosaari, S. Maniscalco, and P. Verrucchi, “Whenever a quantum environment emerges as a classical system, it behaves like a measuring apparatus,” Quantum, vol. 3, p. 179, Aug. 2019.
  * [27] J. Nokkala, F. Arzani, F. Galve, R. Zambrini, S. Maniscalco, J. Piilo, N. Treps, and V. Parigi, “Reconfigurable optical implementation of quantum complex networks,” New Journal of Physics, vol. 20, no. 5, p. 053024, 2018\.
  * [28] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. Kimble, and D. E. Chang, “Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays,” Physical Review X, vol. 7, no. 3, p. 031024, 2017.
  * [29] M. Bello, G. Platero, J. I. Cirac, and A. González-Tudela, “Unconventional quantum optics in topological waveguide qed,” Science advances, vol. 5, no. 7, p. eaaw0297, 2019.
  * [30] Y. Cai, J. Roslund, G. Ferrini, F. Arzani, X. Xu, C. Fabre, and N. Treps, “Multimode entanglement in reconfigurable graph states using optical frequency combs,” Nature communications, vol. 8, no. 1, pp. 1–9, 2017\.
  * [31] J. Roslund, R. M. De Araujo, S. Jiang, C. Fabre, and N. Treps, “Wavelength-multiplexed quantum networks with ultrafast frequency combs,” Nature Photonics, vol. 8, no. 2, pp. 109–112, 2014.
  * [32] R. Vasile, F. Galve, and R. Zambrini, “Spectral origin of non-markovian open-system dynamics: A finite harmonic model without approximations,” Physical Review A, vol. 89, no. 2, p. 022109, 2014.
  * [33] Á. Rivas, S. F. Huelga, and M. B. Plenio, “Entanglement and non-markovianity of quantum evolutions,” Physical review letters, vol. 105, no. 5, p. 050403, 2010.
  * [34] Á. Rivas, S. F. Huelga, and M. B. Plenio, “Quantum non-markovianity: characterization, quantification and detection,” Reports on Progress in Physics, vol. 77, no. 9, p. 094001, 2014.
  * [35] H.-P. Breuer, E.-M. Laine, and J. Piilo, “Measure for the degree of non-markovian behavior of quantum processes in open systems,” Physical review letters, vol. 103, no. 21, p. 210401, 2009.
  * [36] L. Li, M. J. Hall, and H. M. Wiseman, “Concepts of quantum non-markovianity: A hierarchy,” Physics Reports, vol. 759, pp. 1–51, 2018.
  * [37] E.-M. Laine, J. Piilo, and H.-P. Breuer, “Measure for the non-markovianity of quantum processes,” Physical Review A, vol. 81, no. 6, p. 062115, 2010\.
  * [38] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, “Colloquium: Non-markovian dynamics in open quantum systems,” Reviews of Modern Physics, vol. 88, no. 2, p. 021002, 2016.
  * [39] B.-H. Liu, L. Li, Y.-F. Huang, C.-F. Li, G.-C. Guo, E.-M. Laine, H.-P. Breuer, and J. Piilo, “Experimental control of the transition from markovian to non-markovian dynamics of open quantum systems,” Nature Physics, vol. 7, no. 12, pp. 931–934, 2011.
  * [40] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, “An open-system quantum simulator with trapped ions,” Nature, vol. 470, pp. 486–491, Feb. 2011\.
  * [41] S. Cialdi, M. A. C. Rossi, C. Benedetti, B. Vacchini, D. Tamascelli, S. Olivares, and M. G. A. Paris, “All-optical quantum simulator of qubit noisy channels,” Applied Physics Letters, vol. 110, no. 8, p. 081107, 2017\.
  * [42] S. Yu, Y.-T. Wang, Z.-J. Ke, W. Liu, Y. Meng, Z.-P. Li, W.-H. Zhang, G. Chen, J.-S. Tang, C.-F. Li, and G.-C. Guo, “Experimental investigation of spectra of dynamical maps and their relation to non-markovianity,” Phys. Rev. Lett., vol. 120, p. 060406, Feb 2018.
  * [43] Z.-D. Liu, H. Lyyra, Y.-N. Sun, B.-H. Liu, C.-F. Li, G.-C. Guo, S. Maniscalco, and J. Piilo, “Experimental implementation of fully controlled dephasing dynamics and synthetic spectral densities,” Nature Communications, vol. 9, p. 3453, Aug. 2018.
  * [44] G. García-Pérez, M. A. C. Rossi, and S. Maniscalco, “IBM Q Experience as a versatile experimental testbed for simulating open quantum systems,” npj Quantum Information, vol. 6, p. 1, Jan. 2020.
  * [45] S. Gröblacher, A. Trubarov, N. Prigge, G. D. Cole, M. Aspelmeyer, and J. Eisert, “Observation of non-Markovian micromechanical Brownian motion,” Nature Communications, vol. 6, p. 7606, July 2015.
  * [46] P. M. Harrington, E. J. Mueller, and K. W. Murch, “Engineered dissipation for quantum information science,” Nature Reviews Physics, vol. 4, no. 10, pp. 660–671, 2022.
  * [47] R. Menu, J. Langbehn, C. P. Koch, and G. Morigi, “Reservoir-engineering shortcuts to adiabaticity,” Phys. Rev. Res., vol. 4, p. 033005, Jul 2022\.
  * [48] S. Davidson, F. A. Pollock, and E. Gauger, “Eliminating radiative losses in long-range exciton transport,” PRX Quantum, vol. 3, p. 020354, Jun 2022\.
  * [49] C.-F. Li, G.-C. Guo, and J. Piilo, “Non-markovian quantum dynamics: What is it good for?,” Europhysics Letters, vol. 128, p. 30001, jan 2020.
  * [50] K. Head-Marsden, S. Krastanov, D. A. Mazziotti, and P. Narang, “Capturing non-markovian dynamics on near-term quantum computers,” Phys. Rev. Res., vol. 3, p. 013182, Feb 2021.
  * [51] M. Gluza, J. a. Sabino, N. H. Ng, G. Vitagliano, M. Pezzutto, Y. Omar, I. Mazets, M. Huber, J. Schmiedmayer, and J. Eisert, “Quantum field thermal machines,” PRX Quantum, vol. 2, p. 030310, Jul 2021.
  * [52] Y. Shirai, K. Hashimoto, R. Tezuka, C. Uchiyama, and N. Hatano, “Non-markovian effect on quantum otto engine: Role of system-reservoir interaction,” Phys. Rev. Res., vol. 3, p. 023078, Apr 2021.
  * [53] K. G. Paulson, E. Panwar, S. Banerjee, and R. Srikanth, “Hierarchy of quantum correlations under non-markovian dynamics,” Quantum Information Processing, vol. 20, no. 4, p. 141, 2021.
  * [54] M. Carrega, L. M. Cangemi, G. De Filippis, V. Cataudella, G. Benenti, and M. Sassetti, “Engineering dynamical couplings for quantum thermodynamic tasks,” PRX Quantum, vol. 3, p. 010323, Feb 2022.
  * [55] G. Spaventa, S. F. Huelga, and M. B. Plenio, “Capacity of non-markovianity to boost the efficiency of molecular switches,” Phys. Rev. A, vol. 105, p. 012420, Jan 2022.
  * [56] J. Nokkala, F. Galve, R. Zambrini, S. Maniscalco, and J. Piilo, “Complex quantum networks as structured environments: engineering and probing,” Scientific Reports, vol. 6, p. 26861, May 2016.
  * [57] A. L. Barabási, Networks science.  Cambridge University Press, 2016.
  * [58] M. E. J. Newman, Networks, second edition.  Oxford University Press, 2018.
  * [59] C. Bloch and A. Messiah, “The canonical form of an antisymmetric tensor and its application to the theory of superconductivity,” Nuclear Physics, vol. 39, pp. 95–106, 1962.
  * [60] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’networks,” nature, vol. 393, no. 6684, pp. 440–442, 1998\.
  * [61] R. Albert and A.-L. Barabási, “Statistical mechanics of complex networks,” Rev. Mod. Phys., vol. 74, pp. 47–97, Jan 2002.
  * [62] C. Gardiner and P. Zoller, Quantum noise.  Springer Berlin, 2004.
  * [63] H.-P. Breuer, F. Petruccione, et al., The theory of open quantum systems.  Oxford University Press on Demand, 2002.
  * [64] U. Weiss, Quantum Dissipative Systems.  World Scientific, Singapore, 2008.
  * [65] F. Mascherpa, A. Smirne, A. D. Somoza, P. Fernández-Acebal, S. Donadi, D. Tamascelli, S. F. Huelga, and M. B. Plenio, “Optimized auxiliary oscillators for the simulation of general open quantum systems,” Physical Review A, vol. 101, no. 5, p. 052108, 2020.
  * [66] R. Vasile, S. Maniscalco, M. G. Paris, H.-P. Breuer, and J. Piilo, “Quantifying non-markovianity of continuous-variable gaussian dynamical maps,” Physical Review A, vol. 84, no. 5, p. 052118, 2011.
  * [67] T. Kouadou, F. Sansavini, A. Matthieu, H. Johan, N. Treps, and V. Parigi, “Spectrally shaped and pulse-by-pulse multiplexed multimode squeezed states of light,” Preprint, vol. arXiv:2209.10678, 2022.
  * [68] Y.-S. Ra, A. Dufour, M. Walschaers, C. Jacquard, T. Michel, C. Fabre, and N. Treps, “Non-Gaussian quantum states of a multimode light field,” Nature Physics, vol. 16, pp. 144–147, Feb. 2020.
  * [69] M. Walschaers, N. Treps, B. Sundar, L. D. Carr, and V. Parigi, “Emergent complex quantum networks in continuous-variables non-gaussian states,” arXiv preprint arXiv:2012.15608, 2020.
  * [70] G. Ferrini, J. Roslund, F. Arzani, Y. Cai, C. Fabre, and N. Treps, “Optimization of networks for measurement-based quantum computation,” Physical Review A, vol. 91, no. 3, p. 032314, 2015.
  * [71] A. Monmayrant, S. Weber, and B. Chatel, “A newcomer’s guide to ultrashort pulse shaping and characterization,” Journal of Physics B: Atomic, Molecular and Optical Physics, vol. 43, no. 10, p. 103001, 2010.
  * [72] P. Abdolghader, A. F. Pegoraro, N. Y. Joly, A. Ridsdale, R. Lausten, F. Légaré, and A. Stolow, “All normal dispersion nonlinear fibre supercontinuum source characterization and application in hyperspectral stimulated raman scattering microscopy,” Optics Express, vol. 28, no. 24, pp. 35997–36008, 2020.
  * [73] F. Arzani, C. Fabre, and N. Treps, “Versatile engineering of multimode squeezed states by optimizing the pump spectral profile in spontaneous parametric down-conversion,” Phys. Rev. A, vol. 97, p. 033808, Mar 2018\.

## 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:

  1. 1.

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

$F=\frac{2}{\sqrt{2(1+\cosh{2r_{1}}\cosh{2r_{2}}-\cos{\varphi_{0}}\sinh{2r_{1}}\sinh{2r_{2}})}}.$
(8)

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.

$F(t)\approx\frac{2}{\sqrt{2(1+\cosh{2r_{1}(t)}\cosh{2r_{2}(t)}-\cos{\varphi_{0}}\sinh{2r_{1}(t)}\sinh{2r_{2}(t)})}}.$
(9)

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,

$F(t)\approx\frac{2}{\sqrt{\cosh^{2}(r_{1}(t)-r_{2}(t))}},\quad
F(t)\approx\frac{2}{\sqrt{\cosh^{2}(r_{1}(t)+r_{2}(t))}}.$ (10)

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

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

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.

Figure 6: Protocol for the determination of the fidelity and consequently the
BLP witness as a function of time. Here, two types of computed spectrum are
sent to the SLM: one related to state $\rho_{1sq}$ (green) and another one
related to state $\rho_{2sq\perp}$ (blue). The tomography of state, related to
a given time $t$, is performed by the homodyne detection and is then
determined the corresponding fidelity.

Figure 7: Analytical simulation of damping kernel as a function of time for
all networks. For each network implemented as environment of the probe, the
resulting damping kernel exhibits the same shape: first, it is observed a
decay of the oscillations amplitude until reaching then a more or less long
plateau where is set $t_{max}$. Let’s notice that for the linear networks the
plateau is more easily distinguishable than for the complex network
environment. The red lines show the position of $t_{max}$ selected for the
simulation of the different total system dynamics.