A unified model for linear responses of physical networks

José M. Ortiz-Tavárez,1 William Stephenson,1 and Xiaoming Mao1, 2

1
Department of Physics, University of Michigan, Ann Arbor, MI 48109-1040, USA
2
Center for Complex Particle Systems (COMPASS), University of Michigan, Ann Arbor, USA
Many physical systems—from mechanical lattices and electrical circuits to biological tissues and
architected metamaterials—can be understood as networks transmitting physical quantities. We
present a unified mathematical framework for describing linear responses of such physical networks
using tools from algebraic graph theory. This approach captures static and dynamic behaviors across
multiple domains, including mechanical, electrical, thermal, and diffusive responses using node and
edge variables (e.g., potentials, flows). Our formalism connects multiscale and multi-domain re-
sponses to the underlying network structure. We demonstrate how this framework enables efficient,
arXiv:2508.04616v1 [[Link]] 6 Aug 2025

generalizable solutions to a wide class of linear response problems, including stress propagation,
charge transport, and wave dynamics, and provide insights into network duality and entropy pro-
duction.

I. INTRODUCTION In this paper, we provide a mathematical framework to

describe linear responses of physical networks, which can
serve as a unifying tool to accelerate discoveries in this
All materials are fundamentally discrete at the atomic
new front. Discrete network-based models arise in many
scale. As one zooms out to larger scales, microscopic het-
contexts in physics and engineering, from Maxwell’s me-
erogeneities often average out, and a homogenized con-
chanical frames and phonons in atomic lattices [24, 25],
tinuum theory emerges as a good description of physical
to electric circuits [26, 27] and finite element modeling of
properties [1–3]. For instance, in crystals, long-range or-
multiphysics models [28, 29], offering a rich foundation
der based on repeating unit cells enables elegant math-
of knowledge. Building on this, we extract universal ele-
ematical formulations, such as band theory. In glasses,
ments from these approaches and propose a streamlined
despite the lack of long-range order, statistical rotational
mathematical framework grounded in algebraic graph
and translational invariance beyond a few atomic lengths
theory [30, 31], capable of representing diverse physical
allows for the construction of effective continuum theories
phenomena, static and dynamic, in a unified language.
grounded in symmetry principles.
As summarized in Table I, our approach maps linear re-
However, many materials exhibit structural complex- sponses in various domains—mechanical, electrical, ther-
ity across multiple length scales and cannot be captured mal, diffusive—to potentials and flows on nodes and edges
by traditional homogenization alone [4–6]. A prime ex- of a network. This mapping, combined with the linear al-
ample is biological tissue, which incorporates hierarchical gebraic structure of graphs that relates nodes and edges,
organization from macromolecules and organelles to cells enables an efficient and extensible framework for ana-
and vascular networks, spanning nanometers to millime- lyzing physical responses, coupling across domains, and
ters [5]. Such multiscale materials are also ubiquitous connecting physical behavior to underlying network ar-
in engineering—fiber-reinforced composites, architected chitecture.
metamaterials, and hierarchical porous media—where This paper is organized as follows. In Sec. II we in-
structural features span from nanometers to centimeters troduce the basic formulation of how linear responses in
to achieve tailored mechanical, thermal, or acoustic prop- multiple domains (Table I) can be written in a unified
erties [4–6]. Modeling such multiscale architectures re- language based on potential-flow problem on networks.
mains a grand challenge, as conventional multiscale mod- In Sec. III we discuss static responses of mechanical net-
eling methods are often tailored to specific materials and works to showcase the utility of this method in efficiently
difficult to generalize across different material classes [7– calculating a diverse set of problems. In Sec. IV we dis-
10]. cuss a set of irreversible transport problems on networks,
In parallel, advances in network science offer a powerful the formulation of which runs in parallel with static prob-
alternative perspective [11, 12]. Can we characterize the lems, albeit with one time derivative. In Sec. V we discuss
architecture of multiscale materials statistically—using dynamics and wave propagation on networks.
concepts such as connectivity, community structure,
geodesic distance, centrality, and spectral properties—
to build efficient and generalizable theoretical models? II. PHYSICAL NETWORKS
Such an approach holds promise for capturing the essen-
tial features of complex materials where structural details A. Notation
vary widely but exhibit statistical regularities. Recent
developments in this direction have yielded exciting new First we define some conventions of notation that we
insights into systems ranging from granular matter and will use in this paper. Because both network indices and
gels to geological materials [13–23]. Cartesian indices are necessary for mechanical networks,
 2

Properties Mechanical Mechanical Electrical Thermal Diffusive

(static) (mobility)
Potential Position Ui Velocity Vi Electric potential Vj Temperature Chemical potential
Node incoming Force Fi Force Fi Incoming current Ij Incoming heat Incoming matter
current
Potential difference Extension enm Extension rate ėnm Voltage vnm Temperature Chemical potential
difference difference
Flow Tension tnm Tension tnm Current inm Heat current Matter current
Transported Momentum Momentum Charge Heat Matter
quantity
Conductivity Stiffness Mechanical Electrical Thermal Diffusivity
impedance admittance conductivity

TABLE I. Analogous potentials and flows on physical networks across different domains.

we adopt a unified convention of notation here. or diffusive flow, and mechanical tension, among others.
Upper indices identify columns while lower indices This is in contrast to non-physical networks, such as so-
identify rows. Repeated indices are summed over only cial networks, where connections capture relationships
when they appear as both lower and upper indices. Lower and are not necessarily concerned with the flow of any
case Latin indices label nodes. They can also label cycles quantity [12].
or faces. Greek indices label spatial directions x, y, z. Up- In this paper we will describe a large set of physical
per case Latin indices label irreversible processes. These problems in terms of potentials and flows. This gen-
are only used in the Sec. IV. eral approach was first developed in the context of con-
We often group indices. For example eij represents trol theory where flows are referred to as the “through”
a quantity on edge ij. Here ij plays the role of an edge variables and potential differences as the “across” vari-
index and is understood to run only over edges that actu- ables [26]. While this literature mainly focused on elec-
ally exist in the network. We adopt this notation of using tric and pipe flow networks as scalar problems and ex-
the indices on the two nodes an edge connects to label plored dynamics of these networks, other studies used
the edge so it naturally denotes a directionality from i to similar graph theory approaches for static mechanical
j, which is useful for the discussion of physical transport. networks as vector potential problems, exploring zero
Consider Unα in a two dimensional network of three modes, states of self stress, and duality [32]. Here we
nodes. The lower index n labels a row while the upper take our present language from previous work on mechan-
index α labels a column. Written explicitly as a matrix ics [32], but aim at providing a comprehensive approach
we have: that combines both fields, and discuss a unified frame to
 x y solve similar types of linear response problems, static and
U1 U1
y dynamic.
Unα ≡ U2x U2  . (1)
To introduce the basic formulation, let us for concrete-
U3x U3y
ness consider the example of transport on an electrical
If we instead write Unα , then nα labels a row: network (Fig. 1a). The flow of charge through edge nm
is the current inm . Since charge is a conserved quantity
U1x
 
the currents must satisfy the discrete continuity equation:
U1y 
U  Q̇j = Ij = Qnm
 
Unα ≡  2x  . (2) j inm + Sj , (3)
U2y 
U 
3x where Qj is the amount of charge at node j, Ij is the net
U3y current flowing into node j, inm is the current flowing
from node m to node n, and Sj is an external source term.
The “cut-set” matrix (see more discussions in Sec. II C):
B. Flows and potentials on networks
Qnm
j ≡ δjn − δjm , (4)
We use the term physical network to refer to structures
where there occurs a transmission or flow of some phys- adds the current on all edges incident on node j (Fig. 1b).
ical quantity through the edges, and said flow can lead (Note that that nm is a single index for an edge here. If
to accumulation on the nodes. Examples of this trans- nm is not an edge, it simply does not appear as a row of
mitted quantity or flow include electric current, thermal the matrix.) This equation is directly analogous to the
 3

continuity equation in continuous space: In the case of that element being a resistor, Eq. (8) is
simply Ohm’s Law, where the admittance is one over
ρ̇ = −∇ · J + s, (5) resistance ynm = −1/rnm . Note that flows and poten-
tial differences are directed quantities (vnm = −vmn ),
where s is a source term, J is the current density and ρ in contrast to the admittance or its analogues which are
the charge density. undirected (ynm = ymn ). For more on sign conventions
see the appendix (App. A). We can also consider compo-
(a) Sn nents on the nodes, this is equivalent to having an edge
from a node to the ground. The constitutive relation on
the nodes is then essentially identical:
i nm n
m Ij = Yj Vj . (9)

With Eqs. (3,6,9,8) we have a complete set of equations

for the electric network, characterizing the transmission
of charges on the network and the responses of the nodes
and edges. We can write them together as:

edge → node: Ij = Qnm

j inm + Sj
(b) node → edge: vnm = Cnmj
Vj + v0,nm
Potential Flow node: Ij = Yj Vj
Nodes Vn In edge: inm = ynm vnm . (10)
  We will use this electric network problem as our first
Potential Flow example to discuss basic algebraic graph theory tools for
Edges vnm inm physical networks.
The same formulation applies to parallel problems in
thermal and diffusive transport problems, as listed in Ta-
FIG. 1. Nodes and edges of a physical network. (a) Net ble I, which will be discussed in details in the following
flow to a node, illustrating to Eq. (3). (b) Potential and flow
sections.
feedback between nodes and edges, illustrating the effect of
the cut-set (Q) and incidence (C) matrices from Eqs. (10). Eqs. (10) form a closed set of equations to solve linear
response problems in physical networks. One can write
them together in terms of an “equation of motion on the
At the same time, we have electric potentials Vj on the
network”:
nodes leading to potential differences across edges vnm :
j
vnm = Cnm Vj + v0,nm , (6) Yj Vj = Dji Vi + Qnm
j ynm v0,nm + Sj , (11)

j
where the matrix Cnm , generally called the incidence where we have introduced the dynamical matrix
matrix of a graph, projects potentials from nodes to
edges (i.e., calculates the potential drops on edges as Dji ≡ Qnm i
j ynm Cnm . (12)
vnm = Vn − Vm ). v0,nm is an externally imposed po-
tential drop or electromotive force (EMF), for example, We consider Vj , potentials on the nodes, as our funda-
from a battery on the edge. The two matrices are related mental degrees of freedom, and this equation solves the
via responses of the network with respect to external drive
on nodes current sources Sj or on edges voltage sources
Q = CT . (7) v0,mn .
(Using our index notation Qnm j
= Cnm .) In this work we In the simple case of ymn = 1 for all edges, the dynam-
j
analyze linear physical networks meaning that the flow ical matrix reduces to the graph Laplacian:
through and edge is linearly related to a potential differ-
ence across the edge via a constitutive relation. Edges Dji → Lij , (13)
can generally be seen as elements, in an electrical con-
text this would be resistances capacitors or inductors. where the graph spectrum determines the dynamics of
The elements on an edge determine its admittance ynm . the physical problem. Here we present a more general
We then have the constitutive relation on the edges: scenario, where the nodes and edge admittance, Yj , ynm
exhibit their own time scales, and can lead to rich dy-
inm = ynm vnm . (8) namics not captured by the graph Laplacian [30, 31].
 4

