Supplementary Information
Supplementary 1 - Energy, net moment and toroid moment
The magnetic energy U can be scaled withU0 =µ0 m20,4πd3

where d is a characteristic lengthscale (e.g. the edge length of the cube), m0a dipole moment scale and µ0 the vacuum permeability. The dimensionlessenergy E per dipole can then be defined asN
1 X mi · mj − 3 (mi · eij ) (mj · eij )
U,=
E=
N U0
N i<j|rij |3(1)where m1 , .., mN describe the dimensionless moments of the N freely orientable dipoles of equal magnitude |mi | = 1 and rij is the dimensionlessrelative position vector between dipole i and j with eij denoting the corresponding unit vector.
The net magnetic moment is defined as
N

M=1 X
mi ,
N i=1

which is a discrete analogue to the magnetization. The normalization of themoments |mi | = 1 guarantees that |M| ≤ 1. The toroid moment is definedwith respect to a given point c in space as
N

T=1 X
pi × mi ,N i=11

where pi is the vector from the point c to the position of dipole i, see e.g. [1].
The magnetic dipole moment m can be thought of as a current loop in thelimit where its area A goes to zero while its current I diverges, keeping theproduct |m| = AI constant. In the same sense, the toroid moment T can beenvisioned as a toroidal inductor coil in the limit where the torus radius R andcross-sectional tube area A go to zero while the current I diverges, keeping
the product |T| = RAI constant. T is useful to describe vortex-like statesof dipole configurations, like e.g. planar head-to-tail ring configurations, forwhich |T| is maximal. For the cube continuum we have M = 0 and T = 0
with respect to the center of the cube.

Supplementary 2 - System of equilibrium equations
Finding all equilibria of a dynamical system with energy E is equivalent
to finding all states of the system where the gradient vanishes. In general,the gradient has as many components as there are degrees of freedom inthe system, in our case 2N : Two angles per dipole to describe its orientation. However, we found that the use of cartesian coordinates to describe
the dipole orientations is advantageous compared to the angle formulation:
Firstly, we avoid the inherent coordinate singularities of spherical coordinates. Secondly, the degree of the resulting system of polynomial equations
is lower (quadratic instead of quartic). Thirdly, the structure of the equations
is highly symmetric.
Let the orientation of dipole i be given as mi = (xi , yi , zi )T , then we can2

define the full orientation vectorΩ := (x1 , y1 , z1 , ..., xN , yN , zN )T ,which contains the orientations of all dipoles. The energy E can now bewritten as a symmetric bilinear form E(Ω, Ω) with the representation
1E = ΩT W Ω23N

1 X i jE=Ω Wi Ωj ,2 i,j=1

orwhere W is the symmetric 3N × 3N interaction matrix which encodes thepositional information of all dipoles. W is constant for a given arrangement(like e.g. the cube). With the definition of the linear combinations
Lik :=3NXj
W3(i−1)+kΩj ,i = 1, 2, ..., N ,k = 1, 2, 3 ,j=1

the zero gradient equations (i.e. equilibrium conditions) for all dipoles canbe expressed in the form of the 3N cyclic equations
Li1 yi = Li2 xi
Li2 zi = Li3 yi

i = 1, 2, ..., N(2)

Li3 xi = Li1 zi .
Since we assume fixed magnitudes for the moments |mi | = 1 the cartesian
description demands N additional “sphere” equations
x2i + yi2 + zi2 = 1 ,i = 1, 2, ..., N .

(3)

Altogether Eqs. (2) and (3) set up a strongly coupled system of 4N quadratic
polynomial equations in the 3N unknown components of Ω.
3

