Local probe for connectivity and coupling strength in quantum complex networks

We develop a local probe to estimate the connectivity of complex quantum networks. Our results show how global properties of different classes of complex networks can be estimated – in quantitative manner with high accuracy – by coupling a probe to a single node of the network. Here, our interest is focused on probing the connectivity, i.e. the degree sequence, and the value of the coupling constant within the complex network. The scheme combines results on classical graph theory with the ability to develop quantum probes for networks of quantum harmonic oscillators. Whilst our results are proof-of-principle type, within the emerging field of quantum complex networks they may have potential applications for example to the efficient transfer of quantum information or energy or possibly to shed light on the connection between network structure and dynamics.

assume that the links are undirected, meaning that the interactions or relations modeled by the links are taken to be symmetric. Our results are general in the sense that the only assumption one must make about a physical network is that we know the number of nodes within the network and it is possible to perform measurements with results that are in a known relationship with the eigenvalues of the Laplace matrix of the corresponding graph. As an example of this type of system, we use a network of identical quantum harmonic oscillators interacting with spring-like couplings of constant magnitude 9,31 .
Earlier work for quantum networks has been done in the case of networks of spins, based on continuous-measurement-based approach of small networks up to 5 nodes with uniform or approximately uniform couplings 32 , as well as for quantum oscillator networks where the mutual information between a node and the rest of the network was shown to be characteristic of the topology when the network is at or near its ground state 33 . In contrast, our approach can in principle be applied to any kind of classical or quantum networks as long as the Laplace eigenvalues can be extracted. In practice, as producing the estimate requires solving a combinatorial problem, the method is best suited for networks of modest size. That being said, depending on the particular instance of a network, extracting all or nearly all eigenvalues might prove to be the most limiting factor.
Our main result is that it is indeed possible to obtain accurate estimate for the degree sequence of the network by using only a single probe weakly coupled to any of the nodes of the complex network. This result is based on exploiting known mathematical relations between the Laplace eigenvalues and the connectivity, and using the possibilities that quantum probing provides. The numerical evidence shows that the scheme works very well for different classes of network structures and is robust to small errors in the probed quantities. We also consider the case where the coupling strength in an oscillator network is uniform but a priori unknown. It turns out that for some classes of networks the coupling strength can always be correctly deduced, and numerical evidence suggests that the estimation succeeds with high probability in the general case.
For the sake of simplicity, we show first -in terms of classical graph theory -how the degree sequence of complex networks can be estimated once the eigenvalues of the Laplace matrix are known. After this, we turn our attention to quantum networks and develop a scheme to probe locally these eigenvalues within the network of quantum harmonic oscillators.

