Entanglement formation in continuous-variable random quantum networks

Entanglement is not only important for understanding the fundamental properties of many-body systems, but also the crucial resource enabling quantum advantages in practical information processing tasks. While previous works on entanglement formation and networking focus on discrete-variable systems, light---as the only travelling carrier of quantum information in a network---is bosonic and thus requires a continuous-variable description in general. In this work, we extend the study to continuous-variable quantum networks. By mapping the ensemble-averaged entanglement dynamics on an arbitrary network to a random-walk process on a graph, we are able to exactly solve the entanglement dynamics and reveal unique phenomena. We identify squeezing as the source of entanglement generation, which triggers a diffusive spread of entanglement with a parabolic light cone. The entanglement distribution is directly connected to the probability distribution of the random walk, while the scrambling time is determined by the mixing time of the random walk. The dynamics of bipartite entanglement is determined by the boundary of the bipartition; An operational witness of multipartite entanglement, based on advantages in sensing tasks, is introduced to characterize the multipartite entanglement growth. A surprising linear superposition law in the entanglement growth is predicted by the theory and numerically verified, when the squeezers are sparse in space-time, despite the nonlinear nature of the entanglement dynamics. We also give exact solution to the equilibrium entanglement distribution (Page curves), including its fluctuations, and found various shapes dependent on the average squeezing density and strength.

The above works, whether on the basic understanding of scrambling or practical design of networking, mainly focus on entanglement in discrete-variable (DV) systems, which is natural for computing. However, as quantum networks inevitably utilize light as the carrier of quantum information in transmission, the bosonic nature of light makes it necessary to consider entanglement in a continuous-variable (CV) description. Moreover, various applications in the photonic or microwave domain, in- * zhuangquntao@email.arizona.edu cluding universal quantum computing based on cluster states [46], quantum illumination [47][48][49], quantum reading [50], distributed sensing [51][52][53][54] and entanglementassisted communication [55,56], require CV entanglement in the form of Gaussian states [57]. In this regard, noiseless linear amplifiers [58] and novel error correction codes [59,60] provide initial tools for CV networking, and an out-of-time-order correlator (OTOC) has revealed a unique squeezing-dependent butterfly-velocity of operator spreading [61].
In this paper, we study quantum information scrambling in CV quantum networks (see Fig. 1) focusing on the entanglement formation dynamics. Inspired by the classical statistical theory of complex networks [62][63][64], we consider random quantum networking protocols to enable analytical solutions, through a mapping to a random-walk process on graph; at the same time, random protocols are expected to reveal typical and universal characteristics. Our results apply to quantum networks on general graphs, therefore provide a foundation for the statistical theory of complex quantum networks.
We provide an analytical formula connecting the entanglement entropy to weights in the passive linear optical transforms, therefore establishing a mapping between ensemble-averaged entanglement dynamics to the probability evolution of a random-walk process on a general graph (see Fig. 1). The change of the entanglement entropy S(L, t) of a subsystem L, similar to the random walker's probability in subsystem L, is determined by the boundary ∂L. Moreover, we also solve the fluctuations of the entanglement entropy.
In discrete space-time, the ensemble-averaged weights dynamics can be described by a Markov chain, with the transition matrix determined by the graph connectivity. Going into continuous space-time, we can derive simple diffusive partial differential equations (PDEs) to describe the entanglement evolution. In both cases, the model is completed with the analytical formula connecting weights to entanglement entropy. Alternatively, a phenomenological model coupling an epidemiology equation with a nonlinear diffusion can directly capture the entanglement dynamics.
To go beyond the bipartite characterization of entanglement through entanglement entropy, we also give an operational witness of multipartite entanglement. This witness is directly connected to entanglement's advantage over classical correlations in distributed sensing protocols [60]. Maximum values of the entanglement witness are achieved towards the late time, therefore verifying the full scrambling of the entire network.
Through the mapping between quantum dynamics and random walk, we also connect the scrambling timethe time it takes for the entire system to be maximally entangled-directly to the mixing time of the random walk. Moreover, at infinite time, the equilibrium entanglement distribution-analog to the Page curve in DV systems [65][66][67]-can be solved analytically from the stationary state of the random walk. Surprisingly, the Page curve is independent of the topology of the network, as long as the graph is connected. And in general it depends on two statistical properties of the quantum network-the squeezer's density and the average squeezing strength. Interestingly, a small subsystem can almost get close to the maximum entanglement entropy, while in DV systems, half of the system size is necessary. While our theory works for general graphs, we also apply to networks respecting 'locality' of interactions-D-dimensional Cartesian graphs where links only exist between nearest neighbors (see Figs. 2 and 4). In this regard, we identify a diffusive entanglement light cone at the early time, which divides the regions with almost no entanglement and regions with substantial entanglement. After the entanglement light cone reaches each node, there is a period of entanglement sudden growth, where the entanglement entropy quickly gets close to its equilibrium value. In the end, there is a long period of saturation, determined by the mixing time ∼ M 2 quadratic in the length M on each dimension.
Our theory framework provides a unique complement to the DV counterparts [9][10][11][12][13]20]. And our results provide insights into not only CV quantum networks being engineered, but also quantum information scrambling in various physical systems, as any form of bosonic radiation is intrinsically CV.
The paper is organized as the following. In Sec. II, we specify the model in details and give a more specific overview of results; in Sec. III, we present the statistical theory for the single-squeezer case based on a mapping to random-walk dynamics; in Sec. IV, we generalize the single-squeezer results to the general case through providing a linear superposition law.

