Destructive influence of interlayer coupling on Heider balance in bilayer networks

We consider the problem of Heider balance in a link multiplex, i.e. a special multiplex where coupling exists only between corresponding links. Numerical simulations and analytical calculations demonstrate that the presence of such interlayer connections hinders the emergence of the Heider balance. The effect is especially pronounced when the interactions between layers are negative, similarly as in antiferromagnetically coupled spin layers. The larger is the network, the narrower is the region of coupling parameters where the Heider balance can exist. If the interlayer couplings are of opposite signs and are strong enough, then the link dynamics can be reduced to the system of weakly coupled harmonic oscillators. For large strongly-coupled networks and randomly chosen initial conditions the probability of attaining the Heider balance decreases with the network size N as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${2}^{-{N}^{2}/2}$$\end{document}2−N2/2. Our finding can explain a lack of the Heider balance in many social systems, where multilayer structures mediate social interactions.


Supplementary Note
Fixed points in a one-triad network The set of equations readsẋ = (1 − x 2 )(yz + β 1 u) (1) The task is to find the fixed points and to determine their stability. The stability condition is that the real parts of all eigenvalues of the Jacobian at the fixed point are strictly negative.
The Jacobian is In the simplest case β 1 = β 2 = 0 there is no coupling between the layers and the only stable points are when xyz > 0 and uvw > 0, what is just the condition of the Heider balance (HB). Also we have three lines of unstable fixed points (x, 0, 0), (0, y, 0) and (0, 0, z), and the same for (u, v, w).
This is a special case of the more general situation, where β 1 = 0 and β 2 = 0. This is an example of the master-slave relation, as the layer (x, y, z) does not depend on the layer (u, v, w), but the opposite dependence persists. The Jacobian has a block-diagonal structure, with three eigenvalues equal to −2xyz; here again, the only stable solution for this layer is HB. For u, v, w we get three equations of motion:u Looking for stable fixed point for the whole system, we can distinguish some sets of equivalent configurations, and there is no need to consider each of them separately; then we can assume that either x = y = z = +1, or x = y = −1, z = 1. Let us consider the first case. Either some of the variables u, v, w = ±1, or the product of some pairs of them is equal to −β 2 . Truly different options are as follows: * u, v, w = ±1. Then the Jacobian is diagonal, and the eigenvalues are −2u(vw +β 2 ), −2v(wu+ β 2 ), −2w(vu + β 2 ). For |β 2 | < 1 the condition of stability is HB again. For Then, β 2 < 0 and the related eigenvalues are −e, −e, 2e, where e = ± √ −β 2 (1 + β 2 ).
Below we can assume that all weights are different from zero. In each one out of six equations (1-6) for the fixed point, one of two factors is zero. In other words, the set of equations (1-6) can be written asẋ = F (x)G(x) etc., where F (x) = 1 − x 2 and, symbolically, G(x) = yz + β 1 u. Then, all cases to be considered can be written in the form F F GGF G, what means: Obviously, we do not need to distinguish between F GF GGG and GF F GGG, and so on.
This case is discussed in the main text. For completeness of this material recall that the Jacobian is diagonal, and once the coupling is small (|β 1 | < 1 and |β 2 | < 1), the condition of stability is HB for both layers: xyz > 0 and uvw > 0. Once β 1 > 1 and |β 2 | < 1, the only stable solution HB (if achieved) in the layer (u, v, w) forces HB also in the layer (x, y, z), as x = u, y = v and z = w there. If β 1 < −1 and |β 2 | < 1, the HB in (u, v, w) forces x = −u, y = −v and z = −w (no HB). For β 1 > 1 and β 2 > 1, the stability condition is only x = u, y = v and z = w, with HB possible but accidental. If β 1 > 1 and β 2 < −1 or the opposite, we expect oscillations of the weights. This effect is more clear when the equations (1-6) are linearised around zero; we geṫ x = β 1 u,u = β 2 x etc., i.e. three mutually independent harmonic oscillators. See following section for details.
The Jacobian is diagonal, except the last row. The eigenvalue related to w is zero, as uv + β 2 z = 0.
*F F F F GG. In four first rows of the Jacobian, only diagonal elements are non-zero. The eigenvalues related to v and w are ±u (1 − v 2 )(1 − w 2 ).
*F F GF F G. In 1-st, 2-nd, 4-th and 5-th rows of the Jacobian, only diagonal elements are non-zero. The eigenvalues related to z and w are ± β 1 β 2 (1 − z 2 )(1 − w 2 ). *F F GGF F . In 1-st, 2-nd, 5-th and 6-th rows of the Jacobian, only diagonal elements are non-zero. As xy + β 1 w = vw + β 2 x = 0, the eigenvalues related to z and u are equal to zero. *F F GGGF . In 1-st, 2-nd and 6-th rows of the Jacobian, only diagonal elements are non-zero. The three eigenvalues related to z, u, v are ±w (1 − u 2 )(1 − v 2 ) and zero. *F F GF GG. In 1-st, 2-nd, and 4-th rows of the Jacobian, only diagonal elements are non-zero. One eigenvalue related to z, v and w is zero; the remaining two are The Jacobian has a block-diagonal structure, and the eigenvalues related to its upper left part, which is diagonal, has been discussed above. The lower right block has all diagonal elements equal to zero; the sum of its eigenvalues is zero; hence some of them are nonnegative. *F F GGGG. Suppose that x = y = +1. As xy + β 1 w = 0, we get w = −1/β 1 . Further, as wu = −β 2 , we have u = β 1 β 2 ; similarly wv = −β 2 gives v = β 1 β 2 . From uv = −β 2 z we get z = −β 2 1 β 2 . As an effect, yz + β 1 u = 0, what means that the first row of the Jacobian contains only zeros. The same construction works if x = ±1 and y = ±1.
*GGGGGG. All the diagonal elements of the Jacobian are equal to zero. Then, its trace is zero, with the consequence as above.
We have checked that the non-generic fixed points are also unstable. To summarize, the only stable fixed points are those where all weights are ±1.