Results
Connectivity estimation. Once the nodes of a simple and connected graph have been labeled, its structure may be encoded into a matrix in many ways. In particular, the Laplace matrix L of the graph has elements where d i is the degree of node i and l ij = 1 if there is a link between nodes i and j and 0 otherwise; notice that l ij = l ji . Given the eigenvalues λ i of the Laplace matrix, the objective is to estimate the degrees d i . This can be done by combining several results from spectral graph theory, which studies the relationship between graphs and the eigenvalues of their matrices. In addition to bounds on minimum and maximum degree by eigenvalues λ i , the following relations must be fulfilled 34,35 The above restrictions are illustrated using a small example in Fig. 1. The sequences d′ of N positive integers that satisfy simultaneously the degree bounds and restrictions (2), (3) and (4) are referred to as solutions. The number of solutions is clearly always finite. It will be seen that the feasibility of finding them depends strongly on both the network topology and size.
For some graphs there is only one solution, and consequently the degree sequence d is unambigiously determined by the eigenvalues. These include simple chains, completely connected graphs and any graphs for which all degrees coincide. For the first two this follows from the squared sum in Eq. (3) attaining its minimimum or maximum value for a given number of nodes, while a direct application of Cauchy-Schwarz inequality shows that only regular graphs have equality in There is also an important class of graphs called threshold  (2), (3) and (4). The length of the line is given by Eq. (2). The number of boxes coincides with the sum of squared degrees given by Eq. (3), and the ticks are given by the bounds appearing in Eq. (4). By filling the line with squares such that all boxes are used, there are as many squares as there are slots, and each square will not exceed its slot, one will find an integer sequence satisfying simultaneously the three equations. Here this is done in three steps for a small graph to illustrate the problem. The bounds on mimimum and maximum degree are not shown.
SCiEntiFiC REPORTS | (2018) 8:13010 | DOI:10.1038/s41598-018-30863-2 graphs 36 that are uniquely determined by their degree sequence and their degree sequence is in turn determined by the Laplace eigenvalues, however the eigenvalues will typically be degenerate.
More generally, Eqs (2) and (3) for a given N fix the mean, variance and root mean square of the bounded solutions and Eq. (4) further refines them by ruling out cases where deviations from d are bunched together. For any solution, the deviations must cancel out because the correct sum of degrees is enforced; similarly also deviations between any element-wise squared solution and the element-wise squared d must cancel out. The degree sequence is always included in the solutions.
To find the solutions, we considered Eq. (3) as an integer partitioning problem, where the sum of squared elements must be partitioned into N integers. The allowed integers are square numbers with bounds determined from the eigenvalues. Taking the element-wise squareroots of each found partition and filtering the results according to equations (2) and (4) will provide the solutions. Alternatively, one could start from Eq. (2) and then filter but we found that this is more wasteful and consequently uses more memory and computation time. Further details about finding the solutions are provided in Methods.
We tested our estimation scheme on random graphs generated by Erdős-Rényi 37 , Barabási-Albert 38 and Watts-Strogatz 39 random graph models as well as tree graphs. An Erdős-Rényi random graph refers to either of two closely related models of generating random graphs with exactly N nodes. The so-called G(N, L) model chooses uniformly among all possible graphs with exactly L links, while the G(N, p) model starts from a completely connected graph and includes each link in the final graph with probability p. Both models are used in this work. A Barabási-Albert random graph G(N, K) is constructed starting from a cyclic graph of three nodes and iteratively adding a new node with K links until the graph has N nodes, connecting the new links randomly but favoring nodes with higher degree. It can be shown that graphs constructed like this have a degree distribution that follows a powerlaw. Watts-Strogatz graphs G(N, k, p) are constructed by starting from a circular graph where each node is connected to up to k-th nearest neighbors. Then each link is rewired with probability p, creating a graph with small world properties. Finally, a tree of N nodes is any connected graph with exactly N − 1 links; this gives them the property that they have no cycles, i.e. closed walks without repetitions of links or nodes other than the starting and ending node. We also considered two real complex networks: a social network of 62 bottlenose dolphins 40 and the macaque visual cortex network 41 .
As a figure of merit of a solution d′ we chose the  1 distance from d normalized by the total degree of the graph, i.e.
This choice is motivated by the fact that this quantity can be interpreted as the average deviation from the real degree per link. We found that, for all considered cases, f(d′) < 1/2. By choosing as final estimate the solution that has the smallest  1 distance from the mean of solutions it is possible to single out a solution particularly close to d, since the deviations, that must cancel out for any particular solution as explained previously, will then be partly averaged out. By mean of solutions, we indicate the sequence where each element is the corresponding mean degree calculated from all solutions. On the other hand, the set of all solutions always contains d while the estimate is typically not a perfect match.
The results, averaged over 1000 realizations of each random graph with size fixed to N = 30, are shown in Fig. 2. Besides the parameter values used here, we have also checked other values and found similar results for these random graph models. In particular, the estimation was not sensitive to size. After a transient for very small graphs of only a handful of nodes the averaged results converge, and there is no appreciative change as size is increased further.
For Erdős-Rényi random graph, we used L = 87. This would be the expected number of links for G(N, p) of same size with p = 1/5. For Barabási-Albert graph, we used K = 2. In the latter case the estimation performs worst, and in particular none of the estimates coincided with the real degree sequence. This is caused by the high variance of d for this class of random graphs: higher variance allows the solutions to deviate more from d and consequently the estimation is less accurate. Compared to the other two graphs which had typically thousands of solutions, Watts-Strogatz graphs and trees had much less solutions, with the former having tens and the latter only a handful with the used parameter values. Consequently a significant fraction of estimates were a perfect match with d. The plots are not smooth, indicating that certain values are much less likely than others, a feature not present for the other two graphs. For the former, we used k = 2 and p = 0.2. Unlike for the other graphs, more than half of the solutions had the same distance from d. We believe this to be because this class of random graphs had the smallest variance of d since they are generated from regular graphs. Trees had the biggest fraction of perfect matches out of all graphs, but this is mostly because the number of solutions was so small to begin with. This is essentially caused by any tree having the smallest possible number of links for a given number of nodes, greatly restricting also the solutions.
The insets in Fig. 2 show the scaling of the averaged execution time of the program that, given the eigenvalues, finds the solutions and the estimate(s). Each data point was averaged over 1000 realizations for the Erdős-Rényi and Barabási-Albert graphs, and over 10000 realizations for the other two. Different amounts of realizations were used partly because the results for the first two graphs converged significantly faster, and partly because the problem instances took much longer for them. For all four graphs the scaling would appear to be exponential in the size N. The other parameter values are as in the main figures except for Erdős-Rényi graph, for which G(N, p) model was used with p = 1/5. Refer to Methods for more details.
The real networks had considerably more solutions than the random graphs due to their bigger size, but for both a single, relatively accurate estimate was found as there was no averaging over realizations. The dolphin social network 40  While d′ close to mean solution are alike, the outliers are different from both them and d. This is because there are many relatively smooth sequences that satisfy the constraints, but only a few jagged ones that pass. Indeed, the estimation works poorly on graphs with jagged degree sequences since the majority of solutions will be much smoother. We stress that choosing an outlier and realizing it as a network will in general not yield the same solutions since degree bounds and restrictions imposed by Eq. (4) can change even between different realizations of a fixed d.
Application to quantum networks. To exploit the previous results for quantum probing and networks, we consider networks of uniformly coupled quantum harmonic oscillators 9 . We will use units as referred to an arbitrary (but fixed) frequency unit and give coupling strengths, times and temperatures in terms of this unit. We will also set ℏ = 1 and k B = 1. The network is composed of N unit mass quantum harmonic oscillators coupled by springs, each having the same bare frequency ω 0 . The couplings between network oscillators are assumed to be uniform with the strength given by g. We can express the network Hamiltonian in a compact way as where p = {p 1 , p 2 , …, p N } T and q = {q 1 , q 2 , …, q N } T are the vectors of momentum and position operators, I is the identity matrix and L is the Laplace matrix of the underlying graph. We will assume that g and N are known, but make no assumptions on L. Since the row sums of any Laplace matrix are zero, the eigenvalues λ i are non-negative. This, together with a positive coupling constant g, ensures the positivity of Hamiltonian H E . Since the network Hamiltonian is quadratic in position and momentum operators for any configuration given by L, it can be diagonalized with an orthogonal transformation. This allows us to move into an equivalent picture of noninteracting eigenmodes of the network. In this picture, , where P i and Q i are the position and momentum operators of the network eigenmodes and Ω i are their frequencies, related to the eigenvalues λ i of the Laplace matrix L as This is the key equation which allows us to use the previously described estimation procedure for the degree sequence. In other words, if we can probe the eigenfrequencies Ω i of the network, this gives us direct information about the eigenvalues of the Laplace matrix and therefore a way to estimate the connectivity of the network. It is also worth mentioning that, since ω 0 coincides with the smallest eigenfrequency, it is not necessary to know it beforehand.
Assuming that the network is in a thermal state of known temperature T, the detection of eigenfrequencies can be done by measuring the mean excitations 〈n(t)〉 of a bosonic probe weakly coupled to a node in the network and doing a frequency sweep across the range that covers the spectrum 9 . The probe is assumed to be a quantum , where p S and q S are its momentum and position operators and ω S is its frequency, while the interaction Hamiltonian is of the form H I = −kq S q j , where k is the strength of the coupling and q j is the position operator of the node interacting with the probe. By fixing the states of the probe and the network, the reduced dynamics of the probe can be determined exactly by diagonalizing the total Hamiltonian, solving the Heisenberg equations of motion for the decoupled oscillators, and returning to old operators. While here we fix the state of the probe and the network to be vacuum and thermal state of temperature T, respectively, the accuracy is largely insensitive to the state of the probe as long as there is an energy difference between the probe and the network 9 .
When coupled strongly to the network, the probe will exchange energy with all eigenmodes and the reduced dynamics depends on the structure of the network in a complex way. On the other hand, with a sufficiently weak coupling the interaction becomes limited to only the few closest modes in the vicinity of system frequency ω S , and this makes the reduced dynamics very sensitive to the resonance condition in the sense that when ω S matches an eigenfrequency, a significantly larger amount of energy can flow between the network and the probe before finite size effects cause the flow to be reversed. An example is shown in Fig. 3, which demonstrates that even a small difference in frequencies can lead to a very different value of 〈n(t)〉, for sufficiently long interaction times and a weak coupling, provided that there is an energy difference between the probe and the network. While this behaviour is universal to finite networks, the number of nodes N is assumed to be known in the probing protocol because otherwise one does not know when all eigenfrequencies have been found.
The probe must interact with an eigenmode to detect its frequency. The spectrum should also be nondegenerate because any degenerate eigenfrequencies are interpreted as a single frequency. This is typically the case, and it can be seen by considering the oscillators in terms of the eigenmodes: any q i can be expressed as a weighted sum of eigenmode position operators where the weights are given by the elements of the ith eigenvector of the matrix g I L ( ) /2 0 2 ω + . For a generic L, all eigenvalues are distinct and the eigenvectors will not have zero elements, which means that the probe will interact with and resolve all eigenmodes from any node. It should be mentioned however that bipartite graphs, an important class of complex networks, do often have degenerate eigenvalues.
In the non-ideal case, there might be some errors in the values of eigenfrequencies or the coupling strengths might be only approximately uniform. We checked the robustness against both for all four classes of networks. For all of them, 1% unbiased error in either eigenfrequencies or coupling strengths will typically not cause any errors in the detected sum of degrees while perturbing the probed sum of squared degrees, degree bounds, and bounds on partial sums very little if at all. With larger errors, the worst case accuracy of results averaged over many realizations deteoriates slowly, but the differences between individual realizations grows. We also found that the number of solutions had a large impact on the robustness of the best case accuracy, as this was affected very quickly for trees and Watts-Strogatz networks while the other two classes of networks were much more resilient. Sometimes the affected bounds on partial sums did not provide any solutions at all for trees or Watts-Strogatz networks, in which case we considered the accuracy of solutions without this restriction.
In the case of nonuniform coupling strengths, the eigenvalues of a weighted Laplace matrix L can be recovered from i 2 0 2 ω Ω − . Now the off-diagonal elements of L are the coupling strengths between the oscillators and the diagonal has the sums of coupling strengths to each oscillator. While other restrictions still apply as before, the eigenfrequencies only upper bound the sum of the squares of diagonal elements of L and conversely, their variance can only be bounded from above, reducing the accuracy of the estimation considerably. The number of possible solutions can still be finite if the coupling strengths in the network are divisible by the same number, for instance if there is a weakest coupling and others are its integer multiples.
Estimation of an unknown coupling constant. If the coupling strengths are known to be uniform but the value of the coupling constant is not known, one can estimate it from the probed eigenfrequencies using the relation λ ω = Ω − g i i 2 0 2 obtained from Eq. (7). The estimation procedure uses general properties of the eigenvalues λ i of an unweighted connected graph. We stress that the success or failure of the estimation depends only on Figure 3. Evolution of mean excitations 〈n(t)〉 for a probe system weakly coupled to a single node in an Erdős-Rényi network of 40 oscillators and 80 couplings with bare frequency ω 0 = 0.2 and coupling constant g = 0.1, with probe frequency coinciding with an eigenfrequency of the network (squares) and just a little above 1% off (circles). The initial state of the probe and the chain were vacuum and thermal state with T = 0.3, respectively, while coupling strength between the probe and the chain was k = 0.0025. The clear difference in the dynamics for longer interaction times makes the detection of eigenfrequencies possible. Once detected, the eigenfrequencies can be used to determine the Laplace eigenvalues.
SCiEntiFiC REPORTS | (2018) 8:13010 | DOI:10.1038/s41598-018-30863-2 the structure of the graph, rather than on a particular value of g. As will be seen below, for generic degree sequences it succeeds.
Because the graph is connected and simple, we know that are even integers, as they must be for a connected graph. This set can be further refined by using results related to regular graphs and the largest eigenvalue λ N . As mentioned before, for with equality iff the graph is regular. This property can be violated for values of the coupling constant larger than g, which can be used to rule them out. On the other hand, values smaller than g can violate the property λ N ≤ N 42 . Typically several values pass these tests, however as we will argue below, they are not equally likely to be correct.
Clearly, if some g′ satisfies the condition that both ∑ d are even, then so does any g′/x where x = 2, 3, 4, …. This suggests that g is more likely to be among the larger values satisfying the constraints. In fact, for trees and regular networks, the largest possible value coincides with g. In the former case this follows directly from the fact that the sum of degrees attains its minimum value, and hence any g′ > g will violate the assumption that the network is connected. In the latter case this can be seen by letting g′ = ag and noticing that then for all a > 1, where Δ is the constant degree of the network. More generally, for some g′ > g to lead to even sum of degrees and squared degrees, it has to be the case that i is the wrong sum of degrees corresponding to g′. While such a g′ might still be ruled out by the other constraints, this implies that without prior knowledge of the structure of the network, g can be determined unambigiously only when no other value passes the tests. We studied how well the estimation works in the case of the Erdős-Rényi random network as a function of connection probability p, as shown in Fig. 4 -here the G(N, p) model is used since for prime values of the total degree the estimation almost always succeeds, and consequently the G(N, L) model would lead to a discontinuous plot. The results confirm that the largest value coincides with g with high probability, success rate improving for larger values of p. Also shown is the fraction of conclusive cases, i.e. when g is the only possible value. The curve shows an interesting behaviour, with a sudden transient from most cases being inconclusive to most being conclusive, between p = 0.2 and p = 0.4. This is essentially because then λ N > N/2, which will rule out any g′ ≤ g/2. This does not guarantee conclusiveness since some g > g′ > g/2 might still pass, but this requires special values of ∑ d