II. QUANTUM NETWORKS: MODELLING AND MAIN RESULTS
Our overall goal is to characterize generic entanglement formation dynamics towards equilibrium in a CV quantum network (see Fig. 1 for a schematic). In general, a quantum network can have complicated topology, which makes the problem difficult. Moreover, each node can possess multiple optical modes, and perform local operations coordinated by classical communication to entangle them. Considering the optical modes, we can reduce a general entanglement generation protocol to a quantum circuit on a graph, as we illustrate in an one-dimensional (1-D) hopping quantum network in Fig. 2 (a). To establish entanglement, each node performs the following protocol repetitively: it receives a light mode from a neigh- case is plotted). We index the M = 2N + 1 bosonic modes by integers x ∈ [−N, N ]. For convenience, we also introduce a scaled coordinatex = x/2N ∈ [−0.5, 0.5]. Here we exemplify the notation: the connection matrix E x,x = δ |x−x |−1 ; neighbors N (x) = {x − 1, x, x + 1}; it is convenient to consider the left part of the system L = [−N, x], and the right part of the system R = [x + 1, N ]. In (b) the cyan square is a singlemode squeezer, while the orange rectangles are the 2-mode beamsplitters (combing phase shifters). The red dashed line represents the light cone from the center. In (a) each empty circle denotes a local mixing operation by beamsplitters and the dashed arrows indicate the transmission of optical modes in the network. The cyan lines indicate the light cone starting from the vertex with the single-mode squeezer.
bor, which gets entangled with a stored mode through a local unitary; then, it sends out a mode to another neighbor and stores one mode locally. For simplicity, the nodes send light to the left and right neighbors alternatively in even and odd steps. If we focus on the dynamics of the optical modes, the above protocol reduces to a 1-D local circuit in Fig. 2 (b), where local gates apply alternatively on the light modes [61].
The transmission links in a quantum network are in general lossy. To cope of loss, error correction [59] can be applied in each link transmission. On the physical layer, this means including additional components that seemingly complicate the analyses. However, on the logical layer, up to some small residual errors from imperfect error correction, the state being protected is identical to the state being generated in a lossless quantum network, as demonstrated in Ref. [60] for sensing purposes. Therefore, we start with the lossless case.
With the mapping between quantum networks and quantum circuits in mind, we specify the set-up of the circuits on an arbitrary (undirected) graph (see Fig. 1(c)). In general, the topology can be described by an undirected graph (G, E), where G denotes the set of all vertices, each described by a coordinate system x [68]. The set of edges E can be described by a generalized connection matrix E x,x . When E x,x = 1, the vertices x, x are connected by an edge xx , zero when not connected. For simplicity, we write the set of vertices that are directly connected to x (neighbors) as N (x). We are interested in the entanglement between a set of vertices L and the rest R = G\L. In Fig. 2(b), we give a 1-D example of the notations.  Figure 3. Heat map of ensemble-averaged entanglement entropy evolution showing the entanglement light cone. x-axis is the re-scaled spatial coordinatex and y-axis is time t. The color indicates the entanglement entropy between the left and right part of the system, divided by the boundary at x. Number of modes M = 201. (a) Single squeezer r = 8 is placed at the center of the system (shown as a cyan box). The green dashed line represents the entanglement light cone (t = T1) and the black dashed line represents the characteristic time-scale T1 +T2 of entanglement-growth. See Sec. III D 2 for details. (b) Multiple squeezers placed at different spacetime (shown as cyan boxes). The first squeezer r1 = 5 is placed atx = 0, t = 0, the second one r2 = 3 is placed at x = −0.35, t = 200 and the third one r3 = 7 is placed at Unitaries are applied on the edges E. We separate the edges into disjoint sets {E k } K k=1 , such that the edges in each set E k do not have common vertices. The dynamics repeat in a period of K steps; in the k-th step of each period, one applies unitaries U t,x,x on each edge xx ∈ E k . The particular separation of the unitaries is not essential to the dynamics and equilibrium. As an example, in a 1-D local circuit, K = 2 and we alternative between gates {U t,k,k+1 } on k odd and even; in a 2-D local circuit, we have K = 4, as shown in Fig. 4.
To produce the Gaussian states that enables various applications in communication, sensing and computing, we consider Gaussian unitaries [57], which are unitaries generated by Hamiltonians that are second order in the quadrature operators (see Appendix A). Gaussian unitaries include squeezing, which creates asymmetry in quadrature noises; and passive linear optics, which includes beamsplitters and phase-shifters.
Squeezing is essential for entanglement generation. However, as an 'active' component, squeezing is relatively difficult to implement. Thus, we consider the gates {U t,x,x } to be passive linear-optics gates. And squeezing operations are added in between in a sparse way. As an example, in Fig. 2(b), to establish entanglement, in this case a single vertex performs a squeezing operation (the cyan box), and then entanglement is generated by passing it around through passive components (the orange boxes), with vacuum on the other input modes.
We expect random protocols to reveal universal characteristics, therefore we choose the passive linear optics gates {U t,x,x } to be Haar random (see Appendix A). This is also justified by the following reasons: (1) In classical complex network theory [62][63][64], various networks can be modeled as random networks with a proper degree distribution. (2) In condensed matter theory, random quantum circuits and Hamiltonian systems [9][10][11][12][13]20] are able to capture the essential quantum information spreading features in generic many-body interacting systems. (3) In real quantum networks, the form of entanglement required can be complicated, depending on the purpose, e.g. the weights of the global parameter of interest in a distributed sensing protocol [51].
We aim to characterize the entanglement dynamics in the above random circuits. The entanglement entropy, measured by von Neumann entropy or Renyi entropy, can be numerically evaluated efficiently (see Appendix B). Two examples of time evolution of von Neumann entropy in 1-D are given in Fig. 3. In Fig. 3(a), we have a single squeezer at the center in the first step, which is identical to the case depicted in Fig. 2. The entanglement entropy grows diffusively from the source of squeezing (see Sec. III D). Different from the DV case, we can identify an entanglement light cone (green lines), which is the boundary between regions with substantial entanglement and regions with almost-zero entanglement. These phenomena can also be found in higher dimensional random local circuits, as shown in Fig. 4 for the 2-D case. When choosing a subsystem L as an individual mode at (x 1 , x 2 ), we can see similar entanglement light cone (green dashed).
As shown in Fig. 1(b), the above entanglement dynamics can be solved by mapping to a random walk on a graph, which gives the exact ensemble-averaged entanglement entropy (see Sec. III A) S(L, t) as a function of the total probability η L,t of having the walker in region L in the corresponding random walk, is the (von Neumann or Renyi) entropy of a thermal state with mean photon number x (see Appendix B) and r is the original squeezing strength. Note that the mapping holds for arbitrary graphs beyond the local Cartesian graphs shown above. The scrambling time-the time that the entire system becomes maximally entangled-can be calculated by the mixing time of the random walk (see Sec. III C). The mapping also gives the Page curves-the late-time equilibrium entanglement entropy S(L, ∞) = S(|L|/|G|) (2) as the CV analog to Page curve (see Sec. III B), while the fluctuation can be solved as ∝ |L||R|. Moreover, when there are multiple squeezers, we can regard the entanglement dynamics as the superposition of all singlesqueezer dynamics, as depicted in Fig. 3(b) and will be detailed in Sec. IV C. When there are multiple squeezers, we find that the Page curve is determined by the average squeezing strength and density of the squeezers (see Sec. IV B). Therefore, combining the results, we have a complete understanding of the entanglement dynamics in a CV quantum network.

