Probing the spectral dimension of quantum network geometries

We consider an environment for an open quantum system described by a"Quantum Network Geometry with Flavor"(QNGF) in which the nodes are coupled quantum oscillators. The geometrical nature of QNGF is reflected in the spectral properties of the Laplacian matrix of the network which display a finite spectral dimension, determining also the frequencies of the normal modes of QNGFs. We show that an a priori unknown spectral dimension can be indirectly estimated by coupling an auxiliary open quantum system to the network and probing the normal mode frequencies in the low frequency regime. We find that the network parameters do not affect the estimate; in this sense it is a property of the network geometry, rather than the values of, e.g., oscillator bare frequencies or the constant coupling strength. Numerical evidence suggests that the estimate is also robust both to small changes in the high frequency cutoff and noisy or missing normal mode frequencies. We propose to couple the auxiliary system to a subset of network nodes with random coupling strengths to reveal and resolve a sufficiently large subset of normal mode frequencies.


Introduction
Networks [1,2] describe discrete topologies that can capture the architecture of complex systems, from the brain to societies. As an increasing number of datasets about complex systems are becoming available, the characterization of real-world network structures has significantly enriched the understanding about the relation between network topology and dynamics. A classic result of Network Theory is that the statistical properties of the complex network topology, including for instance the degree distribution, strongly affect the properties of dynamical processes such as percolation,epidemic spreading and spin models [3].
However only recently the scientists have realized that a large variety of network datasets have an intrinsic geometrical nature that affects their dynamical properties.
The term "network geometry" [4,5] refers to discrete topological spaces with notable geometrical features. These structures can be modelled by using simplicial complexes [4,6] which are discrete structures not only formed by nodes and links but formed also by triangles, tetrahedra and so on. A general feature of network geometries is that they display a finite spectral dimension [7,8], i.e. the random walks defined on these structures relax to its stationary state only algebraically like in finite dimensional Euclidean lattices. The spectral dimension of an d-dimensional Euclidean lattice is d, however it is well known that also fractals have a finite spectral dimension [7], and only recently it has been shown that several real network datasets also display this important geometrical property [9].
Recently a general theoretical framework called "Network Geometry with Flavor" (NGFs) [10,11] has been proposed to model network geometries using simplicial complexes.
The NGFs depend on the dimension d of their building blocks and on another parameter s called flavor. For any dimension d and any flavor s the NGFs capture the main characteristics of complex networks including modularity, small-world properties and heterogeneous degree distribution. Moreover these topologies display also a finite spectral dimension which is the signature of their geometrical nature. It is also worth noting that the spectral dimension can change significantly the dynamical properties of classical and quantum random walks [8,12,13], and the Kuramoto model leading to an anomalous frustrated synchronization [14,15].
During the recent years, Network Theory has also opened significant possibilities within the framework of quantum physics [4,16]. Indeed, the concept and use of quantum complex networks are becoming increasingly common. Here, the interest is focused, e.g., on quantum correlations and phase transitions [17], quantum communication networks and internet [18], quantum walks [19], and generalizing the classical concepts of Network Theory to quantum domain [20,21].
Our current interest lies on using quantum complex networks within the framework of open quantum systems, where the environment -that an open system interacts with -consists of a network of interacting quantum entities. Intriguing possibilities are now provided by using an open system to probe the properties of the quantum complex network. Here, one assumes that the only available information about the quantum network is the information that can be inferred by studying the properties of the probe system only. This approach has been fully developed in the case in which the quantum network is formed by a network of harmonic quantum oscillators coupled to an open system oscillator that acts as the probe [22,23]. Moreover, the experimental implementation of the theoretical framework and the probing schemes, have been recently proposed using a multi-mode optical setup [24].
Building on these results, in this paper we investigate the QNGF which is the environment of a quantum open system (the probe) and is formed by a network of quantum harmonic oscillators coupled according to the topology of the NGF.
Quantum algorithms to investigate the topological properties of simplicial complexes have been proposed in previous works [25] however not using the theory of open quantum systems.
We are interested in how to obtain crucial information about the geometrical properties of the QNGF without having any a priori knowledge of the generating mechanism of the NGF or of its structure. In particular our main goal is to infer the value of the spectral dimension of the underlying topology of the QNGF. We provide a connection between the eigenfrequencies of the QNGF and the eigenvalues of the Laplacian matrix of the NGF. Using the theory of open quantum systems combined with data science we are able to probe the value of the spectral dimension of the QNGF. Moreover, we also consider cases where the probe has limited capacity to reveal the full set of eigenfrequencies of the oscillator network, and also when the obtained data is noisy, i.e., influenced by random fluctuations of the eigenfrequencies. The results indeed demonstrate that the probing scheme for the spectral dimension is rather robust when not having fully ideal setting at hand.
The paper is organized as follows. We introduce the considered quantum network model in Sections 2.1 and 2.2, focusing on the topology and physical model, respectively, and introduce the open quantum system used to probe its properties as well as the interaction term in Section 2.3. Specifically, we consider networks of identical interacting quantum harmonic oscillators where the network structure is generated with Network Geometry with Flavor (NGF). The relation between the eigenvalues of the unweighted Laplacian matrix and the frequencies of the normal modes of the oscillator networks is given and used to connect the spectral properties of the Laplacian matrix with those of the physical network. In Section 3 we show how the spectral dimension of the eigenvalues can be estimated from the normal mode frequencies. In Section 4, we consider the impact of small errors and missing values on the estimate and find it to be robust to both. While the normal mode frequencies are possible to probe with an auxiliary system, the normal mode structure can make it difficult to reveal and resolve enough of them; we propose to couple to multiple network nodes with random couplings to avoid this problem. We conclude in Section 5.