Derivation of oscillator equation
We assume that the coupling coefficients are β 1 > 1 and β 2 < −1. It means that β 1 β 2 < −1. Assuming, β 1 , |β 2 | 1, then model equations (equation 2 in the main text) can be approximated by a system where each link is directly dependent only on its replica from the other layer: Therefore, a system comprises pairs of coupled variables. Each pair can be transformed into a one-variable derivative equation of second order: Further, if β 1 ≈ |β 2 |, then (17) becomes the non-linear oscillator: On the other hand, if it is assumed that |x 1 ij |, |x 2 ij | 1, then the pair of equations from 16 is simplified to: which transformed to the second order equation of one variable gives equation (8) from the main text.

Details of numerical calculations
Numerical implementation of a model presented in the main text in (1) needs one adjustment: x α ij ∈ [−1 + , 1 − ], where = 10 −5 . The reason to introduce the parameter is to evade the situation when due to finite computation precision the derivativeẋ α ij becomes exactly zero at the fixed point x α ij = ±1. Such a link would be unable to change its weight ever again regardless of the other changes in its "neighbourhood". With the parameter the links' weights that are very close to but not at ±1 still can change their value.
Numerical analysis of a dynamical system requires the proper halting conditions which are directly connected to the issue of the solution stability. The solution stability can be analysed analytically as it is shown in previous subsection. Such an approach describes an expected outcome of the numerical solution, but does not give information about the time needed to reach the solution. This is a nuisance in numerical analysis and that is why the halting conditions are necessary. The conditions used in this paper consider two kinds of system dynamics. The system may achieve a fixed point. A solution is believed to be numerically stable when for each pair of nodes for each layer, two following inequalities are simultaneously fulfilled: The first condition implies that an edge weight is close to the possible stable value (±1) and the second inequality ensures that the weight is changing towards this value. There is also a possibility that the system will never achieve such a point. Instead for example, the trajectory can be trapped around a marginally stable fixed point or even a saddle. The application of the parameter ε can make the trapping time shorter, but it is not a universal cure. Hence, applying a rather long numerical time t max of dynamics is necessary. Summing up, in this paper a single run of calculations is halted when the inequalities (20) are fulfilled or when the numerical time of the variables [x α ij (t)] exceeds the arbitrary chosen value of t max = 1 000.

Heider balance probability and interlayer link order in symmetric systems
Figures 3a and 3b in the main text show the relationship of the probability P HB versus mean interlayer link order S for networks with symmetric coupling β for different network sizes N . Presented in Supplementary Fig. S1, P HB vs S provides no new information but only better perspective of data corresponding to P HB S (β) . For all network sizes the probability P HB ≈ 1 when | S | 1 but the range of high P HB is not symmetrical with respect to S = 0. This range is skew to the left and is narrower for larger N , which corresponds to a narrow peak observed in Fig.  3a for larger N . Moreover, the larger is the system the smaller is the range when P HB > 0, and the sooner P HB decays as the function of | S |. For N = 3 there are such values of β that the probability P HB = 1 and the complete ferromagnetic ILO ( S = 1) are observed simultaneously.

Relaxation time of system dynamics
Relaxation time of the whole system is the time needed to achieve a stationary or a balanced solution. It is expected to be inversely proportional to the value of the r.h.s. of (2) from the main paper; yet the latter changes in time. We have analysed relaxation time for networks networks with symmetric (i.e. β 1 = β 2 ) and antisymmetric coupling (i.e. β 1 = −β 2 ). The results are presented in Supplementary Fig. S2. This time is the longest if oscillations are possible or when the interlayer and intralayer terms in (2) from the main paper are of opposite signs and similar absolute values. Hence, in the case of systems with symmetric coupling this time increases with |β| for |β| 1, and decreases with |β| for |β| 1.

Distribution of group sizes in balanced networks
Nodes from each layer in the HB can be divided into two groups A and B. Internally links in those groups are +1 and all the links connecting these groups are −1.  3,5,9,15. If in any layer a stationary state with the HB was reached, we divided nodes in this layer into groups (two different divisions were possible to occur). In order to divide nodes into groups: a node n = 1 was always assumed to belong to the group A. We tested the hypothesis that nodes are divided into groups with equal probability. Then, the theoretical group A size distribution is related to the binomial distribution: where B(N, s) is a binomial distribution of having s successes in N trials with the probability 0.5. The above way of dividing nodes skews the group A size distribution towards larger sizes, because the chosen node is more probable to belong to the larger group. Moreover, group A can't be empty, but it is possible to have all the nodes belonging to A (paradise). Using Pearson's chi-square test we estimated p-values for each N , for each layer and for each measured coupling coefficient pair when the stationary state with the HB is achieved for at least 20% of the system runs. Afterwards, using Holm-Bonferroni 1 correction with significance level of 0.05 only in one case (N = 3, β 1 = −1.5, β 2 = −1.1, α = 1) the hypothesis should be rejected.
We also grouped all the results over coupling coefficient pairs for different N and different layers ( Supplementary Fig. S3). The differences between the expected and real distributions are small and using Holm-Bonferroni correction with significance level of 0.05 the hypothesis that group A size distribution is given by equation (21) could not be rejected (Supplementary Table 1). Thus, the process of dividing nodes into groups is random and does not depend on the system parameters.