III. STATISTICAL THEORY OF RANDOM CV QUANTUM NETWORKS
In this section, we present a statistical theory of the entanglement growth. We will focus on the single squeezer case in Fig. 3(a), while the extension to multiple squeezers is presented in Sec. IV. We introduce the mapping between random unitary circuits and random walk on a graph in Sec. III A, which allows us to solve the Page curves in Sec. III B and scrambling time in Sec. III C for general graphs. Explicit closed-form solutions can be obtained for local Cartesian graphs of an arbitrary dimension in Sec. III D. Finally, we present the entanglement witness for multipartite entanglement in Sec. III E.
A. Mapping to random walk on graphs Consider the entire unitary evolution U (t) of the random circuit. In the single-squeezer case, the mode annihilation operator a x,t at vertex x ∈ G experiences a passive transform, which in general can be expresses as where mode a SV is in a squeezed-vacuum (SV) state with strength r and 'vac' denotes all vacuum terms that complete the commutation relation. Here the phase θ x,t is entirely random, and the positive weights w t ≡ {w x,t } x∈G describe the overall energy splitting of the single SV among all modes. For any subsystem L with a density operator ρ L (t), we can design a passive linear optics circuit U L,t such that U L,t ρ L (t)U † L,t concentrates all the squeezing parts to a single mode with the total transmissivity and all other modes are in vacuum required by energy conservation. Because unitary operations preserve entropy, the entanglement entropy of L can be calculated from the entropy of mode a L,t as where S(η L,t ) is defined in Eq. (1). We will focus on von Neumann entropy, but all of our results can be adapted to Renyi entropy easily. At the large squeezing limit of η L,t (1 − η L,t )e r 1, for von Neumann entropy we have Eq. (7) to the leading order. When η L,t = 0, 1, subsystem L has zero or entire portion of the SV, indeed from Eq. (6) we have S (L, t) = 0. When η L,t = 1/2, we have the maximum entropy S 0 (r) = g(sinh 2 (r/2)) log 2 (e r /4). When one has large number of modes, this should agree with the result of max L S (L, t) at any time.
So far we have the exact result S (η L,t ) of the entanglement entropy of an arbitrary subsystem L, given the weights w t (which determines η L,t ) obtained in each random circuit realization. Due to the self-averaging in the random circuit, we expect up to corrections that decay with the system size. Thus, we have reduced the problem of solving the ensemble-averaged entanglement dynamics to solving the ensemble-averaged dynamics of the weights w t . We start by focusing on a single gate U t,x,x on the modes at x and x . By considering the Haar random ensemble averaging, we can derive the exact equation of motion of the weights as (see Appendix C) The overall dynamics alternatives in K steps, in the k-th step the transition of Eq. (9) on all edges in E k is applied. An immediate observation from Eq. 9 is that the change of the entanglement entropy of L, determined by the total weights η L,t , is related only to boundary ∂L (schematic in Fig. 1(c)), in the sense that which equals the net flow of the weights from the vertices on the outer boundary ∂L + towards L and the weights from the inner boundary ∂L − out from L. Another observation is that the weights update rule in Eq. (9) also describes the probability evolution of a lazy symmetric random-walk step, where the walker have half probability of staying and half probability of taking a step along xx (see Fig. 5). Combining the K steps, the underlying transition matrix for the weights where I is the identity matrix and E k,x,x describes the adjacency matrix for the corresponding graph (G, E k ) (an isolated mode is regarded as a vertex with a loop). Eq. 11 describes a modified symmetric random walker on the graph (see Fig. 5), with K steps combined to implement a single random-walk step from the current position x to all neighbors N (x) (including x) with equal probability. Utilizing the K-step transition matrix, the ensembleaveraged weights can be solved at any time t as with the initial condition w 0 = δ x0 as the Kronecker delta at the squeezer position x 0 . Thus, one can obtain η L,t and the exact result of S( η L,t ) from Eq. (6) on any graph.
We give examples of the random walk in Fig. 5 in 1-D and 2-D Cartesian graphs, whose entanglement evolution can be found in Figs. 2 and 4. For the later use, we also introduce a general D dimensional Cartesian lattice G D , with the coordinates x = (x 1 , · · · , x D ) on a grid (x d ∈ [−N, N ]). The total number of modes is |G D | = M D , with M = 2N + 1 modes on each dimension.
In the following, we will consider the equilibrium and the dynamics. Some results hold for general graphs, while  In this section, we focus on the Page curves-the equilibrium entanglement distribution at infinite time. In order to share entanglement, squeezers are applied, which are then followed up by the random beamsplitters and phase shifters. As the layers of gates increases, the overall passive linear transform will approach the Haar measure (see Appendix A). Therefore, we can regard the equilibrium entanglement distribution as the CV analog to Page curves.
Considering the mapping from the circuit to the random walk, the equilibration of the entanglement also corresponds to the full mixing of the random walk on the graph. Assuming the full connectivity of the graph, due to the special transform matrix in Eq. (11), the equilibrium (stationary) state of weights is uniform among all vertices, i.e., where |G| is the total number of vertices, despite how one arranges the set of edges E k . Note that this is different from normal random walks on a graph, where the stationary state has weights proportional to the degree of the vertex |N (x)| [70]. Therefore, the total transmissivity (i.e., total weights) In fact, assuming fully random weights from a Haar random unitary, one can obtain the probability density of total weights as (see Appendix D) From Eq. (6), the Page curve is therefore given by The maximum is achieved at |L|/|G| = 1/2, which equals S 0 (r) introduced following Eq. (6). The second equality is the leading order result similar to Eq. (7).
Note that in terms of Page curves, the graph topology is irrelevant as the entire dynamics is equivalent to a single passive global Haar unitary; therefore, we can simply stretch the coordinates x of a general graph G to a single coordinate x in 1-D. It then suffices to verify our theory of Page curves in G 1 , which allows simple visualization. In the 1-D system, when subsystem L contains the left side of mode x, we have η L,∞ = |L|/|G| = (x + N )/(2N + 1) x + 1/2. Therefore, it is convenient to choose the parameterizationx. Fig. 6(a) plots S(x, ∞) for various squeezing values of r and system sizes of M , where we see perfect overlapping among curves with identical r for different system sizes. And they all agree with Eq. (16) very well, as shown in Fig. 6 (b). As a by-product, the maximum entanglement-the maximum height of the Page curve max x S(x, ∞) -agrees with the theory prediction S 0 (r) following Eq. (6), as shown in Fig. 6 Furthermore, we can also consider the fluctuations around the Page curves, as exemplified in Fig. 7 (a) and   21), (22) and (24). Orange scatter dots show the numerical mixing time for corresponding random walk. is set to be 1×10 −7 . Note: In 2-D grid system, we estimate the conductance without full optimization.
(c). Since the entanglement entropy is determined by η L,t = x∈L w x,t , it takes some time for {w x,t } x∈L to entirely change its values; thus, there will be correlations in the entanglement entropy at different times. We consider the equilibrium auto-correlation We expect the decay of this auto-correlation should have a similar time-scale with the mixing time of the entire system (which will be detailed in Sec. III C), as demonstrated in Fig. 7. When ∆t = 0, the auto-correlation goes to the variance where we have used the chain rule of variance and the variance of η L,∞ can be obtained in Appendix D. Numerical results in 1-D and 2-D Cartesian graphs G D agree well with Eq. (20), as shown in Fig. 8.
On the other hand, if we look at the change of subsystem entropy within a short period of time, e.g., a single-step, similar to Eq. (10) the short-time fluctuating (S (L, t + 1) − S (L, t)) 2 will mainly come from the boundary. This interplay of short time fluctuation related to the boundary, while long time fluctuation related to the bulk manifests the rich entanglement dynamics in random quantum networks.