C. Kirchhoff ’s laws and compatibility law, which is to say they add up to zero along cycles:
X
In this section we review Kirchhoff’s laws [33, 34] and 0= vnm ∀ cycles of G. (17)
their relation to basic concepts from graph theory, cycles nm∈cycle
and cut sets, and the fundamental theorem of linear al-
gebra [35]. This discussion is the basis for the universal This is Kirchhoff’s voltage law. Cycles make a linear
applicability of Kirchhoff’s laws beyond electrical circuits subspace which means any linear combination of cycles
and also gives a deeper and more expansive understand- is also a cycle, in the sense that it satisfies Eq. (17).
ing of them. One can always find a set of linearly independent cycles
We first start from the condition of current conserva- such that all cycles are linear combinations of them, a
tion, Kirchhoff’s current law, meaning that for all sub- complete cycle basis. A matrix such that its rows form
graphs the flow coming in equals the flow coming out. a cycle basis is called the cycle or circuit matrix B [32].
Note that this is assuming there are no sources and no In terms of this matrix Kirchhoff’s voltage law is written
accumulation on the nodes which in the context of Eq. (3) as:
means Si = 0, Ii = 0. To make this statement mathemat-
0 = Bfnm vnm . (18)
ically precise we need to define cut sets [30, 31]. A cut
set is a set of edges that if removed divide the graph into
The space of all vectors which are potential differences
two subgraphs now mutually disconnected. A flow inm
is then the null space of B, which is also the orthogo-
is conserved if the sum of the flows along every cut set is
nal complement of the cycle space. Since this space is
zero:
also the column space of the incidence matrix, one can
X show via the fundamental theorem of linear algebra [35]
0= inm . (14)
that the cycle space and the cut space are orthogonal
nm∈ cut set
complements. Therefore, these spaces are equivalent to
This is Kirchhoff’s current law, also called the flow law. the space of flows that satisfy the flow law and potential
Cut sets can be represented as vectors and form a linear differences respectively (Fig. 2).
subspace called the cut space. A basis for this space is
called a cut basis. If the flows satisfy the flow law for
all elements of a cut basis they satisfy it for any cut
set, since any cut set is a linear combination of the basis
vectors (shown in App. A). A matrix such that its rows
make a complete cut basis is called a cut-set matrix. A
convenient choice for the cut-set matrix is Q = C T , since
each cut set in the basis is simply the set of edges incident
on a node, as we used in Eq. 4. In terms of this matrix
Kirchhoff’s current law reads:

0 = Qnm
j inm . (15)

Here the label j formally labels a cut set but it also labels
a node since for our chosen basis each cut set is identi- FIG. 2. A Venn diagram representation of vectors on RN e ,
fied with one node. This equation implies that all flows where Ne is the number of edges. The space of vectors that
that satisfy the flow law are orthogonal to the cut space. satisfies the voltage law is the orthogonal complement of the
Furthermore, the space of flows that satisfy the flow law space of vectors that satisfy the current law and are therefore
is the orthogonal complement of the cut space [30, 36]. shown as sharing only one point, the zero vector. All the de-
scriptions within the circles are equivalent. The area outside
Given a potential Vi for each node we can compute
of both circles but in the box represent vectors that contain
potential differences vij = Vi − Vj across each edge. The components from both subspace.
incidence matrix previously introduced performs this op-
eration:
i
vnm = Cnm Vi . (16)
D. Vector flow problems
Note that this matrix has as many rows as there are edges
and as many columns as there are nodes. In general a Potentials and flows need not be one dimensional or
directed quantity on the edges that can be obtained as scalar quantities. An example of a case where these quan-
potential differences are called compatible. It is easy to tities are vectors is mechanics. For a spring network in d
show that for a given set of directed quantities vnm there dimensions the potentials can be the d dimensional po-
exist a potential such that this flows are obtained as the sitions or the displacements and the flows are generally
potential differences if and only if they satisfy the cycle the tensions of the edges.
 5

Following the notation we defined in Sec. II A, Latin in-

ℓα 37 7
4...
dices label network elements like nodes, edges and cycles, 3
and Greek indices label spatial directions. For example, ℓα 67
2 ℓα 34
potential Viα , which should be interpreted as the α com-
ponent of the potential on node i. Flows on the edges 4 Face 1
α 3
will be denoted as fij , which is the α component of flow 6
ℓ α 5 α
f to node i through edge ij. 45 ℓ 56
Kirchhoff’s laws discussed previously for scalar quanti- 1
ties apply to the vector quantities component by compo- 2
nent. First, the flow law, or Kirchhoff’s current law is the
ℬ f ij
force balance condition. Let tα mn represent the tensions, ℓαij
the force balance condition is then written as: 3,4 4,5 6,7 3,7
ℓx 34 ℓy 34
0 = Qmn α
j tmn . (19)
ℓx 45 ℓy 45
If we take the node positions as the potentials, then the 1 1 1 1 -1 0 0
geometry of the edges described by a vectors ℓα α
ij = Xi − ℓx 67 ℓy 67
α
Xj is a potential difference. The incidence matrix relates ℓx 37 ℓy 37
potentials to potential differences:
ℓα k α
ij = Cij Xk . (20) FIG. 3. A mechanical network with a closed loop of edges
Kirchhoff’s voltage law (potential differences add to zero has a non-trivial cycle basis, and therefore a non-trivial cycle
matrix (Bfij ). For a given face (defined by the cycle of edges
along cycles) correspond to:
that surround it), the elements of its row in the cycle matrix
0 = Bfij ℓα will be +1 (-1) in a column if that column corresponds to
ij . (21)
an edge that points in the direction of (opposite of) a loop
Here, the index f labels a cycle, as illustrated in Fig. (3) around the face, all other elements will be 0. The cycle matrix
along with a diagrammatic representation of Eq. (21). multiplied by the edge length matrix, for a configuration that
is embedded in the flat space (i.e., compatible), will result in
For ℓα α α
ij = Xi − Xj from a given geometry, this voltage an all-zero matrix, per Kirchhoff’s voltage law (Eq. (21)).
law is automatically satisfied. In Sec. II D 2 we discuss
the more general implications of this rule.
For two-dimensional (2D) planar graphs the set of in- 1. Mechanics in the linear regime
ner faces makes a complete cycle basis. If one chooses
this basis then f labels a face. This is particularly rel-
In the linear regime, we are concerned with small, for-
evant for the study of mechanical duality and graphic
mally infinitesimal, node displacements from some given
statics [32, 37, 38].
reference positions Xkα → Xkα +Ukα . Using this in Eq. (22)
So far the “vector problem” of mechanics seems a triv-
and expanding to first order in U , we obtain:
ial extension of scalar problems, such as electricity. The
β
crucial difference comes from the constitutive relation- tα α
ij ≃ −kij ℓ̂0,ij (Uiβ − Ujβ )ℓ̂0,ij , (23)
ship, Hook’s law, which is not linear. We can write
Hooke’s law as: where ℓ̂β0,ij (Uiβ − Ujβ ) = eij is the length extension of
q  edge ij (which is eij = ℓij −ℓ0,ij , to first order). Note that
β
tα = −k ℓ ℓ − ℓ ℓ̂α
ij , (22) ℓ̂α
0,ij is not the orientation of an edge’s rest length, but
ij ij ij,β ij 0,ij
the orientation of the unperturbed edge. To first order,
where kij is the spring stiffness, ℓ0,ij is the rest length ℓ̂α α
0,ij is equivalent to ℓ̂ij . For the linear regime, it is then
(distinct from ℓij , the realized length of an edge) and ℓ̂ is convenient to identify these displacements Uiβ , instead
the unit vector that points along the edge, from node j of the node positions, as the node potentials, while the
to node i. Note that the relationship between potential flows are still the edge tensions.
difference ℓ, and flow t is not linear in general. It reduces The relationship between node displacements and edge
kα
to linear exactly for one dimensional spring networks and elongations is given by the compatibility matrix Cij :
when the rest lengths are zero. The relationship is effec- kα
tively linear in the regime of small deformations, aligning eij = Cij Ukα , (24)
the tension to the direction of the edge at the same time. which is obtained from the incidence matrix via:
In the following subsections we will first present the basic
kα k α
concepts pertaining to this linear regime and the analogy Cij = Cij ℓ̂0,ij . (25)
to the scalar version of the problem, taking node displace- Forces on the nodes are related to tensions via the equi-
ments U as potentials, and then discuss the general case librium matrix:
of geometric compatibility without assuming the linear
regime. fkα = Qij
kα tij , (26)
 6

which is obtained from the cut-set matrix Q = C T via: ( a) (b) ( c)

Qij
kα = Qij ij
k ℓ̂0,α . (27)
Note that we have written the tensions as scalar quan-
tities. It is clear that Q = C T , a relationship inherited
from the cut-set and incidence matrices. This implies
then that the column space of C and the null space of Q
are orthogonal complements and correspond to compati- FIG. 4. Examples of a floppy mode, where no edge is stretched
ble elongations and force balanced tensions respectively, (a), an isogonal mode, where no hinge rotation occurs at any
a relation used in studies of kinematic and static deter- node (b), and a state of self stress, where edges carry tensions
minacy of mechanical frames [39–42]. that are force-balanced on the nodes (c).
The analogous condition to Kirchhoff’s current law
(Eq. (14)) is force balance. Dropping the indices we can
write it as: The dynamical matrix in index notation is:
kβ kβ
0 = Qt, (28) Dpα = Qij
pα kij Cij . (32)

which is equivalent to Eq. (19) using the fact that ten- To write this without the indices we define diagonal ma-
nm nm
sions are along the edges (Eq. (23)). The sets of tensions trix Kij = kij δij . Then we rewrite the previous ex-
that satisfy this equation are the states of self stress. See pression as:
Fig. (4) for an example.
The analogous condition to Kirchhoff’s voltage law D = QKC. (33)
Eq. (17) in this linear regime is:
Se = 0 , ∀ e ∈ col(C), (29) 2. Geometry, isometries and floppy modes

where S is a matrix such that its rows form a basis for

We have previously remarked that in the context of
null(Q), which are the states of self stress [40].
mechanics, that the potential describes only geometry.
Analogously to the scalar case Eq. (10) we have the
We now expand upon this, in a more general way that
four equations:
does not use the linear expansion (Eq. (23))
edge → node: Fαj = Qnm ext
jα tnm + Fjα A network embedded in space can generally be de-
kα scribed by specifying the graph G (the connectivity), and
node → edge: enm = Cnm Ukα + e0,nm
the position of the nodes X, which together form (G, X),
node: Fjα = Mj ∂t2 Ujα a realization of G. A realization is simply a way to draw
edge: tnm = −knm enm . (30) the graph in a given space with all edges as straight lines.
This description is common in the field of combinatoric
Here, forces map to currents and displacements map to rigidity [49]. It is clear that the geometry of the edges ℓα
ij
voltages. We see that net force on the nodes F takes specifies a realization up to a global translation (note that
the role of net current on the nodes I. Correspondingly ℓα
ij means all x, y components of each edge are specified).
the external force F ext is the current source term. Iden- Therefore the space of edge geometries ℓα ij that satisfy
tifying force as the flow implicitly identifies momentum Eq. (21) is equivalent to the space of realizations of G.
as the quantity that is “flowing,” analogously to charge In other words, each element of this space corresponds
in electric transport. The potential here describes only to a realization and vice versa, up to global translation.
geometry. The quantity e0,nm corresponds to a small This is a vector space which implies that we can add and
change in the rest length of an edge. The node consti- subtract realizations to each other, and that the result
tutive relation is simply Newton’s second law while the would still be a valid realization. We see here a very
edge constitutive relation is Hooke’s law. This type of basic example of how connectivity constrains geometry.
formulation has been central in the study of topological Now, consider we are given a graph and the direction
mechanics in Maxwell lattices and networks [24, 25, 43– of the edges (G, ℓ̂), but not the length. We will now
48]. introduce a new matrix, which we shall name the “closure
More generally, the last two equations can include more matrix,” B:
diverse relations, such as
ij
Bcα = Bcij ℓ̂ij
α. (34)
node: Fjα = Mj ∂t2 Ujα + γj ∂t Ujα + Kj Ujα ,
edge: tnm = −knm enm − ηnm ∂t enm , (31) The null space of the closure matrix is then:
where γj , Kj are the friction coefficient and the pinning ij
Bcα ℓij = 0, (35)
spring constant of the node relative to the substrate/-
matrix, and ηnm represent the “loss” of the edge as a (for all c and α) which gives all possible sets of lengths ℓij
viscoelastic spring. that along with (G, ℓ̂) yield a valid realization. In other
 7

words, all realizations of G such that its edges are parallel which leads to:
to the directions ℓ̂ij
α . The difference between any two
such realizations is also an element of this space. These Flows that satisfy current law of G = null(Q(G))
deformations are called isogonal modes, they change the ↕
length of the edges but not their orientation in space such Flows that satisfy voltage law of G̃ = null(B(G̃)). (41)
that they are angle preserving (Fig. 4, b). These modes
have been discussed in the context of the mechanics of A straightforward extension of this duality relation
epithelial tissues [50]. in the vector problem is the construction of Maxwell-
Although the isogonal modes defined above do not in- Cremona reciprocal diagrams [39, 51–56], where the dual
volve linear approximation, they are related to another graph is embedded in the same space, (G̃, X̃) with each
type of special deformations—zero modes—which are de- edge in the reciprocal diagram parallel to its correspond-
fined in the linear regime. Zero modes are node displace- ing edge in the original. Note that it is more common
ments that do not elongate any edges (Fig. 4, a). The to rotate the reciprocal diagram by 90◦ such as the cor-
space of zero modes is the null space of the compatibility responding edges are perpendicular (as shown in the ex-
matrix. In other words all zero modes correspond to sets ample Fig. 5), but it is more convenient to keep them
of node displacements that satisfy: parallel for this discussion.
The extension of Eq. (40) is then:
0 = CU. (36)
Q(G, X) = B(G̃, X̃), (42)
The term floppy mode is used to denote “non-trivial”
zero modes, the trivial ones being rigid translations and giving
rotations. Floppy modes, along with global rotations can
be characterized in terms of the rotation of the edges. For States of self stress of G = null(Q(G))
two dimensional networks we can obtain the perpendic- ↕
ular component of the elongation as:
Isogonal modes of G̃ = null(B(G̃)), (43)
β
e⊥ α α


