Heat flux localization and abnormal size effect induced by multi-body vibration in complex networks

Heat conduction in real physical networks such as nanotube/nanowire networks has been attracting more and more attention, but its theoretical understanding is far behind. To open a way to this problem, we present a multi-body vibration model to study heat conduction in complex networks, where nodes’ degrees satisfy a random distribution, and links consist of 1D atom chains with nonlinear springs. Based on this model, we find two interesting phenomenons: (1) The main heat fluxes of a network always localize in a skeleton subnetwork, which may have potential applications in thermal management and thermal concentrators, and (2) there exists an abnormal size effect of heat conduction in complex networks, i.e., the total heat flux of a network will enlarge with the increase of atoms on links, which is in contrast to the previous result on a 1D chain. Furthermore, we introduce a transmission diagram to characterize the skeleton of localized heat fluxes and then discover a jumping transition of total heat flux in the process of removing links, implying that the control of heat flux can be effective only when the change in a network topology focuses on the links within the skeleton. A brief theory is introduced to explain the abnormal size effect.


Introduction
Complex networks are ubiquitous in both natural and artificial systems and have been well studied in the last two decades [1][2][3]. It is revealed that most of the complex networks have the features of a small world and heterogeneity, which are the basis for tremendous applications such as in the aspects of communication on the Internet, synchronization in the brain, and epidemic spreading in social networks [4][5][6][7][8]. However, little attention has been paid to the aspect of heat conduction in complex networks, except for a few works [9][10][11][12][13][14][15]. In fact, this problem is becoming more and more important because the fast development of massive integration methods boost the interconnected intact structure of the nanonetwork extending from microscale to macroscale, and provide practical applicability in device designs [16,17]. On the other hand, more and more experiments on nanonetworks have demonstrated remarkable properties for many applications, such as nanofiltration [18], energy storage [19], electronic skin [20], sensors [21], neuromorphic net-works [22], thermal management [23], optical elements [24], and electrocatalysis [25]. A characteristic feature of these real networks is that their links can be considered as 1D atoms chains [13], implying that the knowledge obtained in 1D and 2D lattices [26][27][28] will be helpful for us to study heat conduction in nanonetworks.
A distinctive difference between complex networks and 1D/2D lattices is that the degrees of nodes are distributed in a range for the former but the same for the latter. This difference will make their heat conduction be fundamentally different. For example, it was shown by an oversimplified network model that the network topology can seriously influence its heat conduction and result in a rectification effect [9], positive degree correlation (assortativity) can enhance thermal transport, but negative degree correlation (disassortativity) can weaken it [12], interface thermal resistance exists at all the nodes of complex networks and their values depend sensitively on the sites of nodes, and the coupling strengths of connected links, which make the theoretical analysis of heat conduction in complex networks be a challenging problem [11,29]. In addition, it was shown that the heat flux can go from a lowtemperature node to a higher-temperature node, and there exists an optimal network structure that displays small thermal conductance but large electrical conductance [14].
However, realistic nanotube/nanowire networks are much more complicated than the oversimplified network model. Take the links as an example, the former can be considered as 1D atom chains [13], while the latter are only simplified as springs. Then, an exciting but challenging problem is how the interaction of different 1D atoms chains at each individual node influences the heat conduction on the whole network. To solve this problem, we here present a multi-body vibration model of heat conduction where nodes' degrees satisfy a random distribution, and links consist of 1D atom chains with nonlinear springs, as shown in Fig. 1.
The remainder of the paper is organized as follows. In Sect. 2, we present a random network to represent the nanotube and nanowire networks, where links are 1D atom chains and their motion satisfies the canonical equations. In Sect. 3, we do a lot of numerical simulations and unfold some fascinating phenomenons including the localization of heat fluxes, a jumping transition of total heat flux in the process of removing links, and skeleton subnetwork of heat conduction in net-works. Moreover, an abnormal size effect of heat conduction in physical networks is revealed, i.e., the total heat flux increases with the increase of links length, in contrast to the case of 1D chain. In Sect. 4, a brief theory analysis referring to interface thermal resistances is introduced to explain the abnormal size effect. Conclusion and discussion are presented in Sect. 5.