C. Entanglement scrambling time
The entanglement scrambling time-the time for the system to be maximally entangled and reach the equilib- rium is an important quantity for many physical problems, especially those related to the black hole [28,33]. In the CV quantum network, this scrambling time can be directly obtained from the mixing time t of random walks on a graph, where the probability measure (weights) get -close to the stationary state in Eq. (13). Formally, we define the mixing time t to be the time when the deviation | w x,t − 1/|G|| ≤ for all x ∈ G. Multiple estimates can be obtained, as we explain bellow [70]. The first estimate relies on the eigenvalues of the transition matrix E x,x in Eq. (11). The stationary state w x,∞ = 1/|G| corresponds to the largest eigenvalue of unity, and the second largest eigenvalue λ < 1 gives the decay of deviations From the above, one can obtain t ∼ K ln (1/ ) / ln(1/λ ).
Here we have taken into account that, one needs K steps to implement the transition in Eq. (11).
For the Cartesian graphs G D , we can obtain the convergence towards w x,∞ = 1/|G| from the expected hitting time in the large t limit as [70] | where M = |G| 1/D is the length of the Cartesian graph in each dimension. This gives another estimate t ∼ ln (1/ ) DM 2 / ln (6).
A third bound can be obtained by calculating the conductance of a graph as where |∇L| denotes the set of boundary edges connecting L to R and |E| is the total number of edges. The minimization will be able to capture the slowest part of the mixing. We can obtain from Ref. [70] as We compare the three estimates of mixing time on G D for D = 1, 2 in Fig. 9, and find a good agreement in the scaling of the scrambling time as t ∝ DM 2 .