Network Geometry with Flavor
NGF [26] is a model of random simplicial complexes, also interpretable as random hyper-graphs. Simplicial complexes are ideal to model discrete network geometry because they are formed by geometrical building blocks such as triangles, tetrahedra and so on. Specifically these building block are called simplices. A d-dimensional simplex include d + 1 vertices forming a fully connected graph with d + 1 nodes. A face of a simplex is a lower dimension simplex formed by any proper subset of its nodes. A simplicial complex is formed by simplices that are connected along their faces. The dimension of a simplicial complex is defined to be the highest dimension of its simplices. The NGF are pure d-dimensional simplicial complexes, i.e. they are formed for a set of d-dimensional simplicies connected along their (d−1)-dimensional faces; a NGF of d = 1 consists of links connected through their nodes, a NGF of d = 2 of triangles connected through their links, a NGF of d = 3 of tetrahedra connected through their triangles and so on.
A d-dimensional NGF of N nodes evolves from a single d-dimensional simplex by choosing at each timestep a (d − 1)-dimensional face to which an additional simplex is added, increasing the number of nodes by one. The probability to choose a particular face is affected by the number of simplices already incident with it, as controlled by flavor s. Let n α be the number of d-dimensional simplices incident with face α minus one. The probability Π s (α) to choose face α is given by where Z s is Z s = α (1 + sn α ) (2) and the sum is over all (d − 1)-dimensional faces currently in the growing NGF, ensuring normalization. A rich variety of simplicial complexes with distinct properties can be generated by adjusting the dimension d, and flavor s. In particular a NGFs with s = 1 evolves thanks to an explicit generalized preferential attachment mechanism while NGFs with flavor s = −1, 0 do not obey an explicit preferential attachment rule. Additionally, a d-dimensional NGF with flavor s = −1 is a discrete manifold; all (d − 1)-dimensional faces are incident to either 1 or 2 d-dimensional simplices.
The NGFs display an emergent hyperbolic geometry [11] as well as a finite spectral dimension d S [27] which is a further indication of their spontaneous geometric character. The spectral dimension is a important spectral property strongly affecting the long time behaviour of the random walk defined on networks. In particular the spectral dimension of a generic (connected) network of N nodes is defined as follows. Let L be a Laplacian matrix of elements where δ ij indicates the Kronecker delta, d i is the degree of node i and a ij is an element of the adjacency matrix of the network. Let {λ i } indicate the eigenvalues of the Laplacian matrix L listed in a non-increasing order, i.e.
If density of eigenvalues ρ(λ) satisfies for λ 1, with the first non-zero eigenvalue (also called the Fiedler eigenvalue) λ 2 vanishing in the large network limit as then d S is called the spectral dimension of the network. Consequently the cumulative distribution ρ c (λ) scales as for λ 1, where θ(x) = 1 if x ≥ 0 and θ(x) = 0 otherwise. In other words, when the Fiedler eigenvalue of the networks λ 2 goes to zero in the large network limit, and ρ(λ) and ρ c (λ) have a power-law scaling, the network has a spectral dimension determined by this scaling at low values of λ. The spectral dimension for a regular lattice reduces to the Hausdorff dimension The panel on the right shows ρc(λ) for two NGFs with larger size (N = 2000) but otherwise same parameters. The power-law scaling for λ 1 is a sign of a finite spectral dimension d S .
of the lattice. In general for networks and fractal geometry, the spectral dimension can be understood as the effective dimension of the networks as probed by a random walker moving on the network. NGFs display a finite spectral dimension for every flavor s ∈ {−1, 0, 1}. Interestingly the presence of a finite spectral dimension is robust to modification of the models including non-integers flavors s = −1/m with m > 1 and the generalization of NGFs to cellcomplexes, i.e. discrete topological structures obtained by gluing regular polytopes others than simplices such as hypercubes, orthoplex and so on.
Here for simplicity we focus exclusively on the case s = −1. In this case the spectral dimension d S is an increasing function of d. Examples of both NGF with d = 2 and d = 3 flavor s = −1 and the scaling of the associated ρ c (λ) are shown in Fig. 1. From this figure it possible to clearly appreciate that the spectral dimension d S is increasing with increasing values of d.

