Reconfigurable optical implementation of quantum
complex networks

J Nokkala, Francesco Arzani, F Galve, R Zambrini, S Maniscalco, J Piilo,
Nicolas Treps, V Parigi

► To cite this version:

J Nokkala, Francesco Arzani, F Galve, R Zambrini, S Maniscalco, et al. Reconfigurable opti- cal implementation of quantum complex networks. New Journal of Physics, 2018, 20, pp.053024. 10.1088/1367-2630/aabc77. hal-01798144

HAL Id: hal-01798144

https://hal.sorbonne-universite.fr/hal-01798144

Submitted on 23 May 2018

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Reconfigurable optical implementation of quantum complex networks

To cite this article: J Nokkala et al 2018 *New J. Phys.* **20** 053024

View the [article online](#) for updates and enhancements.

## Related content

- *An Introduction to the Formalism of Quantum Information with Continuous Variables: Quantum information with continuous variables*
C Navarrete-Benlloch

- *Spacetime replication of continuous variable quantum information*
Patrick Hayden, Sepehr Nezami, Grant Salton et al.

- *Coupled harmonic systems as quantum buses in thermal environments*
F Nicacio and F L Semião
New Journal of Physics

The open access journal at the forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

28 December 2017

REVISED

25 March 2018

ACCEPTED FOR PUBLICATION

9 April 2018

PUBLISHED

9 May 2018

Original content from this work may be used under the terms of the [Creative Commons Attribution 3.0 licence](https://creativecommons.org/licenses/by/3.0/).

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Reconfigurable optical implementation of quantum complex networks

J Nokkala<fnref n="1" />, F Arzani<fnref n="2" />, F Galve<fnref n="3" />, R Zambrini<fnref n="3" />, S Maniscalco<fnref n="1" />, J Piilo<fnref n="1" />, N Treps<fnref n="2" /> and V Parigi<fnref n="2" />

E-mail: valentina.parigi@lkb.upmc.fr

**Keywords:** optical frequency comb, quantum networks, arbitrary complex topology

Abstract

Network theory has played a dominant role in understanding the structure of complex systems and their dynamics. Recently, quantum complex networks, i.e. collections of quantum systems arranged in a non-regular topology, have been theoretically explored leading to significant progress in a multitude of diverse contexts including, e.g., quantum transport, open quantum systems, quantum communication, extreme violation of local realism, and quantum gravity theories. Despite important progress in several quantum platforms, the implementation of complex networks with arbitrary topology in quantum experiments is still a demanding task, especially if we require both a significant size of the network and the capability of generating arbitrary topology—from regular to any kind of non-trivial structure—in a single setup. Here we propose an all optical and reconfigurable implementation of quantum complex networks. The experimental proposal is based on optical frequency combs, parametric processes, pulse shaping and multimode measurements allowing the arbitrary control of the number of the nodes (optical modes) and topology of the links (interactions between the modes) within the network. Moreover, we also show how to simulate quantum dynamics within the network combined with the ability to address its individual nodes. To demonstrate the versatility of these features, we discuss the implementation of two recently proposed probing techniques for quantum complex networks and structured environments.

1. Introduction

During the last twenty years, network theory has experienced remarkable progress and revolutionized the research in diverse disciplines ranging, e.g., from technology to social sciences and biology [1–4]. The discovery of new types of networks, such as having scale-free [5] or small-world properties [6], and development of new tools, e.g. community detection [7], has led to the invaluable role that network theory currently has in empirical studies of many real-world complex systems.

By character, the combination of network theory and its applications are multidisciplinary and during the recent years a new area applying network theory and complex networks to quantum physical systems has emerged [8, 9]. Furthermore, specific features of quantum networks with no classical equivalent have been reported [10–13].

Indeed, due to their ubiquitous nature, complex networks have found applications also in diverse topics in theoretical physics. The properties and topology of the underlying complex network influence localization properties of coherent excitations [14], phase transition of light [15], Bose–Einstein condensation of non-interacting bosons [16], and quantum walks with the associated quantum transport [17]. With technological advances for quantum communication and development of quantum internet [18], the future information networks have to handle not only their computational complexity but also their complex geometrical structure, and there are already results on how the network topology influences entanglement percolation [10, 19].

<footnote n="1">¹ Turku Centre for Quantum Physics, Department of Physics and Astronomy, University of Turku, FI-20014, Turun Yliopisto, Finland</footnote>

<footnote n="2">² Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-PSL Research University, Collège de France; 4 place Jussieu, F-75252 Paris, France</footnote>

<footnote n="3">³ IFISC (UIB-CSIC), Instituto de Física Interdisciplinar y Sistemas Complejos, UIB Campus, E-07122 Palma de Mallorca, Spain</footnote>
Moreover, complex graphs of higher dimensions are strongly connected with special quantum algorithms [20] and networks have also been used for studies in quantum gravity [21]. In general, earlier theoretical discussions on quantum complex networks include topics from the control, probing and engineering of the network [22–26], to the generation and protection of quantum resources [25, 27–29], as well as the study of collective phenomena [30, 31], like quantum synchronization [32, 33].

However, despite the theoretical progress, a large gap remains hindering the practical developments: is it actually possible to create and control quantum complex networks with arbitrary topology in a single setup?

We give a positive answer for this question by combining earlier theoretical work on complex networks of interacting harmonic oscillators [32, 34] with recent experimental advances in creating multipartite entangled states in a multimode optical system [35]. Interacting optical modes represent an interesting option because they are highly resilient to noise, highly controllable with classical instruments and efficiently detectable. Our proposal is based on a compact optical setup including optical frequency combs and parametric processes which, along with mode selective and multimode homodyne measurements, allows for the implementation of networks with reconfigurable coupling and topology, with sufficient size and diversity to be relevant in the context of complex networks. Here, the modes represent the nodes in the network and the interactions between them the links. We want also to stress here that although the efficient implementation of quantum networks has been achieved in other platforms, e.g. superconducting qubit [36], atoms in optical lattices [37] and single photons [38], the topologies which have been explored are mainly regular ones (lattice, circulant graphs, triangular graphs) and they mainly concern logical encoding rather than physical interactions between the quantum systems. On the contrary, our platform allows for: the deterministic implementation of the network, as the mapping is based on continuous variables, and the implementation of arbitrary complex topology within the limits of the experimentally achievable size. Moreover, the reconfigurability, i.e. varying the topology for the network, does not require to change the optical setup.

The quantum nature of the networks is provided by the ability to initialize the oscillators in quantum states and/or generate entanglement connections between nodes. Part of this strategy has already demonstrated the fabrication of multipartite entangled states [39] and cluster states [35, 40]. Here we address a more general scenario, with additional tools and a specific mapping, for the implementation of quantum complex networks.

Several features make the scheme truly appealing for the general study of quantum complex networks. We can control the number of nodes in the network and in principle any network—whether complex or not—can be created from a given set of nodes. Moreover, the system allows to simulate quantum dynamics within the network by mapping the dynamical results of [34] for optimized experimental parameters of the optical multimode set-up. It is also important to note that each node of the network can be individually addressed which opens significant possibilities to probe the global properties of the network by detecting the local properties, as proposed in [34, 41]. The proposal also opens the possibility to design quantum simulators for continuous variable open quantum systems. This is important because open system dynamics and memory effects have been very actively studied in the last ten years both theoretically [42–44] and experimentally (see, e.g., [45]). On the other hand, memory effects have been predicted considering microscopic bath realization of coupled oscillators [46] and can be realized with the present set-up. However, most of the experimental work in this context has dealt with conceptually rather simple single qubit (photon) open systems and no experiments on controlled continuous variable open system dynamics with memory effects exist yet, to the best of our knowledge. Our results also benefit research on complex quantum communication and information networks [13, 47], and in energy transport and harvesting within biological systems [11, 48]. In the following sections, we first describe the dynamics within the complex network of coupled harmonic oscillators and then show in detail how this can be mapped to the optical platform. We then continue by demonstrating the implementation of two probing techniques for complex networks, and discuss the requirements of scalability and reconfigurability of the system.

## 2. Network dynamics

We consider bosonic quantum complex networks composed by an ensemble of quantum harmonic oscillators, with frequencies $\omega_i$, linked by spring-like couplings according to a specific topology. This is defined by the adjacency matrix V, which contains the coupling terms $v_{ij}$ for any couple of nodes. Notice that since the interactions between the oscillators are symmetric, the adjacency matrix must be symmetric as well, meaning that the networks will be undirected. This system can provide a quantum model of a heat bath sustained by an electromagnetic field, or it can be interpreted as a generic bosonic structure which can be exploited, e.g. for transporting energy or information. The network Hamiltonian is $H_E = \mathbf{p}'^\mathrm{T}\mathbf{p}'/2 + \mathbf{q}'^\mathrm{T}\mathbf{A}\mathbf{q}'$ where we have vectors of momentum and position operators $\mathbf{p}'^\mathrm{T} = \{p_1', ..., p_N'\}$, $\mathbf{q}'^\mathrm{T} = \{q_1', ..., q_N'\}$ for N harmonic oscillators. The matrix A has diagonal elements which can be written by introducing the effective frequencies $\tilde{\omega}_i$
as $A_{ii} = \tilde{\omega}_i^2/2 = \omega_i^2/2 + \sum_j v_{ij}/2$, and off-diagonal elements $A_{i \neq j} = -v_{ij}/2$. Note that A can be expressed in terms of the adjacency matrix V of the network and a diagonal matrix containing the effective frequencies $\Delta_{\tilde{\omega}} = \text{diag}\{\tilde{\omega}_1^2, \dots, \tilde{\omega}_N^2\}$, as $A = \Delta_{\tilde{\omega}}/2 - V/2$.

As the bosonic oscillators are intended to describe electromagnetic fields we use a more common formalism in quantum optics by substituting $\mathbf{q}'$ and $\mathbf{p}'$ with the renormalized quadratures $\mathbf{q}$ and $\mathbf{p}$ by defining $\mathbf{q}^T = \mathbf{q}'^T \sqrt{\Delta_{\tilde{\omega}}}$ and $\mathbf{p}^T = \mathbf{p}'^T \sqrt{\Delta_{\tilde{\omega}}^{-1}}$ with $\Delta_{\tilde{\omega}} = \text{diag}\{\omega_1, \dots, \omega_N\}$. So the network Hamiltonian is rewritten as

$$H_E = \frac{\mathbf{p}^T \Delta_{\tilde{\omega}} \mathbf{p}}{2} + \mathbf{q}^T \sqrt{\Delta_{\tilde{\omega}}^{-1}} \mathbf{A} \sqrt{\Delta_{\tilde{\omega}}^{-1}} \mathbf{q}. \quad (1)$$

By diagonalizing the symmetric matrix A with an orthogonal matrix K, such as $K(\Delta_{\tilde{\Omega}}/\sqrt{2})^2 K^T = A$, where $\Delta_{\tilde{\Omega}}$ is a diagonal matrix = diag{$\Omega_1, \dots, \Omega_N$}, we can define decoupled oscillators, also named normal modes, with quadratures

$$ \begin{aligned} \mathbf{Q} &= \sqrt{\Delta_{\tilde{\Omega}}} \mathbf{K}^T \sqrt{\Delta_{\tilde{\omega}}^{-1}} \mathbf{q}, \\ \mathbf{P} &= \sqrt{\Delta_{\tilde{\Omega}}^{-1}} \mathbf{K}^T \sqrt{\Delta_{\tilde{\omega}}} \mathbf{p}. \end{aligned} \quad (2) $$

Plugging equation (2) in (1) the network Hamiltonian becomes

$$H_E = (\mathbf{P}^T \Delta_{\tilde{\Omega}} \mathbf{P} + \mathbf{Q}^T \Delta_{\tilde{\Omega}} \mathbf{Q})/2 = \sum_{i}^{N} \Omega_i (P_i^2 + Q_i^2)/2. \quad (3)$$

From the free temporal evolution of the decoupled quadratures {$\mathbf{Q}, \mathbf{P}$} and by inverting the relations (2) we can recover the temporal evolution of the network oscillators:

$$\begin{pmatrix} \mathbf{q}(t) \\ \mathbf{p}(t) \end{pmatrix} = \begin{pmatrix} T_1 D_{\cos}^{\Omega} & T_1 D_{\sin}^{\Omega} \\ -T_2 D_{\sin}^{\Omega} & T_2 D_{\cos}^{\Omega} \end{pmatrix} \begin{pmatrix} \mathbf{Q}(0) \\ \mathbf{P}(0) \end{pmatrix} \quad (4)$$

and then

$$\begin{pmatrix} \mathbf{q}(t) \\ \mathbf{p}(t) \end{pmatrix} = \begin{pmatrix} T_1 D_{\cos}^{\Omega} T_1^{-1} & T_1 D_{\sin}^{\Omega} T_2^{-1} \\ -T_2 D_{\sin}^{\Omega} T_1^{-1} & T_2 D_{\cos}^{\Omega} T_2^{-1} \end{pmatrix} \begin{pmatrix} \mathbf{q}(0) \\ \mathbf{p}(0) \end{pmatrix}, \quad (5)$$

where $T_1 = \sqrt{\Delta_{\tilde{\omega}}} \mathbf{K} \sqrt{\Delta_{\tilde{\Omega}}^{-1}}$ and $T_2 = \sqrt{\Delta_{\tilde{\omega}}^{-1}} \mathbf{K} \sqrt{\Delta_{\tilde{\Omega}}}$ and we have introduced the diagonal matrices $D_{\sinii}^{\Omega} = \sin(\Omega_i t)$ and $D_{\cosii}^{\Omega} = \cos(\Omega_i t)$. The equations can be given in the compact form

$$ \begin{aligned} \mathbf{x}(t) &= S'_{\mathbf{V}} \mathbf{X}(0), \\ \mathbf{x}(t) &= S_{\mathbf{V}} \mathbf{x}(0), \end{aligned} \quad (6) $$

where $\mathbf{X}^T = \{\mathbf{Q}^T, \mathbf{P}^T\}$ and $\mathbf{x}^T = \{\mathbf{q}^T, \mathbf{p}^T\}$. The optical implementation of the evolution described by (4) or (5) requires a multimode system, consisting of a number of modes equal to the number of nodes in the network. The optical modes need to be prepared in the initial state described by the quadratures {$\mathbf{Q}(0)^T, \mathbf{P}(0)^T$} in the normal modes picture or, equivalently, by {$\mathbf{q}(0)^T, \mathbf{p}(0)^T$} in the bare oscillators picture. Finally the implementation involves reconfigurable optical processes, able to drive the quadratures evolution $S'_{\mathbf{V}}$ or $S_{\mathbf{V}}$ at different times, which are intimately connected with the network structure V, the bare frequencies of the network oscillators and the frequencies of the normal modes. The evolutions described by equations (4) and (5) being equivalent, we can choose one of the two pictures according to which is more convenient for evaluating the desired expectation values or which one between $S'_{\mathbf{V}}$ and $S_{\mathbf{V}}$ could be more efficiently implemented in the experimental setup. The matrices $S'_{\mathbf{V}}$ and $S_{\mathbf{V}}$ in equation (6) are symplectic and they can then be factored via the Bloch–Messiah (B–M) decomposition [49], for example:

$$S_{\mathbf{V}} = R_{1}^{T} \Delta^{\mathrm{sq}} R_{2}. \quad (7)$$

$\Delta^{\mathrm{sq}}$ being a diagonal matrix and $R_1, R_2$ symplectic and orthogonal matrices. In figure 1 three different networks composed by $N = 51$ (figure 1 A) or $N = 50$ (figures 1(B), (C)) nodes are shown: a periodic chain, a linear chain with shortcuts and a random network, the bare frequencies of the network oscillators are set to $\omega_i = \omega_0 = 0.25$ for the three different topologies. The network matrices $S_{\mathbf{V}}$ calculated at different times from the evolution in equation (5) are shown in figure 2, while the B–M decomposition at time $t = 50$ for the three network is shown in figure 3.

Note that the complexity of B–M decomposition is of the order of $O(N^3)$ [50]. It takes around 0.03 s to implement the B–M decomposition for networks of 50 nodes with a standard laptop.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Graphical structure of the investigated networks. (A) is a periodical chain with 51 nodes, the coupling strength $v_{ij} = 0.1$, is reduced to $v_{ij} = 0.06$ every third link. (B) is a linear chain, with $v_{ij} = 0.1$, where weak shortcuts with strength $v_{ij}/50$ are added. In the centre the shortcuts structure is shown. (C) is an Erdős–Rényi random network of 50 nodes with $v_{ij} = 0.01$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Matrix plot of the calculated evolution $S_v$, i.e., the matrix of equation (5), for the three networks of figure 1, at different times, all the oscillators frequencies are set to $\omega_j = \omega_0 = 0.25$. The matrices contain 100 (or 102) rows as they concern the evolution at the given time of the 50 (51) q quadratures plus the 50 (51) p quadratures of the network oscillators. A temperature colour map is used: positive values are in red colours, negative in blue while white colour indicates the zero value, as shown by the legend in the plot at $t = 999$.</figcaption></figure>

### 3. Mapping the network to a multimode experimental platform

The matrices $\Delta^{\text{sq}}$, $R_1$ and $R_2$ in equation (7) have a well-known physical interpretation in optics: the first is a squeezing operation while each of the two others corresponds to the action of a linear multiport interferometer, i.e. a basis change. The implementation of $S_v$ can then be decomposed in a combination of squeezing and linear optical operations in a multimode scenario. In the following we are going to present a particular compact and completely reconfigurable experimental implementation of such strategy.

The experimental implementation proposed here is based on parametric processes pumped by optical frequency combs. A very large network of entangled photons is produced [35, 51] under a quadratic Hamiltonian, that leads to a symplectic evolution. The mapping we present here aims at matching this evolution
<figure><img src="image_6.png" /></figure>

to $S_V$ using the experimentally controllable degree of freedoms: the pump spectral amplitude and the measurement basis.

In practice, the laser sources considered here are mode-locked lasers emitting pulses from 30 to 150 fs with repetition rates from 70 to 150 MHz at a wavelength centred around $\lambda = 800$ nm. The spectrum of these lasers is constituted of hundreds thousands of frequency components. A second harmonic generation process is used for generating a frequency-doubled comb which serves as pump of the parametric process. The pump frequencies are written as $\omega_p = \omega_{p0} + p \cdot \omega_{RR}$, where $\omega_{RR}$ is the repetition rate of the original laser source, $\omega_{p0}$ is the pump central frequency and p an integer. The Hamiltonian describing the parametric process, which is parametric down conversion in a $\chi^{(2)}$ crystal, is $H = ig \sum_{m,n} \mathcal{L}_{m,n} a_m^\dagger a_n^\dagger + \text{h.c.}$; it couples pairs of frequencies $\omega_m$ and $\omega_n$ via the term $\mathcal{L}_{m,n}$ which is the product of crystal's phase matching function $f_{m,n}$ and of the complex pump amplitude $\alpha_p$ at frequency $\omega_p = \omega_m + \omega_n$. $\mathcal{L}_{m,n}$ describes the probability amplitude that a photon at frequency $\omega_p$ is converted into two daughter photons at frequencies $\omega_n$ and $\omega_m$. Given the number of frequency components in the pump, a number of down-converted frequencies, which is of the same order of magnitude that the one in the original spectrum, is combined in a non-trivial multipartite entangled structure [39]. This Hamiltonian governs the system evolution that we express for the optical quadrature operators $\{q_1, q_2, \dots, q_N, p_1, p_2, \dots, p_N\} = \{q, p\} = X_b^T$ in a given mode basis b, with the convention $q_j = (1/\sqrt{2})(a_j + a_j^\dagger)$ and $p_j = (i/\sqrt{2})(a_j^\dagger - a_j)$. A natural choice for the basis would be the frequency components of the comb spectrum. The parametric Hamiltonian leads to a symplectic transformation $X_f = S_{ex}X_i$ which can be decomposed via B-M in

$$S_{ex} = R_{1ex}^T \Delta_{ex}^{sq} R_{2ex}. \qquad (8)$$

If the input modes $X_i$ are in the vacuum state the $R_{2ex}$ term can be neglected

$$X_f = R_{1ex} \Delta_{ex}^{sq} X_i. \qquad (9)$$

The squeezing parameters in $\Delta_{ex}^{sq}$ can be derived from the eigenvalues of $\mathcal{L}$ and $R_{1ex}$ can be reinterpreted as change of basis from the basis of the squeezed modes, named 'supermodes', towards the measurement basis. In the case of a pump and phase matching function with Gaussian spectral profile, the supermodes basis obtained from $\mathcal{L}$ coincide with Hermite-Gauss modes. If the symplectic transformation $S_{ex}$ can be experimentally controlled in order to coincide with $S_V$ we are obtaining an optical implementation of the complex network described by the equation (5). This means that by addressing the optical modes we can measure the properties of the oscillators in the network.
<figure><img src="image_7.png" /><figcaption>**Table 1.** The mapping of the quantum networks to the experimental platform. Refer to main text for more details.</figcaption></figure>

<table><thead><tr><th></th><th>Quantum network</th><th>Experimental implementation</th></tr></thead><tbody><tr><td>Node</td><td>Quantum harmonic oscillator</td><td>Optical mode</td></tr><tr><td>Link</td><td>Spring-like coupling term</td><td>Nonlinear mode coupling</td></tr><tr><td>t</td><td>Evolution time</td><td>Parameter controlling the symplectic matrix S<sub>ex</sub></td></tr><tr><td>Frequency</td><td>Oscillator frequencies</td><td>Squeezing of optical modes</td></tr><tr><td>Decoupled mode</td><td>Normal mode of the network</td><td>Supermode, a non-correlated squeezed optical mode</td></tr><tr><td>Addressing a node</td><td>Local measurement</td><td>Pulse shaping and mode-selective measurement of the supermodes</td></tr></tbody></table>

We would like to stress here the special features of our mapping: to emulate the network dynamics we only need to know its symplectic evolution $S_V$ at time $t$, which can be experimentally emulated by $S_{ex}$. We do not need to know the physical phenomenon that generates the network's dynamics or the nature of the oscillators in the network (they could be mechanical or optical oscillators) or even the specific picture that can be chosen to describe the Hamiltonian evolution. In the case described by equation (5), for example, the symplectic evolution is derived from the full Hamiltonian also containing the non-interacting terms, while the experimental process is described in the interaction picture. As a consequence, the frequencies of the network oscillators will not match the frequencies of the correspondent modes, but instead will be mapped to squeezing parameters in the optical setup. By focusing on the emulation of the symplectic evolution, we can then use the same experimental resources to map different network parameters according to the features of the phenomenon we are studying.

The mapping of the quantum networks to the experimental platform is summarized in table 1.

In the following we explain how the symplectic transformation can be controlled by manipulating the terms of its B–M decomposition. The squeezing operations in $\Delta_{ex}^{sq}$ can be controlled by tailoring the spectral shape of the pump in the parametric process. This can be done experimentally via a pulse shaper based on a spatial light modulator (SLM) in a 4-f configuration (see figure 4) which provides control over amplitude and phase of the spectral components by exploiting a 2D geometry. The phase mask of the SLM can be set according to the results provided by an optimization procedure, based on an evolutionary algorithm which directly relates experimental parameters with the diagonal terms of the matrix $\Delta_{ex}^{sq}$ [52]. Note that the optimization of $\Delta_{ex}^{sq}$ also influences the spectral shape of the supermodes, and thus the squeezed modes basis.

The linear transformation $R_{1ex}$ can be tuned adequately choosing the measurement basis, which is the equivalent of choosing the basis on which evolution (9) is calculated. Experimentally, this is implemented using homodyne detection where the measured mode is governed by the spectral shape of the local oscillator (LO) [40]. The squeezing matrix $\Delta_{ex}^{sq}$ acts on the supermode basis, then the quadratures of the $m$th network oscillator can be addressed by applying the basis change defined by the $m$th row of the matrix $R_{i}^{T}$. This means that we can measure the quadratures $q_{m}(t)$ and $p_{m}(t)$ of the $m$th network oscillator by performing homodyne detection with a LO shape defined by $u_{LO-m}(\omega) = \sum_{j=1}^{N} R_{1mj}^{T} SP_{j}(\omega) + i \sum_{j=N+1}^{2N} R_{1mj}^{T} SP_{j-N}(\omega)$, where $SP_{i}(\omega)$ are the supermodes. Figures 4(A) and (B) show the scheme for the experimental control of $\Delta_{ex}^{sq}$ via the control of the
<figure><img src="image_8.png" /><figcaption>Figure 5. Simulation of the periodic and the random network evolution by the experimental setup. 1, 2, 3 and 4 stand for $t = 50, 150, 500$ and $999$, respectively. Upper part of each box: diagonal terms of the matrix $\Delta^{\text{sq}}$ in dB, the red are the theoretical values coming from the B–M decomposition, the blue are the optimized values which can be obtained from an optimized pump shape with a BiBO crystal of length 2.5 mm, in the case of the periodic network, and of 1.5 mm, in the case of random network. Lower part of each box: spectral shape of mode 26 representing the 26th oscillator obtained by applying the linear transformation $R_1^T$ on the calculated supermodes. The amplitude is represented by a plot with rainbow colours while the grey line is the spectral phase. Units are arbitrary. This corresponds to the spectral shape which has to be set in the LO in order to address the mode.</figcaption></figure>

pump spectral shape, after the optimization, and the control of $R_{1ex}$ by choosing the measurement basis via pulse shaping of the LO.

A second point has to be stressed here: even if the B–M reduction of general symplectic transformations and its physical interpretation is well known, the full control allowing for complete reconfigurability of the squeezing operations and basis changes in a single optical setup, enabling from trivial to complex topologies, has never been addressed before.

The control of $R_{ex}$ and $\Delta_{ex}^{\text{sq}}$ is sufficient for the simulation of the dynamics if the particular feature we are interested in is independent on the initial quantum state of the network oscillators. In that case we can consider the oscillators being in the vacuum state at time $t = 0$, by choosing $R_{ex} = R_1^T \Delta_{ex}^{\text{sq}} = \Delta^{\text{sq}}$, and discarding $R_2$ the network simulation can be performed. Otherwise we need to initialize the optical modes in the desired quantum states. The oscillators can be prepared in an initial pure Gaussian state, by considering a first symplectic transformation $S_0$ which prepares the modes in the desired state, followed by network evolution $S_V$. The effective total evolution to be simulated is

$$x(t) = S_V S_0 x(0) = S_{ef} x(0), \qquad (10)$$

where the initial states are vacua. Then via B–M decomposition we can set

$$S_{ef} = R_{1ef} \Delta_{ef}^{\text{sq}} R_{2ef}, \qquad (11)$$

experimentally implement $\Delta_{ex}^{\text{sq}} = \Delta_{ef}^{\text{sq}}$ and $R_{ex} = R_{1ef}$ and finally discard $R_{2ef}$ as it is acting on vacuum states. If the complex networks we are implementing are intended to simulate environments of open quantum systems, a typical requirement we may have to fulfil is to set a finite temperature for the network. This can be achieved by initializing the oscillators in an optical thermal state, which can be built, for instance, as the reduced state of an entangled state. Indeed, a maximally entangled state of 2 qubits has marginal distributions which are maximally mixed (equivalent to $T \to \infty$), and the same is true in continuous variable Gaussian states. Consider e.g., a two-mode squeezed state of two bosonic modes with squeezing parameter $r$. The reduced covariance matrix of only one of the modes reads $\sigma = \text{diag}\{\hbar\omega \cosh(2r)/2, \hbar\omega \cosh(2r)/2\omega\}$, which describes a thermal state with $\cosh(2r) \to \coth(\hbar\omega/2T)$. For $T/\omega > 1$ the required squeezing grows as $r \approx \text{acosh}(2T)/2$ and this scales like $\log T$, e.g. $r \sim 2$ for $T/\omega \sim 20$, which is reasonable experimentally. We may want to initialize one oscillator of the network in a thermal state. This is then done by having two probes initially in a two-mode squeezed state, which is obtained by a $\pi/4$ rotation in the $2 \times 2$ subspace, individual squeezings, and a rotation back. If instead we need to initialize the entire network into a thermal state we should duplicate it and impose a two-mode squeezed state between each pair of node copies. Network dynamics then follows by applying the time evolution only on half of the network units, and the identity on the other half. This scheme increases the overhead of the protocol but is feasible in principle.

The figure 5 shows the results of the simulation for the periodic network and the random network based on experimental parameters. The optimization procedure for the pump shape is run in order to minimize the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Simulation of the periodic (A) and the random (C) network in the eigenmode picture at $t = 50$ with a crystal length of 1.5 mm. The squeezing values here are the diagonal terms of $\Delta^{\text{sq}}$ from the B–M decomposition of $S'_V$ in equation (6). In the centre: the pulse shape for the 26th node. In the lower part the calculated mean photon number of each oscillator (brown curve) and the simulated one (black curve) from the experimentally achievable squeezing values.</figcaption></figure>

distance between the theoretical values (the diagonal of $\Delta^{\text{sq}}$) and the experimentally realizable ones ($\Delta_{\text{ex}}^{\text{sq}}$) by varying the spectral shape of the pump. The simulated squeezing values are close to the theoretical ones: the relative distance of the simulated curves to the theoretical ones are in between 6.8% and 10.4% for the periodic network and in between 3.5% and 5.3% for the random network. The spectral shape of the mode corresponding to one oscillator of the network is given: 0 nm corresponds to the central wavelength of the laser source (800 nm). The relative errors on the experimentally simulated $S_V$ are exactly the same of $\Delta^{\text{sq}}$ if the exact required basis change $R_1$ is performed. This necessitates that any of the 50 modes simulating the oscillators can be addressed by choosing the right shape of the LO in homodyne detection. As shown in figure 5, this requires to handle laser light of large bandwidth, which can be generated by coherent broadening of the main laser source. The results of figure 5 are obtained by setting two different crystal lengths: in the case of the shorter crystal (1.5 mm) the experimentally achievable squeezing matrix ($\Delta_{\text{ex}}^{\text{sq}}$) is closer to the theoretical one ($\Delta^{\text{sq}}$) than in the case of 2.5 mm crystal length. However the case of shorter crystal generally requires larger bandwidth for the LO as the resulting supermodes have a larger spectrum. In figure 6 we show one observable, the mean photon number, for all the oscillators in the network. We evaluate both its theoretical value with the exact evolution and the simulated one with the squeezing values derived from the optimization. For this calculation we choose the vacuum state as the initial state. The relative errors on the calculated excitation numbers is 5.3% for the periodic and 5.0% for the random network. The simulation is performed with a crystal length of 1.5 mm, also in this case by admitting larger relative errors (of the order of 10% in the excitation numbers) it is possible to reduce the bandwidth of the LO (to half the value of the shapes shown in figure 6) with a crystal length of 2.5 mm.

In order to show that the setup has the potential of replicating the behaviour of complex networks with any topology by simulating $S_V$ or $S'_V$ at different times, we repeat the procedure for other networks with more complex topology, which are shown in figure 7. We study the Barabási–Albert model [53], characterized by a power law of the degree distribution, and the Watts–Strogatz model [54], which describes networks with small-world properties as clustering. Additionally, we consider two graphical structures derived from real natural networks, studied by Scannell et al [55] and Lusseau [56], that we call Scannel and Lusseau networks. The simulations are shown in figure 8 and they show the same performances of the simulations in figure 5 with a similar trade-off between relative errors limitation and reduction of the LO bandwidth.

## 4. Interaction with additional oscillators

In a more general scenario we would like to describe the network interacting with additional oscillators: these would, for example, depict the principal quantum systems which interact with the bath described by the network or the systems in which the information we want to convey through the network is initially encoded. We take *M* oscillators with frequencies {$\omega_{S1}, \omega_{S2}, \dots, \omega_{SM}$}, each of them interacting with one node according to the Hamiltonian $H_{\text{lr}} = -kq_{S\text{r}}q_i/\sqrt{\omega_{S\text{r}}\omega_i}$. The interaction couples the supplementary oscillators with any of the
<figure><img src="image_10.png" /></figure>

normal mode, as we can see from the interaction Hamiltonian in the decoupled quadrature picture: $H_{Ir} = -kq_{S,r} \sum_j K_{ij}Q_j / \sqrt{\omega_{S,r}\Omega_j}$. The new total Hamiltonian again can be written in diagonal form by orthogonal diagonalization of $B = O(\Delta_f / \sqrt{2})^2 O^T$, the matrix analogue to A containing the extra terms given by the interaction. $\Delta_f$ is a diagonal matrix diag{$f_1, \dots, f_{N+M}$}.

$$H_{S+E+I} = \sum_{j}^{N+M} f_j (P_j^2 + Q_j^2)/2. \qquad (12)$$

It is convenient to describe the dynamics of the network plus the new oscillators in terms of the evolution of the quadratures for the network eigenmodes and the supplementary systems $\dot{X}^T = \{Q, q_S, P, p_S\} = \{Q, q_{S1}, \dots, q_{SM}, P, p_{S1}, \dots, p_{SM}\}$

$$\begin{pmatrix} Q(t) \\ q_S(t) \\ P(t) \\ p_S(t) \end{pmatrix} = \begin{pmatrix} O_1 D_{\cos} O_1^{-1} & O_1 D_{\sin} O_2^{-1} \\ -O_2 D_{\sin} O_1^{-1} & O_2 D_{\cos} O_2^{-1} \end{pmatrix} \begin{pmatrix} Q(0) \\ q_S(0) \\ P(0) \\ p_S(0) \end{pmatrix}, \qquad (13)$$

where we introduced the diagonal matrices $D_{\cos i} = \cos(f_i t)$ and $D_{\sin i} = \sin(f_i t)$, and the matrices $O_1 = \sqrt{\Delta_{\Omega, \omega_S}} O \sqrt{\Delta_f^{-1}}$ and $O_2 = \sqrt{\Delta_{\Omega, \omega_S}^{-1}} O \sqrt{\Delta_f}$, being $\Delta_{\Omega, \omega_S} = \text{diag}\{\Omega_1, \dots, \Omega_N, \omega_{S1} \dots \omega_{SM}\}$. Equation (13) can be again written in a more compact form as $X(t) = SX(0)$ with S a function of time. The simulation of the dynamics of the network interacting with the system can then be performed by experimentally addressing a collection of oscillators X and implementing the transformation S at any time.
<figure><img src="image_11.png" /></figure>

## 5. Probing the spectral density of a structured environment

The additional oscillators introduced above can also be considered as probes for the network in the framework of the theory of open quantum systems [57, 58]. Here, the probes would be considered the open quantum systems and the network their environment, and the objective would be to deduce some properties of the environment based on knowing the reduced dynamics of the probes. In the following we consider the case of a single probe and we describe the simulation, based on experimentally realizable operations, of the probing technique for the network spectral density as proposed in [34].

The spectral density $J(\omega)$, a key quantity in open quantum systems theory, is determined by the environment and interaction Hamiltonians $H_E$ and $H_I$, and gives the density of environmental modes as well as coupling strengths to them as a function of frequency [46, 58]. Nothing else about $H_E$ and $H_I$ needs to be known to determine the open system dynamics. Conversely, knowing the open system dynamics reveals information about $J(\omega)$. Indeed, provided that the coupling between the probe and the environment is sufficiently weak, measuring the photon number operator $\langle n \rangle$ of the probe having frequency $\omega_S$ allows to approximate $J(\omega_S)$ from

$$J(\omega_S) = \frac{\omega_S}{t} \ln \left( \frac{N(\omega_S) - \langle n(0) \rangle}{N(\omega_S) - \langle n(t) \rangle} \right), \quad (14)$$

where $N(\omega_S) = (e^{(\omega_S/T)} - 1)^{-1}$ is the thermal average boson number and $t$ is a large enough interaction time. In the case at hand, one may couple the probe weakly to any subset of network oscillators and use equation (14) to probe $J(\omega_S)$, however different subsets correspond to different spectral densities as the function depends also on $H_I$.

The mean photon number here can be accessed from the homodyne measurement of the probe as $\langle n(t) \rangle = 1/2 \langle q_S(t)^2 \rangle + 1/2 \langle p_S(t)^2 \rangle - 1/2$. The quadratures of the probe $\{q_S(t), p_S(t)\}$ appear in the left hand side of the equation (13), which can be used to derive their expectation value at any time. We have to stress here
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Purple: diagonal elements of the target $\Delta_{ef}^{sq}$; green: experimental achievable squeezing spectrum, i.e. diagonal elements of Kex in the case of a Gaussian pump spectrum for the parametric process with a crystal length of 1.5 mm and with no pump optimization. Brown and red: experimental achievable squeezing spectrum, when optimization on the pump shape is performed with a crystal length of 0.5 and 1.5 mm. In the three cases (green, red and brown) the first diagonal element is exactly the same as the target one (purple). The pictures in the circles show the spectral shape of the first and the third supermode in the three possible experimental configuration.</figcaption></figure>

that when the initial mode basis $X(0)$ is defined, the B–M decomposition of $S$ (the matrix appearing on the left side of equation (13)) completely identifies, via the basis changes, the final mode corresponding to $\{q_S(t), p_S(t)\}$, which can be simply addressed by pulse-shaped homodyne detection. In this case we cannot initialize all the network nodes and the probe in the vacuum state, because no energy exchange between the network and the system would be mapped in $\langle n \rangle$. We then will use a first operation able to set a non-zero value for the mean photon number of the probe $\langle n(0) \rangle$, as for example a squeezing operation $S_0 = \Delta'$, and then we will implement the B–M decomposition of $S_{\text{eff}} = S\Delta'$.

The protocol is the following one: given a network structure, a probe with frequency $\omega_S$ is chosen and the evolution of the total system is computed at an appropriate time $t$, the experimental parameters are set according to the calculated B–M giving $R_{1\text{ef}}$, $\Delta_{ef}^{\text{sq}}$. The probe mode is measured by pulse shaping the LO and one value of $J(\omega_S)$ is obtained. The procedure is repeated for other $\omega_S$ in order to retrieve the full shape of the spectral density.

We now look at how the network and probe parameters are mapped in the model and particularly in the B–M decomposition of the symplectic matrix in equation (13). The network is totally identified by the matrix $A$, which contains the geometrical structure, the frequency of the oscillators and the strength of the couplings between them. The probe is identified by its frequency $\omega_S$ and the strength of its coupling $k$ to one of the network node which is supposed to be weak in the probing scheme. For simplicity we set a constant value for the frequencies and the couplings in the network, i.e. $\omega_i = \omega$ and $v_{ij} = g$. If we fix the set of values $\{\omega, g, \omega_S, k\}$ we find almost the same matrix $\Delta^{\text{sq}}$ in the B–M decomposition of $S$ for any network at any time, while the different network geometries lead to completely different basis changes $R_1, R_2$. The frequency terms $\{\omega, \omega_S\}$ in fact enter in the renormalized quadrature definitions and they contribute in the Hamiltonian as squeezing terms, while the parameters $\{g, k\}$ are the ones driving the quadratic Hamiltonians which contain squeezing operations. On the other hand the basis change transformations have almost the same role as in cluster state generation, by defining the graphical structure, and they are hence intimately connected with the network structure. When the probe is in a very weak coupling regime $k < 0.01$, condition which can be always chosen without affecting the results of the probing technique, the influence of $\omega_S$ on the squeezing matrix $\Delta^{\text{sq}}$ is very small. This simplifies the experimental procedure of the probing protocol because scanning $\omega_S$ requires almost no changes in the squeezing structure (i.e. in pump shaping) but only appropriate pulse shaping of the LO mode, according the updated basis change. Figure 9 shows in purple the target squeezing, i.e. the diagonal elements of $\Delta_{ef}^{\text{sq}}$, in the case of a periodic chain, with parameters set as in figure 1; ten different purple lists are superposed corresponding to
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The spectral density $J(\omega)$ for the periodic (A), the shortcuts (B) and the random network (C). Blue line: analytical calculation, violet circles: exact probing method. Green dots: calculation of $J(\omega)$ by replacing $\Delta_{ef}^{sq}$ with the $\Delta_{ex}^{sq}$ which can be obtained in the experimental setup with no pump optimization. The red and the brown dots are calculated with the $\Delta_{ex}^{sq}$ obtained from pump optimization with a crystal of 1.5 and 0.5 mm. Top: the calculated spectral density is plotted after normalization to its maximal value. Down: the results of direct calculation are shown.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** The spectral density $J(\omega)$ for the Barabási–Albert (D), the Watts–Strogatz (E), the Scannell (F) and the Lisseau (G) networks. Shapes and colours have the same meaning of figure 10.</figcaption></figure>

ten different values of $\omega_S$. The configuration in which we probe with an initially squeezed oscillator requires only one significantly squeezed mode in the matrix $\Delta_{ef}^{sq}$, this is the reason why all the 10 purple lists almost coincide. The diagonal elements of three experimentally achievable $\Delta_{ex}^{sq}$ are shown: the green points correspond to a Gaussian profile for the pump of the parametric process (no pump shaping), while the brown and red points are obtained from the optimization procedures run considering a crystal length of 0.5 and 1.5 mm. In the three cases the first mode can be set to have exactly the required value of squeezing of the first target mode. The insets show the shape of the first and the third supermodes corresponding to the experimentally achievable $\Delta_{ex}^{sq}$. In figures 10 and 11 we compare the results for the analytically calculated $J(\omega)$ (blue lines), the one derived from the exact probing (violet circles), and the ones calculated by replacing the exact $\Delta_{ef}^{sq}$ with the experimentally achievable $\Delta_{ex}^{sq}$ (green, red and brown dots). We show the normalized value of $J(\omega)$ (top) and the raw values after applying equation (14). The experimentally achievable configurations allow for the retrieval of the peculiar shape $J(\omega)$ for the seven kinds of network. In the case of no optimization of the pump shape, even if the intrinsic structure of the spectral density is mainly conserved, the numerical values of the calculated $J(\omega)$ are quite low, while larger values are obtained when optimization via pump shaping is included.

We finally report, in figure 12, the first 10 non-zero values of the $N+1$th row of $R_{1f}$ at different frequencies $\omega_S$. The $N+1$th row identifies the shape of the mode encoding the probe in term of the supermodes basis, the resulting shape has to be set as the shape of the LO in homodyne detection in order to measure $\langle n(t) \rangle$. Examples of the LO shapes are shown in figures 13–15. As in the case of simulation of the networks dynamics, longer
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Calculated first 10 elements of the $N+1$th row of $R_{1f}$, defining the shape of the probe mode as a linear combination of the initial supermodes, for (A): periodic (B): shortcuts and (C): random network. They define the shape of the LO as shown in figure 4(B).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Shape of the LO which has to be set for measuring the probe at three different frequencies $\omega_S$ in the case of a Gaussian pump spectrum for three networks. (A), (B) and (C) stand for periodic, shortcuts and random network.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Shape of the LO which has to be set for measuring the probe at three different frequencies $\omega_S$ when pump shaping with a 0.5 mm crystal is considered.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Shape of the LO which has to be set for measuring the probe at three different frequencies $\omega_S$ when pump shaping with a 1.5 mm crystal is considered for Barabási–Albert (D), the Watts–Strogatz (E), the Scannell (F) and the Lisseau (G) networks.</figcaption></figure>

crystals allow for reduced LO bandwidth. In any case all the LO shapes calculated in the figures 13–15 are experimentally doable.

## 6. Probing network from the entropy of the probe

An interesting concept, developed in [41] is whether a simple, resource friendly, measurement protocol can distinguish different classes of complex networks. It is indeed possible to distinguish among Random Regular, Erdös–Renyi and Scale-Free bosonic networks from features displayed by the averaged single node entropy $\langle S_h \rangle$, where the average is taken over all nodes of the same connectivity $h$. Interestingly, also a set-up where only the entropy of a probe is measured allows to assess such different classes of networks [41]. In this set-up, the probe is sequentially attached to more and more network nodes (increasing its connectivity $h$) and measuring the probe's entropy as a function of $h$ suffices to discriminate network classes. The required von Neumann entropy is easy to access through the covariance matrix of a single node (or the probe) $\sigma = \{\{\alpha, \gamma\}, \{\gamma, \beta\}\}$, and reads $S = (\mu + 1/2)\log(\mu + 1/2) - (\mu - 1/2)\log(\mu - 1/2)$ and $\mu = \sqrt{\alpha\beta - \gamma^2}$ [59]. The protocol is valid for a network in the ground state, but is robust for small temperatures, so the implementation devised above can be easily applied.

## 7. Conclusions and outlook

In conclusion we have shown a protocol for implementing reconfigurable experimental quantum networks with a complex topology in a multimode quantum optical setup. The temporal evolution of networks of various topology can be implemented by tailoring multimode squeezing operations and multimode measurements. This can be done via an optimization protocol on the pump shape for the parametric process and mode selective homodyne measurements. Moreover the probing with one additional oscillator can be simulated: the recovered behaviours of the spectral density are able to univocally identify the corresponding network topology. Also probing based on entropy measurements is realizable. The proposed protocol establishes a connection between the field of continuous variables experimental quantum optics and quantum complex network theory and it shows, for the first time, the full reconfigurability of an optical setup in implementing arbitrary B–M decompositions associated with topologies which range from regular to complex.

We have considered networks both from well-established theoretical models and from real biological and community examples. In order to assess the complexity of the generated networks we are considering geometries involving at least of 40–50 nodes, this is for example the requirement for univocally identifying the reconstructed $J(\omega)$ with a non-trivial graphical structure. The experimental setup already demonstrated the capability of addressing networks made of 12 nodes [35]. While the spectral shape of the required LO for the probing of $J(\omega)$ is within the experimental reach, the required shapes for following the network dynamics will require more experimental effort, depending on the precision we want to reach for the simulated quantities. The extension of involved modes numbers up to 50 is within the capability of the present setup by improving the pulse shaped homodyne detection procedure via the coherent broadening of LO spectrum in order to address the squeezed high-order modes. Moreover the production of multimode quantum states involving a large number of frequency- or temporal-modes has already been demonstrated in other experimental setups [60, 61].
Hence the merging of techniques based on spectral and temporal multiplexing can allow the implementation of larger quantum complex networks.

The quantum character of the networks is provided by the initialization of the nodes in quantum states or establishing entanglement correlations between different nodes. This condition already enables experimental study of control, engineering and probing of the network [25, 34, 41], generation and protection of quantum resources [28, 29], quantum synchronization [32], memory effects and dynamics of open quantum system [42–44, 46]. Also, as in the case of classical networks [62], the resilience of quantum task over network topology can be investigated. In particular quantum communication tasks, like entanglement percolation [10, 19], and their extension to multiparty protocols, like secret sharing, [63] can be studied. The subject of this work is the topological complexity but the present setup with a supplementary tool, i.e. mode-dependent coherent single photon subtractor, can be exploited in future works to demonstrate possible computational complexity, i.e. a form of quantum advantage. The mode-dependent coherent single photon subtraction, which has been recently experimentally demonstrated by some of the authors [64], allows for the introduction of non-Gaussian processes in the network [65] and a specific circuit including multimode squeezing and single photon operations recently demonstrated to be hard to sample [66]. The strategy can be also exploited to test proposed models of quantum enhanced transport [11, 67] via the propagation of non-Gaussian excitation over complex graphs.

A second supplementary tool which can be added is multi-pixel homodyne detector which enables simultaneous measurement of several network components. Here after mixing the signal with the LO in the 50:50 BS the spectral components of the light are spatially dispersed and focussed on two linear array of detectors. Subtracting signals from the detectors corresponding to the same colours allows for simultaneous homodyne detection at each band. Pulse shaping on the LO also allows to independently choose the phase of the quadrature to be measured by independently modifying the phase of each spectral band. Finally computer post-processing provides an additional basis change [68]. The simultaneous multimode analysis can be exploited for implementing a multi-probes scenario and moreover to test measurement based continuous variables quantum information protocols over complex graphs.

## Acknowledgments

JN, FG, RZ, SM and JP acknowledge financial support from the Horizon 2020 EU collaborative projects QuProCS (Grant Agreement No. 641277). FG and RZ acknowledge MINECO/AEI/FEDER through projects NoMaQ FIS2014-60343-P, QuStruct FIS2015-66860-P and EPheQuCS FIS2016-78010-P. NT is members of the Institut Universitaire de France. FA, VP, and NT acknowledge financial support from the European Union Grant QCUMbER (No. 665148). VP acknowledges financial support from Emergence Sorbonne Universités project PiCQuNet.

## ORCID iDs

J Nokkala https://orcid.org/0000-0002-5052-9813

V Parigi https://orcid.org/0000-0002-0779-1971

## References

[1] Albert R and Barabási A-L 2002 Statistical mechanics of complex networks *Rev. Mod. Phys.* **74** 47–97

[2] Newman M, Barabási A-L and Watts D J 2006 *The Structure and Dynamics of Networks* (Princeton, NJ: Princeton University Press)

[3] Newman M 2010 *Networks, An Introduction* (Oxford: Oxford University Press)

[4] Trabesinger A 2012 Complexity *Nat. Phys.* **8** 13 and following articles in this insight issue

[5] Goh K-I, Oh E, Jeong H, Kahng B and Kim D 2002 Classification of scale-free networks *Proc. Natl Acad. Sci.* **99** 12583–8

[6] Amaral L A N, Scala A, Barthélémy M and Stanley H E 2000 Classes of small-world networks *Proc. Natl Acad. Sci.* **97** 11149–52

[7] Rosvall M and Bergstrom C T 2008 Maps of random walks on complex networks reveal community structure *Proc. Natl Acad. Sci.* **105** 1118–23

[8] Bianconi G 2015 Interdisciplinary and physics challenges of network theory *Europhys. Lett.* **111** 56001

[9] Biamonte J, Faccin M and De Domenico M 2017 Complex networks: from classical to quantum arXiv:1702.08459

[10] Acin A, Cirac J I and Lewenstein M 2007 Entanglement percolation in quantum networks *Nat. Phys.* **3** 256–9

[11] Faccin M, Migdal P, Johnson T H, Bergholm V and Biamonte J D 2014 Community detection in quantum complex networks *Phys. Rev. X* **4** 041012

[12] Paparo G D, Muller M, Comellas F and Martin-Delgado M A 2013 Quantum google in a complex network *Sci. Rep.* **3** 2773

[13] Perseguers S, Lewenstein M, Acin A and Cirac J I 2010 Quantum random networks *Nat. Phys.* **6** 539–43

[14] Jahnke L, Kantelhardt J W, Berkovits R and Havlin S 2008 Wave localization in complex networks with high clustering *Phys. Rev. Lett.* **101** 175702

[15] Halu A, Garnerone S, Vezzani A and Bianconi G 2013 Phase transition of light on complex quantum networks *Phys. Rev. E* **87** 022104

[16] Burioni R, Cassi D, Rasetti M, Sodano P and Vezzani A 2001 Bose–Einstein condensation on inhomogeneous complex networks *J. Phys. B At. Mol. Opt. Phys.* **34** 4697
[17] Mulken O and Blumen A 2011 Continuous-time quantum walks: model for coherent transport on complex networks *Phys. Rep.* **502** 37–87

[18] Kimble H J 2008 The quantum internet *Nature* **453** 1023–30

[19] Cuquet M and Calsamiglia J 2009 Entanglement percolation in quantum complex networks *Phys. Rev. Lett.* **103** 240503

[20] Rossi M, Bruß D and Macchiavello C 2014 Hypergraph states in Grover’s quantum search algorithm *Phys. Scr.* **2014** 014036

[21] Bianconi G and Rahmede C 2016 Network geometry with flavor: from complexity to quantum geometry *Phys. Rev. E* **93** 032315

[22] Kato Y and Yamamoto N 2014 Structure identification and state initialization of spin networks with limited access *New J. Phys.* **16** 023024

[23] Jurcevic P, Lanyon P, Hauke B P, Hempel C, Zoller P, Blatt R and Roos C F 2014 Quasiparticle engineering and entanglement propagation in a quantum many-body system *Nature* **511** 202–5

[24] Elliott T J and Mekhov I B 2016 Engineering many-body dynamics with quantum light potentials and measurements *Phys. Rev. A* **94** 013614

[25] Moore D W, Tufarelli T, Paternostro M and Ferraro A 2016 Quantum state reconstruction of an oscillator network in an optomechanical setting *Phys. Rev. A* **94** 053811

[26] Gokler C, Lloyd S, Shor P and Thompson K 2017 Efficiently controllable graphs *Phys. Rev. Lett.* **118** 260501

[27] Yung M-H and Bose S 2005 Perfect state transfer, effective gates, and entanglement generation in engineered bosonic and fermionic networks *Phys. Rev. A* **71** 032310

[28] Kostak V, Nikolopoulos G M and Jex I 2007 Perfect state transfer in networks of arbitrary topology and coupling configuration *Phys. Rev. A* **75** 042319

[29] Manzano G, Galve F and Zambrini R 2013 Avoiding dissipation in a system of three quantum harmonic oscillators *Phys. Rev. A* **87** 032114

[30] Arenas A, Díaz-Guilera A, Kurths J, Moreno Y and Zhou C 2008 Synchronization in complex networks *Phys. Rep.* **469** 93–153

[31] Dorfler F, Chertkov M and Bullo F 2013 Synchronization in complex oscillator networks and smart grids *Proc. Natl Acad. Sci.* **110** 2005–10

[32] Manzano G, Galve F, Giorgi G L, Hernández-García E and Zambrini R 2013 Synchronization, quantum correlations and entanglement in oscillator networks *Sci. Rep.* **3** 1439

[33] Bellomo B, Giorgi G L, Palma G M and Zambrini R 2017 Quantum synchronization as a local signature of super- and subradiance *Phys. Rev. A* **95** 043807

[34] Nokkala J, Galve F, Zambrini R, Maniscalco S and Piilo J 2016 Complex quantum networks as structured environments: engineering and probing *Sci. Rep.* **6** 26861

[35] Roslund J, Medeiros de Araújo R, Jiang S, Fabre C and Treps N 2014 Wavelength-multiplexed quantum networks with ultrafast frequency combs *Nat. Photon.* **8** 109

[36] Sameti M, Potočnik A, Browne D E, Wallraff A and Hartmann M J 2017 Superconducting quantum simulator for topological order and the toric code *Phys. Rev. A* **95** 042330

[37] Robens C, Alt W, Meschede D, Emary C and Alberti A 2015 Ideal negative measurements in quantum walks disprove theories based on classical trajectories *Phys. Rev. X* **5** 011003

[38] Qiang X, Loke T, Montanaro A, Aungskunsiri K, Zhou X, O’Brien J L, Wang J B and Matthews J C F 2016 Efficient quantum walk on a quantum processor *Nat. Commun.* **7** 11511

[39] Gerke S, Sperling J, Vogel W, Cai Y, Roslund J, Treps N and Fabre C 2015 Full multipartite entanglement of frequency-comb Gaussian states *Phys. Rev. Lett.* **114** 050501

[40] Cai Y, Roslund J, Ferrini G, Arzani F, Xu X, Fabre C and Treps N 2017 Multimode entanglement in reconfigurable graph states using optical frequency combs *Nat. Commun.* **8** 15645

[41] Cardillo A, Galve F, Zueco D and Gómez-Gardeñes J 2013 Information sharing in quantum complex networks *Phys. Rev. A* **87** 052312

[42] Rivas A, Huelga S F and Plenio M B 2014 Quantum non-Markovianity: characterization, quantification and detection *Rep. Prog. Phys.* **77** 094001

[43] Breuer H-P, Laine E-M, Piilo J and Vacchini B 2016 Colloquium *Rev. Mod. Phys.* **88** 021002

[44] de Vega I and Alonso D 2017 Dynamics of non-Markovian open quantum systems *Rev. Mod. Phys.* **89** 015001

[45] Liu B-H, Li L, Huang Y-F, Li C-F, Guo G-C, Laine E-M, Breuer H-P and Piilo J 2011 Experimental control of the transition from Markovian to non-Markovian dynamics of open quantum systems *Nat. Phys.* **7** 931–4

[46] Vasile R, Galve F and Zambrini R 2014 Spectral origin of non-Markovian open-system dynamics: a finite harmonic model without approximations *Phys. Rev. A* **89** 022109

[47] Paparo G D and Martin-Delgado M A 2012 Google in a quantum network *Sci. Rep.* **2** 444

[48] Walschaers M, Schlawin F, Wellens T and Buchleitner A 2016 Quantum transport on disordered and noisy networks: an interplay of structural complexity and uncertainty *Annu. Rev. Condens. Matter Phys.* **7** 223–48

[49] Braunstein S L 2005 Squeezing as an irreducible resource *Phys. Rev. A* **71** 055801

[50] De Gosson M A 2006 *Symplectic Geometry and Quantum Mechanics* (Basel: Birkhäuser) Proposition 2.13 and its proof (https://doi.org/10.1007/3-7643-7575-2)

[51] Medeiros de Araújo R, Roslund J, Cai Y, Ferrini G, Fabre C and Treps N 2014 Full characterization of a highly multimode entangled state embedded in an optical frequency comb using pulse shaping *Phys. Rev. A* **89** 053828

[52] Arzani F, Fabre C and Treps N 2018 Versatile engineering of multimode squeezed states by optimizing the pump spectral profile in spontaneous parametric down-conversion *Phys. Rev. A* **97** 033808

[53] Barabási A-L and Albert R 1999 Emergence of scaling in random networks *Science* **286** 509–12

[54] Watts D J and Strogatz S H 1998 Collective dynamics of ‘small-world’ networks *Nature* **393** 440–2

[55] Scannell J, Burns G, Hilgetag C, O’Neil M and Young M 1999 The connectional organization of the cortico-thalamic system of the cat *Cerebral Cortex* **9** 277–99

[56] Lusseau D, Schneider K, Boisseau O J, Haase P, Slooten E and Dawson S M 2003 The bottlenose dolphin community of doubtful sound features a large proportion of long-lasting associations *Behav. Ecol. Sociobiol.* **54** 396–405

[57] Davies E 1976 *Quantum Theory of Open Systems* (New York: Academic)

[58] Breuer H-P and Petruccione F 2007 *The Theory of Open Quantum Systems* (Oxford: Oxford University Press)

[59] Agarwal G S 1971 Entropy, the Wigner distribution function, and the approach to equilibrium of a system of coupled harmonic oscillators *Phys. Rev. A* **3** 828–31

[60] Yokoyama S, Ukai R, Armstrong S C, Sornphiphatphong C, Kaji T, Suzuki S, Yoshikawa J i, Yonezawa H, Menicucci N C and Furusawa A 2013 Ultra-large-scale continuous-variable cluster states multiplexed in the time domain *Nat. Photon.* **7** 982
[61] Chen M, Menicucci N C and Pfister O 2014 Experimental realization of multipartite entanglement of 60 modes of a quantum optical frequency comb *Phys. Rev. Lett.* **112** 120505

[62] Gao J, Barzel B and Barabási A-L 2016 Universal resilience patterns in complex networks *Nature* **530** 307

[63] Hillery M, Bužek V and Berthiaume A 1999 Quantum secret sharing *Phys. Rev. A* **59** 1829–34

[64] Ra Y-S, Jacquard C, Dufour A, Fabre C and Treps N 2017 Tomography of a mode-tunable coherent single-photon subtractor *Phys. Rev. X* **7** 031012

[65] Walschaers M, Fabre C, Parigi V and Treps N 2017 Entanglement and Wigner function negativity of multimode non-Gaussian states *Phys. Rev. Lett.* **119** 183601

[66] Chabaud U, Douce T, Markham D, Van Loock P, Kashefi E and Ferrini G 2017 Continuous-variable sampling from photon-added or photon-subtracted squeezed states *Phys. Rev. A* **96** 062307

[67] Walschaers M, Diaz J F-d-C, Mulet R and Buchleitner A 2013 Optimally designed quantum transport across disordered networks *Phys. Rev. Lett.* **111** 180601

[68] Ferrini G, Gazeau J P, Coudreau T, Fabre C and Treps N 2013 Compact Gaussian quantum computation by multi-pixel homodyne detection *New J. Phys.* **15** 093015