D. Closed form solutions to Cartesian graphs
With the understanding of the equilibrium Page curves, we now proceed to characterize the dynamical evolution towards the equilibrium. We will focus on the Cartesian graphs G D , which allows closed-from solutions by analog to random walkers. In the continuum limit, the dynamics can be well-described by a diffusive PDE (Sec. III D 1). In particular, we identify a unique parabolic entanglement light cone, followed by an entanglement sudden growth phenomenon, in CV networks (Sec. III D 2), as has already been shown in Figs. 3 and 4. In Cartesian graphs, the random walk analog to Eq. (9) can be understood as independent along each dimension. Thus, we modify the Pascal's triangle from a usual random walk to obtain the solution with ] and a b as the binomial factor of a-choose-b.
With the weights in hand, we can calculate the entanglement entropy of an arbitrary subsystem L. For example, we can consider L = {x |x d < x d , 1 ≤ d ≤ D}, i.e., the system is cut into two parts by a high-dimensional plane. Then we can obtain the ensemble averaged total transmissivity where F is the hypergeometric function. We expect this to hold up to rounding errors from the integers. In the following, we compare the exact solution of weights w t from numerically solving Eq. (12) and the binomial solution of weights w t (Bi) in Eq. (25). As shown in Fig. 10, we see a very good agreement in the 1-D case, up to rounding errors before boundary effects comes in. When there is a finite boundary, standard techniques like image source methods can give more precise solutions. In fact, the continuum limit of Gaussian solutions, as we will present in Sec. III D 1, also agrees well with the above results. The perfect agreement of the weights directly indicates the validity of the solution for the entanglement entropy, as demonstrated in Fig. 11. Note that due to the symmetry among different dimensions, it suffices to consider the 1-D case; however, results in higher dimension reveals more interesting dynamics.  Fig. 14(a1)-(a4). We see a gradual saturation to the equilibrium, where entanglement entropy along (x 1 + N )(x 2 + N ) = constant are about equal, as the total weights η L,t are equal along this line. Second, we can choose L = [−|x 1 |, |x 1 |] × [−|x 2 |, |x 2 |] as squares centered at the origin. In Figs. 14(b1)-(b4) We can also see gradual saturation to equilibrium, where the entanglement entropy along |x 1 ||x 2 | = constant are equal, as the total weights η L,t are equal along this line. To enable comparison in 2-D, we consider cross sections with x 1 = x 2 . As shown in Fig. 12, in both choices of the region L, good agreement between the circuit results and the random-walk results can be seen.
In the above, we see the entanglement entropy from exact ensemble-averaged evolution of Eq. (12) (combined with Eq. (6)) and the actual results from numerical solving the entropy agree well, therefore verifying the underlying random-walk model. Below, we further address the continuum limit and the entanglement light cone.

Continuum limit
We can take the continuum limit of Eq.
In Fig. 10, we see a good agreement of the above Gaussian approximation with the other solutions. The above solutions can be written as a function of x = x/M andt = t/M 2 up to normalization, which is well-defined in the continuum limit of M → ∞. We verify the continuum limit in 1-D by calculating the deviation measured by the relative 1-norm between the entanglement entropy S(x, t) and the static value S(x, ∞) as in the dynamic (short as 'dy') process (see Fig. 13 (a)), where f (x) 1 = x |f (x)| sums over the spatial coordinates. We see the re-scaled curves overlap well for systems with different sizes. We can also directly verify in Fig. 13 (b)(c) that the entanglement entropy agrees well after re-scaling. It is also worthy to point out that the continuum limit identified for the above single-squeezer case also holds for the multiple-squeezer cases (see Appendix E). The continuum limit in Eq. (27) naturally brings a PDE that descrbes the dynamics The above holds for Cartesian graphs G D , in general one needs to adopt the '∇ 2 ' operator to a graph. Nevertheless, we can obtain observations from G D . First, for any region L, as the continuum limit of Eq. (5), we have η L,t =´L dxw x,t , therefore the time derivative becomes a loop integral on the boundary ∂L of the 'flow' ∇w x,t along the normal directionn x . As the continuum limit of Eq. (10), it shows that the entanglement entropy's dynamics is governed by the boundary. The above PDEs (30) and (31) describe the evolution of weights, one relies on Eq. (6) to connect to the entanglement evolution. Alternatively, one can directly focus on the entanglement entropy and design a coupled nonlinear diffusive epidemiology model to describe the entanglement dynamics, as we present in Appendix F. We note that both nonlinear diffusion equations [9] and epidemiology models [78] have been separately used in modeling quantum information scrambling. This phenomenological model shows an interesting combination of both to describe a unique CV entanglement growth process.

Entanglement light cone and sudden growth
We now proceed to identify unique phenomena for the entanglement evolution. One of such an evolution is depicted in Fig. 3(a), where we see entanglement diffusively spreads from the source at the origin. We can introduce an entanglement light cone (green lines) and a wave front of the entanglement sudden growth (black lines), as will be explained in the following paragraphs.
To further understand the dynamics, we focus on particular modes and consider S(x, t) as a function of t (see Fig. 15(a) for examples). We observe a three-stage evolution: (1) In the first period 0 ≤ t < T 1 , the entanglement S(x, t) is almost zero. This is the time period before the entanglement light cone reaches location x. (2) In the second period T 1 ≤ t < T 1 + T 2 , the entanglement light cone reaches the spot and causes a rapid increase in S(x, t) , after which S(x, t) gets close to S(x, ∞) .
(3) In the last period T 1 + T 2 ≤ t, S(x, t) gradually saturates towards S(x, ∞) . We numerically investigate the scaling of T 1 , T 2 with respect to the distance ∆x to the initial squeezer, which reveals the physics of entanglement growth.
In the first period t < T 1 , the entanglement entropy at each spot is negligible. As we see in Fig. 3(a), the green curves show the threshold T 1 for spots at different distances ∆x-a parabolic entanglement light cone much slower than the usual linear light cone of operator spreading. The second period T 1 ≤ t < T 1 + T 2 describes the wave-front of entanglement rapid growth. As we see in Fig. 3(a), the black curve depicts the threshold T 1 +T 2 for spots at different distances ∆x. The parabolic shape again indicates a diffusion behavior.
This parabolic light cone can be explained by our statistical theory. We want a constant fraction of the maximum entanglement S 0 (r) = S(η L,t ) in Eq. (6), combining with Eq. (26) or Eq. (28) we can solve T 1 , T 1 + T 2 precisely, despite the analytical formula being lengthy, one immediately recognizes the scaling with some function f of the squeezing strength r. Indeed, as shown by Fig. 15(b), the green curve (entanglement light cone) and the black curve (entanglement sudden growth) both agrees well with the quadratic fitting. This is indeed consistent with the OTOC diffusion identified in [61], revealing an unique universal behavior intrinsic to CV quantum networks and absent in DV circuits [9,10].