ij = ℓ̂ij ϵβα Uj − Ui , (37)
which further equals non-translational zero modes of G in
where the antisymmetric tensor ϵ performs a 90 degree◦ the linear regime. This duality shows up in various con-
rotation. We can write this in matrix terms as: texts such as force network ensembles of granular mat-
ter [57–59] and topological floppy modes in disordered
e⊥ = C ⊥ U. (38) networks [38].

When described in this manner, floppy modes and ro-

tations can be obtained from the closure matrix. The
perpendicular component of the elongation corresponds
to a zero mode if and only if

0 = Be⊥ . (39)

Therefore, in two dimensions, there is a one-to-one

mapping between isogonal modes and zero modes (ex-
cluding translations). Note that the row space of B
(row(B)) in general is not the same as null(Q). The
two spaces are equivalent only in the one dimensional
case. Instead, the null(Q) is equivalent to the row space
of row(S) defined in the linear regime (Eq. (29)).

3. Dual graphs and reciprocal diagrams

It is well established in graph theory that cycles and

cut sets of a planar graph G map to cut sets and cycles FIG. 5. A graph with A-D, faces α-γ, and edges 1-6 has a
of its dual graph G̃, which is constructed by mapping dual graph with faces A-D, nodes α-γ, and edges still 1-6, but
faces, nodes, and edges of G to nodes, faces, and edges of now lying perpendicular to their original direction. Exterior
G̃ [30, 31]. Using the formulation we defined in Sec. II, edges in the original graph become attached to a pin in the
and take cycles in the B matrix to be faces of the graph, dual. The states of self-stress in a graph are mapped to floppy
this can be written as modes (as in this example) or overall rotations in the graph’s
dual.
Q(G) = B(G̃), (40)
 8

III. STATIC RESPONSES OF MECHANICAL A. Response to forces on nodes

NETWORKS
When a force F (a vector in the RN d space of node
In this section we study the linear response of mechan- degrees of freedom, where N, d are the number of nodes
ical networks to static forces and displacements. This and spatial dimension) is exerted on the nodes, the net-
implies that all of our discussion here is only valid in the work responds with node displacements which give rise
“linear regime” which means close to mechanical equilib- to edge tensions that balance the imposed force. We have
rium (Eq.(23)). We first discuss the case where this me- then the following problem. We must minimize |w| given
chanical equilibrium is stress-free (all edges are at their the constraint that it balances the external force on the
rest lengths), and then generalize to the stressed equi- nodes,
librium case in Sec. III F. This set of tools can be useful
min |w|2 subject to F = Q′ w. (47)
in computational studies of mechanical networks such as w
jamming and rigidity percolation [60–66] and metamate-
rials [67–69].
The methods described here are also applicable to
physical networks which are not mechanical. In a later
section we will discuss how this applies to networks of
irreversible flows near thermodynamic equilibrium.
Mechanical linear response problems are in general
elastic energy minimization problems where there is some
imposed quantity that acts as a constraint, displacement,
force or edge swelling to give some examples. We start
with a spring mass network in an unstressed configura-
tion. For convenience we will define the following matrix:
√
C′ = KC, (44)
√
where K is a diagonal matrix where each element of the FIG. 6. The minimization problem in Eq. (47) over w re-
diagonal is the square root of the corresponding spring quires F = Q′ w in order to balance the applied external force
stiffness (assumed to be all positive). We can define F. Therefore, any possible tension response for a given ex-
Q′ = C ′T . These matrices are analogous to the com- ternal force must equal [Q′ ]+ F when projected onto col(C ′ ).
patibility and equilibrium matrices (C, Q), just with the In this schematic, the valid domain over which we minimize
|w|2 is therefore any point collinear to the horizontal dashed
spring constants incorporated. We will describe the state
gray line. Shown in blue is an arbitrary hypothetical |w| that
of the√edges, both elongations and tensions, with a vector meets this condition.
w = Ke that is related to displacement and forces via:
Any w can be written as a sum of two components—
w = C ′ U, (45) its projection along the compatible space col(C ′ ) and the
self stress space null(Q′ ) (Fig. 6). Here we assume that
before applying the force the network is unstressed. The
F = Q′ w. (46) network can only realize tensions via node displacements,
wmin = C ′ U for some U. As a result, wmin ∈ col(C ′ ).
It should be clear that the total elastic energy is E = To see that this wmin is unique note that the difference
|w|2 /2. For this reason we will solve all the following between any two w that result in the same force on nodes
linear response problems by writing an expression for and has to be an element of null(Q′ ). This would imply that
then minimizing |w| given the constraints. We will do for any w that balances the external forces F its projec-
this by exploiting the fact that the space of compatible tion on to col(C ′ ) is wmin . To obtain this projection one
elongations col(C ′ ) and the space of self stresses null(Q′ ) can use the Moore-Penrose pseudo inverse [70], which we
are orthogonal complements. In what follows we will here denote by +:
often use the projection operator S ′ S ′T where S ′ is a +
matrix such that its columns make an orthonormal basis [Q′ ] F = wmin , (48)
for null(Q′ ) (similar to Eq. (29) but for Q′ ). At other as in Fig. (6). We can also solve for node displacements
times we will use SS T where S is a matrix so that its such that wmin = C ′ U by solving:
columns make an orthonormal basis for null(Q).
F = DU, (49)
It is worth noting that one could always solve the prob-
′ ′
lems via the equations of motion without using this min- where D is the dynamical matrix D = Q C = QKC. An
imization approach. Nevertheless we find that this min- example of this computation is shown in Fig. (7). This
imization approach is offers an elegant perspective that can also be done by using the pseudo inverse of D or by
deepens and unifies the understanding of the subject. the use of other numerical methods.
 9

FIG. 8. The decomposition into compatible and force bal-

anced components of an imposed edge swelling w0 , repre-
FIG. 7. Static forces (green dashed arrows) are applied to senting changes in the rest length.
a few nodes in the corners of a disordered triangular lattice
generated by randomly displacing nodes on a regular trian-
gular lattice by small amounts. The resulting response of the
system sees many edges extended (red) and some edges com-
pressed (blue) (thickness proportional to magnitude of ten-
sion) in order to for every node to be force-balanced.

B. Response to edge swelling

There are circumstances where one can impose a

change in the state of the edges by means other than node
displacement. One straight forward example is a change
in rest-length the springs via some active element, which
we call “swelling” for simplicity [71], but the actual cause
can be diverse, such as piezoelectric coupling, chemical
reactions, growth, etc.
In all of these cases we are producing some change
in the edges by accessing some degrees of freedom be- FIG. 9. The rest length of an edge in a disordered triangular
yond node displacements, which are normally available lattice is increased (black arrows). The system is incapable
for a mechanical network. We can describe the imposed of reducing potential energy to zero in response to this rest-
length change, so some edges carry tension (magnitude pro-
change as some w0 , the system will then try to release
portional to edge thickness) from being compressed (blue) or
the resulting stress via node displacements to minimize extended (red).
the energy. The final state of the system is described by

min |w|2 subject to w = w0 − C ′ U. (50) In terms of elongations and tensions this is:
U
√ √
This is a well known minimization problem in the context t = KS ′ S ′T Ke0 , (52)
of the method of least squares [35].
where e0 are the imposed elongations. An example of the
Note the vector w0 can be decomposed as shown in computation of tensions in response to edge swelling is
Fig. (8). Only the compatible component can be relaxed shown in Fig. (9). A set of displacements that minimize
via node displacements. Once all of this is relaxed only |w| can be obtained by use of the pseudo inverse:
the force balanced component remains. Therefore, the
minimal w can be obtained by projecting the imposed +
U = [C ′ ] w0 . (53)
change into null(Q′ ). As we defined, columns of S ′ form
an orthonormal basis for null(Q′ ), then minimal w is The set of displacements obtained via pseudo inverse will
given by be orthogonal to all zero modes.
It is worth noting how Eq. (52) looks when we don’t
w = S ′ S ′T w0 . (51) do the change of variables C → C ′ . Let S be a matrix
 10

such that its rows are an orthonormal basis for null(Q).

Note that for one dimensional mechanical networks, and
scalar flow problems null(Q) is just the cycle space. Then
one can show,
−1
t = S T SK −1 S T

Se0 . (54)

This is a more complicated expression than Eq. (52) but

it is also a more general expression since it does not de-
pend on K being a positive diagonal matrix. We will see
that analogous expressions are applicable to the case of
prestressed networks and dynamics. In all of these cases
the tension, or its analogue, only depends on the incom-
patible part of the applied edge elongations, or its ana-
logues, a consequence of Kirchhoff’s laws. Similar results
have showed up in different contexts in Ref. [24, 25, 72].
For one dimensional spring networks this S is an or-
thonormal basis for the cycle space.

C. Response to prescribed displacements on nodes

Consider a spring mass network where some agent se- FIG. 10. The response of a mechanical network to prescribed
lects a set of nodes—let us label the set as A and all displacements (nodes within shaded gray regions). Displace-
ments of all nodes shown by green arrows. The resulting ten-
other nodes as B—and produces displacements uA on
sions (magnitude proportional to edge thickness) are colored
them without directly interacting with any other nodes. red (blue) if extended (contracted).
The response of the network is obtained by solving the
minimization problem:

min |w|2 ′
subject to w = CA ′
UA + C B UB . (55) Once UB is known, we can solve for FA using the top row
UB of Eq. (58). An example of the computation of tensions
and displacements in response to applied displacements
Here we have split the C ′ matrix in two. CA
′
contains only is shown in Fig. (10). The solution for the force is unique,
the columns identified with nodes in A while the rest of but the solution for the displacements in B may not be
′
the columns are in CB . This problem is very similar in unique due to zero modes. These modes, by definition,
form to the edge swelling problem of the previous section do not alter the energy cost of a deformation. This type
in Eq. (50). We therefore solve it analogously: of calculation, where the dynamical matrix is separated
′ + ′
UB = −[CB ] C A UA , (56) into controlled and free nodes, has appeared in Ref. [19].

T ′
wmin = SB SB CA UA , (57) D. Condensing nodes

where SB is an orthonormal basis for null(QB ) and QB = There are a number of situations in which one would
T
CB . Therefore, the solution for wmin will produce no net wish to constrain a group of nodes to move with the
forces on the nodes of set B, which is an obvious physical same displacement, or in other words have the same po-
requirement for the static response of any system. tential. This can be useful for cases such as defining
We can alternatively solve the problem with the dy- coarse-grained degrees of freedom which represent a set
namical matrix. We first split said dynamical matrix of nodes in a network. Note that we call this proce-
into blocks for the subspaces A, B: dure “condensing” the nodes, but we are not shrinking
  
FA DAA DAB

UA
 the physical distances between them. We simply require
− = . (58) that they move together as a rigid body.
0 DBA DBB UB
The compatibility matrix has one column for each de-
Here, DBA is the submatrix made up of the rows that gree of freedom. Condensing a group of nodes into a
map to nodes in B and the columns that map to nodes super node supposes diminishing the degrees of freedom
in A, the other blocks are similarly defined. To solve and so we must change the compatibility matrix. The
this problem we first look for the displacements of B by super node will have d columns in the new compatibility
solving the bottom row, matrix. The column for super node displacement along
the x direction is the sum of the x direction columns of
−DBA UA = DBB UB . (59) each of the nodes condensed into the super nodes, rep-
 11

resenting a combined constraint on the super node from such that:

all the individual constraints in the original network. We
do the same procedure for all the spatial directions. Uiα = Γµiα Vµ . (63)
Note that the elements Vµ are Vx , Vy , Vθ .
Now we group together all the columns of the original
compatibility matrix that map to nodes in the super node
into a matrix we call CA , where the rest of the columns
we group into CB , the elongations are given by:
 