Quantum Network Geometry with Flavor
In this work we investigate the properties of a quantum network geometry called QNGF. The QNGFs are formed by a set of N quantum harmonic oscillators interacting through the topology described by NGF. Each quantum harmonic oscillator has the same frequency ω 0 ≥ 0, and each pair of connected oscillators is coupled with the same interaction strength g > 0. The resulting Hamiltonian reads where p and q are vectors indicating the momentum and position operators of each of the N nodes of the NGF. The matrix A present in the Hamiltonian H can be expressed in term of the Laplacian L as Since the Hamiltonian (8) is quadratic in position and momentum operators and A is positive definite, we may find a basis of non-interacting normal modes [28]. To this end let us now define the matrix U whose columns are the eigenvectors of A, then the diagonal matrix D whose diagonal elements are the eigenvalues of A, is given by It follows that by defining new operators for the normal modes as Q = U q and P = U p the Hamiltonian takes a diagonal form where the modes are clearly decoupled, i.e.
and their frequencies ω i are given by D ii = 1 2 ω 2 i . We may now relate the eigenvalues of the NGF to the normal mode frequencies of the associated QNGF. In fact from Eq. (9) it is clearly evident that the matrix A and the Laplacian can be diagonalized on the same basis, i.e. the eigenvectors of A are eigenvectors of L and vice versa. Since the eigenvalues of A are given by D ii = ω 2 i /2 Eq. (9) implies that the eigenvalues λ i of the Laplacian can be expressed in terms of the frequencies ω i as Therefore for quantum network geometries displaying a finite spectral dimension (like QNGF), we can consider the cumulative distribution of the normal mode frequencies P c (ω), defined as The cumulative distribution of normal model freqquencies P c (ω) is related to the cumulative density of the eigenvalues ρ c (λ) of the Laplacian by In presence of a finite spectral dimension d S the cumulative density of eigenvalues ρ c (λ) of the NGF Laplacian obeys Eq. (7). Therefore it follows that P c (ω) scales like This relation is key for our work as it implies that the knowledge about the normal mode frequencies reveals information about the spectral dimension d S of the QNGF.