E. Growth of multi-partite entanglement measured by distributed sensing
So far we have focused on bipartite entanglement between a subsystem L and its complement R. In quantum networks, many applications often require multipartite entanglement, which is in general difficult to characterize [72]. Here we take an operational approach from a quantum sensing perspective. An important application of the entanglement generated in such a random quantum network is distributed sensing [51][52][53][54], where multi-partite entanglement enables an improvement in the measurement sensitivity. In the case of measuring uniform real displacements of amplitude α on all modes, one can prove that considering a total mean photon number |G|N S , the optimum |G|-mode separable state can only offer a variance in estimating the displacement α (the standard quantum limit). Therefore beating the above precision limit is an evidence of entanglement. In fact, one can show that the optimum precision attainable by all entangled state is which possesses the Heisenberg scaling of V E ∼ 1/|G| 2 that is only possible with multipartite entanglement. Therefore, we can define an entanglement witness for any |G|-mode state ρ with total energy |G|N S as where V (ρ) is variance achievable by performing a localoperations (LO) and classical communications (CC) on the input ρ (hence LOCC ), while keeping the total energy conserved (see Fig. 16). We maximize over all such LOCC schemes. We have E(ρ) = 0 for all separable states [75], and E(ρ) ≤ log 2 (V C /V E ) ∼ log 2 |G| for all states.
For the state ρ(t) generated in the single-squeezer random network at time step t, we can design an estimator to obtain a lower bound of E(ρ(t)). Given the weight {w x,t } on distributing the SV state to |G| modes, we can design the following measurement protocol. First, we perform a phase rotation on each mode such that the displacements act on the corresponding squeezed quadrature; then we perform homodyne measurements on the corresponding quadratures to obtain the results {α x,t } x∈G . The estimatorα = x∈G √ w x,tαx,t / x∈G √ w x,t , which gives the variance where the effective number of modes |G| = ( x∈G √ w x,t ) 2 ∈ [1, |G|] that are entangled provides the advantage.
Combing the weights in Eqs. (25) and (27), we can obtain the effective entangled mode number for the D dimensional Cartesian graph which leads to the entanglement witness before the boundary effect comes in, when the effective modes become comparable to |G| at t ∼ |G| 2 .

IV. MULTIPLE SQUEEZERS: SPARSE LIMIT
In Sec. III, we focus on random networks with a single squeezer and present a thorough theory for the entanglement dynamics and equilibrium, via an exact mapping to random walk on a graph. Quantum networks are likely to have multiple squeezers; therefore, we extend our analyses to random networks with multiple squeezers in this section. A surprising linear superposition law is numerically observed and theoretically explained. We begin with an intuitive example of three squeezers in 1-D Cartesian graph in Fig. 3(b). The overall evolution of the entanglement entropy looks like a linear superposition of three independent squeezers, despite the nonlinearity of the entanglement dynamics. Following this observation, we consider a random circuit C with N q squeezers at different space-time coordinates {ξ k = (x k , t k )} is a simple sum of the ensemble averages S k (L, t) . Here S k (L, t) is generated from a random circuit C k with a single squeezer of strength r k at ξ k , therefore can be calculated by the random-walk mapping in Eqs. (12) and (6). Note that the random beamsplitters in all N q singlesqueezer circuits {C i } Nq k=1 and the original circuit C are independent. To test linear superposition, we numerically calculate the deviation per mode To evaluate the relative deviation, we can also rescale the deviation relative to the steady state value, as δ spp (t) = ∆S spp (t)/( S(L, ∞) 1 /|G|). Both deviations are system-size independent in the continuum limit.  Figure 18. Schematic of the multi-mode transforms. Solid lines with color red, green and blue represent energy flow from three squeezed mode at different space-time. After the random circuit, we apply a set of unitaries UL s ,t to concentrate the squeezed modes into each single modes aL s;s . The dashed lines represent residues from other squeezed modes after passing through the unitary transform UL s,t . The residue of other squeezing modes can be neglected in the sparse-squeezing limit.
In Fig. 17(a), we evaluate the deviation for the 1-D three-squeezer case considered in Fig. 3(b). We see that the relative deviation δ spp (t) is small (< 2%) through the entire dynamical evolution. To be more explicit, in Figs. 17(b) and (c), we directly plot S(x, t) (blue) at various times, which agrees well with the superposition result S spp (x, t) (red).
Following the above observation, Sec. IV A provides a theory at the sparse squeezers limit, which predicts linear superposition for both the equilibrium Page curves and the dynamical evolution; These two aspects are then investigated in Sec. IV B and Sec. IV C.
A. Theory of the sparse squeezers limit Inspired by the above numerical findings in 1-D, we present the following theory to explain the linear superposition law. For the circuit C, similar to Eq. (3), each mode at x ∈ G and time t can be written as where we use the simplified notation ξ = (x, t). Since the squeezers' locations ξ k are different, all fully random phases {θ ξ;k } are independent. Each set of positive weights {w ξ;ξ k } x∈G describes the energy-splitting among the modes of each single-mode SV a SV;k . From these weights, it immediately follows that a total portion η L,t;k = x∈L w ξ;ξ k of each single mode squeezer end up in subsystem L.
To evaluate the entropy of L, similar to Sec III A, we consider a set of passive linear optics unitaries {U Ls,t } Nq s=1 to manipulate the power distribution of the SV within L (see Fig. 18). The first transform U L1,t acts on the entire system L 1 = L to concentrate the first SV part to the n q = 1, r = 5 n q = 0.5, r = 3 n q = 0.5, r = 5 n q = 0.1, r = 5 Figure 19. (a) Schematic of the circuit generating the CV Page curves. Because the passive transform U K acting on vacuum still produces vacuum, the Euler decomposition (bottom) on vacuum input is equivalent to applying the passive transform U K after a layer of squeezers with strength {r k } mode at an arbitrary mode y 1 ∈ L 1 as where η L1,t;1 = x∈L1 w ξ;ξ k is the total portion of a SV;1 in the subsystem L 1 . For the first transform, we have η L1,t;1 = η L,t;1 . We can denote the remaining (M − 1) modes after the transform as L 2 , with the weights {w (1) ξ;ξ k } x∈L for squeezers {a SV;k } Nq k=2 . Note that a w (1) (y1,t);ξ k portion of each squeezer a SV;k (k ≥ 2) is also mixed in the mode a L1;1 , and not in the new subsystem L 2 .
After the first (s − 1) transforms (N q ≥ s ≥ 2), we have subsystem L s with |G| − s + 1 modes . The s-th transform U Ls,t acts on the subsystem L s to concentrate the power of the SV a SV;s to the mode at y s ∈ L s as a Ls;s = √ η Ls,t;s a SV;s + Nq k=s+1 e iθ (s) First, we assume the weights {w ξ;ξ k } x∈G are independent between different squeezers (indexed by k). This assumption is true in two scenarios: (1) when any two squeezers are far away from each other in terms of the shortest path on the graph and/or the time difference. In this case, the independence of weights hold at any time t; (2) as long as the number of squeezers N q |L| and the squeezers do not act on top of each other, the correlations between the weights will be small at late time. This is especially true for the Page curve S(L, ∞) . The independence of the weights means that after the s-th passive linear unitary, the new weights {w (s) ξ;ξ k } x∈L for the remaining squeezers are still randomly distributed across the entire system (not concentrated on any mode).
Second, we assume that the total number of squeezers is small, i.e., N q |G|. This allows us to approximate the weights concentrated from subsystem L s as the global total weights, i.e., η Ls,t;s η L,t;s , 1 ≤ s ≤ N q . Therefore, we conclude that after the set of passive unitaries {U Ls,t } Nq s=1 , we can approximately obtain a set of modes which is a product of lossy single-mode SVs. From additivity of entropy, we arrive at linear superposition S(L, t) S spp (L, t) , note that each single-squeezer term can also be written out explicitly as S(η L,t;s ) from Eq. (1) by replacing r with r k .
Note that linear superposition law holds not only for the dynamical evolution at any t (as systematically explored in Section IV C), but also for the equilibrium Page curves at t = ∞ (see Section IV B).