UB
e = CB , CA Γ . (64)
V
By grouping the two matrices CB , CA Γ together, we have
formed a new matrix that plays the role of the usual com-
patibility matrix from before, but that now acts on the
degrees of freedom of the rigid cluster instead of the de-
grees of freedom of its components. This “new” compat-
ibility matrix’s transpose, the new equilibrium matrix,
obeys a similar expected relationship:
   
FB QB
= t. (65)
FV ΓT QA
Note that FV has 3 components that correspond to the
FIG. 11. The response node displacements (arrows) and edge forces along each degree of freedom of the super node.
tensions (proportional to edge thickness, blue for compressed, In this case the rotational component would be the net
red for extended) of a mechanical network with a number of torque on the super node around X0 . The procedure
nodes condensed into a rigid body (nodes within the gray cir- described here seems quite generic, but this is simply
cle) under shear (imposed displacements within shaded rect- because the information regarding how nodes are con-
angles). As a rigid body, the displacements of the nodes
strained is encapsulated entirely within the matrix Γ.
within the gray circle must be a combination of an overall
x, y translation and θ rotation. Edges between two nodes in
With any matrix Γ that is suitable for a given system,
the rigid body cannot carry tension, as the distance between one can solve linear response problems of the previous
the two nodes they connect does not change, thus each node section on that system by simply using the new com-
within the rigid body is not necessarily force-balanced from patibility and equilibrium matrices in place of the “old”
the edges alone. However, the rigid body can “transmit” ten- compatibility and equilibrium matrices used previously.
sion across itself and must be force-balanced as a whole. An example of the computation of tensions and displace-
ments in response to applied displacements in a system
More generally one can constrain the relative motion with condensed nodes is shown in Fig. (11).
of a subset of nodes to be within a prescribed linear sub- One common scenario where nodes are condensed is to
space. One natural example in mechanics would be to make a ground node. In mechanics, the resulting network
make the group move together as a rigid body, allowing is called a pinned graph, and the grounded nodes are
rigid rotations. We will now describe the method for two called pins (since they are “pinned” to the ground). To
dimensions. Let the 3 degrees of freedom of a super node ground a group of nodes we follow the general procedure
be Vα , Vθ , where α goes over x, y directions and Vθ is the described above but make Γ = 0. The new compatibility
rotation in radians. The displacement of a node included matrix is simply CB .
in this super node is therefore:

Ui = V + Vθ ẑ × (Xi − X0 ) , (60) E. Response to global strains and stresses

where Xi is the position of the node and X0 is the center

Most of the preceding sections offer a microscopic per-
of rotation which can be arbitrarily chosen. In index
spective. By this we mean that the forces and displace-
notation, this is:
ments are specified at the level of individual nodes or
  edges. Materials science is most often concerned with
Uiα = Vα + Vθ ϵαβ Xiβ − X0β . (61) macroscopic quantities, stresses strain and elastic mod-
uli. In this section we describe the response of a network
Now we can define a matrix Γ that maps the degrees of to macroscopic stresses and strains and compute elas-
freedom of the rigid cluster (V) to the displacements of tic moduli in two approaches: a node-based approach
individual nodes (U): through the control of boundary nodes, and an edge-
  based approach where a macroscopic strain is treated as
Γµiα = δαµ + δθµ ϵαβ Xiβ − X0β , (62) analogous to changes in rest lengths.
 12

1. Node-based approach or “hard wall” boundary conditions but the edge-based

method can be applied in more scenarios such as periodic
Consider placing a spring network, 2D for this exam- boundary conditions where there is no obvious choice of
ple, inside a square box. We fix some of its nodes to the “boundary nodes.” Two examples, one of the node-based
boundary box, creating a super node with the method approach, one of the edge-based approach under hard
described in the previous section (Fig. 12, a). We then wall boundary conditions, are shown in Fig. (12). Note
allow this boundary to deform as prescribed by an affine that the tension response of these two examples is exactly
strain ϵ. In 2D the strain tensor has three independent equal.
components (ϵxx , ϵyy , ϵxy ) therefore we have reduced the The edge based approach is mathematically very
degrees of freedom of the boundary nodes to three. We similar to the edge swelling problem described earlier
can find a matrix Γ such that: Eq. (50), where the imposed elongations are:
µν
UA = Γϵ, (66) e0,ij = Cmacro,ij ϵµν . (71)

where we represent the strain as a 3 × 1 vector, known This is the same as the “affine extension” (eaff ) used in
as the Voigt notation [73]. Here uA is the displacements Refs. [24, 25]. It is worth noting that although mathe-
of the nodes in the boundary, where external forces are matically we treat this problem using the edge swelling
applied. As in the previous section, we arrive at equation: formulation, the physical meaning is different. In the
edge swelling problem, we consider the change of the
e = CB UB + CA Γϵ, (67) rest length of some edges due to internal mechanisms
in the edge, whereas here we consider externally im-
where CB and UB are the compatibility matrix and the posed macroscopic strain. Following the methods in
displacements of the nodes in the bulk of the material Refs. [24, 25], we first find the “affine extension” using
respectively. We then have a new compatibility matrix: Eq. (71), and then project it to the space of states of self


C = CB , CA Γ , (68) stress (obtained under appropriate boundary conditions),
as we discuss below.
which also gives us a new Q = C T and the following The matrix Cmacro is obtained from the relationship
relationships: between strains and elongations,
µν

UB
 Cmacro,ij = ℓij ℓ̂µij ℓ̂νij (72)
ϵ 
e = C  xx  , (69)
ϵyy The vectors ℓij , ℓ̂µij are respectively the magnitude and
ϵxy direction of all edges ℓα α α
ij = Xi − Xj . Here the pair of
Greek indices act as one index in Voigt notation. As
defined above the matrix has four columns (in 2D) one
FB
 
for every component of the strain tensor xx, yy, xy, yx.
σxx  Since the strain tensor is symmetric the xy and yx can be
σ  = Qt. (70)
yy replaced with a single column which is the addition of the
σxy two. The resulting Cmacro will then have three columns
and as many rows as there are edges. Keeping the four
Here, FB is the force by the network on the bulk nodes
columns introduces a zero mode which corresponds to a
and the σµν is the components of the stress tensor. We
global rotation. We can obtain an extended compatibility
can verify that these are indeed the components of the
matrix by concatenation Cext = [C, Cmacro ] such that
stress tensor by taking the derivative of the elastic energy
with respect to the strain. This now allows us to do lin-
U

ear response problems where we prescribe a macroscopic
ϵ 
strain or stress and find the edge extensions and node e = Cext  xx  . (73)
ϵyy
displacements in the bulk. These problems are solved
ϵxy
in the same way we solved the problems of prescribed
displacements (now prescribed strains) and prescribed Fig. (13) gives a diagramatic representation of Cext , CB ,
forces (now imposed stresses) (as shown in Fig. 12a). and their transposes.
For an imposed strain it does not matter if we used the
node method or the edge method, both problems reduce
2. Edge-based approach
to the edge swelling problem:

We have just described a node-based approach to pre- e = CB UB + e0 . (74)

scribing a global affine strain. Now we describe an edge-
based approach. The two methods are equivalent for For the edge method, e0 = Cmacro ϵ, while for the node
the purposes of calculating elastic moduli under “fixed” method, e0 = CA Γϵ. These two vectors are not the same
 13

(a)

FIG. 13. Venn diagram representation including all vectors

in RNe with the exception of the zero vector 0. We show here
the relationship between the relevant linear subspaces. Note
that the intersection corresponds to responses of the bulk to
forces or displacements of the boundary.

(b)

FIG. 12. Response of a mechanical network to a global strain.