Discussion
Connectivity is an important structural property of complex networks. We considered simple connected graphs and showed how connectivity can be estimated from the eigenvalues of the Laplace matrix. Our estimation scheme is applicable to any network, quantum or classical, amenable to the extraction of Laplace eigenvalues from measurement results. While the accuracy is best for networks with a degree sequence having small variance, the estimation performs well also for, e.g., networks where the degrees follow a powerlaw. Even if the size of the network prevents finding the estimate in a reasonable amount of time, the connectivity can still be classified as the mean and variance as well as bounds are easily determined from the eigenvalues. We stress that in actual situations the main limiting factor might instead be the extraction of the eigenvalues, as this becomes increasingly difficult for large networks due to the growing density of the eigenvalues. For actual networks of uniform link strength, the eigenvalues are in principle extractable locally whenever a basis of non-interacting normal modes exist. Hence in these cases the probing scheme requires only access to any single node, and once in weak contact, only local operations on the probe are needed. While this is an advantage in any situation, it is also the only option should the access to the network be restricted to only one node. As an example, we applied our results to networks of identical uniformly coupled quantum harmonic oscillators and showed how not only the connectivity but also the uniform coupling strength can be estimated with local probing of any of the oscillators in the network.
In this work, we have demonstrated how even in the quantum case, graph theory can be highly useful in eludicating the properties of coupled many-body systems. While here we used information extractable from a quantum network with minimal access, it would be interesting to study the case where a small subset of nodes could be accessed, or investigate how knowing also some of the eigenvectors of the graph could be exploited.