QNGF as a quantum environment of an open quantum system
We consider the QNGF as a quantum environment interacting weakly with a single probe formed by an open quantum system. By only considering local observables of the probe our general problem is to study which properties of the QNGF can be inferred.
In previous studies of quantum complex networks it has been shown that a similar approach is able to infer several important properties of the quantum networks, e.g., the network spectral density [22] or degree sequence [23]. Here we focus specifically on the problem of inferring the spectral dimension of a given QNGF from the normal mode frequencies {ω i }. In this Section the total Hamiltonian is defined and the principle behind probing the normal mode frequencies is explained. Although here we assume that the network geometry is generated with NGF, our approach is very general and can be applied to any quantum network geometry displaying a finite spectral dimension, provided that the Hamiltonian is of the form given by Eq. (8).
The open system is assumed to be an additional quantum harmonic oscillator with a Hamiltonian H s , given by where ω s is the frequency of the probe and q s and p s are the position and momentum operators respectively. The interaction Hamiltonian is assumed to be of the form where the vector k has elements k i ≥ 0 indicating the interaction strength between the probe and network oscillator i. The total Hamiltonian of the system including both the probe and the environment is therefore given by In the bases of the normal modes of the QNGF the total Hamiltonian H tot can be expressed as where From this expression of the total Hamiltonian it can be appreciated how the eigenvalues of the Laplacian control the physical frequencies affecting the open system while the eigenvectors affect how strong is the open system's coupling to each eigenmode. Indeed Eq. (20) relates the eigenvectors of the NGF that constitute matrix U to coupling strengths to the network normal modes. An important special case is the one in which the probe is coupled to a single node (oscillator). In this case k has only one non-vanishing element and g given by Eq. (20) is proportional to a single left eigenvector of the Laplacian L. Since according to Eq. (20) g is linear in k, and since k can be written as where e i is the vector of elements e i j = δ ij , it follows that in the general case the vector g can be interpreted as a linear combination of the interaction strength resulting from coupling to single nodes of the network.
We assume that the network is initially in vacuum, i.e. its ground state, while the probe is in some single-mode Gaussian state. Since additionally the Hamiltonian is quadratic, the dynamics can be both solved and simulated efficiently [29]; for an explicit example, see, e.g., [30].
The considered local observable is the expected value of the probe excitations, a † a , where a † and a are the probe creation and annihilation operators, respectively. This quantity is proportional to the energy of the probe.
Let a(0) † a(0) be the expected value of the probe initial excitations, and suppose it evolves with the network according to total Hamiltonian H tot for time t, taking probe excitations to a(t) † a(t) . The probe is then decoupled from the network and its excitations are measured, allowing the determination of the quantity ∆n = | a(t) † a(t) − a(0) † a(0) |. It is well-known that an open system weakly interacting with a bosonic heat bath must be in resonance with the bath to exchange information and excitations with it. In the present case this means that for sufficiently weak interaction, ∆n 0 implies that ω S ≈ ω i for some i ∈ {1, 2, . . . , N }. The principle is then to repeat the protocol for many different values of ω S ; a large ∆n indicates that a normal mode frequency is in the vicinity of the used ω S .
In the following Section we focus on the estimation of d S assuming that we have acquired the full set of normal mode frequencies {ω i } from probing them. We will discuss only later, in section 4, the more realistic situation in which we are able to obtain only some approximate values of normal mode frequencies {ω i } and briefly consider general strategies to infer {ω i } from ∆n.

