Disorder and dephasing as control knobs for light transport in optical fiber cavity networks

Transport phenomena represent a very interdisciplinary topic with applications in many fields of science, such as physics, chemistry, and biology. In this context, the possibility to design a perfectly controllable experimental setup, where to tune and optimize its dynamics parameters, is a challenging but very relevant task to emulate, for instance, the transmission of energy in light harvesting processes. Here, we experimentally build a scalable and controllable transport emulator based on optical fiber cavity networks where the system noise parameters can be finely tuned while maximizing the transfer efficiency. In particular, we demonstrate that disorder and dephasing noise are two control knobs allowing one to play with constructive and destructive interference to optimize the transport paths towards an exit site. These optical setups, on one side, mimic the transport dynamics in natural photosynthetic organisms and, on the other, are very promising platforms to artificially design optimal nanoscale structures for novel, more efficient, clean energy technologies.

Very recently several theoretical and experimental studies have shown that the remarkably high efficiency (almost 100%) of the excitation energy transfer through a network of chromophores in photosynthetic systems seems to be the result of an intricate interplay between quantum coherence and noise [1][2][3][4][5][6][7][8][9][10][11][12][13] . Indeed, the presence of coherence leads to a very fast delocalization of the excitation that can hence exploit several paths to the target site (named also as sink or reaction center). However, since destructive interference among different pathways and the energy gaps between different sites are obstacles to the transmission of energy, this regime is not optimal by itself. In fact, the additional and unavoidable presence of disorder and noise, which is usually assumed to be deleterious for the transport properties, here instead positively affects the transmission efficiency. This can be explained in terms of the inhibition of destructive interference and the opening of additional pathways for excitation transfer 5,14,15 . This mechanism, known as Noise-Assisted Transport (NAT), has been recently observed in some physical platforms [16][17][18][19] , but a deep analysis of the underlying contributions, such as interference, disorder and dephasing, is still missing from the experimental point of view.
More in general, the possibility to experimentally realize simple test platforms being able to mimic transport on complex networks and to reproduce NAT effects, would be of high interest from several perspectives. First of all, it would clarify the role of the different disorder, interference and noise contributions in the transport behavior of natural photosynthetic complexes. Secondly, it might represent a model of different complex networks where the role of topology can be further investigated. Finally, one could play with this platform in order to engineer new artificial molecular structures where all these control knobs are optimized to achieve some desired tasks, such as the maximization of the transferred energy or its temporary storage in some part of the network. Such a simple setup would have a remarkable advantage with respect to real biological samples or expensive artificial systems, which are very difficult to manipulate both in their geometry and in the system parameters, since these aspects are often governed by specific bio-chemical laws.
For these reasons, we have recently realized an optical platform, entirely based on fiber-optic components, which has allowed us to provide the first experimental observation of the typical NAT peak in the dependence of the network transmission rate as a function of the amount of noise introduced into the system 18 . In such a simple and scalable kind of setup, which has variable and controllable parameters, the coherent propagation of excitons in a N-site quantum network is simulated by the propagation of photons in a network of N coupled optical cavities. The optical cavities are based on Fiber Bragg Grating (FBG) resonators 20,21 , which have been mainly exploited for sensing applications so far 22,23 , also because they are characterized by a straightforward alignment and easy tunability by tiny deformations of the fiber section within the Bragg mirrors. Each j-site of the network (corresponding to each chromophore in a photosynthetic complex) is represented by a FBG resonator, while the local site excitation energy is associated to the resonance frequency of the cavity (ω j ).
Here we provide a complete theoretical and experimental analysis of NAT effects for all the different configurations of our optical model system, demonstrating that constructive and destructive interference, static disorder and dephasing noise can be successfully exploited as feasible control knobs to manipulate the transport behavior of coupled structures and, for example, to optimize the transmission rate.