The model
The multi-body vibration model of heat conduction on physical networks is shown in Fig. 1, where the "green" points represent the nodes, and the "blue" points denote the atoms on links, the nodes with the "red" and "black" circles contact heat baths with high temperature T h and low temperature T , respectively, and n i j represents the number of atoms on the link i j connecting the two nodes i and j. It is clear that each atom on a link is connected only to its two nearest neighbors, while the atom at node i is connected to its k i nearest neighbors. To study heat conduction of this model, we first construct a random network by the approach in Ref. [1] with size N and average degree k . Then, we put atoms on the network, represented by all the points in Fig. 1. We assume n i j to be a random number from a uniform distribution and satisfies n 1 ≤ n i j ≤ n 2 . We here set n 1 = 3 and n 2 = 11, if without specific illustration. Following Ref. [30], we make both the nodes and the atoms on links be the same FPU-β oscillator with Hamiltonian for the atoms on a link i j with amount n i j , they have a potential for the atoms at nodes, their potential satisfy where the "green" points represent the nodes, the "blue" points denote the atoms on links, the "red" and "yellow" points represent two heat source nodes contacting two thermostats with the high temperature T h (the red circle) and the low temperature T (the yellow circle), respectively, two shadowlike marks connected to i 0 and j 0 represent the fixed boundaries, and n i j represents the number of atoms on link i j between the two connected nodes i and j where x i represents the displacement from the equilibrium position of the ith atom, C i, j is the coupling strength between the node i and its neighboring atom j, C i,i+1 is the coupling strength between two neighboring atoms i and i + 1 on a link and equal to C i, j for the same link, which is 1.0 in this paper unless otherwise stated, k i is the number of links connecting the ith atom, the sum is for all the nearest neighboring atoms j of node i, and g 4 is the normalized nonlinear constant, which is 0.1 in this paper unless otherwise stated. In our simulations, fixed boundary conditions are used, which means two source nodes i 0 and j 0 are connected to two fixed boundaries, as shown in Fig. 1. The heat bath attached to the two source nodes is stochastic Langevin heat bath and the coupling strength determined by the friction coefficient γ in Langevin dynamics is 5 in this paper. This value is within the range of γ ∈ (1, 100) that is recommended by Chen et al. [31] so that a meaningful physics can be obtained. The motion of all nodes and atoms satisfies the canonical equations To discuss heat conduction on the constructed network, we randomly choose two nodes i 0 and j 0 as the source nodes to contact two heat baths with temperature T h and T , respectively (see Fig. 1). Therefore, they satisfy where Γ h,l is the Gaussian white noises with where k B is the Boltzmann constant, and we adopt the dimensionless unit by setting k B = 1.
As T h > T , there will be heat fluxes continuously pumped from node i 0 to node j 0 through other nodes and links in the network. Specifically, the heat fluxes will be firstly transmitted from node i 0 to other nodes through all the links of node i 0 , and then gradually merge and finally go to node j 0 through the links of node j 0 . It should be noted that each link is considered as a one-dimensional FPU-β chain. For simplicity, only vibration in the longitudinal direction is considered, the transverse vibration is neglected. The fourth-order Runge-Kutta method is used to integrate the equations of motion in the simulation; after integration of 10 6 dimensionless time units with a time step of 0.01, the network will reach a stationary state. The local temperature at each atom i can be defined as [26,27] and the local flux J i j on each link i j can be calculated by [26,[32][33][34]] where · · · is the time average.