Probing the spectral dimension
Given the normal mode frequencies {ω i } of a QNGF, the objective is to estimate its spectral dimension d S . Our strategy is to find the least-squares fit of {ω i } into a model with d S as a fitted variable. We begin by recasting the power-law scaling in relation (15) to a linear one by considering the logarithm of both sides, namely The fitting is now amenable to ordinary least squares and the best fit can be found applying the Moore-Penrose pseudoinverse [31]. While fitting to a linear model requires to know the bare frequency ω 0 it is not necessary to know it a priori since it will coincide with the smallest normal mode frequency; this follows from the well-known property min({λ i }) = 0 of Laplace eigenvalues and from Eq. (12). At this point we face two additional problems. First, the normal mode frequencies are expected to behave according to relations (15) and (22) only up to some upper limit in frequencies; therefore to estimate d S we should have a way to truncate the frequencies appropriately. Second, we would like to be able to say something about the quality of the estimated value and in particular be able to compare different estimates resulting from different points of truncation. To address the first point we consider the goodness-of-fit between the linear model and actual data derived from {ω i } for different points of truncation and choose the optimal value; specifically, we consider the coefficient of determination R 2 [32]. For the second, we consider the 95% confidence intervals. Further details are given in the Appendix Appendix A.
Remarkably, the estimate is independent of the values of ω 0 and g. On the one hand this is because the value of ω 0 can in principle be probed exactly which facilitates the use of units where it becomes irrelevant. On the other hand g can only scale the L.H.S. of Eqs. (15) and (22); this does not affect the goodness-of-fit and therefore not the values of R 2 nor the confidence intervals as a function of the index i which determines the cutoff frequency ω i .
A specific example of this procedure is shown in Fig. 2 where we consider a QNGF whose underlying topology is NGF of N = 2000 nodes, with d = 2 and s = −1. As explained before, the particular values of ω 0 = 0.25 and g = 0.1 do not affect the estimate. In the top panel the logarithm of P c (ω) is shown as a function of log(ω 2 − ω 2 0 ) together with R 2 as a function of the point of truncation; higher values indicate a better  agreement between the fit and the data. Here d S is proportional to the slope of the linear model. Below the corresponding confidence intervals for d S are shown. The optimal value of R 2 , showing the optimal point to truncate the frequencies, is shown by the vertical dashed line. There is clear correlation between the behaviour of R 2 and the confidence intervals, with higher values leading to a smaller interval which we interpret as a more reliable estimate for d S .
We compare the confidence intervals of Fig. 2 with those of two other networks which have otherwise the same parameters but have dimensions d = 3 and d = 4, shown as a function of the truncation frequency in Fig. 3. Here it can be better appreciated that the estimates are not particularly sensitive to small changes in the truncation frequency. On the other hand, the results suggest that this robustness decreases for increasing values of the spectral dimension d S .
Before concluding the Section, we point out the possibility to fit Eq. (15) directly to a nonlinear model. Similar arguments can be employed to show that such an estimate is also independent of the values of ω 0 and g, however in general finding a nonlinear leastsquares fit is more challenging and in particular the available algorithms, such as Levenberg-Marquardt method [33], nonlinear conjugate gradient method [34,35] and quasi-Newton methods [36], do not necessarily converge to a global optimum. This can lead to, e.g., jumps from one local optimum to another as the cutoff frequency is varied. We also remark that there are many other goodness-of-fit measures that could be employed; in the present case numerical experiments suggest that they all lead to similar results.