Results
Network of fiber-optic resonators. We build up an experimental fiber-optic setup to reproduce energy transport phenomena, including the NAT effect, in the 4-site network shown in the inset (a) of Fig. 1 18 . A detailed scheme of the experimental apparatus is shown in Fig. 1. It can be schematically described by the Mach-Zehnder setup of inset (b), which presents two main differences with respect to a standard Mach-Zehnder interferometer: (i) the insertion of a FBG resonator in each path of the interferometer, and (ii) the presence of two additional mirrors (M1 and M2) at the normally unused input and output port of the interferometer. The resonators FBG1 and FBG2 represent the sites 1 and 2 with variable local excitation energy ω 1 and ω 2 . The role of the other two sites 0 and 3 is played by two fiber optic couplers (C), which represent two sites with fixed local excitation energy resonant with the energy of the propagating excitation ω S (ω 0 = ω 3 = ω S ). The presence of the two additional mirrors M1 and M2 makes it possible to couple sites 1 and 2, since the light reflected by the resonators is partially re-inserted into the network by M1, while transmitted light is partially recycled by M2. The setup is entirely based on single-mode fiber-optic components at telecom wavelength (1550 nm). The choice of fiber components presents several advantages. First, it completely removes issues related to matching the transverse spatial mode of the fields and considerably simplifies the alignment of sources, cavities, and detectors, thus allowing one to easily adjust the network size and topology. Second, working at telecom wavelengths guarantees low optical losses and a low cost of the fiber components. These are two other important factors for the scalability of the apparatus.
The FBG resonators have a finesse of about 100 and a Free Spectral Range of 1 GHz (corresponding to a cavity length of about 10 cm), which result in a full-width at half-maximum linewidth FWHM cav of about 10 MHz. The Light is injected in the network by a continuous-wave diode laser, transmitted over coupled optical cavities, and then irreversibly absorbed by a detector measuring the transmission rate. O.I., Optical Isolator; M1 and M2, mirrors; P.C., polarization controller; C, 50 × 50 fiber coupler; C1 and C2, 90 × 10 fiber couplers; FBG1 and FBG2, Fiber Bragg Grating resonators; PZT piezoelectric transducer; D1, D2 and D3, detectors. The PZT can be driven to introduce a difference between the resonance frequencies of FBG1 and FBG2, which can be either constant in time (V 0 ) or variable (V(t)), in order to insert disorder and/or dephasing into the network. Inset (a): scheme of the 4-site network mimicked by the optical setup. Inset (b): simplified scheme of our optical platform. reflectivity peak of the Bragg gratings is around 1550.13 nm, with a width of 0.7 nm. Each resonator is inserted in a home-made mounting to isolate it from environmental noise and allowing the piece of fiber containing the cavity to be stressed and relaxed in a controlled way by the contact with a piezoelectric transducer (PZT). In such a way the length of each cavity and, consequently, its resonance frequency, can be finely tuned. The input of each resonator is equipped with a fiber polarization controller to allow only one polarization mode to be confined inside the cavity. The light source is a continuous-wave, single-frequency, external-cavity laser (Thorlabs SFL1550S), emitting at 1550 nm, with a linewidth of about 50 kHz, narrower than the 10 MHz cavity linewidth. The laser source injects light of frequency ω S into one input port of the first (50:50) fiber coupler C. Light exiting the two cavities passes through two more polarization controllers before being coupled by a second (50:50) fiber coupler C. Finally, one of the interferometer outputs is measured by detector D3.
All the parameters characterizing the network and that are described in the following section are expressed in terms of the cavity detuning parameter Δ x. This parameter is defined as the difference (ω 2 − ω 1 ) between the resonance frequencies of the two cavities in units of their linewidth FWHM cav . Two additional fiber couplers (C1 and C2) are used to split a small portion (about 10%) of the light in each interferometer arm in order to measure the transmission peaks of each single cavity before interference. The transmission signals are measured by detectors D1 and D2 while scanning the laser frequency ω S over an interval including a single longitudinal mode of both cavities; from the difference between the frequency positions of the two peaks it is possible to infer Δ x. The three detectors D1, D2 and D3 are biased fiber optic InGaAs photodiodes.
Network parameters. The network configuration can be completely described by 3 parameters: the initial conditions related to global interference, the disorder (static disorder) and the dephasing (dynamical disorder). All these parameters are defined in terms of the previously described cavity detuning parameter Δ x. We theoretically and experimentally characterize the network transmission as a function of these parameters.
Initial Conditions for Interference. Our interferometric apparatus involves no active stabilization. Consequently, the system is intrinsically unstable and the network response will be time dependent (on the time scale of the order of hundreds of ms). This intrinsic instability can be used to establish different initial conditions of global interference for our network. The system throughput when the two cavities are resonant (Δ x = 0) and without any kind of noise will vary between a minimum and a maximum value in correspondence of global destructive or constructive interference. The output signal measured by detector D3 in this case will thus set the initial conditions of global interference of the network.
Disorder or Static Disorder. A network is said to be disordered if the local excitation energies of different sites are unequal (ω j ≠ ω k , with j ≠ k) but constant in time. In our setup, the disorder (or static disorder) of the network can be quantified with the cavity detuning parameter Δ x. The case of an ordered system is illustrated in Fig. 2a, where Δ x = 0 and constant. The case of a disordered system with a static disorder Δ x = Δ 0 ≠ 0 is shown in Fig. 2b. We assume that a system presents a medium level of disorder if Δ x ~ 1, i.e. (ω 2 − ω 1 ) ∼ FWHM cav , and a high level of Dephasing or Dynamical Disorder. Dephasing (or dynamical disorder) introduces a random phase perturbation in one or more sites of the network, thus resulting in temporal fluctuations of the corresponding resonance frequencies ω j 's around their stationary values. We can introduce it into our network by slightly changing the value of ω 2 during measurement by means of the PZT. In such a way, Δ x is not time constant but can vary within an interval of ± δx/2 around Δ x. The amount of dephasing can be quantified by the amplitude δx of this interval. The case of a network with disorder Δ 0 and a dephasing δ 0 is illustrated in Fig. 2c, where Δ x is variable in the interval Δ 0 − δ 0 /2 and Δ 0 + δ 0 /2. We assume that a system presents a medium level of dephasing if δx ~ 1, and a high level of dephasing if δx > 1.
Network Transmission. We define the transmission of the network as the output signal measured by detector D3 in correspondence of the transmission peak of cavity 1 (ω = ω 1 ) -see Methods for technical details on the acquisition procedure. The value of the transmission is normalized to the measured value of the output signal in conditions of constructive interference, without disorder and without dephasing (i.e., for Δ x = δx = 0).