(a) The node-based approach to finding response to a global
strain (here, global dilation), is accomplished by prescribing FIG. 14. A mechanical network subject to uniform edge rest
affine displacements to nodes on the boundary (within the length decreases under periodic boundary conditions (with
shaded region). Resulting displacements for the bulk (ar- the boundary denoted by dashed lines). Edge tensions (pro-
rows) and tensions (proportional to edge thickness) can be portional to edge thickness) show that almost all edges are
calculated using Eq. (59) with the new compatibility matrix extended (red). Node displacements are denoted by green ar-
from Eq. (68). (b) The edge based analog to global strain rows.
requires imposing (negative) edge elongations according to a
given affine strain (again, contraction; arrows on edges, from
Eq. (71), while pinning the boundary (squares, within shaded This is due to the fact that both operators move the
region). The resulting tensions (proportional to edge thick- nodes in the boundary in the same way, according to
ness, red for extension) are identical to the node based ap- strain ϵ, but CA Γϵ does not result in bulk displacements.
proach if both systems are under the same global strain. This means that on can go from one to the other via
displacements of the nodes in the bulk.
Note that the notion is “fixed boundary condition”
in general. When we are dealing with fixed boundary is associated with the edge swelling formulation we use
conditions, their difference is compatible: here. Translating to the problem of imposed strain, this
correspond to rigidly deforming the boundary, which is
C macro ϵ − CA Γϵ ∈ col(CB ). (75) the same as in the node-based approach. Another way,
 14

Note, from Eq. (70):

 
σxx
σyy  = QA t, (79)
σxy
and
√ √
S T KCA Γϵ = S T Ke0 . (80)
By substituting these expressions into Eq. (76), we ob-
tain a linear relationship between stress and strain as in
Eq. (78) and identify the stiffness tensor as:
√ √
K = ΓT QA KSS T KCA Γ. (81)

F. Prestressed networks

Previously, we assumed that the equilibrium config-

uration we find the system in is stress free. Here we
consider spring networks whose initial configuration is
stressed and in equilibrium. The solutions of linear re-
sponse problems will still minimize the elastic energy, but
FIG. 15. A mechanical network subject to pure shear (via
edge elongations) under periodic boundary conditions (with representing them as minimization problems in the same
the boundary denoted by dashed lines). Edge tensions (pro- way as before is not as straightforward. We will discuss
portional to edge thickness) show that some edges are com- this in the context of two dimensional networks. The
pressed (blue) and extended (red), resulting in the displace- procedure is easily generalized to networks of higher di-
ment of nodes (green arrows). This response if found via mensions by adding more ⊥ components in a similar way
Eq. (71), where the strain tensor is anti-diagonal with ϵ12 = the one that we will describe.
ϵ21 = −0.01. It has been shown in Ref. [72] that in two dimensions
the change in elastic energy due to a small node displace-
ments U around equilibrium is given by:
as we mentioned above, is to use periodic boundary con-

K∥ 0
  ∥
ditions across the box. In this case, the imposed elon- 1 C
gations are the same as given in Eq. (71), but the states δE = U Q∥ Q⊥ U, (82)
2 0 K⊥ C⊥
of self stress are slightly different (due to different CB ),
which will cause small differences in the response, but where Q∥ , C ∥ are the usual equilibrium and compatibility
these differences will in general vanish for large systems. matrix and Q⊥ , C ⊥ are the equilibrium and compatibility
Example of the computation of tensions and displace- matrix of the network if it where rotated by 90◦ . We
ments in response to global strains through edge elonga- previously defined C ⊥ in Sec. II. The K ∥ matrix is equal
tions are shown in Figs. (14, 15). to the usual diagonal matrix of spring stiffness (as in
The solution to the edge swelling problem tells us that, Eq. (33)), while
√ √
t = KSS T Ke0 , (76) ⊥nm tij nm
Kij =− δ , (83)
ℓ0,ij ij
where the columns of S are an orthonormal basis √ for the
states of self stress, the null space of Q′B = QB K. The coming from the second order expansion of Eq. (22), as
node and edge methods are equivalent since shown in Ref. [72]. Here, tij is the tension of edge ij at
√ √ √ equilibrium and ℓ0,ij is the rest length of the spring. We
KCbdry ϵ − KCA Γϵ ∈ col( KCB ). (77) consider tij to be positive when it pushes the nodes apart
This difference would then be annihilated by S T . and negative when it pulls them together.
One can then represent edge elongations as:
 ∥  ∥
e C
3. Finding the elastic moduli e= = U. (84)
e⊥ C⊥
As discussed above, given an imposed strain ϵ, we can Note, now we have a vector for each edge. The change
compute an effective stress tensor σ as a response. We in tension due to the displacements is given by:
can then determine the effective stiffness tensor K,  ∥  ∥ ∥ 
δt K C
δt = = U. (85)
σµν = Kαβ
µν ϵαβ . (78) δt⊥ K ⊥C ⊥
 15

We can define Q = Q∥ Q⊥ , C = QT and K =

Irreversible transport phenomena can be described
K∥ 0 through the language of thermodynamic forces and flows.
 
and write equations analogous to the non This and the other relevant thermodynamic concepts in
0 K⊥
prestressed case. Linear response problems here are bet- this section are explained in Ref. [74]. The entropy pro-
ter solved by use of the dynamical matrix D = QKC. duction per unit volume ṡ, is given by the sum of the
We have previously considered the response to perturb- products of thermodynamic forces which are gradients of
ing rest lengths, in prestressed systems we can also per- potentials ϕK and flows which are current densities JK ,
turb the stiffness and obtain non trivial response. We
can differentiate Hooke’s law to obtain: ṡ = −∇ϕK · JK . (90)

∥ Above, the index K is summed over and labels the differ-

δtij = −kij e∥ − (ℓij − ℓ0,ij )δkij + kij δℓ0,ij , (86)
ent types of transport—electrical, thermal, diffusive, etc.
Note that ṡ is not the rate of change of the entropy but
tij ⊥ the entropy production rate, the difference is that the
δt⊥
ij = − e . (87) second quantity does not take into account exchange of
ℓ0,ij ij
entropy with the surroundings. When the system is close
We can now solve the general linear response problem for to equilibrium, in the “linear regime” the forces have a
prestressed networks with equation: linear relationship with the flows.:

JK = −LJK ∇ϕJ . (91)

 
U
0 = − D, −QK ∥ , QΛ δℓ0  + F.


(88)
δk By the second law of thermodynamics, ṡ ≥ 0, where
the equality only being satisfied when all currents are
where the Λ is a diagonal matrix whose elements are zero. This implies that the matrix L is positive definite.
(ℓij − ℓ0,ij ). This problem is solved in an analogous fash- Therefore, the matrix is invertible. It is also known that
ion to the prescribed displacement problem. this matrix is symmetric—the Onsager reciprocal rela-
The solution to the edge swelling problem in pre- tions [75, 76].
stressed networks is given by Now let us apply this to networks. Consider a system
consisting of a network of pipes, wires or in general any
−1 conducting elements which come together at joints. We
δt = S T SK −1 S T

Sδℓ0 (89)
identify the joints with the nodes and the edges with the
which is essentially the same as Eq. (54) where S is such conducting elements. We make the assumption that the
that its rows are an orthonormal basis for Q as defined volume of the joints is negligible compared to the volume
in this section and K as it is defined in this section. of the conducting elements. By integrating over volume,
we can “discretize” Eqs. (90,91) on this network as:

IV. STEADY STATE RESPONSE OF Ṡ = (ϕm n K

K − ϕK ) inm , (92)
IRREVERSIBLE TRANSPORT PROBLEMS
anm K J
iK L ϕm − ϕJn ,


nm = (93)
1. Entropy production in irreversible transport ℓnm J

Many if not most transport phenomena are irre- where iKnm is the current of type K going through con-
versible, meaning that they dissipate energy and produce ducting element nm and anm , ℓnm are the cross sectional
entropy. This applies to electrical, diffusive and thermal area and the length of the conducting element respec-
transport to give a few examples. We will see that we can tively. We can define a generalized resistance matrix for
treat these problems with the same tools we have devel- each edge as:
oped in previous sections. We will describe the physics ℓnm  −1 JK
JK
in terms of graph potentials and flows and reduce lin- Rnm = L . (94)
ear response problems to minimization problems solvable anm
through the use of the fundamental theorem of linear al- Using this relationship, we can calculate the response on
gebra. networks with coupled conductivity, a small example of
First, the minimized quantity will be the entropy pro- which is shown in Fig. (16).
duction rate, often referred to as the entropy production. This resistance matrix is symmetric and positive def-
This has been called the law of minimum entropy pro- inite, it inherits this quantities from L. To make the
duction, which was introduced by Prigogine and stated analogy with Ohm’s law obviously we can label potential
as: “In the linear regime, the total entropy production K
differences as “voltages” (vnm = ϕK K
n − ϕm ) to write the
rate in a system subject to flow of energy and matter, generalized Ohm’s law:
reaches a minimum value at the nonequilibrium station-
J JK
ary state” [74]. vnm = −Rnm inm,K . (95)
 16

Heat Sink linear regimes we will show that entropy production plays
the role of elastic energy.
T1=1 Currents are analogous to tensions, voltages to elon-
4
0.

q 13
gations and potentials to displacements, as we listed in
=- μ1=3.3

=-
Tab. I. The stiffness is analogous to conductance which


7

m1
2
0.

0.
q1

4
3
 12= is the inverse of resistance. We can write some basic

=-
m

0.
relationships in terms of the incidence matrix C and its

7
0.23 0.13 transpose Q,
LJ K = 0.13 0.78

Mass Source
Mass Sink

T2=3.6 T3=2.4 K p
vnm = Cnm ϕK
p , (97)
μ2=2 μ3=4
m 24 =-1.4
q24 =0.0 IpK = Qnm K
p inm , (98)


7
m2

0.
4

 34= and the stiffness matrix:

=-

4
0.
q 24

m
0.

T4=5 -
7

 −1 J nm
=-

= nmJ
KpqK = Rpq δ . (99)
q 34 K pq
0 .4

μ4=2.7
Note that this matrix is symmetric and positive definite.
In mechanics the constitutive relation is what couples
Heat Source
the different spatial components of potentials and flows,
analogously here it is the constitutive relation that cou-
FIG. 16. A four-node system where heat and mass are both
ples the different types of irreversible transport to each
able to flow between nodes. nodes 1, 4 are connected to a heat
sink/source (fixed held temperatures in green, T1 = 1, T4 = 5)
other.
while nodes 2, 3 are connected to a mass sink/source (fixed We can also build a compatibility matrix:
held chemical potentials in green, µ2 = 2, µ3 = 4). The nJ n J
edges of this system are made from some hypothetical ma- CpqK = Cpq δK , (100)
terial where heat and mass transport are coupled such that
the Onsager coefficients are: Lheat heat
heat = 0.23, Lmass = 0.13,
where the equilibrium matrix would be its transpose.
mass
Lmass = 0.78. The fact that the off-diagonal terms are non- Now the total entropy production can be written as:
zero (Lheat
mass = 0.13) indicates the coupling between heat and
mass in this material. The calculated response of this system Ṡ = ϕ C T KC ϕ, (101)
is shown by edge and node properties that are not in green,
where arrows for edge properties point in the direction of pos- which is of the same form as that of the elastic energy of
itive flow. a spring network. Note that the generalized dynamical
matrix D = C T KC is semi positive definite. We can
perform the same change of variables as in the mechanics
Note the negative sign indicating that there is a voltage case. Because K is positive definite it can be decomposed
drop as we go from node m to node n given that the as
current goes from m to n. √ T√
The equation for the entropy production in terms of K= K K. (102)
currents and resistances is now:
We can include the conductances (“stiffnesses”) in the
Ṡ = i nm,J K
Rnm,J inm,K , (96) compatibility matrix by defining:
√
C ′ = KC, (103)
Where Ṡ is the total entropy production rate. For elec-
trical transport this reduces to the power dissipated ri2 . and describe currents and voltages with a single vector:
All of this terms are dissipative in nature. For this rea-
son the “law of minimum entropy production” can be √ √ T
−1
understood as a law of least dissipation. w= Kv = K i. (104)

The entropy production is then the norm squared of this

vector:
A. Analogy between static mechanics and steady
state in irreversible flow
Ṡ = |w|2 . (105)

Here we will write the equations for irreversible flow in This allows us to solve linear response problems in a
the same form as the ones we used for mechanics. In this standardized manner just like it did for mechanical
way we show that there is an analogy between mechanical networks.
equilibrium and thermal equilibrium. In their respective
 17

There are two basic response problems, imposing po- A. General analogy
tential differences or imposing current sources. Current
sources on the nodes are analogous to forces on the nodes We start from the general physical network relations
on the mechanical case and imposing potentials is analo- we wrote in Eq. (10). In the case of electric networks,
gous to prescribing displacements. We have solved these Yn , ynm represent node and edge admittance. Using the
problems in previous sections. These solutions apply frequency domain notation, familiar circuit elements have
identically to irreversible flow problems. admittance
Imposing potential differences on the edges is analo-
gous to changing the rest lengths of the springs. For Ycapacitor = iωC, (107)
an electrical problem, rest-lengths are voltage sources, Yresistor = 1/R, (108)
batteries [77]. The analogy can be extended to voltage Yinductor = 1/(iωL). (109)
sources of other types, not just electrical. We are using
the word voltage here to refer to any potential difference. The admittance of an edge or a node can be any combi-
What these voltage sources are physically will depend on nations of these elements in parallel or series, which can
the context. be written as a complex number in general.
It is useful to write the mechanical problem in an anal-
In our search to find realizations of these voltage
ogous form, so we have a universal language for the dy-
sources we might take a clue from the everyday battery,
namical response of networks. To do this, we adopt an
which is a chemical battery. Chemical reactions can be
analogy developed by Firestone in 1933 [78]. Here, in-
sources, or sinks, of heat, chemical species and charge.
stead of node displacement, we take node velocity as po-
This naturally relates to thermal, diffusive and electrical
tentials, while the flows remain forces. This allows a more
transport respectively. In fact, chemical reactions gen-
natural analogy where the admittance of the mechanical
erally involve irreversibility and can be written in terms
problem correspond to, as shown in Fig. (17):
of thermodynamic forces and flows, included in Eq. (90)
and coupled with the other irreversible processes via the Ymass = iωm, (110)
matrix LKJ . And so, it seems there is a great richness of
Ydashpot = η, (111)
phenomena within the linear regime that could be stud-
ied with the framework presented here. Yspring = k/(iω). (112)
For chemical reactions, the thermodynamic force is the We describe other analogies between mechanical and
affinity A, which depends linearly on the chemical po- electrical networks in the Appendix.
tentials, but cannot generally be expressed as a potential
difference. This is in contrast to transport phenomena.
The thermodynamic flow is the velocity of reaction, also B. Timescales of physical networks
called “rate of conversion,” ξ. ˙ The entropy production
due to chemical reactions is: Beyond the static and steady state responses we dis-
cussed in Secs. III,IV, which are governed by the null
AK ˙K space of the Q, C (or Q, C) matrices, we need to utilize
Ṡchem = ξ , (106) the full spectrum of these matrices to characterize the
T
dynamical response, where forces don’t balance on the
where the index K labels the different chemical reactions nodes and currents are not conserved, representing the
that may be taking place. conversion among different forms of energy.
Following the discussion in Sec. II, the general equa-
tion of motion in the frequency domain, using matrix
notation, takes the form

V. DYNAMICS AND WAVES Y (ω)V (ω) = A(ω)V (ω) + S(ω) + Qy(ω)w, (113)
where
In this section we present a more general analogy of
dynamical processes on networks that involve both en- Aij ≡ Qnm i
j ynm Cnm . (114)
ergy conserving, dissipative, and active components. The Note that we base our discussion here on the scalar prob-
mechanical equilibrium case presented in Sec. III and the lem for simplicity. The vector problem can be formulated
thermal equilibrium case presented in Sec. IV can be seen using the scheme presented in Sec. II D.
as special cases of this more general analogy, where the Here we first discuss timescales of the problem, and the
simple frequency dependence of the edge and node ad- general response to a dynamic signal will be discussed in
mittance permits a neat scheme based on minimization. the next subsection (Sec. V C). To do this, we study the
Here, due to the more general dynamical rules, we will characteristic equation of these problems
directly write dynamical equations on nodes and edges,
instead of constructing a total energy or entropy. Y (ω)V (ω) = A(ω)V (ω), (115)
 18

edges y k
mass iωm dashpot η spring iω
nodes Y
ω ′′ = spec(M −1 A) ω ′ = spec(M −1 K)
p
mass iωM δ = 2π

dashpot H ω ′′ = spec(M−1 H) δ=π ω ′′ = spec(H −1 K)

ω ′ = spec(M−1 K) ω ′′ = spec(KA−1 )
K
p
spring iω
δ=0

TABLE II. Mechanical network dynamical timescales for different combinations of edge and node elements.

problem as a simple example. In this case,

(a)
Y (ω) = iωM, (116)
Uext A(ω) = A, (117)
where M is the constant matrix containing all node
masses and A = Qnm i
j ηnm Cnm is the constant matrix
η i containing all edge dashpot coefficients (i.e., all edges
K
are characterized by dashpots containing viscous fluids
where force is proportional to velocity with coefficient
inert. ηnm ). The characteristic equation is then

M iωV (ω) = M −1 AV (ω), (118)

determining that the homogeneous solution is a linear
combination of exponential decay modes with decay rate
(b ) determined by the eigenvalues of M −1 A,
Vext iω = spec(M −1 A), (119)

where spec(M −1 A) denote the set of eigenvalues of

(M −1 A). This enables us to write the general homo-
C
R L geneous solution as
X
V (t) = e−λα t , (120)
α

where λα are eigenvalues of M −1 A, labeled by α.

One can similarly write this homogeneous solutions for
C other cases. For the dashpot-on-nodes and spring-on-
edges problem, we have
Y (ω) = H, (121)
FIG. 17. Analogous components (same color/position) be- 1
tween mechanical (a) and electrical (b) networks as described A(ω) = K, (122)
by Firestone’s analogy [78]. iω
giving decay rates
iω = spec(H −1 K). (123)
which correspond to the homogeneous solution to
Eq. (113). Note that this is not yet a conventional eigen- The cases of dashpot-on-nodes and mass-on-edges, as
value problem, as Y (ω), A(ω) depend on ω in the general well as spring-on-nodes and dashpot-on-edges, follow
case, as we discussed above. The general homogeneous similar analysis, giving pure decay solutions.
solution involves a nonlinear equation of ω. These cases are characterized by exponential decay be-
cause the equation of motion involves only the first order
Nevertheless, it is instructive to consider a few simple time derivative. We get oscillatory solutions when there
cases, where we have only one type of elements for the are second order time derivatives, which shows up in two
nodes and one type of elements for the edges. This en- cases in this context. The first one is the mass-on-nodes
ables us to pull out the ω dependence, and turn the char- and spring-on-edges problem,
acteristic equation into a conventional eigenvalue prob-
lem. Let’s take the mass-on-nodes and dashpot-on-edges Y (ω) = iωM (124)
 19

1 which is typically evaluated using contour integrals

A(ω) = K, (125)
iω (adding a semicircle at Ω′′ → −∞ for retarded Green’s
giving oscillatory frequency scales functions (t > 0)), and thus governed by the poles of the
p frequency domain solution Eq. (131).
ω = spec(M −1 K), (126) Similarly, for physical networks, the calculation is sim-
1
which is a generalization of the familiar harmonic oscil- ilar, with the simple factor mΩ2 −iγΩ−k replaced by the
network green’s function (Y (Ω)I −A(Ω))−1 , but the gen-
p
lator frequency k/m.
A similar case arises when we put masses on edges and eral formulation follows similarly. The most nontrivial el-
springs on nodes, ement here, in some sense, is the complex Ω dependence
in both Y (Ω) and A(Ω), causing nonlinear equations de-
1 termining the poles.
Y (ω) = K (127)
iω We show some examples of this computation in
A(ω) = iωM, (128) Fig. (18) where a mechanical network is driven at given
giving oscillatory frequency scales frequencies at chosen boundary nodes.
p
ω = spec(M−1 K), (129)
which shares a form that appears similar to Eq. (126) but VI. CONCLUSION AND OUTLOOK
has important differences, as M describes edge masses
and K nodes spring constants. We have developed a comprehensive and unified frame-
Finally, we should also include the “diagonal” cases work for analyzing the linear response of physical
where both nodes and edges have only mass, only dash- networks—ranging from mechanical lattices to electri-
pot, and only springs. These cases exhibit no timescales cal circuits, thermal transport systems, and beyond—
and only trivial solution of V = 0 for the homogeneous grounded in algebraic graph theory. This approach of-
equation. A summary of these cases can be found in fers a mathematically rigorous yet versatile language that
Tab. V B. Other types of problems on physical networks, captures both static and dynamic behaviors, scalar and
such as electrical, heat and mass transport, can be cal- vectorial flows, and conservative as well as dissipative
culated in similar ways. processes. By identifying flows and potentials on net-
It is also worth noting that such eigenvalues analysis work edges and nodes, and systematically relating them
is enabled by assuming simple elements on the physical via incidence and cut-set matrices, we derived a set of
network. The general case where Y (ω) and A(ω) depend universal equations governing physical responses across
on ω in general ways can not be fully captured in such domains.
simple analysis. One will need to solve the full equation. Our formulation reveals deep analogies between seem-
ingly disparate systems. For instance, mechanical force
balance and thermal current conservation share algebraic
C. Finite-frequency response of physical networks structure through Kirchhoff’s laws, while entropy pro-
duction in irreversible systems mirrors elastic energy in
With these homogeneous solutions in place, we can static networks. Furthermore, our use of linear alge-
now discuss the inhomogeneous solution to the general bra—particularly the interplay of column spaces, null
equation (113). The solution can be formally written as spaces, and dual graphs—unifies these analogies and sug-
gests that many known results (such as self-stress states,
V (Ω) = (Y (Ω)I − A(Ω))−1 S̃(Ω), (130) floppy modes, and reciprocal diagrams) are manifesta-
tions of a more general structure.
where S̃(Ω) = S(Ω) + Qy(Ω)w is the effective node
sources including the edge driving, and I is the iden- Looking forward, this framework opens multiple av-
tity matrix. We use uppercase Ω to signify the signal enues for exploration. One immediate extension is multi-
frequency. physics [29] on networks: The algebraic structure nat-
It is helpful to think of this as a generalization of a urally accommodates multiple physical variables (e.g.,
driven damped harmonic oscillator problem, where the stress, heat, mass, charge) and their couplings via gen-
frequency domain response is eralized admittance or resistance matrices, providing a
principled path to study cross-phenomena effects like
f (Ω) thermoelectric or chemo-mechanical responses.
u(Ω) = , (131)
mΩ2 − iγΩ − k While our work focuses on linear responses, the under-
lying network formalism is amenable to generalization.
where m, γ, k are the mass, drag (dashpot) coefficient,
Incorporating nonlinear elements (e.g., threshold conduc-
and spring constant of the oscillator. The time domain
tances, active feedback) or time-dependent control (e.g.,
response is then
Z space-time modulation) could extend this method to ac-
f (Ω) tive matter, neuromorphic systems, and non-Hermitian
u(t) = dΩ e−iΩt (132) physics.
mΩ2 − iγΩ − k
 20

Springs Only Dashpots Only Springs and Dashpots

ω = 0.05
ω = 0.5

FIG. 18. The dynamic response of spring only (left), dashpot only (middle), and spring with dashpot (right) networks to
imposed dynamic displacements at frequencies ω = 0.05 (top) and ω = 0.5 (bottom). Applied displacements (dark green) are
periodic horizontal oscillations to nodes at the top and bottom of the network. The osculations to the top of the network are
exactly out of phase with those at the bottom of the network, shown by dots indicating the position of the node at an equal
time (light green). The response displacements of nodes in the bulk (red), which in general are elliptical orbits (direction and
equal time positions shown by orange arrow, orange dot when no eccentricity), penetrates further into the bulk of system at
lower frequencies, more closely resembling the static displacement problem. Edges (black lines) have thickness proportional to
the magnitude of their complex tension response.

Furthermore, the ability to condense nodes and de- bul Chakraborty, Mark Newman, and Anthony Grbic for
scribe macroscopic strain or stress from microscopic de- helpful discussions. This work was supported in part
tail hints at a promising direction for renormalization- by the National Science Foundation Center for Complex
like procedures on networks [79]. This could lead to sys- Particle Systems (Award #2243104), National Science
tematic coarse-graining tools that preserve key physical Foundation Award #2026825, and the Office of Naval
responses across scales. Research (MURI N00014-20-1-2479).
Beyond physics and engineering, our formulation res-
onates with recent developments in machine learning and
artificial intelligence, particularly in graph-based models Appendix A: Graph Theory and Linear Algebra
and network inference. The algebraic structure underpin-
ning physical networks parallels the architecture of graph Here we will review some basic concepts of linear al-
neural networks (GNNs) [80], where information propa- gebra and graph theory as they apply to our present
gates along edges and is aggregated at nodes, opening a work. There appears to be much variation in the lan-
route for data-driven discovery of effective models and guage and conventions previously used when discussing
materials. At the same time, these tools also provide this topic. Here we make our own language and con-
a convenient basis for the study of physical neural net- ventions clear. First we discuss the concept of cycles
works, where intrinsic physical dynamics on networks are and cut sets [27, 30, 36, 87]. Then we discuss the fun-
utilized for information processing [81–86]. damental theorem of linear algebra (FTLA) [35] and its
Acknowledgments.—The authors thank Nick Kotov, Bul- consequences for physical networks.
 21

1. Cycles and Cut Sets A cut set is a subset of the edges that when removed
divides the graph into two disconnected parts. More
When studying transport on a network it is useful to formally a cut set is a set of edges such that there ex-
focus on the edges, since this is where the transport takes ists a partition of a graph into two sets of nodes such
place. Certain sets of edges are of fundamental impor- that all edges in the set have one node in each element
tance. Here we will discuss cycles and cut sets. A cut set of the partition. Now, consider G our example graph
is a set of edges such that if eliminated splits the graph (1, 2), (1, 3), (2, 3), (2, 4). Note G is connected, if we elim-
into two connected components. A cycle is a set of edges inate edges (1, 2) and (2, 3) we disconnect it into two com-
that form a closed path. In what follows we will describe ponents Ga and Gb whose nodes are {1, 3} and {2, 4}. We
cycles and cut sets of a graph in linear algebra terms as have found a cut-set. Like in the case of cycles, cut sets
vectors in RNe where Ne is the number of edges in the have to be written with a consistent orientation of the
graph. edges. Cut sets should be written all pointing from com-
Firstly we start with a undirected graph G and choose ponent Ga to Gb (or all vice-versa). Our cut set is then
some ordered list of its edges to represent it. We call (1, 2), (3, 2), which maps to the vector ⟨1, 0, −1, 0⟩. The
this ordered list the edge list. For example, let G be linear subspace spanned by the ct sets is called the cut
the graph with edge list (1, 2), (1, 3), (2, 3), (2, 4). This space and a basis for it is called a cut basis.
chosen ordering of the graph, for both edges and nodes One can easily verify that our example cycle is orthog-
within edges, is a convention and should have no effect onal to our example cut set. This is in fact always the
on physics. Note that this ordering allows us to refer to case. For a given graph all of its cut sets are orthogo-
edges by their index. For example edge number 3 is (2, 3). nal to all of its cycles. Since the cut space and the cycle
In terms of vectors we write (2, 3) as ⟨0, 0, 1⟩ and (3, 2) as space are orthogonal complements. This is a well
⟨0, 0, −1⟩. In this manner we can represent any subset of established result and is related to the fundamental the-
edges and the orientation in which they are given. For ex- orem of linear algebra. We will not provide a proof. One
ample (1, 2), (3, 1) would correspond to vector ⟨1, −1, 0⟩. way to intuitively see why this is the case is to consider
To be more explicit the vector corresponding to a given a cycle and a cut set which are not disjoint. The cycle
edge subset G′ with is own edge list, is obtained in the must intersect the cut set an even number of times. Say
following way: if the nth edge of G is not on G′ the nth the cut set separates the graph into components Ga and
element of the vector is 0, it is 1 if it is on G′ and listed Gb . The cycle must exit Ga the same number of times
the same way and −1 if it is there but listed in reverse it enters Ga . In each exit the cycle is oriented along the
order. cut set but on each entry its oriented opposite to the cut
set. These naturally makes the two vectors orthogonal.
In other presentations of this topic the graph is treated
Now, one can relate a number of matrices with a graph
as a directed graph. Here we have chosen to use this
G. It is common to use the adjacency matrix as a repre-
object we call the edge list and continue to refer to G
sentation of a graph. One inconvenience of this approach
as undirected for a number of reasons. First there is a
is that the adjacency matrix goes from the nodes to the
conceptual reason, we want to make clear that the ori-
nodes. Since we are interested in studying what goes on
entation induced on the graph by this ordered list will
in edges we will need different matrices. Let us intro-
not correspond to any inherent property of our object of
duce the incidence matrix C which in our convention has
study, a physical network. The orientation is just a con-
as many rows as there are edges and as many columns
vention, similar to choosing a reference frame. There is
as there are nodes. Here we will also make use of an or-
also a practical computational aspect, which is that this
dered edge list for the graph. Let EdgeList(n) refer to
ordered edge list allows us to build vectors and matrices
the ordered pair that is the nth element of the edge list.
and perform linear algebraic operations in a clear stan-
We then write:
dardized way. This is has proven crucial in the process 
of writing code to numerically perform all the methods 1,
 if (i, j) = EdgeList(n), ∃j
we will discuss. Cn,i = −1, if (j, i) = EdgeList(n), ∃j (A1)
A cycle in a graph can be represented by a list of edges 
0, otherwise.
for example (1, 2), (2, 3), (3, 1), note how we organize the
edges so that the second node of one edge is the first node The column space of this matrix is the cut set of the
of the next one in the list and the first and last nodes are graph.
the same. In this way we represent that we have taken a Another relevant matrix is the cycle matrix B, (also
closed walk visiting the nodes in sequence 1, 2, 3, 1. The circuit matrix) which for us is any matrix such that its
cycle, a subgraph G′ of our previously defined G , de- rows form a complete cycle basis. Traditionally this ma-
scribed as an ordered list (1, 2), (2, 3), (3, 1) maps to the trix is more narrowly defined by requiring that all rows
vector w = ⟨1, −1, 1, 0⟩. We will use the word cycle to re- correspond to cycles as we have described them before.
fer to the closed trail, the subgraph, and its vector repre- This results in all entries being 0,1 or -1. Since the cycle
sentation interchangeably. The linear subspace spanned space and cut space are orthogonal complements:
by the cycles is called the cycle space and a basis for
this space is called a cycle basis. 0 = BCx , ∀x. (A2)
 22

Further, since the two spaces are complements, any vec- 3. General Least Squares Minimization Problem
tor w in the space of edges can be represented as a sum
of a cut set component and a cycle component: Consider the following problem. Given a real valued
w = Cx + B x̃. T
(A3) m × n matrix C and an m × 1 vector w find w − Cx such
that |w − Cx| is minimized . We will see that we can cast
a variety of physical network problems into this form and
2. Fundamental Theorem of Linear Algebra solve them by the same generalized approach.
This problem can be solved by using the Moore Pen-
rose pseudo inverse C + ,
Here we will discuss the fundamental theorem of linear
algebra and explain what it tells us about cycles and xmin = C + w. (A9)
cut sets of a graph. The fundamental theorem of linear
algebra states [35]: Given a real valued m × n matrix A, The matrix C + is obtained from the singular value de-
composition. The solution to our problem is then:
1. dim col(A) = dim col(AT ) and, dim col(A) +
w − Cxmin = 1 − CC + w.


dim null(AT ) = m, (A10)

2. null(AT ) and col(A) are orthogonal complements. Let us explore the problem from a different perspective.
It is clear that w − Cx = 0 only has solution if w is in
Let us apply the FTLA to the incidence matrix. We the column space of C. From a geometric point of view
know that col(C) is the cut space. The FTLA then tells the closest Cx can get to w is the projection of w along
us that null(C T ) is the orthogonal complement of the cut the column space of C. The minimum difference is then
space: the cycle space, given by the projection into the orthogonal complement
of C (See Fig. (8)):
null(C T ) ≡ col(B T ). (A4)
In other words: w − Cxmin = B̂ T B̂w, (A11)

0 = C T B T x̃ , ∀x̃. (A5) where B̂ is a matrix such that its rows form an orthonor-
mal basis for null(C T ). We can alternatively define a
Part one of the FTLA theorem is also known as the matrix Cˆ such that its columns are an orthonormal basis
rank nullity theorem from which one can easily show that for col(C). One suitable choice for this basis are the left
for any real valued m × n matrix , say C: singular vectors of C. We can then write:
dim null(C T ) − dim null(C) = m − n. (A6)
 
w − Cxmin = 1 − CˆCˆT w. (A12)
One might call this a fundamental counting rule. If C is
the incidence matrix of a given graph, then:
Appendix B: Other Electrical-Mechanical Analogies
• dim null(C T ) = Ncycles : The dimension of the cycle
space.
As mentioned in the main text, there are a number of
• dim null(C) = Nc.c. : The number of connected ways one can map electrical circuits to mechanical ones
components. and vice versa. In addition to Firestone’s analogy [78]
discussed in the main text, here we will discuss the more
• m = Nedges .
conventional impedance analogy,sometimes attributed to
• n = Nnodes . Maxwell [88], and our own original analogy.
For the discussion of these analogies it is convenient to
The counting rule then gives us the following relation- talk in terms of impedance z which is the reciprocal of the
ship: admittance z = 1/y. “Mechanical impedance” is defined
Ncycles − Nc.c. = Nedges − Nnodes , (A7) differently for each analogy. We will present a number of
equations of motion all analogous to Eq. (113).
note that Nc.c. (the number of connected components)
is 1 for a completely connected graph. We will later
see that this is a case of what is often referred to as 1. Impedance Analogy: force as potential, velocity
Maxwell’s counting [39]. as current

This counting relationship could also be written in The impedance analogy is the oldest and more con-
terms of the cut space by making use of the fact that be- ventional analogy [78]. It maps electrical impedance to
cause it is the orthogonal complement to the cycle space mechanical impedance defined as the frequency depen-
its dimension Ncut satisfies dent complex quantity that relates velocity to force as:
Nedges = Ncycles + Ncut . (A8) f = z u̇, (B1)
 23

where u̇ is the node velocity and f the force on a node. (a)

Similarly, edge impedance is the coefficient that relates
the rate of change of the elongation and edge tension:
t = z ė, (B2)
η Fext
while electrical impedance is the coefficient that relates M
voltage and current: 1
V = zI. (B3) K 2
The Impedance Analogy is based on the analogy be-
tween Eqs. (B2,B3). Voltage maps to tension, electric
potential to force, the rate of elongation maps to edge
currents and node velocity maps to current sources. This
mapping is opposite to the electrical-mechanical analogy
we used in previous sections for the static case. The (b )
R
impedance analogy has a serious drawback which is that 1
the mechanical circuit analogous to a given mechanical
circuit does not have the same topology (as shown in
Fig. 19). Elements that are in parallel in one are in se-
ries in the other and vice versa. More generally under the L
impedance analogy cycles in the electrical circuit map to
cut sets in the mechanical and vice versa. This means Vext 2
that the graphs underlying the two circuits are the dual
of each other.
We can see this duality at play when we consider a 1D C
spring network initially unstressed and then driven at a
frequency ω. The elongation rates will be compatible just
like the elongations themselves meaning that they add to
zero along cycles:
FIG. 19. Analogous components (same color/position) be-
X
0= ėij . (B4)
tween mechanical (a) and electrical (b) networks as described
ij∈cycle
by the Impedance Analogy. Note that the topology of these
For the electric analogous electric circuit we will have the two networks is different. In the mechanical analogy, the dash-
current conservation law: pot and spring are in paralell, but their analogous components
X in the electrical network (resistor and capacitor) are in series.
0= Iij . (B5)
ij∈cut set
where C̃ is the incidence matrix (or compatibility ma-
In order for these two conditions map to each other under
trix) for the dual graph. Note that this would be the
the impedance analogy cycles must map to cut sets and
more natural way of writing this equation according to
vice versa.
the impedance analogy where force plays the role of po-
Our general equation Eq. (113) for mechanical systems,
tential.
with impedance defined as z = t/ė would read
Λ−1 u̇ = −QZC u̇ + QZ ė0 + f0 (B6)
2. Firestone’s Analogy: velocity as potential, force
Here Z and Λ are the diagonal matrices of mechanical as current
impedances of edges and nodes respectively. In the con-
text of Eq. (113) Y = Λ−1 , A = QZC. Note that al-
though in the impedance analogy u̇ is formally analogous In 1933, Firestone, noticing the problem of the
to current it appears in the equation in the same way as impedance analogy not being topology preserving on a
electrical potential. From an algebraic perspective one network, proposed a different analogy [78]. This anal-
would more naturally arrive at a different analogy where ogy is sometimes called the admittance analogy since it
velocity is potential. This is precisely Mobility analogy maps electrical elements to mechanical ones in such a
we will discuss next. way that the electrical impedance matches the mechani-
One could alternatively write the equation of motion cal admittance (reciprocal of impedance). For clarity, we
as: refer to this analogy as the “Firestone” analogy or “Fire-
stone’s” analogy, which we discussed significantly in the
Λf = −Q̃Z −1 C̃f + Q̃Z −1 t0 + u0 , (B7) main text. Firestone’s analogy preserves topology and is
 24

Mechanical Elements Mechanical Impedance Electrical Elements Electrical Impedance

Mass/Inertia iωM inductor iωL

Dashpot η resistor R

Spring K/iω capacitor i/iωC

TABLE III. Impedances of mechanical and electrical elements in the Impedance Analogy.

more in line with the theory we have presented in this

work. In his original paper, Firestone went as far as to (a)
strongly propose that mechanical impedance be redefined
to u̇/f . For this reason Firestone called the mechanical Uext
admittance the bar impedance z̄ a notation we will make
use of here. Mechanical admittance is also called mobility
and so the analogy is also sometimes called the mobility
analogy. η
Mobility z̄ is the coefficient that relates forces and dis- K f
placements as:
fixer
u̇ = z̄f, (B8)
M
where u̇ is the node velocity and f the force on a node.
Similarly, edge mobility is the coefficient that relates the
rate of change of the elongation and the tension:
(b )
ė = z̄ t. (B9) Vext
When we map the mobility to electrical impedance we
find that tension maps to current, force on nodes to cur-
rent sources, elongation rate to voltage and velocity to C
electric potential. Now elongation rates that sum to zero R L
along cycles map to voltages that sum to zero along cy-
cles.
In this analogy, charge maps to momentum, therefore
charge flow (current) maps to momentum flow (force).
The analogies between the different elements can be ex-
plained conceptually through this idea. For example, in C∝ω
the same way a capacitor stores charge a mass or a inertia
stores momentum.
In terms of bar impedance we can write the equation
FIG. 20. Analogous components (same color/position) be-
of motion as: tween mechanical (a) and electrical (b) networks as described
by the Static Consistent Analogy. Note that the node-
Λ̄u̇ = −QZ̄ −1 C u̇ + QZ̄ −1 ė0 + f0 . (B10) capacitors (analogous to mass) must have a capacitance that
is proportional to frequency.
Note how this equation is essentially identical to the elec-
trical case. In the context of Eq. (113) Y = Λ̄, A =
QZ̄ −1 C.
fact a third analogy. The usefulness of the “displace-
ment as potential” analogy lies in the fact that it holds
3. “Static Consistent” Analogy: displacement as for both static and dynamic problems where as the other
potential, force as current analogies do not apply in the static case. For this reason,
we call this newly introduced analogy the “Static Con-
In this work we previously presented a mechanical- sistent” Analogy, as it is consistent between the static
electrical analogy between static equilibrium and a con- and dynamic cases. Fig. (20) shows the analogous ele-
stant current steady state. We identified the node dis- ments of this analogy, note that springs and resistors are
placement with the electric potential. This does not agree analogous here, even with time-dependence.
with either the impedance or mobility analogy, it is, in We define our stand in for the impedance as z̃ which
 25

Mechanical Elements Mechanical Mobility Electrical Elements Electrical Impedance

Mass/Inertia 1/iωM Capacitor 1/iωC

Dashpot 1/η Resistor R

Spring iω/K Inductor iωL

TABLE IV. Mobilities and impedances of mechanical and electrical elements in the Firestone Analogy.

satisfies: length. The corresponding equation for this element is:

u = z̃f (B11) ṫ = −η e. (B13)

It is worth noting that this is a mechanical element with

for nodes, and memory. The value of tension depends on the past states
of the edge:
e = z̃ t (B12) Z T
tT = t0 + −η e dT. (B14)
for edges. The quantity z̃ is a generalized “compliance” 0
(reciprocal of stiffness) in the same way impedance is a
One peculiar aspect of this analogy is that it maps con-
generalized resistance.
servative electrical elements to non conservative mechan-
An edge mass/inertia maps to a capacitor such that ical elements. For example, capacitors (conservative) are
the capacitance depends linearly on the frequency. Such mapped to dashpots (dissipative). This analogy, when
an element is not so unrealistic since the electrical per- applied to dynamics, gives us a way to map a family of
mittivity is frequency dependent, so all real capacitors active mechanical systems to passive electrical ones.
that make use of dielectrics have frequency dependent The equation of motion for this analogy now written
capacitances. A capacitor with capacitance c = c0 + cω ω as:
can be represented as a dashpot and an edge mass/inertia
connected in parallel. Λ̃u = −QZ̃ −1 Cu + QZ̃ −1 e0 + f0 . (B15)
The “fixer” is an active element that keeps increasing
or decreasing the tension as long as it is not at its rest In the context of Eq. (113) Y = Λ̃, A = QZ̃ −1 C.

[1] D. Cioranescu and P. Donato, An introduction to homog- Society A: Mathematical, Physical and Engineering Sci-
enization (Oxford university press, 1999). ences 371, 20120375 (2013).
[2] S. Torquato et al., Random heterogeneous materi- [12] M. Newman, Networks (Oxford university press, 2018).
als: microstructure and macroscopic properties, Vol. 16 [13] D. S. Bassett, E. T. Owens, K. E. Daniels, and M. A.
(Springer, 2002). Porter, Physical Review E—Statistical, Nonlinear, and
[3] G. W. Milton and A. Sawicki, Appl. Mech. Rev. 56, B27 Soft Matter Physics 86, 041306 (2012).
(2003). [14] L. Zhang, D. Z. Rocklin, L. M. Sander, and X. Mao,
[4] R. Lakes, Nature 361, 511 (1993). Physical Review Materials 1, 052602 (2017).
[5] J. Park and R. S. Lakes, Biomaterials: an introduction [15] L. Papadopoulos, M. A. Porter, K. E. Daniels, and D. S.
(Springer Science & Business Media, 2007). Bassett, Journal of Complex Networks 6, 485 (2018).
[6] X. Mao and N. Kotov, MRS Bulletin 49, 352 (2024). [16] J. E. Kollmer and K. E. Daniels, Soft matter 15, 1793
[7] D. Lukkassen and G. W. Milton, in Proceedings of the (2019).
Conference on Function Spaces, Interpolation Theory [17] N. Dehmamy, S. Milanlouei, and A.-L. Barabási, Nature
and Related Topics in Honour of Jaak Peetre on his 65th 563, 676 (2018).
Birthday (2000) pp. 311–324. [18] S. Zhang, L. Zhang, M. Bouzid, D. Z. Rocklin,
[8] J.-L. Bouvard, D. K. Ward, D. Hossain, S. Nouranian, E. Del Gado, and X. Mao, Physical review letters 123,
E. B. Marin, and M. F. Horstemeyer, (2009). 058001 (2019).
[9] E. Weinan, Principles of multiscale modeling (Cambridge [19] D. Z. Rocklin, L. Hsiao, M. Szakasits, M. J. Solomon,
University Press, 2011). and X. Mao, Soft Matter 17, 6929 (2021).
[10] A. Ramı́rez-Torres, R. Penta, R. Rodrı́guez-Ramos, [20] S. Machlus, S. Zhang, and X. Mao, Physical Review E
A. Grillo, L. Preziosi, J. Merodio, R. Guinovart-Dı́az, 103, 012104 (2021).
and J. Bravo-Castillero, Computing and Visualization in [21] D. A. Vecchio, S. H. Mahler, M. D. Hammig, and N. A.
Science 20, 85 (2019). Kotov, ACS nano 15, 12847 (2021).
[11] A.-L. Barabási, Philosophical Transactions of the Royal [22] R. Yang, K. Bernardino, X. Xiao, W. R. Gomes, D. A.
 26

Mechanical Elements Mechanical Compliance (z̃) Electrical Elements Electrical Impedance

2
Mass/Inertia 1/ω M Capacitor* 1/iωC, C ∝ ω

Dashpot 1/iωη Capacitor 1/iωC

Spring 1/K Resistor R

“Fixer” iωf Inductor iωL

TABLE V. Mobilities and impedances of mechanical and electrical elements in the Static Consistent Analogy.

Mattoso, N. A. Kotov, P. Bogdan, and A. F. de Moura, communications 8, 14201 (2017).

Advanced Science 11, 2402464 (2024). [47] H. Xiu, H. Liu, A. Poli, G. Wan, K. Sun, E. M. Ar-
[23] G. Pete, Á. Timár, S. Ö. Stefánsson, I. Bonamassa, and ruda, X. Mao, and Z. Chen, Proceedings of the National
M. Pósfai, nature communications 15, 4882 (2024). Academy of Sciences 119, e2211725119 (2022).
[24] T. Lubensky, C. Kane, X. Mao, A. Souslov, and K. Sun, [48] H. Xiu, I. Frankel, H. Liu, K. Qian, S. Sarkar, B. Mac-
Reports on Progress in Physics 78, 073901 (2015). Nider, Z. Chen, N. Boechler, and X. Mao, Proceedings
[25] X. Mao and T. C. Lubensky, Annual Review of Con- of the National Academy of Sciences 120, e2217928120
densed Matter Physics 9, 413 (2018). (2023).
[26] B. W.A., Mathematical Modeling of Physical Networks [49] L. Asimow and B. Roth, Transactions of the American
(The Macmillan Company, New York, 1968). Mathematical Society 245, 279 (1978).
[27] B. Bollobás, “Electrical networks,” in Modern Graph [50] N. Noll, M. Mani, I. Heemskerk, S. J. Streichan, and
Theory (Springer New York, New York, NY, 1998) pp. B. I. Shraiman, Nature Physics 13, 1221 (2017).
39–66. [51] J. C. Maxwell, The London, Edinburgh, and Dublin
[28] B. Szabó and I. Babuška, (2021). Philosophical Magazine and Journal of Science 27, 250
[29] W. B. Zimmerman, Multiphysics modeling with finite ele- (1864).
ment methods, Vol. 18 (World Scientific Publishing Com- [52] J. C. Maxwell, Transactions of the Royal Society of Ed-
pany, 2006). inburgh 26, 1–40 (1870).
[30] M. Swamy and K. Thulasiraman, Graphs, Networks, [53] L. Cremona, Graphical statics: two treatises on the
and Algorithms, Wiley-Interscience publication (Wiley, graphical calculus and reciprocal figures in graphical stat-
1981). ics (Clarendon Press, 1890).
[31] N. Biggs, Algebraic graph theory, 67 (Cambridge univer- [54] H. Crapo and W. Whiteley, Beitrage zur Algebra und
sity press, 1993). Geometrie 35, 259 (1994).
[32] O. Shai, Mechanism and Machine Theory 36, 343 (2001). [55] T. Mitchell, W. Baker, A. McRobie, and A. Mazurek,
[33] G. Kirchhoff, Annalen der Physik 148, 497 (1847). International Journal of space structures 31, 85 (2016).
[34] G. Kirchhoff, IRE Transactions on Circuit Theory 5, 4 [56] T. Mitchell, W. Baker, and A. Mcrobie, in Proceedings
(1958). of IASS Annual Symposia, Vol. 2015 (International As-
[35] G. Strang, The American Mathe- sociation for Shell and Spatial Structures (IASS), 2015)
matical Monthly 100, 848 (1993), pp. 1–12.
[Link] [57] J. H. Snoeijer, T. J. Vlugt, M. van Hecke, and W. van
[36] E. K. Lloyd, J. A. Bondy, and U. S. R. Murty (1978). Saarloos, Physical review letters 92, 054302 (2004).
[37] D. Zhou, L. Zhang, and X. Mao, Phys. Rev. Lett. 120, [58] B. P. Tighe, J. H. Snoeijer, T. J. Vlugt, and M. van
068003 (2018). Hecke, Soft Matter 6, 2908 (2010).
[38] D. Zhou, L. Zhang, and X. Mao, Phys. Rev. X 9, 021054 [59] S. Henkes and B. Chakraborty, Physical Review
(2019). E—Statistical, Nonlinear, and Soft Matter Physics 79,
[39] C. Calladine, International Journal of Solids and Struc- 061301 (2009).
tures 14, 161 (1978). [60] D. J. Jacobs and M. F. Thorpe, Physical review letters
[40] S. Pellegrino and C. R. Calladine, International Journal 75, 4051 (1995).
of Solids and Structures 22, 409 (1986). [61] C. S. O’hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,
[41] C. Calladine and S. Pellegrino, International Journal of Physical Review E 68, 011306 (2003).
Solids and Structures 27, 505 (1991). [62] W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. vAN
[42] S. D. Guest and J. W. Hutchinson, Journal of the Me- SAARLoos, Physical review letters 97, 258001 (2006).
chanics and Physics of Solids 51, 383 (2003). [63] W. G. Ellenbroek, V. F. Hagh, A. Kumar, M. Thorpe,
[43] K. Sun, A. Souslov, X. Mao, and T. C. Lubensky, Pro- and M. Van Hecke, Physical review letters 114, 135501
ceedings of the National Academy of Sciences 109, 12369 (2015).
(2012). [64] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Physical
[44] C. L. Kane and T. C. Lubensky, Nature Physics 10, 39 review letters 109, 095704 (2012).
(2014). [65] D. B. Liarte, X. Mao, O. Stenull, and T. Lubensky,
[45] D. M. Sussman, O. Stenull, and T. Lubensky, Soft Mat- Physical review letters 122, 128006 (2019).
ter 12, 6079 (2016). [66] E. Lerner, G. Düring, and E. Bouchbinder, Physical re-
[46] D. Z. Rocklin, S. Zhou, K. Sun, and X. Mao, Nature view letters 117, 035501 (2016).
 27

[67] K. Bertoldi, V. Vitelli, J. Christensen, and (1933), [Link]

M. Van Hecke, Nature Reviews Materials 2, 1 (2017). pdf/4/3/249/18722594/249 1 [Link].
[68] J. Ortiz-Tavarez, E. Stanifer, and X. Mao, Physical Re- [79] A. Gabrielli, D. Garlaschelli, S. P. Patil, and M. Á. Ser-
view Research 5, L042001 (2023). rano, Nature Reviews Physics , 1 (2025).
[69] T. Castle, Y. Cho, X. Gong, E. Jung, D. M. Sussman, [80] M. Gori, G. Monfardini, and F. Scarselli, in Proceed-
S. Yang, and R. D. Kamien, Phys. Rev. Lett. 113, ings. 2005 IEEE International Joint Conference on Neu-
245502 (2014). ral Networks, 2005., Vol. 2 (2005) pp. 729–734 vol. 2.
[70] R. Penrose, Mathematical Proceedings of the Cambridge [81] M. Stern, D. Hexner, J. W. Rocks, and A. J. Liu, Phys-
Philosophical Society 52, 17–19 (1956). ical Review X 11, 021045 (2021).
[71] D. Z. Rocklin, New Journal of Physics 19, 065004 (2017). [82] J. W. Rocks, H. Ronellenfitsch, A. J. Liu, S. R. Nagel,
[72] S. Zhang, E. Stanifer, V. V. Vasisht, L. Zhang, and E. Katifori, Proceedings of the National Academy of
E. Del Gado, and X. Mao, Phys. Rev. Res. 4, 043181 Sciences 116, 2506 (2019).
(2022). [83] J. W. Rocks, A. J. Liu, and E. Katifori, Physical Review
[73] W. Voigt, Lehrbuch der kristallphysik:(mit ausschluss der Research 2, 033234 (2020).
kristalloptik), Vol. 34 (BG Teubner, 1910). [84] M. Stern and A. Murugan, Annual Review of Condensed
[74] D. Kondepudi and I. Prigogine, “Modern thermodynam- Matter Physics 14, 417 (2023).
ics: from heat engines to dissipative structures,” (John [85] S. Li and X. Mao, Nature communications 15, 10528
Wiley & Sons Inc., Chichester, West Sussex; Hoboken, (2024).
New Jersey, 2015) Chap. 17, second edition. ed. [86] S. Li and X. Mao, arXiv preprint arXiv:2503.07796
[75] L. Onsager, Phys. Rev. 37, 405 (1931). (2025).
[76] L. Onsager, Phys. Rev. 38, 2265 (1931). [87] G. Strang, “Introduction to linear algebra,” (Wellesley,
[77] J. M. Ortiz-Tavárez, Z. Yang, N. Kotov, and X. Mao, 2016) Chap. 10.
Phys. Rev. Lett. 134, 147401 (2025). [88] I. J. Busch-Vishniac, “Electromechanical sensors and ac-
[78] F. A. Firestone, The Journal of the tuators,” (Springer, New York, 1999) p. 20.
Acoustical Society of America 4, 249