Robustness to missing or noisy normal mode frequencies
As outlined in Section 2.3, normal mode frequencies {ω i } can be inferred from the dynamics of an open quantum system, or probe, weakly interacting with the QNGF, by performing a frequency sweep in the low frequency regime and tracking , e.g., the change in probe excitations. In practice, only a finite set of values for ω S can be considered and consequently the probed values of {ω i } will not be exact. On the other hand, the cumulative distribution P c (ω) can be expected to be very robust to small errors. Provided that we are able to resolve different normal mode frequencies these errors would be smaller than the difference between adjacent normal mode frequencies, leaving P c (ω) almost unaffected. We have checked that even a very pessimistic case of i.i.d. relative error up to 5% leaves all results shown in Figs. 2 and 3 virtually the same. Rather than finding accurate values for {ω i }, the challenge is to find values for them at all. In the extreme case where g i = 0 for some index i in Eq. (20) it is impossible to learn the corresponding normal mode  Top: comparison of confidence intervals derived from complete (solid lines) and partial (dashed lines) data. In the latter case a random sample of 10% of the eigenfrequencies is discarded before estimation. The dots show the fitted value of d S at the optimal value of the goodness-of-fit measure R 2 . Bottom: the estimated value of d S (middle line) and confidence intervals (outer lines) when the optimal cutoff frequency is used as a function of the probability p of missing data. The network parameters are as in Fig. 2. frequency ω i since the probe is decoupled from this normal mode. Since g is a linear combination of the left eigenvectors of matrix L of the NGF, this is highly unlikely to occur as long as k = 0. What is to be expected, however, is that normal modes with frequencies in close proximity to each other interact with the probe with very different strengths. In such a case the impact on probe dynamics is dominated by the normal modes with a stronger coupling strength, making it almost impossible to reveal the other nearby frequencies. It is also difficult to resolve frequencies that are close and interact with the probe with similar strengths, in which case many frequencies might be mistakenly counted as a single frequency.
We consider the case where each normal mode frequency is equally likely to be missing in Fig. 4. In the top panel we fix the probability to miss a normal mode frequency to p = 0.1 and compare the confidence intervals when full and partial set of normal mode frequencies is available, indicated by solid and dashed lines, respectively. The fitted value of d S at the optimal cutoff frequency is shown by the dot. In the bottom panel we show the estimated value of d S which maximizes the goodness-of-fit against probability p. It can be seen that neither the confidence intervals nor the optimal estimate for d S are particularly sensitive to uniformly missing frequencies. As p increases the trend seems to be that the optimal estimate decreases, however the result is still reasonably close to original even at p = 0.3.
In a probing scenario missing frequencies might not be uniformly distributed, while in some cases a frequency might be mistakenly interpreted as a normal mode frequency. To further investigate these challenges we simulated a simple proof-of-concept scheme where the QNGF spectrum is swept using equidistant values of probe frequency and the change in probe excitations ∆ n = | n(0) − n(t) | from initial to final time t is considered as a function of probe frequency ω S . Prominent local maxima indicate resonance, i.e. that a normal mode frequency is close. A general purpose non-adaptive algorithm is employed to find these maxima and d S is estimated from the associated values of probe frequency. The considered network is a small QNGF with N = 50, d = 2, and s = −1. Further details about the dynamics, the frequency sweeps and the algorithm used to find the peaks are given in Appendices Appendix B, Appendix C and Appendix D, respectively.
The results are shown in Fig. 5. In the top panel an example of the behaviour of ∆ n is shown.
Here the shown values are averaged over 10 frequency sweeps. For each sweep the probe is coupled to 8 randomly chosen network nodes. The vertical lines indicate the values of {ω i } and the dots are their probed values, found by identifying the peaks. As can be seen, the peaks align quite well with {ω i }. Even though the majority of normal mode frequencies are found, some are missed, including rather obvious cases in the high frequency regime. While indeed there is a local maximum in the vicinity of every normal mode frequency it will be seen that avoiding false positives is important for probing d S . Consequently very cautious settings for the peak finding algorithm are used in the entire spectrum; further improvements could be expected by using a custom algorithm optimized for the task at hand.
In the bottom panel we compare d S estimated from original {ω i } (dashed vertical line) to values estimated from probed normal mode frequencies as a function of the number of performed frequency sweeps and when a different amount of randomly selected network nodes is directly interacting with the probe. The results are not sensitive to cutoff frequency in the vicinity of its optimal value, which is used throughout. It can be seen that the accuracy tends to improve with the number of sweeps. This is to be expected since each set of network nodes corresponds to a different linear combination of the associated left eigenvectors of L. Different sweeps therefore change which normal modes are interacting too weakly to be distinguished. While the trend is not as strong, it also seems to be the case that it is better to couple to multiple network nodes than just one. This is actually connected to the number of false positives which tends to be highest when coupling to a small number of network nodes; in terms of correctly identified {ω i } alone there is no evident difference between the three cases, which all saturate to approximately 80% of {ω i } with the used algorithm and its settings.
In general, the challenges revolve around resolving a sufficiently large fraction of {ω i } while avoiding false positives, which in practice requires large differences of ∆ n when in and out of resonance. The difficulty increases with the density of {ω i } and therefore with N and d, while it tends to decrease with longer interaction times, weaker interaction and larger difference between probe and network initial energy. Multiple frequency sweeps and coupling to multiple network nodes can be expected to remain beneficial. False positives in particular seem to cause a loss in the rise and fall behaviour of the goodness-of-fit measure, moving the maximum value to high frequency regime. This could potentially be used to eliminate them.