Heat flux localization
In numerical simulations, we first randomly choose two nodes as the source nodes-i 0 and j 0 , respectively, and set T h = 1 and T = 0.1 in this work unless otherwise stated. Then, we generate a network with size N = 80, average degree k = 4, and total links m = 100, and calculate the heat fluxes J i j for all the links in this network by Eq. 11. Finally, we schematically illustrate transmission diagram of heat fluxes in Fig. 2a. For clarity, we only plot all the nodes and their links in In order to further understand the influence of different links on heat conduction, one of the richest tool is removing links [35]. Hence, we probe the influence of removing links in the network on the total heat flux J of the network, which is the sum of heat fluxes of all the paths from the high-temperature source node i 0 to its neighbors, i.e., J = Taking the network corresponding to Fig. 2a as an example, without any change of other parameters, we remove the links in the network one by one according to the order of its heat flux value from small to large, and present the dependence of total heat flux J on the amount of removing links m r , as shown in Fig. 3a. We can see clearly that the total heat flux J will gradually increase to a optimal value with the increase of m r at first, but surprisingly appear a sharp jump after m r increases to 60. According to the parameters of the generated network corresponding to Fig. 2a For visualization of the skeleton subnetwork, we schematically illustrate the transmission diagram of heat fluxes after graying the removed 60 links, i.e., the remaining structure corresponding to the point in a red circle in Fig. 3a, as shown in Fig. 2b. Moreover, we further illustrate the heat flux transmission diagram when m r = 84, i.e., the remaining structure corresponding to the point with a yellow circle in Fig. 3a, as shown in Fig. 2c. Figure 2c clearly shows that after the removed 84 links are grayed out, only one heat conduction path remains from the high-temperature source node in the red circle to the low-temperature source node in the black circle, which means m r = 84 is the maximum amount of removed links that the network can endure, and leads to the minimum total heat flux J in Fig. 3a. In addition, we are surprised to find that the skeleton subnetwork is closely related to the topological path between two source nodes i 0 and j 0 . To clarify this point, we reorder the transmission diagram of Fig. 2a into a hierarchical structure. Specifically, we first arrange the high-temperature source node i 0 at the top and then arrange the nodes of its nearest neighbor and the next nearest neighbor until all the nodes are arranged, as shown in Fig. 2d. The hierarchical structure clearly shows that the links with large heat flux are almost located on the short paths between two source nodes i 0 and j 0 , such as two dark green paths in Fig. 2d. In other words, the shortest topological paths between the two extremities i 0 and j 0 are important elements of the skeleton subnetwork.
In order to further understand the influence of skeleton subnetwork on network heat conduction, we regulate the coupling strength C i, j of Eq. 3 inside and outside the skeleton, respectively. Firstly, we maintain the coupling strength of the part outside the skeleton subnetwork as the constant value 1, only regulate the coupling strength C in i, j inside the skeleton subnetwork, and present the dependence of the total heat flux J in the network on the coupling strength C in i, j , as shown by the green stars in Fig. 3a. Then, on the contrary, we maintain the coupling strength of the part inside the skeleton subnetwork as the same value 1, only regulate the coupling strength C out i, j outside the skeleton subnetwork, and present the dependence of the total heat flux J in the network on the coupling strength C out i, j , as shown by the solid pink circles in Fig. 3a. We interestingly find that the total heat flux J will monotonically increase with the increase in the coupling strength C in i, j inside the skeleton subnetwork, while it will weakly decrease with the increase in the coupling strength C out i, j outside the skeleton subnetwork, indicating that the heat conduction in the original network can be approximately replaced by that in the skeleton subnetwork. Moreover, we have checked the results of other networks for a case of N = 50, k = 6, total links m = 150, and another case of N = 80, k = 4, total links m = 160, as shown in Fig. 3c, d,  that the control of heat flux will be effective only when the controlled links are chosen from the skeleton subnetwork. In this sense, maybe we can make the heat fluxes of networks be pumped into the skeleton from outside the skeleton by the way of regulating the coupling strength inside the skeleton subnetwork, implying potential applications in controlling heat diffusion and thermal concentrators [37,38].

Abnormal size effect
Now, we turn to the size effect of n i j . For this purpose, we consider a specific situation with n 1 = n 2 = n 0 + 2 and focus on the dependence of total heat flux J on n 0 , where the other parameters are the same as the network in Fig. 2a. Figure 4a shows the dependence of total heat flux J on n 0 , where we can see that J increases monotonically with n 0 , in contrast to the 1D case where J decreases with the increase in the system size [26][27][28], indicating an abnormal size effect. To figure out the answer, we randomly choose 5 links i j from the network corresponds to Fig. 4a, and calculate their thermal resistances R i j = |T i − T j |/J i j under different link length n i j , as shown in Fig. 4b. We see that thermal resistance R i j decrease monotonically with an increase of n 0 , which means the total thermal resistance of a network will be weaken by increasing atoms, thus making total heat flux J increase with n 0 . Further, we checked the results of other network topologies with  Fig. 5a, b. The results in Fig. 5a, b consistently show that the total heat flux of the network J increases monotonically with n 0 , which confirm that the abnormal size effect is a general phenomenon in physical networks. To note that, for each n 0 in Fig. 5a, b, the total heat flux J is an average on 50 realizations. In order to better understand this phenomenon, more explanations should be given to the effect of interface thermal resistance. For fixed temperature of Langevin heat bath, J is determined by the thermal resistance of links in a network. For each link i j , there is always interface thermal resistances at its two ends [29]. It should be emphasized that the simulation results in Fig. 4a are obtained for the case of a small value of normalized nonlinear constant g 4 = 0.1, which means the phonon transport in the links (1D atom chains) is approximate to ballistic transport; thus, the total thermal resistance of the network is mainly dominated by the interface thermal resistance at nodes. Moreover, we added the theoretical result for f = 0.226 calculated by Eq. (13) in Sect. 4, which is equal to the total heat flux J of the network for the case of without interface thermal resistance, as shown by the red guide line in Fig. 4a. We can see clearly that the simulation results of total heat flux J increase monotonically with the increase of n 0 from 1 to 15, and gradually approaches the red guide line in Fig. 4a. Based on the above analysis, we deduce that the interface thermal resistance can be weakened by increasing atoms, thus causing the decrease of R i j and the abnormal size effect.