Theoretical model and numerical results.
To describe the dynamics of our optical cavity network, we consider the following Hamiltonian 24 where â i and ˆ † a i are the usual bosonic field operators, that annihilate and create an excitation (a photon in our case) in the i-th site of the network, ω i the corresponding resonance frequency, and g ij are the coupling constants between all the connected sites. The first term in (1) describes the energy structure of the system, while the second one is related to the hopping process among the network nodes. In the following we will refer to the (random) energy level spacings of the network's sites as static disorder and it will be obtained by tuning the frequencies ω i of the network sites. In addition, the presence of dephasing noise, randomizing the photon phase during the dynamics, will be introduced in terms of dynamical disorder, i.e. time-dependent random variation of the site energies. We denote with γ i the dephasing rate for the site i. The latter will be mathematically described by a so-called Scientific RepoRts | 6:37791 | DOI: 10.1038/srep37791 , where ρ is the density matrix describing the system state that evolves according to the following differential Lindblad (Markovian) master equation: in Eq. (2) describe two distinct irreversible transfer processes, respectively, from the light source to the network (energy injection) and from the exit site to an external sink (energy detection) -see Methods for technical details. The transferred excitation energy, reaching the sink at time t, is defined as where ⋅ Tr[ ] is the trace operation, Γ det is the rate at which the photons reach irreversibly the output detector, and â k refers to the effective absorption of photons from the site k corresponding to the output site. To compare the theoretical results with the experimental data where the light is continuously injected into the network (with rate Γ 0 ) and absorbed from the detector, we need to define the network transmission as the steady-state rate for the photons in the detector, i.e.
. Let us point out that, if one repeats a single-photon experiment many times, one obtains the same statistics corresponding to an injected coherent state 25 , but this holds just because nonlinear processes are not present in our setup. Moreover, note that the model is able to take into account also the slight asymmetry that is present in the experimental setup by imposing different coupling rate in the two paths of the networks, i.e. g 01 ≠ g 02 and g 13 ≠ g 23 . Such asymmetry is mainly caused by different loss rates in the two resonators. In Fig. 3 different time behaviours of the transferred energy are shown for different initial conditions of global constructive (3a) and destructive (3b) interference, and different values of dephasing and static disorder (see Methods for details). Since the photon injection is continuous in time, the energy in the steady-state condition increases monotonically in time with an asymptotic linear behavior whose slope is indeed the transmission rate. An exponential behaviour is instead observed for the initial temporal regime, as show in the inset of Fig. 3. Then we find that, for constructive interference both disorder and dephasing individually reduce the transferred energy (no NAT). However, in presence of some disorder inhibiting the path to constructive interference, dephasing slightly assists the transport by opening additional pathways. This can be also intuitively explained by the fact that the two cavities are not energetically on resonance (because of the disorder) but the line broadening effect, induced by dephasing, allows again the hopping between them. On the other side, destructive interference leads to very small transferred energies since the two transmission paths over the two cavities cancel each other (opposite phase). In this case, dephasing and disorder can hinder such a perfect cancellation, thus partially restoring transport, i.e. NAT behaviour.
Experimental results. The network transmission has been investigated for different initial conditions of global constructive or destructive interference, as a function of both disorder and dephasing on the experimental 4-site network of fiber-optic resonators in Fig. 1. In particular, the noise effects on the amount of transferred energy will be shown from two points of view: we analyze the network transmission as a function of dephasing for different values of the disorder, and then as a function of disorder for different values of the dephasing. Notice that our experimental fiber-optic setup has been designed as a simple, scalable and low-cost platform with not perfectly identical FBG resonators and small but significant losses within the network. Let us stress that, although it is challenging to estimate the model parameters to reproduce the experimental imperfections, as site energies, coupling rates, relative phases between the two arms, and asymmetric cavity losses, we have however achieved a sufficiently high level of control by tuning the global interference, static and dynamical disorder. Besides, in the model the static disorder is added by tuning the cavity frequency of site 2, ω 2 , while dephasing is given by γ 2 . For both the experimental data and numerical results, the network transmission is normalized to the value of the output signal in the condition of constructive interference without disorder nor dephasing.
Constructive Interference. We start by investigating the behavior of the network transmission for initial conditions of global constructive interference. In Fig. 4 the experimental and theoretical network transmission are shown as a function of dephasing for three different configurations of disorder. Without disorder (Δ x = 0), the only effect of dephasing is to reduce the transferred energy, i.e. no NAT is observed because the different pathways do already constructively interfere. However, if the system's energy landscape presents some disorder (Δx > 0.4), NAT effects can be detected and, in particular conditions of disorder, one finds the typical bell-shaped NAT behavior with an optimal value of dephasing that maximizes the network transmission. In other terms, dephasing enhances the transport efficiency if the additional presence of disorder inhibits the otherwise fast constructive-interference path. A similar behavior has been found for the theoretical model, although a discrepancy with the experimental data for small values of dephasing is obtained. This difference is mainly due to the fact that the experimental realization of constructive interference is not perfect and the initial negative impact of dephasing is reduced. A similar reasoning will explain in the following the discrepancy for the destructive interference case.
In Fig. 5 the transmission is shown as a function of disorder for two different dephasing configurations. Here, as expected, we find that disorder has a generally negative impact on the transport performance since it always leads to the suppression of the initial constructive interference. However, while a little bit of disorder quickly deteriorates transport in the case without dephasing, the presence of some dephasing noise that broadens the  resonances has the effect of making the system more robust against static disorder, with an evidently smoother decay in both the theoretical and experimental cases.
Destructive Interference. Repeating the analysis above for the case of initial destructive interference, it turns out that both dephasing and disorder independently assist transport, i.e. NAT behavior, since they reduce the amount of interference that prevents the transmission of energy -see Fig. 6. When the network is in a regime of high-disorder, again the typical bell-like NAT shape is recovered. However, if a lower amount of disorder is included, the additional contribution of dephasing does not further improve the network transmission that instead shows a minimum value for a dephasing value of δx ≈ 2.2. Finally, in Fig. 7 we show the role of disorder in enhancing the transmission rate with its peak moving to higher value of disorder for increasing dephasing values.