Methods
Finding the solutions. For each problem instance, the integer sequences bounded by d min and d max and satisfying equations (2) and (3) were found with the algorithm outlined below.
The output was further filtered according to eq. (4). The function IntegerPartitions is the one from Mathematica (Wolfram Research, Champaign, IL) version 7 and onwards, and it implements a multiply restricted integer partitioning algorithm that in this case finds all ordered sequences of exactly N integers adding up to A using only the elements appearing in the list s defined in the lines 1-4 above. The vast majority of time is taken by the call to IntegerPartitions. For the two real networks which had millions of solutions, IntegerPartitions was run in block mode where a subset of partitions was found with each call until all partititions had been found. While it is conceivable that more direct and efficient methods might exist to find the solutions, such discussion is outside the scope of this work.
Time complexity analysis. The data for time complexity analysis was gathered on a 64-bit PC running Microsoft Windows 7 Ultimate 6.1.7601 Service Pack 1 Build 7601. The processor was Intel Xeon CPU e5-1650 v3 @ 3.50 GHz, 3501 Mhz, with 6 cores and 12 logical processors, and the installed physical memory (RAM) was 48 GB.
For each data point, 1000 problem instances were generated and solved with Mathematica 11.2.0 using the method described above. The CPU time spent in the Wolfram Language kernel was recorded starting from having the graph eigenvalues to finding all solutions and choosing the estimate. In particular, this does not include the generation of the random graphs or finding the eigendecomposition of the graph Laplace matrix.