B. Page curves with multiple squeezers
In terms of Page curves, as we explained in Sec. III B, the graph topology is irrelevant as the entire dynamics is equivalent to a passive, Haar-random, global unitary; therefore, we can simply consider the 1-D case. Furthermore, a Gaussian unitary can be decomposed into a layer of squeezers concatenated by passive unitaries (Euler decomposition, see Appendix A), we can effectively push all squeezers to the first step and then apply a global Haar random passive linear transform (see Fig. 19(a)). Therefore, each Page curve is characterized by a list of squeezing strength {r k } Nq k=1 . In this case, as long as N q M , we can regard the squeezers as sparse. We denote the average squeezing strength as r = Nq =1 r /N q and the density of squeezers n q = N q /M . When N q |L|, superposition S(L, ∞) S spp (L, ∞) holds, with each single-squeezer result given in Eq. (16). In the large squeezing limit, we can further utilize Eq. (17) to obtain We see a dependence on statistical quantities r, n q , while the shape of the curve is invariant in the bulk at large squeezing limit, as can be seen in Fig. 20 (b) and (c). We numerically examine the validity of the linear superposition in Page curves via the relative deviations δ spp (t) in Fig. 20. Indeed when the squeezer density is low, we see a good agreement, as shown in Fig. 20 (b); while when the squeezers are dense, substantial deviation can be found, as shown in Fig. 20 (c). The transition is captured by the relative deviation δ spp (∞) in Fig. 20 (a), where δ spp (∞) increases linearly with the squeezer density n q .
Although when n q is not small, superposition does not hold, we can further numerically explore the Page curves' dependence on the parameters. To consider the typical case, we randomly generate {r k } Nq k=1 and plot the normalized Page curves in Fig. 19(b) for different system sizes. We find that all Page curves coincide as long as r, n q are the same; and independent of M when M is large Similar to the single-squeezer case. Therefore, we can plot the case with identical squeezer r as a benchmark. To demonstrate the coincidence, we consider random cases with r i ∈ [0.8r, 1.2r], 1 ≤ i ≤ N q , while guaranteeing the average equaling r [74]. Indeed, all the random points lie right on top of the benchmark of uniform squeezing with the same r, n q . Since r and n q are both scale-free statistical quantities, this indicates a well-defined continuum limit of CV Page curves. Fig. 17 already confirms superposition through the entire evolution for a simple case, to verify it in a more general setting, we consider circuits with squeezers randomly distributed in space-time. Guided by the theory in Sec. IV A, we expect the minimum distance between the squeezers to be the dominating factor of the devia- tion ∆S spp (t); However, simple random sampling methods inevitably lead to the appearance of clusters, where squeezers can be close to each other. In order to tune the minimum distance d, while preserving the random nature of the circuit set-up so that the results apply in general, we adopt the random Poisson sampling method [79] to control the space-time distances between the squeezers. With the squeezers randomly chosen, we evaluate the entanglement dynamics for circuits with random squeezers of different minimum distances. Two examples are given in Fig. 21: When the squeezers are sparse (d = 60 is large) as indicated by the cyan dots in Fig. 21(c), the superposition results S spp (x, t) (dashed lines) agree well with the true values S(x, t) (solid lines) at various time steps in Fig. 21(d); when the squeezers are dense (d is small), substantial deviations from the linear superposition can be observed, as demonstrated in Fig. 21(a)(b) for the case of d = 20.

C. Dynamics with linear superposition
To systematically examine the transition from 'dense' to 'sparse' squeezers, we calculate the deviation ∆ spp (t) for distances 0 ≤ d ≤ 500 at various time steps in Fig. 22, where the deviation decreases monotonically with d up to numerical precision. Although for a smaller d, the deviation is larger, the relative deviation is still below 10% of the static value S(x, ∞) 1 /M (orange). Therefore, we conclude that in generic CV quantum networks with sparse squeezers, linear superposition of entangle- ment growth holds.