Conclusions
Finite spectral dimension is characteristic of network geometries. A very general mathematical framework for generating a large variety of network geometries is the "Network Geometry with Flavor" model. The simplicial complexes generated by this model obey simple stochastic rules yet from these simple rules network geometries emerge spontaneously. Here we have proposed the model QNGF which is formed by a set of quantum harmonic oscillators interacting through a hyperbolic simplicial network geometry given by NGF. Using the relation between the normal mode frequencies of QNGFs and the eigenvalues of the Laplacian of NGFs we show that the spectral dimension of NGF can be probed with an open quantum system.
The obtained estimate of the spectral dimension is independent of the bare frequency of the network oscillators and the strength of the couplings between them and is not sensitive for the choice of the cutoff frequency, especially when the spectral dimension is small. Our estimate of the spectral dimension is formed from the normal mode frequencies of the network, which are in principle possible to deduce from the reduced dynamics of the probe. In practice some deviation from the exact values and incompleteness of the set of frequencies can be expected. Our results show that the estimate remains robust to small deviations and uniformly missing frequencies. As a further proof-of-concept we simulated probing of the normal mode frequencies and the spectral dimension for a small network. The results reveal that misinterpreting a frequency as a normal mode frequency can be just as harmful for the estimation as not finding one, and in particular the two effects do not tend to cancel each other out. To reduce both errors simultaneously we propose to perform multiple frequency sweeps where the probe is coupled to multiple randomly selected network oscillators. Indeed, in this way the estimate is close to that of the ideal case. We expect the proposed method to be useful also for probing just the normal mode frequencies.
The QNGF provides a very flexible benchmark model to investigate quantum network geometry. In particular, this allows us to explore the robustness of the results in higher dimension, by comparing the results obtained on NGF evolving according to the same rules but having different dimension d of their building blocks and different spectral dimension d S . We show that probing the spectral dimension can be performed on QNGF formed by simplices of different dimension d although we achieved higher accuracy of the results for lower dimensions. Interestingly the approach proposed in this work and tested on QNGF can be applied to arbitrary quantum network geometries displaying a finite spectral dimension.
Acknowledgments J.N. thanks Turku Centre for Quantum Physics for the hospitality.
Appendix A. Used goodness-of-fit measure and determination of parameter confidence intervals The coefficient of determination R 2 is a measure of how well the predictions of the fit approximate the data points, defined as where f i are predicted values andȳ is the mean value of the data points. Notice that the denominator is proportional to the variance of the data. Normally 0 ≤ R 2 ≤ 1, where R 2 = 0 indicates that the variance in the data cannot be explained at all by the fit while R 2 = 1 indicates that the variance can be explained perfectly.
Besides R 2 , we also considered adjusted R 2 [37], Bayesian information criterion [38] and Akaike information criterion with [39] and without [40] correction for small sample sizes. Unsurprisingly, adjusted R 2 behaves very similarly to R 2 , while the latter three have very similar behaviour with each other but differ from that of R 2 and adjusted R 2 . Specifically, for them the optimal value seems to appear at somewhat higher frequencies than for R 2 and adjusted R 2 , which in turn tends to correspond to a bit smaller estimate for d S .
The considered parameter confidence intervals are defined to be intervals of values such that they include the true value 95% of the time. To determine the intervals it is assumed that the error in the parameters is normally distributed. Let α = 1−0.95, n the number of data points and p = 2 the number of parameters. Then the intervals are where d S is the estimated spectral dimension, t(1 − α/2, n − p) is the 100(1 − α/2) percentile of Student's The Gaussian states considered in this work are a paradigmatic class of states in continuous-variable quantum information science, and may be defined as the states determined completely by the second and first moments of the momentum and position operators [41,42]. Such states are conveniently described by their covariance matrix σ. Consider the Hamiltonian H tot in the basis of network normal modes, given by Eq. (19). Let X = {Q 1 , Q 2 , . . . , Q N , q s , P 1 , P 2 , . . . , P N , p s }. In this basis the covariance matrix is where the angle brackets denote an expectation value over the state and [X i , X j ] + = X i X j + X j X i is the anticommutator. As a Hamiltonian quadratic in momentum and position operators, H tot preserves the Gaussian character of the state. Consequently the evolution it induces is completely captured by the evolution of the covariance matrix. Let σ(0) be the initial form of the covariance matrix of Eq. (B.1) and let the probe and network interact for time t, taking the covariance matrix to σ(t). Then the initial and final forms of the total covariance matrix are related as where the matrix S(t) is induced by the Hamiltonian H tot and the interaction time t.
The explicit form of S(t) is easily found by using the analytic solution for the equations of motion of noninteracting oscillators. To apply it, the total system of network and probe is diagonalized, propagated in the diagonal basis, and taken back to the original basis. This corresponds to decomposing S(t) as where O is the orthogonal matrix diagonalizing H tot , akin to U in Eq. (10), and S diag (t) propagates the covariance matrix of non-interacting oscillators. Let us define the following diagonal matrices: In the simulations the initial state of the network is assumed to be a stationary state of the network Hamiltonian H, i.e. a thermal state of some temperature T . The corresponding covariance matrix is diagonal in the basis of network normal modes and has elements The initial state of the probe is taken to be a squeezed state, i.e. a state where the second moment in one quadrature is lowered below that of vacuum, at the expense of increasing the second moment in the other quadrature. In the simulations we consider squeezing of the momentum. The corresponding covariance matrix is also diagonal, with elements