Discussion
Noise-assisted transport phenomena occur in several physical systems where noise can open additional transport pathways and suppress the ineffective slow ones. In the last years, this scheme has been applied to better understand energy transport in photosynthetic light-harvesting complexes where dephasing noise remarkably enhances the transmission of an electronic excitation from the antenna complex to the reaction center where such energy is further processed 5 . The basic underlying mechanisms of such behavior are mainly due to line-broadening effects and the suppression of destructive interference, hence the interplay of quantum coherence and noise is responsible for the observed very high transport efficiency. While several theoretical studies have been performed, it is very challenging to actually test these ideas either in the real photosynthetic pigment-protein complex or in artificial ones, since their structure and dynamical properties cannot be controlled or even indirectly measured 26,27 with the required resolution. Moreover, these samples are usually quite expensive or difficult to synthesize. For these reasons, it is very convenient to reproduce such transport phenomena in a controlled system where one can tune the parameters and measure the corresponding dynamical behaviors, while also playing with the underlying network geometry. This will allow one to better understand the underlying problems and to start engineering new molecular/nano-structures for more efficient and feasible technological applications.
In this work, we have experimentally investigated how to design optimized accessible noise features for energy transport in a simple and scalable test optical platform, represented by a network of fiber-optic resonators. Although the experiment could be fully described by classical optics, the role of the quantum coherence in energy transport processes, as observed for photosynthetic systems, can be mimicked by the propagation and interference of light in our experimental setup in correspondence of different global interference conditions. In the proposed interferometric apparatus, the effective controllability of parameters is reached by exploiting the intrinsic instability of the setup, where no phase locking is used for the stabilization of the system. Therefore, our setup not only allows one to observe the NAT peak in the network transmission as a function of the amount of noise in the system dynamics, but, more importantly, can be exploited to monitor and control the different noise sources to show their individual contribution to the transport dynamics in presence of different types of global interference. Then, we find that, when constructive interference already provides a very fast path to the exit site, dephasing has a detrimental effect and reduces the amount of transferred energy. However, if some disorder is present in the system energy levels, hence blocking the constructive interference path, then dephasing represents a recovery tool to achieve again higher transport efficiency, thanks to NAT effects. On the other side, in presence of destructive interference, both dephasing and disorder are able to speed up the energy transport, with dephasing often providing a faster NAT mechanism. Intuitively, the role of noise in increasing the transferred energy can be explained by considering how the pathways of energy transfer are modified, by destroying the inefficient ones (i.e. inhibition of destructive interferences in the network) or by giving access to more efficient network hubs.
Let us further stress the particular usefulness of the fine controllability of the system parameters and topology in such a simple experimental model compared to the control 28,29 of the molecular geometries in natural and artificial light-harvesting molecules. This is very hard or even unfeasible in presence of just a few chromophores, with the problem becoming even harder when many chromophores are taken into account. Furthermore, by working with real molecules it is not possible to tune at will their parameters, such as the exciton energies, couplings, noise strength, etc., since the corresponding structure has to be synthesized again each time and just a few molecular configurations are usually available. Our scalable and tunable optical model does instead provide a very useful benchmark platform where all system parameters and geometrical configurations can be very easily changed and individually controlled, even including tens or more of sites. The power of such optical models is also demonstrated by the fact that they have allowed us to observe, for the first time to the best of our knowledge, the only-theoretically predicted bell-like behavior of the network transmission as a function of the dephasing noise strength 18 . Then, increasing the network size, these results are promising in order to observe even more complex transport behaviors that can be very hardly simulated on a computer (see Methods). In particular, optimizing a network for energy transport will correspond to determine the optimal values of the parameters referring to accessible sources of static and dynamical disorder, given the topological features of the network. Notice that in this context the network topology or connectivity plays also a crucial role that can be further investigated by these scalable optical setups.
Therefore, this platform gives us the crucial possibility of testing new complex structures and transport regimes that would not be accessible by studying single molecular systems, but that can be triggered by the observations on these optical models and then synthesized later on artificial molecular devices 30 . Indeed, our theoretical and experimental results show that disorder and dephasing can be exploited as control knobs to manipulate the transport performance of complex networks, e.g. optimizing the final transmission rate or the capability of temporarily storing energy/information in the system. In particular, noise-enabled information transfer phenomena 10,31,32 exhibit optimal features when both static and dynamical disorder are tuned to specific values, so that the information transfer within complex networks can be considerably enhanced. For these reasons, these achievements are expected to find promising bio-inspired applications for practical and much more efficient light-harvesting and energy/information 33,34 transport nano-technologies. Methods Experimental data acquisition methods. The network transmission is measured by monitoring the output signal acquired by detector D3 while the laser frequency ω S is scanned over an interval comprising a single mode of both cavities. The laser frequency is tuned by applying a voltage ramp, with a frequency of 30 kHz, to the driving current by means of a modulation input port (tuning coefficient about 90 MHz/mA). In order to scan over the cavity modes a variation of less than 1 mA is needed, which yields a residual amplitude modulation smaller than 1% of a typical transmission peak amplitude. We define the network transmission as the D3 signal for ω S = ω 1 , where ω 1 can be obtained by the peak of the signal measured by D1. A control software drives the data acquisition chain according to the following steps: 1. ω 2 is initially adjusted by driving the corresponding FBG2 PZT voltage in order to bring the transmission peaks of the two cavities within a few linewidths from each other, thus reducing their detuning Δ x. We monitor the D1 and D2 signals until Δ x < 0.1 (i.e., until the distance of their transmission peaks is smaller than one tenth of their linewidths) and use this condition to trigger the rest of the acquisition; 2. the value of the D3 signal at this time is recorded and it defines the initial conditions of global constructive or destructive interference for the network transmission; 3. from this condition, a 1 kHz sinusoidal modulation is applied to the FBG2 PZT so that ω 2 is varied in a range of about ± 5 linewidths around its initial value. The three signals from D1, D2 and D3 are then acquired continuously for a time of 1 ms so that, for a 30 kHz laser-driving ramp, the signal from D3 allows us to retrieve 30 different values of the network transmission. Since an acquisition time of 1 ms is so short that the entire setup is sufficiently stable to maintain the initial conditions of global interference unchanged, all these acquisitions correspond to the same initial global conditions of interference but to different values of the static disorder parameter Δ x; 4. after a time acquisition of 1 ms, the acquisition is stopped and the procedure is repeated several times.
In such a way, for each initial condition of constructive or destructive interference, a statistically significant data-set of network transmissions (about 14000) has been obtained for different values of Δ x.
From this data-set, the transmission for a network with static disorder Δ 0 is simply obtained by averaging all the network responses with the same Δ x = Δ 0 . An amount δx of dephasing (dynamical disorder) is introduced by averaging the network responses corresponding to different values of Δ x included in an interval ± δx/2 around Δ x. In particular, the transmission of a network with static disorder Δ 0 and dephasing δ 0 , is calculated by averaging all the network transmissions with Δ x included in the interval Δ 0 ± δ 0 /2. Statistical uncertainties resulting from this procedure are reported as error bars in all the data points for the experimental network transmissions of Figs 4, 5, 6, and 7.
Scientific RepoRts | 6:37791 | DOI: 10.1038/srep37791 Computational complexity. Here we evaluate the computational time for the single realization of the (numerical) dynamical evolution of networks of increasing size (i.e., number of sites/cavities N) with the model as in the Main. The computational complexity does increase exponentially with N, as shown in Fig. 8. In other words, already adding a few cavities to our model would take months to theoretically simulate the corresponding dynamics for a given set of parameters, and the computation becomes unfeasible if one wants to reconstruct the figures in our manuscript where thousands of simulations are needed to take into account dynamical disorder, dephasing, etc. In particular, the critical number of sites above which it becomes very hard to reproduce the theoretical data is around 8 (corresponding to around six months of simulations for one thousand realizations) -see Fig. 8. Let us notice that, while the experimental scheme has been realized with coherent states of light that in principle allow its classical simulation time to scale polynomially with the number of optical elements, if one instead considered a full quantum regime (for example, with several single-photon walkers), the computation complexity would have indeed scaled exponentially, as observed above. On the other side, the experimental complexity is not so affected by the network size and, at most, linearly increases in terms of both the cost of the optical components (cavities, beam-splitters, etc.) and the practical realization and observation time of the stationary behavior of the optical system. The use of single-mode telecom optical fiber components is an essential ingredient in this context, as it dramatically reduces costs and losses, besides canceling all issues related to the alignment and spatial mode-matching procedures that should be faced when adding other cavities to the system or changing its topology. However, even if the increase of the network size is not a significant limit for the present setup employing a classical coherent source, the same would not be true when operating in a quantum regime, with one or multiple single-photon sources and coincidence detection. In such a case, the increased losses in a larger network (mainly due to the presence of several FBG resonators) could substantially reduce the efficiency of the setup and imply much longer acquisition times. Computational time (in hours) as a function of the simulated network size (i.e. N number of cavities) for a single realization of the dynamical evolution described in the Main, with only one choice of the system parameters. An exponential behaviour is observed (green dashed line), hence a scheme with more than ten sites/cavities (red dashed region in the inset) becomes very hard to be simulated by a powerful workstation, corresponding to at least one month of simulation. Our experimental setup corresponds to the case of four cavities (blue dashed region), but can be easily extended to more cavities.