Reconfigurable optical implementation of quantum complex networks

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 in a non-regular topology, have been 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 geometries. However, the question on how to produce and control general quantum complex networks in experimental laboratory has remained open. 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. Overall, our general method for implementing quantum complex networks with reconfigurable set-up has potential to define an experimental playground for designing and controlling complex networks -- and dynamics therein -- for several quantum physical frameworks.

During the last twenty years, network theory has experienced remarkable progress and revolutionised the research in diverse disciplines ranging, e.g., from technology to social sciences and biology [1][2][3][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][11][12][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]. 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][23][24][25][26], to the generation and protection of quantum resources [25,[27][28][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, higly 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 Ref. [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][43][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.

NETWORK DYNAMICS
We consider bosonic quantum complex networks composed by an ensemble of quantum harmonic oscillators, with frequencies ω 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 sys- tem 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 = p T p /2 + q T Aq where we have vectors of momentum and position operators p T = {p 1 , ..., p N }, q T = {q 1 , ..., q N } for N harmonic oscillators. The matrix A has diagonal elements which can be written by introducing the effective frequencies ω i as A ii =ω 2 i /2 = ω 2 i /2 + j v ij /2, and off-diagonal elements A i =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 ∆ω = diag{ω 2 1 , ...,ω 2 N }, as A = ∆ω/2 − V/2. As the bosonic oscillators are intended to describe electromagnetic fields we use a more common formalism in quantum optics by substituting q and p with the renormalized quadratures q and p by defining q T = q T √ ∆ ω and p T = p T ∆ −1 ω with ∆ ω = diag{ω 1 , ..., ω N }. So the network Hamiltonian is rewritten as By diagonalizing the symmetric matrix A with an orthogonal matrix K, such as K(∆ Ω / √ 2) 2 K T = A, where ∆ Ω is a diagonal matrix = diag{Ω 1 , ..., Ω N }, we can define decoupled oscillators, also named normal modes, ∆ sq being a diagonal matrix and R 1 , R 2 symplectic and orthogonal matrices. In Fig. 1 three different networks composed by N=51 ( Fig.1 A) or N=50 ( Fig.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 ω i = ω 0 = 0.25 for the three different topologies. The network matrices S V calculated at different times from the evolution in Eq. 5 are shown in Fig. 2, while the B-M decomposition at time t = 50 for the three network is shown in Fig. 3. Note that the complexity of B-M decomposition is of the order of O(N 3 ) [50]. It takes around 0.03 second to implement the B-M decomposition for networks of 50 nodes with a standard laptop.

MAPPING THE NETWORK TO A MULTIMODE EXPERIMENTAL PLATFORM
The matrices ∆ sq , R 1 and R 2 in eq. (7) have a wellknown 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 multi-mode 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 fre-  quency 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 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 modelocked lasers emitting pulses from 30 fs to 150 fs with repetition rates from 70 MHz to 150 MHz at a wavelength centred around λ = 800nm. 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 ω p = ω p0 + p · ω RR , where ω RR is the repetition rate of the original laser source, ω p0 is the pump central frequency and p an integer. The Hamiltonian describing the parametric process, which is parametric down conversion in a χ (2) crystal, is H = ıg m,n L m,n a † m a † n + h.c.; it couples pairs of frequencies ω m and ω n via the term L m,n which is the product of crystal's phase matching function f m,n and of the complex pump amplitude α p at frequency ω p = ω m + ω n . L m,n describes the probability amplitude that a photon at frequency ω p is converted into two daughter photons at frequencies ω n and ω 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 opera- . A natural choice for the basis would be the frequency components of the comb spectrum. The parametric Hamiltonian leads to a symplectic transformation If the input modes X i are in the vacuum state the R 2ex term can be neglected The squeezing parameters in ∆ sq ex can be derived from the eigenvalues of 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 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.
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 don't need to know the physical phenomenon that generates the networks dynamics or the nature of the oscillators in the network (they could be mechanical or optical oscillators) . The first is controlled by shaping the pump of the parametric process (here in cavity as [35]), the experimental parameters of the pulse shaper are found by running an optimization procedure based on an evolutionary algorithm. The basis change defines the modes which carry the information on the oscillators evolution as a combination of the supermodes shape.
or even the specific picture that can be chosen to describe the Hamiltonian evolution. In the case described by eq. (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 ∆ sq ex 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 2-D 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 ∆ sq ex [52]. Note that the optimization of ∆ sq ex also influences the spectral shape of the supermodes, and thus the squeezed modes basis.
The linear transformation R 1ex can be tuned ade-quately 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 ∆ 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 T 1 . 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 Figures 4 A and B show the scheme for the experimental control of ∆ sq ex via the control of the 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 ∆ sq ex 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 T 1 ∆ sq ex = ∆ 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 where the initial states are vacua. Then via B-M decomposition we can set experimentally implement ∆ sq ex = ∆ sq ef 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 fulfill 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 → ∞), 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 σ = diag{ ω cosh(2r)/2, cosh(2r)/2ω}, which describes a thermal state with cosh(2r) → coth( ω/2T ). , characterized by a power law distribution of the degree with connection parameter set to 3. E is a 50-nodes net work derived from the WattsStrogatz model [54], characterized by short average path lengths and high clustering, the rewiring probability is set to be 0.2. F and G are networks retrieved from a public repository of real-world complex networks. F is a 52 nodes-network which represents the connectional organization of the cortico-thalamic system of the cat, as studied by J.W. Scannell et al. [55]. F is a 62-nodes network representing the social network of dolphins, derived by D. Lusseau [56]. The coupling strength between the nodes for the four networks is set to be vij = 0.05. acosh(2T )/2 and this scales like log T , e.g. r ∼ 2 for T /ω ∼ 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 π/4 rotation in the 2 × 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 distance between the theoretical values (the diagonal of ∆ sq ) and the experimentally realisable ones (∆ sq ex ) 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 ∆ 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 Fig. 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 Fig. 5 are obtained by setting two different crystal lengths: in the case of the shorter crystal (1.5 mm) the experimentally achievable squeezing matrix (∆ sq ex ) is closer to the theoretical one (∆ 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 Fig. 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 Fig. ) 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 Fig. 7. We study the Barabási-Albert model [53], characterized by a power law of the degree distribution, and the WattsStrogatz model [54], which describes networks with small-world properties as clustering. Additionally, we consider two graphical structures derived from real natural networks, studied by J.W. Scannell et al. [55] and D. Lusseau [56], that we call Scannel and Lusseau networks. The simulations are shown in Fig. 8 and they show the same performances of the simulations in Fig. 5 with a similar trade-off be-tween relative errors limitation and reduction of the LO bandwidth.

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 {ω S1 , ω S2 , ..ω SM }, each of them interacting with one node according to the Hamiltonian H Ir = −kq Sr q i / √ ω Sr ω i . The interaction couples the supplementary oscillators with any of the normal mode, as we can see from the interaction Hamiltonian in the decoupled quadrature picture: 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 X T = {Q, q S , P, p S } = {Q, q S1 , .., q SM , P, p S1 , .., (13) where we introduced the diagonal matrices D cos ii = cos(f i t) and D sin ii = sin(f i t), and the matrices O 1 = (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.

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 realisable operations, of the probing technique for the network spectral density as proposed in [34].
The spectral density J(ω), 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(ω). Indeed, provided that the coupling between the probe and the environment is sufficiently weak, measuring the photon number operator n of the probe having frequency ω S allows to approximate J(ω S ) from where N (ω S ) = (e (ω 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 Eq. (14) to probe J(ω 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 n(t) = 1/2 q S (t) 2 + 1/2 p S (t) 2 − 1/2. The quadratures of the probe {q S (t), p S (t)} appear in the left hand side of the Eq. (13), which can be used to derive their expectation value at any time. We have to stress here that when the initial mode basis X(0) is defined, the B-M decomposition of S (the matrix appearing on the left side of Eq. (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 n . We then will use a first operation able to set a non-zero value for the mean photon number of the probe n(0) , as for example a squeezing operation S 0 = ∆ , and then we will implement the B-M decomposition of S ef f = S∆ The protocol is the following one: given a network structure, a probe with frequency ω S is chosen and the evolution of the total system is computed at an appropiate time t, the experimental parameters are set according to the calculated B-M giving R 1ef , ∆ sq ef . The probe mode is measured by pulse shaping the LO and one value of mode index ii FIG. 9. Purple: diagonal elements of the target ∆ sq ef ; 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 mm 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.
J(ω S ) is obtained. The procedure is repeated for other ω S in order to retreive 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 Eq. (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 ω 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. ω i = ω and v ij = g. If we fix the set of values {ω, g, ω S , k} we find almost the same matrix ∆ 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 {ω, ω 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 graph- ical 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 ω S on the squeezing matrix ∆ sq is very small. This simplifies the experimental procedure of the probing protocol because scanning ω 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. Fig. 9 shows in purple the target squeezing, i.e. the diagonal elements of ∆ sq ef , in the case of a periodic chain, with parameters set as in Fig. 1; ten different purple lists are superposed corresponding to ten different values of ω S . The configuration in which we probe with an initially squeezed oscillator requires only one significantly squeezed mode in the matrix ∆ sq ef , this is the reason why all the 10 purple lists almost coincide. The diagonal elements of three experimentally achievable ∆ sq ex 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 ∆ sq ex . In Fig. 10 and 11 we compare the results for the analytically calculated J(ω) (blue lines), the one derived from the exact probing (violet circles), and the ones calculated by replacing the exact ∆ sq ef with the experimentally achievable ∆ sq ex (green, red and brown dots). We show the normalized value of J(ω) (top) and the raw values after applying Eq. FIG. 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 Fig. 4 B. (14). The experimentally achievable configurations allow for the retrieval of the peculiar shape J(ω) 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(ω) are quite low, while larger values are obtained when optimization via pump shaping is included.
We finally report, in Fig.12, the first 10 non-zero values of the N+1-th row of R 1f at different frequencies ω 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 n(t) . Examples of the LO shapes are shown in Fig. 13, 14  tures displayed by the averaged single node entropy S h , 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) σ = {{α, γ}, {γ, β}}, and reads S = (µ + 1/2) log(µ + 1/2) − (µ − 1/2) log(µ − 1/2) and µ = αβ − γ 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.

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 reconfiguratibility 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 theroretical 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(ω) 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(ω) 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][43][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.