V. DISCUSSION
In this paper, we reveal a mapping between entanglement formation dynamics in random CV networks to random walk on general graphs. This mapping allows analytical solutions of the entanglement entropy dynamics, Page curves and scrambling time for an arbitrary network topology. On networks respecting locality, the solution enables the understanding of three unique features of entanglement formation dynamics-an entanglement light cone, an entanglement sudden-growth period and parameter-dependent Page curves. Our results have implications in quantum network protocol design, e.g., the entanglement light cone will place bounds on the latency in entanglement distribution, and also on the fundamental understanding of many-body systems. Lastly, let us point out some future directions: it will be interesting to extend the entanglement light cone to long-range interacting systems; extending the multipartite entanglement witness to more general input states will bring further insights of quantum entanglement; exploration of the connection between entanglement dynamics with statistical properties of random networks such as connectivity distribution will lead to a full statistical theory of CV quantum networks. It is worthy to comment that the diffusive growth of entanglement identified in Sec. III D 2 does not contradict the known result that good approximations (t-design) of the global Haar random unitary only requires a random circuit with the number of layers linear in system size [76]. Mathematically, in a CV system it is the covariance matrix V that is transformed as V → SV S T under the Haar random matrix S. Both the von Neumann and Renyi entanglement entropies are given by the symplectic eigenvalues of V , which is not a simple polynomial of S. This contrasts with the case in DV systems, where the random matrices act on the density matrix directly and Renyi entropies are simple polynomials of them.
Appendix C: Details of the random-walk mapping In terms of the weights, consider two modes x, x on a connected edge A general passive gate U (t, x, x ), with transmissivity τ , on modes x, x would lead to the evolution where the angles θ = (θ x ,t − θ x,t ), θ x,t+1 , θ x ,t+1 are uniform in [0, 2π) and we left out the vacuum terms. Therefore we can write the new modes in the same form as Eq. (3) with and the angles θ x,t+1 , θ x ,t+1 uniformly random in [0, 2π).
We can obtain the ensemble averaged (over the gates U (t, x, x ), with transmissivities τ uniform in [0, 1)) evolution as This describes a symmetric random walk along the edge xx .
Appendix D: Calculation of the variance of ηL,∞ at equilibrium The weights {w x,t } determine the entanglement, essentially they can be considered as the amplitude squared of a random complex number α x,t , i.e. w x,t = |α x,t | 2 . Normalization requires x∈G |α x,t | 2 = 1. Due to the Haar randomness at t = ∞, we assume that this is the only constraint, therefore α x,t 's are random complex numbers on a high dimension sphere. We can calculate the distribution of weights η L,∞ through Here we utilized the N -dimensional sphere area formula where we have used δ(x 2 −a 2 ) = (δ(x−a)+δ(x+a))/(2a).
One can easily verify that Eq. (D4) is normalized and η L,∞ = |L|/|G|. Moreover, the fluctuation can be obtained as var(η L,∞ ) = |L||R| |G| 2 (|G| + 1) . (D7) Appendix E: Continuum limit of the multiple-squeezer case In this section, we show that the same continuum limit, identified in Sec. III D 1 (see Fig. 13) for a single-squeezer, still holds when there are multiple squeezers. In Fig. 23, one layer of N q = M squeezers are applied in the first time step, and we still see a good agreement between the relative 1-norm ( Fig. 23(a)) at different times. Snapshots of entanglement entropy curve at various times also overlap entirely (Fig. 23(b), (c)), after a re-scaling of spacetime. Note that in this case, even with a dense squeezer distribution where the superposition principle does not hold anymore, this continuum limit still holds. (c) Figure 23. Relative 1-norm deviation δ dy (t) for dynamic entropy curves in 1-D Cartesian graph. One layer of Nq = M squeezers is applied at t = 0 with an identical strength r = 5.
All curves are re-scaled and shown in the same way as in Fig. 13. Although the dynamics of S(x, t) can be captured by the diffusion of weights and Eq. (6) that connects weights to entanglement entropy. Alternative models are possible to capture major phenomena in Sec. III D 2 for the 1-D Cartesian graph. We notice that in our CV circuits, squeezing behaves as the source of entanglement generation [77]; and while it diffusively spreads out, its strength at each mode also decays due to the effective loss during interactions with other modes; therefore we introduce a field G(x, t) to model this diffusive source that triggers entanglement growth and device the following set of coupled diffusion-growth equations to give a theory prediction S T (x, t) for the ensemble averaged entanglement entropy S (x, t) where A, c, D, f are four constants. Eq. (F1b) describes a nonlinear diffusion process of the 'source', where the nonlinear term f G(x, t) 2 corrects the early stage dynamics. Eq. (F1a) describes an 'infectious' saturation process with a growth rate proportional to 'uninfected population', c − S T (x, t), triggered by the source G(x, t). We note that both nonlinear diffusion equations [9] and epidemiology models [78] have been separately used in modeling quantum information scrambling. This model shows an interesting combination of both to describe a unique CV entanglement growth process.
The initial condition for the source field G (x, t) is a delta-function at the squeezer's position, and uniform zero for the entanglement S T (x, t). We choose G (x, t) obeying the von Neumann boundary condition (reflection boundary conditions) such that ∂ x G(±N, t) = 0. It is straightforward to check that S T (x, t) = c is a steady state solution of Eqs. (F1). Considering the steady state of the real entanglement entropy in Fig. 6, we expect Eqs. (F1) to only describe the evolution of entanglement in the bulk. Indeed, when we numerically solve Eqs. (F1) in a system with M = 501 modes and fix the steady state c to be the maximal height S 0 (r), we find good agreement before boundary effects come in at late time (see Fig. 24(a)). In particular, this model does not guarantee zero entanglement entropy on the boundary. Therefore, Eqs. (F1) will be able to capture the entanglement growth of the bulk, before boundary effects become important. More precise models can be constructed by further fine-tuning Eqs. (F1) to different boundaries. As an example, we can implement a position-dependent c = S (x, ∞) , such that steady state is consistent with the Page curve. Then, we numerically solve Eqs. (F1) and find the theory prediction in a nice agreement with the real data from random circuit simulations in Fig. 24(b).