Reconfigurable optical implementation of quantum complex networks

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

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 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}={\bf{p}}{{\prime} }^{T}{\bf{p}}^{\prime} /2+{\bf{q}}{{\prime} }^{T}{\bf{Aq}}^{\prime}$ where we have vectors of momentum and position operators ${\bf{p}}{{\prime} }^{T}=\{{p}_{1}^{{\prime} },\,\ldots ,\,{p}_{N}^{{\prime} }\}$, ${\bf{q}}{{\prime} }^{T}=\{{q}_{1}^{{\prime} },\,\ldots ,\,{q}_{N}^{{\prime} }\}$ for N harmonic oscillators. The matrix A has diagonal elements which can be written by introducing the effective frequencies ${\tilde{\omega }}_{i}$ as ${{\bf{A}}}_{{ii}}={\tilde{\omega }}_{i}^{2}/2={\omega }_{i}^{2}/2+{\sum }_{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 ${{\boldsymbol{\Delta }}}_{\tilde{\omega }}={\rm{diag}}\{{\tilde{\omega }}_{1}^{2},\,\ldots ,\,{\tilde{\omega }}_{N}^{2}\}$, as ${\bf{A}}={{\boldsymbol{\Delta }}}_{\tilde{\omega }}/2-{\bf{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 ${{\bf{q}}}^{T}={\bf{q}}{{\prime} }^{T}\sqrt{{{\boldsymbol{\Delta }}}_{\omega }}$ and ${{\bf{p}}}^{T}={\bf{p}}{{\prime} }^{T}\sqrt{{{\boldsymbol{\Delta }}}_{\omega }^{-1}}$ with ${{\boldsymbol{\Delta }}}_{\omega }={\rm{diag}}\{{\omega }_{1},\,\ldots ,\,{\omega }_{N}\}$. So the network Hamiltonian is rewritten as

$$\begin{eqnarray}&&{H}_{E}=\displaystyle \frac{{{\bf{p}}}^{T}{{\boldsymbol{\Delta }}}_{\omega }{\bf{p}}}{2}+{{\bf{q}}}^{T}\sqrt{{{\boldsymbol{\Delta }}}_{\omega }^{-1}}{\bf{A}}\sqrt{{{\boldsymbol{\Delta }}}_{\omega }^{-1}}{\bf{q}}.\end{eqnarray} \tag{ 1 }$$

By diagonalizing the symmetric matrix A with an orthogonal matrix K, such as ${\bf{K}}{({{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}/\sqrt{2})}^{2}{{\bf{K}}}^{T}={\bf{A}}$, where Δ_Ω is a diagonal matrix $=\,{\rm{diag}}\{{{\rm{\Omega }}}_{1},\,\ldots ,\,{{\rm{\Omega }}}_{N}\}$, we can define decoupled oscillators, also named normal modes, with quadratures

$$\begin{eqnarray}{\bf{Q}} & = & \sqrt{{{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}}{{\bf{K}}}^{T}\sqrt{{{\boldsymbol{\Delta }}}_{\omega }^{-1}}{\bf{q}},\\ {\bf{P}} & = & \sqrt{{{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}^{-1}}{{\bf{K}}}^{T}\sqrt{{{\boldsymbol{\Delta }}}_{\omega }}{\bf{p}}.\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&{H}_{E}=({{\bf{P}}}^{T}{{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}{\bf{P}}+{{\bf{Q}}}^{T}{{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}{\bf{Q}})/2=\displaystyle \sum _{i}^{N}{{\rm{\Omega }}}_{i}({P}_{i}^{2}+{Q}_{i}^{2})/2.\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}&&\left(\begin{array}{c}{\bf{q}}(t)\\ {\bf{p}}(t)\end{array}\right)=\left(\begin{array}{cc}{T}_{1}{D}_{{\rm{\cos }}}^{{\rm{\Omega }}} & {T}_{1}{D}_{\sin }^{{\rm{\Omega }}}\\ -{T}_{2}{D}_{\sin }^{{\rm{\Omega }}} & {T}_{2}{D}_{{\rm{\cos }}}^{{\rm{\Omega }}}\end{array}\right)\left(\begin{array}{c}{\bf{Q}}(0)\\ {\bf{P}}(0)\end{array}\right)\end{eqnarray} \tag{ 4 }$$

and then

$$\begin{eqnarray}&&\left(\begin{array}{c}{\bf{q}}(t)\\ {\bf{p}}(t)\end{array}\right)=\left(\begin{array}{cc}{T}_{1}{D}_{{\rm{\cos }}}^{{\rm{\Omega }}}{T}_{1}^{-1} & {T}_{1}{D}_{\sin }^{{\rm{\Omega }}}{T}_{2}^{-1}\\ -{T}_{2}{D}_{\sin }^{{\rm{\Omega }}}{T}_{1}^{-1} & {T}_{2}{D}_{{\rm{\cos }}}^{{\rm{\Omega }}}{T}_{2}^{-1}\end{array}\right)\left(\begin{array}{c}{\bf{q}}(0)\\ {\bf{p}}(0)\end{array}\right),\end{eqnarray} \tag{ 5 }$$

where ${T}_{1}=\sqrt{{{\boldsymbol{\Delta }}}_{\omega }}{\bf{K}}\sqrt{{{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}^{-1}}$ and ${T}_{2}=\sqrt{{{\boldsymbol{\Delta }}}_{\omega }^{-1}}{\bf{K}}\sqrt{{{\boldsymbol{\Delta }}}_{{\rm{\Omega }}}}$ and we have introduced the diagonal matrices ${D}_{{\rm{\sin }}{ii}}^{{\rm{\Omega }}}={\rm{\sin }}({{\rm{\Omega }}}_{i}t)$ and ${D}_{{\rm{\cos }}{ii}}^{{\rm{\Omega }}}={\rm{\cos }}({{\rm{\Omega }}}_{i}t)$. The equations can be given in the compact form

$$\begin{eqnarray}{\bf{x}}(t) & = & S{{\prime} }_{{\bf{V}}}{\bf{X}}(0),\\ {\bf{x}}(t) & = & {S}_{{\bf{V}}}{\bf{x}}(0),\end{eqnarray} \tag{ 6 }$$

where ${{\bf{X}}}^{T}=\{{{\bf{Q}}}^{T},{{\bf{P}}}^{T}\}$ and x^T = {q^T, 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 {Q(0)^T, P(0)^T} in the normal modes picture or, equivalently, by {q(0)^T, p(0)^T} in the bare oscillators picture. Finally the implementation involves reconfigurable optical processes, able to drive the quadratures evolution S'_V or S_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'_V and S_V could be more efficiently implemented in the experimental setup. The matrices S'_V and S_V in equation (6) are symplectic and they can then be factored via the Bloch–Messiah (B–M) decomposition [49], for example:

$$\begin{eqnarray}&&{S}_{{\bf{V}}}={R}_{1}^{T}{{\rm{\Delta }}}^{{\rm{sq}}}{R}_{2}.\end{eqnarray} \tag{ 7 }$$

Δ^sq being a diagonal matrix and R₁, R₂ 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 ω_i = ω₀ = 0.25 for the three different topologies. The network matrices S_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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the complexity of B–M decomposition is of the order of O(N³) [50]. It takes around 0.03 s to implement the B–M decomposition for networks of 50 nodes with a standard laptop.

The matrices ${{\rm{\Delta }}}^{{\rm{sq}}}$, R₁ and R₂ 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 to ${S}_{{\bf{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 λ = 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 ω_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 χ⁽²⁾ crystal, is $H={\imath }g{\sum }_{m,n}{{ \mathcal L }}_{m,n}{a}_{m}^{\dagger }{a}_{n}^{\dagger }+{\rm{h}}.{\rm{c}}.;$ it couples pairs of frequencies ω_m and ω_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 α_p at frequency ω_p = ω_m + ω_n. ${{ \mathcal 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 operators $\{{q}_{1},{q}_{2},\,..{q}_{N},{p}_{1},{p}_{2},\,\ldots {p}_{N}\}=\{{\bf{q}},{\bf{p}}\}={{\bf{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_exX_i which can be decomposed via B–M in

$$\begin{eqnarray}&&{S}_{{\rm{ex}}}={R}_{1{\rm{ex}}}^{T}{{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}{R}_{2{\rm{ex}}}.\end{eqnarray} \tag{ 8 }$$

If the input modes X_i are in the vacuum state the R_2ex term can be neglected

$$\begin{eqnarray}&&{{\bf{X}}}_{f}={R}_{1{\rm{ex}}}{{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}{{\bf{X}}}_{i}.\end{eqnarray} \tag{ 9 }$$

The squeezing parameters in ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{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.

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 networks 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 ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{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 ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}$ [52]. Note that the optimization of ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}$ also influences the spectral shape of the supermodes, and thus the squeezed modes basis.

Zoom In Zoom Out Reset image size

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 Δ^sq acts on the supermode basis, then the quadratures of the mth network oscillator can be addressed by applying the basis change defined by the mth row of the matrix ${R}_{1}^{T}$. This means that we can measure the quadratures q_m(t) and p_m(t) of the mth network oscillator by performing homodyne detection with a LO shape defined by ${u}_{{\rm{LO}}-m}(\omega )={\sum }_{j=1}^{N}{R}_{1{mj}}^{T}{{SP}}_{j}(\omega )+{\imath }{\sum }_{j=N+1}^{2N}{R}_{1{mj}}^{T}{{SP}}_{j-N}(\omega )$, where SP_i(ω) are the supermodes. Figures 4(A) and (B) show the scheme for the experimental control of ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}$ 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 ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{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}_{{\rm{ex}}}={R}_{1}^{T}$ ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}={{\rm{\Delta }}}^{{\rm{sq}}}$, and discarding R₂ 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₀ which prepares the modes in the desired state, followed by network evolution S_V. The effective total evolution to be simulated is

$$\begin{eqnarray}&&{\bf{x}}(t)={S}_{{\bf{V}}}{S}_{0}{\bf{x}}(0)={S}_{{\rm{ef}}}{\bf{x}}(0),\end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray}&&{S}_{{\rm{ef}}}={R}_{1{\rm{ef}}}{{\rm{\Delta }}}_{{\rm{ef}}}^{{\rm{sq}}}{R}_{2{\rm{ef}}},\end{eqnarray} \tag{ 11 }$$

experimentally implement ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}={{\rm{\Delta }}}_{{\rm{ef}}}^{{\rm{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 ={\rm{diag}}\{{\rm{\hslash }}\omega {\rm{\cosh }}(2r)/2,{\rm{\hslash }}{\rm{\cosh }}(2r)/2\omega \}$, which describes a thermal state with $\cosh (2r)\to \coth ({\rm{\hslash }}\omega /2T)$. For T/ω > 1 the required squeezing grows as $r\simeq {\rm{acosh}}(2T)/2$ and this scales like $\mathrm{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 realizable ones (${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{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 Δ^sq if the exact required basis change R₁ 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 (${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}$) 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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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},\,..{\omega }_{{SM}}\}$, each of them interacting with one node according to the Hamiltonian ${H}_{{Ir}}=-{{kq}}_{{Sr}}{q}_{i}/\sqrt{{\omega }_{{Sr}}{\omega }_{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: ${H}_{{Ir}}=-{{kq}}_{{Sr}}{\sum }_{j}{{\bf{K}}}_{{ij}}{Q}_{j}/\sqrt{{\omega }_{{Sr}}{{\rm{\Omega }}}_{j}}$. The new total Hamiltonian again can be written in diagonal form by orthogonal diagonalization of ${\bf{B}}={\bf{O}}{({{\boldsymbol{\Delta }}}_{f}/\sqrt{2})}^{2}{{\bf{O}}}^{T}$, the matrix analogue to A containing the extra terms given by the interaction. Δ_f is a diagonal matrix ${\rm{diag}}\{{f}_{1},\,\ldots ,\,{f}_{N+M}\}$.

$$\begin{eqnarray}&&{H}_{S+E+I}=\displaystyle \sum _{j}^{N+M}{f}_{j}({{ \mathcal P }}_{j}^{2}+{{ \mathcal Q }}_{j}^{2})/2.\end{eqnarray} \tag{ 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 ${{\bf{X}}}^{T}=\{{\bf{Q}},{{\bf{q}}}_{S},{\bf{P}},{{\bf{p}}}_{S}\}=\{{\bf{Q}},{q}_{S1},\,..\,,{q}_{{SM}},{\bf{P}},{p}_{S1},\,..,{p}_{{SM}}\}$

$$\begin{eqnarray}&&\left(\begin{array}{c}{\bf{Q}}(t)\\ {{\bf{q}}}_{S}(t)\\ {\bf{P}}(t)\\ {{\bf{p}}}_{S}(t)\end{array}\right)=\left(\begin{array}{cc}{O}_{1}{D}_{{\rm{\cos }}}{O}_{1}^{-1} & {O}_{1}{D}_{{\rm{\sin }}}{O}_{2}^{-1}\\ -{O}_{2}{D}_{{\rm{\sin }}}{O}_{1}^{-1} & {O}_{2}{D}_{{\rm{\cos }}}{O}_{2}^{-1}\end{array}\right)\left(\begin{array}{c}{\bf{Q}}(0)\\ {{\bf{q}}}_{S}(0)\\ {\bf{P}}(0)\\ {{\bf{p}}}_{S}(0)\end{array}\right),\end{eqnarray} \tag{ 13 }$$

where we introduced the diagonal matrices ${D}_{{\rm{\cos }}{ii}}={\rm{\cos }}({f}_{i}t)$ and ${D}_{{\rm{\sin }}{ii}}={\rm{\sin }}({f}_{i}t)$, and the matrices ${O}_{1}=\sqrt{{{\boldsymbol{\Delta }}}_{{\rm{\Omega }},{\omega }_{S}}}{\bf{O}}\sqrt{{{\boldsymbol{\Delta }}}_{f}^{-1}}$ and ${O}_{2}=\sqrt{{{\boldsymbol{\Delta }}}_{{\rm{\Omega }},{\omega }_{S}}^{-1}}{\bf{O}}\sqrt{{{\boldsymbol{\Delta }}}_{f}}$, being ${{\boldsymbol{\Delta }}}_{{\rm{\Omega }},{\omega }_{S}}={\rm{diag}}\{{{\rm{\Omega }}}_{1},\,\ldots ,\,{{\rm{\Omega }}}_{N},{\omega }_{S1}...{\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.

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(ω), 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 $\langle n\rangle$ of the probe having frequency ω_S allows to approximate J(ω_S) from

$$\begin{eqnarray}&&J({\omega }_{S})=\displaystyle \frac{{\omega }_{S}}{t}\mathrm{ln}\left(\displaystyle \frac{N({\omega }_{S})-\langle n(0)\rangle }{N({\omega }_{S})-\langle n(t)\rangle }\right),\end{eqnarray} \tag{ 14 }$$

where $N({\omega }_{S})={({{\rm{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(ω_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 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₀ = Δ', and then we will implement the B–M decomposition of S_eff = SΔ'.

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{}}{\rm{ef}}},{{\rm{\Delta }}}_{{\rm{ef}}}^{{\rm{sq}}}$. The probe mode is measured by pulse shaping the LO and one value of J(ω_S) is obtained. The procedure is repeated for other ω_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 ω_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₁, R₂. 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 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 ω_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. Figure 9 shows in purple the target squeezing, i.e. the diagonal elements of ${{\rm{\Delta }}}_{{\rm{ef}}}^{{\rm{sq}}}$, in the case of a periodic chain, with parameters set as in figure 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 ${{\rm{\Delta }}}_{{\rm{ef}}}^{{\rm{sq}}}$, this is the reason why all the 10 purple lists almost coincide. The diagonal elements of three experimentally achievable ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{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 ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}$. In figures 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 ${{\rm{\Delta }}}_{{\rm{ef}}}^{{\rm{sq}}}$ with the experimentally achievable ${{\rm{\Delta }}}_{{\rm{ex}}}^{{\rm{sq}}}$ (green, red and brown dots). We show the normalized value of J(ω) (top) and the raw values after applying equation (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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We finally report, in figure 12, the first 10 non-zero values of the N+1th row of ${R}_{1f}$ at different frequencies ω_S. The N+1th 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 crystals allow for reduced LO bandwidth. In any case all the LO shapes calculated in the figures 13–15 are experimentally doable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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)\mathrm{log}(\mu +1/2)-(\mu -1/2)\mathrm{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.

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

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.