Appendix C. Frequency sweeps
In Sec. 4 the quantity ∆n = | a(t) † a(t) − a(0) † a(0) | is considered when frequency sweeps are made. The evaluation of a † a is covered in Appendix Appendix B, while here the protocol for frequency sweeps is explained. A single value for ∆n is acquired as follows.
(i) the value for probe frequency ω S is chosen (ii) the probe is prepared to the initial state with excitations a(0) † a(0) (iii) the interaction Hamiltonian H I is switched on for a time t (iv) the probe excitations are measured (v) the process is repeated as necessary to get an accurate value for a(t) † a(t) In the simulations the state of the network resets with the state of the probe, but in general the results can be expected to be similar as long as there is an energy difference between the probe and the network.
To perform a frequency sweep, the above protocol is repeated for many different values of ω S , resulting in a list of values for ∆n. In Sec. 4 both single and multiple frequency sweeps are considered. In the latter case a frequency sweep is repeated for many different realizations of the interaction Hamiltonian H I , resulting in an array of values for ∆n. This array is reduced to a list of values by averaging over different H I . When choosing network nodes to couple the probe to, each node is selected with equal probability. The selections are independent, consequently overlap between sets for different sweeps may happen.

Appendix D. Normal mode frequency probing
Here it is explained how the results shown in Fig. 5 were made.
Details about the determination of the open system dynamics and the value of the quantity a † (t)a(t) are given in Appendix Appendix B while the protocol for frequency sweeps is covered in Appendix Appendix C.
The considered network consists of identical quantum harmonic oscillators with a bare frequency of ω 0 = 0.25 interacting with uniform coupling strengths g = 0.1. The probe couples to the network in such a way that i k i = 0.05g while all non-vanishing elements k i have the same magnitude. At t = 0 the state of the probe is squeezed vacuum, defined by Eq. (B.6), with squeezing parameter r = 2.5, while the network is in vacuum, defined by Eq. (B.5) with T = 0.
For all shown results, the interaction time in ∆ n = | a † (t)a(t) − a † (0)a(0) | is t = 40000ω s , while 500 equidistant values in the closed interval [0.9ω 0 , 1.1ω N ] are used for ω s .
Here ω N is the largest normal mode frequency. As explained in Appendix Appendix C, the frequency sweeps result in 500 values for ∆ n .
To find the normal mode frequencies from the values of ∆ n , we employ the following algorithm. We perform a Gaussian blurring on the data up to scale σ = 0.55. Out of all of the local maxima of the blurred data, we select any that have a minimum sharpness of 1 and a value at least 0.1. In other words, the data is first convolved with a Gaussian kernel of standard deviation σ. The surviving maxima are chosen if they have a negative second derivative greater in magnitude than 1 and a value at least 0.1. The second derivative is estimated by padding the data by a single repetition of both the first and last element and convolving it with the kernel {1, −2, 1}. The role of blurring is to smooth out weak peaks, and sharpness and minimum value further restrict which surviving maxima are selected.
The spectral dimension is estimated from the found normal mode frequencies exactly the same way as in Sec. 3.