Supplementary 3 - Upper bound for the number of isolated equilibriaTo get an upper bound for the number of possible isolated solutions, ouroverdetermined system Eqs. (2) and (3) with 4N equations for 3N unknowns has to be reformulated. We “randomize down” to a square system:
The system is replaced by 3N random linear combinations of the original 4N
equations, cf. [2]. Then, a variant of Bertini’s theorem (see e.g. [3] p.163)
assures that the set of isolated solutions for the randomized system is a superset of the original system. Therefore, every upper bound (for the numberof isolated solutions) for the randomized system is one for the original system
as well.
Now we use the basic version of Bézout’s theorem: A polynomial system with n variables and n equations of total degree d1 , .., dn has at most
Smax := d1 d2 ...dn isolated solutions. The randomized system with 3N variables consists of 3N quadratic polynomial equations (i.e. the total degree isalways 2). Therefore we have Smax = 23N possible isolated solutions, whichis also an upper bound for the original system (see above). Compared to theanalogous analysis for the angle formulation (Smax = 26N ) this is a significantly better bound.

Supplementary 4 - Numerical solution method
The system of polynomial Eqs. (2) and (3) is solved with the open-source
software Bertini™: Software for Numerical Algebraic Geometry [4]. This nu-4

merical software framework is developed for industrial problems including
mechanism and robot kinematics, chemistry, and computer-aided design, toname a few. However, it can be used to solve any system of polynomial equations. The manageable system sizes for numerical programs like Bertini™ aremuch greater compared to symbolic manipulation programs which utilizeGröbner bases and related methods, see e.g. [5].
The numerical algorithm in Bertini™ is based on homotopy continuation
and uses modern methods like randomization techniques and regenerative
cascades, see [6, 3]. The general idea of homotopy continuation in order tosolve a system Σ, is to consider another system Σ? with known solutions anddeform Σ? to Σ. The deformation is achieved through a parameterizationof the system equations, such that they reduce to Σ and Σ? for different
parameter values.
A parallel version of the software is available based on the Message Passing Interface (MPI). An overview of the theoretical background in algebraic
geometry, a detailed description of the numerical algorithm and a comprehensive user manual for Bertini™ can be found in book [3].
The software provides information about the number of solution components and their dimensions. Further, for 0D (isolated) components it gives
the numerical values of the solutions up to any required precision. For higher
dimensional components any number of sampling points can be generated.

5Supplementary 5 - Stability
We now prove the formula connecting the trace of the Hessian matrix H (thesum of all eigenvalues λk ) with the magnetic energy E2NX

λk = Tr (H) = −4E.

(4)

k=1

Here it is convenient to use spherical coordinates with the polar and azimuthal angles θi and φi to describe the orientation of dipole i xi
cos φi sin θimi =  yi  =  sin φi sin θi  .
zi
cos θi(5)

The energy E in Eq. (1) consists of sums of component-wise products of thedipole moments. If we define i to collect all terms of E containing θi and φiwe can write with Eq. (5)i := Ai cos φi sin θi + Bi sin φi sin θi + Ci cos θi ,(6)where Ai , Bi , Ci depend on all other angles but not on θi and φi . Note thatsince all terms of E consist of pairs of different moments, every term appears
in two different i and we haveN1X
i .
E=2 i=1(7)

We need the 2N diagonal entries of H, which are the second covariant derivatives of E with respect to all angles θ1 , φ1 , θ2 , φ2 , ..., θN , φN . We may replace
E with the respective i when differentiating with respect to θi and φi and6

get the diagonal entries∂ 2 1,∂θ12∂ 2 2H33 =,∂θ22..
.

H11 =1 ∂ 2 1
∂1+ cot,22∂θ1
sin θ1 ∂φ11 ∂ 2 2
∂2
H44 =+ cot,22
∂θ2
sin θ2 ∂φ2..
.
H22 =(8)The extra “cot . . . ” term when differentiating with respect to φ is a direct
result of the second covariant derivative in spherical coordinates. Inserting
Eq. (6) into (8) we find
H11 = H22 = −1 ,H33 = H44 = −2 ,...

and therefore
Tr (H) =2NX

Hkk = −2

NX

i .

i=1

k=1