Theoretical analysis
The existence of interface thermal resistances makes the theoretical derivation of heat conduction on complex networks be a big challenge [9]. However, it is maybe possible to make a brief theoretical derivation for the influence of network topology on heat conduction when the interface thermal resistances is ignored. In this situation, we ignore all the information on the detailed dynamics of atoms, but only concern how the temperatures and heat fluxes are distributed on a network based on the topological structure of a given network and two chosen source nodes. Thus, the physical network of Fig. 1 can be simplified as a matrix A = (a i j ) with a i j = 1 if i and j are connected and 0 otherwise. Specifically, we let a ii = 0. Making a transformation of A, we have L = I − D −1 A, with I being the identity matrix and D being the diagonal degree matrix. We let M be the temperature vector of N nodes. For convenience, we reorganize the network and let the source nodes with T h and T be the nodes 1 and 2, respectively, which results in M = (T h , T , T 3 , · · · , T N ) T . M can be divided into two blocks with block-1 M 1 = (T h , T ) T and block-2 Assuming the Fourier's law J = −κ∇T is valid to network, we have J(r) = −κ∇M(r). Further, we obtain ∇ · J(r) = −κ∇ 2 M(r).
The discrete Laplace operator L can be considered as an analog of −∇ 2 [36]. Substituting it into Eq. (12)  Substituting it into L 11 H 1 + L 12 H 2 = f ,we can obtain f. That is, both H 2 and f can be analytically solved. As L contains all the information of network topology, H 2 and f will be seriously influenced by the network.
To confirm the theoretical derivation, we take the network corresponds to Fig. 4a as an example. Firstly, we calculate the theoretical temperature distribution by Eq. (13), see the "red stars" in Fig. 6, where the same results are given twice in (a) and (b) are for clear comparison with the different simulation results. Then, we present the temperature distributions of simulation results for two specific cases of n 0 = 0 (without atoms on links) and n 0 = 7 (with 7 atoms on each link), respectively, as shown in Fig. 6, where "i" of the x-axis represents the serial number of the node, the "T i " of the y-axis represents the temperature of the corresponding node, the "black solid circles" in Fig. 6a are for the case of n 0 = 7, and "blue open circles" in Fig. 6b are for the case of n 0 = 0. Comparing the two simulation results with the theoretical result in Fig. 6a, b, we can see clearly that the case of n 0 = 7 is much closer to the theoretical result than that of n 0 = 0 not only by the distance but also by the correlation. It should be emphasized here that the theoretical result is obtained by Eq. (13) under the situation of ignoring the interface thermal resistance, but it exists in numerical simulation [9]. These results can confirm again that the added atoms on links can significantly reduce interface thermal resistances at nodes and thus lead to the abnormal size effect. This finding shows that the heat conduction in physical network is more complicated than our intuition, but the influence of network topology can still be traced.

Discussions and conclusions
In summary, we have presented a multi-body vibration model to study heat conduction in complex networks, and numerically find that their main heat fluxes are always localized in a subset for a variety of complex network topologies, indicating a localization phenomenon of heat conduction. To understand its mechanism, we introduce a transmission diagram to clearly show the distribution of heat fluxes on complex networks. Then, we investigate the different roles of the links in complex networks by removing them one by one, and surprisingly discovery a skeleton subnetwork. Furthermore, we find that there is a jumping transition where the total heat flux of a network has no noticeable change when those links out of the skeleton are gradually removed but an enormous change when the links in the skeleton are removed. More importantly, by regulating the coupling strength between atoms inside and outside the skeleton, we confirm that the part inside the skeleton can promote the heat conduction in the network, while the part outside the skeleton can weaken the heat conduction in the network. In this sense, maybe we can make the heat fluxes of networks be pumped into the skeleton from outside the skeleton by the way of regulating the parameters of links inside the skeleton subnetwork, thus forming a new type of thermal concentrator [37,38], implying potential applications in thermal management and thermal metamaterials. After that, we study the influence of link length and surprisingly find an abnormal size effect of heat conduction, where the total heat flux of network can be enlarged with the increase of link length, in contrast to the previous results on 1D chains. Moreover, a brief theory ignoring interface thermal resistances is introduced to explain the abnormal size effect. These findings may have potential applications in the aspects of optimizing integration methods for nanotubes/nanowires and thermal management in nanonetworks, etc.