This together with Eq. (7) proves relation (4). Note that the only condition
on E for Eq. (4) to hold is that E consists only of component-wise productsof two moments. There may be other types of interactions to which thisrelation is applicable.

Supplementary 6 - Robustness of the continuous groundstate
To demonstrate the existence of the continuous ground state as well as therobustness and practicability of the magnetic coupling, we actually designed
and built two different constructions of the cubic dipole cluster, see Figure 1.
In both experiments the theoretical finding that the ground state continuum
7Figure 1: Practical realizations of the cubic dipole cluster. Left: A construction which allows 8 Geomag™ toy magnets to rotate freely in their gimbalmounts. All white parts are produced with a 3D printer, the black rods are
made of carbon fiber and the bearings are realized with non-magnetic screws.
The centerpoints of the dipoles form a cube with an edge length of exactly
40 mm, see [7]. Right: A Teflon™ block with holes drilled in the corners,
allowing 8 spherical magnets to glide easily and therefore rotate freely. The
resulting magnetic forces are attractive and automatically fix the spheres in
the corners (i.e. the spheres are attracted towards the center). With a color
mark on the spheres (in red) one can see the motion of all spheres along thecontinuum if one rotates one of them around its respective space diagonal(done here with a glass tube with a rubber end).
is the only stable state is confirmed: Any other initial states relaxes immediately back to the ground state. The motion along the continuum is smoothwith no noticeable force variations already for these first two prototype constructions.
To quantify the influence of possible imperfections in the constructionof the magnetic coupling, we perform an analysis of the continuum for aperturbed system. For that we displace one of the dipoles in a random di8

0.004

relative energy variation

0.003

0.002

0.001

0.000
0.000

0.002

0.004

0.006

0.008

0.010

perturbation radiusFigure 2: Smoothness of the appearing soft mode for the perturbed systemcharacterized by the relative energy variation. The variation is shown as
a function of the perturbation radius by which one dipole is displaced in arandom direction. The radius is given relative to the cube’s edge length. For
each radius, 1000 random directions are computed. The largest variation ofabout 0.37% (for a radius of 1%) occurs for displacement directions perpendicular to the rotation axis, i.e. the space diagonal of the cube. For directions
along the rotation axis the variations are smallest. The upper limit for thevariation depends linearly on the perturbation radius.
rection with a given perturbation radius away from its original position. In
the considered range for the perturbation radius (it does not exceed 1% of
the cube’s edge length) the perturbed ground state remains the only stableconfiguration. Still the system can not jump to a different stable equilibrium.
The former continuum becomes a soft mode in the perturbed state. The important quantity characterizing how flat this soft mode will be, is the relative9energy variation, given as: The magnitude of the sinusoidal-like magnetic energy variation along the mode divided by the mode’s mean magnetic energy.
The dependence of the relative energy variation on the perturbation radius
is shown in Figure 2. The maximum variation depends linearly on the perturbation radius, showing that no catastrophic alterations in the dynamics
occur, the system is structurally stable.

References[1] Melenev, P, Rusakov, V, & Raikher, Y. (2006) Journal of Magnetismand Magnetic Materials 300, e187 – e190.
[2] Sommese, A. J & Verschelde, J. (2000) Journal of Complexity 16, 572 –602.
[3] Bates, D. J, Hauenstein, J. D, Sommese, A. J, & Wampler, C. W. (2013)
Numerically solving polynomial systems with Bertini. (SIAM).
[4] (2014) (https://bertini.nd.edu/).
[5] Decker, W, Greuel, G.-M, Pfister, G, & Schönemann, H. (2012) Singular 3-1-6 — A computer algebra system for polynomial computations(http://www.singular.uni-kl.de).
[6] Hauenstein, J. D, Sommese, A. J, & Wampler, C. W. (2011) Applied
Mathematics and Computation 218, 1240 – 1246.
[7] Grunwald, M. (2014) (http://audiots.wordpress.com/).
10