Spin-orbit coupling in a hexagonal ring of pendula

We consider the mechanical motion of a system of six macroscopic pendula which are connected with springs and arranged in a hexagonal geometry. When the springs are pre-tensioned, the coupling between neighbouring pendula along the longitudinal (L) and the transverse (T) directions are different: identifying the motion along the L and T directions as a spin-like degree of freedom, we theoretically and experimentally verify that the pre-tensioned springs result in a tunable spin-orbit coupling. We elucidate the structure of such a spin-orbit coupling in the extended two-dimensional honeycomb lattice, making connections to physics of graphene. The experimental frequencies and the oscillation patterns of the eigenmodes for the hexagonal ring of pendula are extracted from a spectral analysis of the motion of the pendula in response to an external excitation and are found to be in good agreement with our theoretical predictions. We anticipate that extending this classical analogue of quantum mechanical spin-orbit coupling to two-dimensional lattices will lead to exciting new topological phenomena in classical mechanics.


Introduction
The topological effects that underlie intriguing quantum mechanical phenomena, such as the quantum Hall effect, are not the prerogative of quantum mechanical systems, but have recently been observed also in classical systems governed by Newton's equations [1]. This has sparked particular research interest as the robust modes that are characteristic of topological systems [2] are especially promising in view of applications. In some types of classical mechanical structures, there can be topologically-protected zero-frequency modes [3,4,5]. These may find key applications in the emerging field of acoustic metamaterials used for controlled stress design, structural engineering and vibration isolation [1,6,7,8]. Whereas in analogue quantum Hall systems, topologically-protected finite-frequency edge states could lead to the implementation of an acoustic isolator, in which sound waves propagate along the edges of a structure without penetrating into the bulk and without being backscattered by system imperfections [9,10,11].
In classical analogues of the integer quantum Hall effect, the mechanical properties of a system should be designed so as to engineer topologically non-trivial phonon bands with non-zero Chern numbers. A first step in this research direction was a theoretical proposal for coupled pendula to simulate the Peierls phase factor of a charged particle hopping in the presence of a nonzero magnetic vector potential [12]. Since then, there have been many important theoretical and experimental works demonstrating how to engineer artificial magnetic fields and topological lattice models for classical systems, such as lattices of pendula [13,14], coupled gyroscopes [15,16] and acoustic crystals [17,18,19].
Alongside these advances in topological classical mechanics, the photonics community has also been strongly active in the theoretical study and the experimental realization of topological lattice models [20,21,22,23,24] and of a spin-orbit coupling for photons [25,26,27,28,29,30]. Spin-orbit coupling in such photonic systems has been predicted to induce various topological phenomena, such as topological phase transitions [30,31]. A detailed theoretical and experimental study of such a spin-orbit coupling for a hexagonal ring of exciton-polariton microcavities was reported in [29], where two polarization states provided the pseudospin degrees of freedom.
In the present work, we show that similar physics can be observed also in systems governed by Newtonian classical mechanics. Towards this goal, we theoretically and experimentally investigate an analogous mechanical model consisting of six pendula arranged in a hexagonal ring structure and coupled by springs. In particular, we observe that the frequency spectrum strongly depends on the ratio between the rest length of the springs and the equilibrium distance between neighbouring pendula. A mismatch between these two quantities results in a finite pre-tensioning of the springs which, as first anticipated by [32], is responsible for different effective spring constants along the longitudinal (L) and the transverse (T) directions with respect to the spring axis. This is analogous to how the coupling between neighbouring sites depends on the two polarization states in the polariton lattice of [29]. On this basis, the mechanical system can also be interpreted as being subject to an effective spin-orbit coupling, as we show in the following.
This article is organized as follows. In section 2 we introduce the mechanical system under consideration and we theoretically review the origin of the spin-orbit coupling term, first in an infinite honeycomb lattice, then in the hexagonal ring of pendula. For this latter case, we discuss in detail the oscillation eigenfrequencies and eigenmodes and we classify them in terms of the symmetry of the oscillation pattern. In section 3 we present the experimental setup and we summarize its geometrical details and physical parameters. Section 4 is dedicated to the presentation of our experimental results, where we show that a comparison with the predictions of a theory for simple pendula gives already a qualitatively good agreement in terms of both the eigenfrequencies and the symmetry of the modes. As we then discuss, this agreement becomes quantitatively excellent when we take into account, for example, the effects of the non-zero radius of the spheres making up the pendula. These results confirm the presence of spin-orbit coupling in a system of pendula coupled by pretensioned springs and show how the strength of this spin-orbit coupling can be tuned by adjusting the amount of pre-tensioning. Conclusions and perspectives for our work are finally discussed in section 5. Figure 1. Left panel: lateral view of a system of two coupled pendula connected by a pre-tensioned spring. The spring is pre-tensioned when its rest length 0 is smaller than the distance D between the hanging points of the two pendula. As a result, in the equilibrium configuration, the total elongation of the spring is such that D ≥ ≥ 0 , and the pendula have a non-zero angle with respect to the vertical direction. Right panel: top view of the pre-tensioned spring, subjected to a further displacement that changes the relative position of the masses respectively by δ L and δ T along the longitudinal (L) and transverse (T) directions with respect to the link direction in the horizontal plane.

Theoretical model
In this section we provide a short theoretical derivation of how the pre-tensioning of the springs gives rise to a spin-orbit coupling term in the equation of motion for the coupled pendula. As a first step, in subsection 2.1 we show how the pre-tensioning splits the transverse and longitudinal oscillation modes of a system of two pendula. Then in subsection 2.2 we will proceed with a review of the theory for the infinitely extended two-dimensional honeycomb lattice of pendula. Finally, in subsection 2.3 we will specialize the theory to a finite geometry with a benzene-like ring of six pendula as considered in the experiment. Throughout this section we will focus on the simplified case of simple pendula consisting of a point-mass m attached to a wire of length L, whose natural oscillation frequency is then ω 0 = g/L. Extension to a slightly more sophisticated model, taking into account the non-zero radius of the spheres making up the pendula used in our experimental setup, will be discussed in section 3.

System of two pendula
We start by considering a system of two pendula coupled with a spring of spring constant κ and rest length 0 smaller than the distance D between the hanging points, 0 < D. In the equilibrium configuration, the elastic force of the spring has to be balanced by the gravitational force and the tension of the wires: as a consequence, the spring is elongated to a length such that so that the pendula have some non-zero angle with respect to the vertical direction, as sketched in the left part of figure 1. In the following, we refer to this feature by saying that the elongated spring in the equilibrium configuration is pre-tensioned.
As was first anticipated in [32] for a system of masses and springs, such a pretensioned spring induces a splitting between the longitudinal (L) and transverse (T) degree of freedom of the two coupled pendula in the horizontal plane. In fact, when one pendulum is displaced from the equilibrium position by δ L and δ T along the two L and T directions, as shown in the right part of figure 1, the elastic energy stored in the spring grows as From this expression, we can see that the small oscillations along the L-and Tdirections around the equilibrium position are characterized by two different effective spring constants, while the cross coupling remains zero: (3) If the equilibrium length of the spring is exactly equal to the rest length = 0 , from (3) we get that the restoring force of the spring is restricted to the longitudinal direction κ L , while the force in the transverse direction is exactly zero, κ T = 0. In the generic = 0 case, the motion along the longitudinal and transverse directions experiences different spring constants κ T = κ L , where the transverse one κ T depends on the initial elongation / 0 and becomes stronger for more pre-tensioned springs.
It is then straightforward to obtain the oscillation frequencies along the longitudinal and transverse directions as: with the coupling frequencies Ω L,T defined as: The small oscillations of the system of two pendula show four eigenmodes. These modes include firstly a pair of eigenmodes where the two pendula oscillate in phase along either the longitudinal or the transverse direction, then a mode where they oscillate out of phase in the longitudinal direction and, finally, a mode where they oscillate out of phase in the transverse direction.
In the rest of this work, the oscillation of a given pendulum along the two directions will be considered as the two components of a pseudo-spin where, as we show in the next subsection, the difference Ω L − Ω T = 0 provides the spin-orbit coupling in a system of many coupled pendula.

Spin-orbit coupling in a honeycomb lattice of pendula
The general idea of a spin-orbit coupling arising from Ω L − Ω T = 0 is best understood in the theoretically simplest case of an infinitely-extended two-dimensional honeycomb lattice of coupled pendula. Due to the presence of the pseudospin degrees of freedom, our model is similar to the p x,y -orbital bands in a honeycomb lattice studied in [33] in the context of ultracold gases. Electrons in solid state graphene would instead correspond to a p z -orbital band model in which there is only one valence-bond orbital Figure 2. Sketch of the honeycomb lattice, whose sites are arranged on the vertices of the hexagons. We also show the indexing of the unit cells, each containing two lattice sites labelled as A and B and coloured in black and grey respectively. The two basis vectors α 1 and α 2 that generate the lattice from the unit cell are also indicated. For each link, orange and purple arrows represent the longitudinalê L i and transverseê T i unit vectors.
per lattice site. The honeycomb lattice is a Bravais lattice with two atoms per unit cell, labelled as A and B and separated by a distance D equal to the lattice spacing. The two generators of the lattice are α 1 = 3D/2, √ 3D/2 and α 2 = 3D/2, − √ 3D/2 , such that the whole lattice can be recovered from one unit cell by a translation of an integer multiple of the generators. In this lattice, pendula are assumed to be coupled to their nearest neighbours through pre-tensioned springs such that the spring restlength is smaller than the lattice spacing 0 < D. In an infinite system, the equilibrium positions of the pendula exactly reproduce the honeycomb geometry of the hanging points, meaning that the pendula hang vertically and that the equilibrium length of the springs matches the lattice spacing D. In a more realistic finite system, this configuration can be attained by applying suitable boundary conditions at the edges of the lattice, e.g. by keeping the position of the outermost pendula fixed.
We use the labelling shown in Fig. 2 and introduce the following unit vectors indicated by coloured arrows in the figure: We denote with a i,j = a x i,j , a y i,j ( b i,j ) the displacement from equilibrium of the pendulum located on an A-site (a B-site) in the unit cell i, j. Newton's equations of motion for the A-site pendula are the following: while the ones for B-site pendula are: The calculation of the normal mode dispersion is made easier by a Fourier transform to the momentum-space variables: where the integrals run over quasi-momenta in the Brillouin zone and V is the total area of the Brillouin zone. The Brillouin zone can be taken in the form of a hexagon and is delimited by the highly-symmetric points K and K . The extra term e iD k·ê L 1 accounts for the intra-cell distance between A and B sites. We project each of the equations in (7) and (8) along the x-and y-direction and The eigenvectors of the dynamical matrix correspond to the normal modes for a given momentum k.
The nature of the spin-orbit coupling becomes clearer in the circularly polarized +/− basis, to which we can transform by means of the unitary matrix: In terms of the transformed vector and where we have introduced J = Ω 2 The matrix in (12) can be expressed in terms of a spin operator acting on the pseudospin of the sublattice and another spin operator acting on the polarization degree of freedom. In our formalism, the ones acting on the sublattice degree of freedom read: while the ones acting on the polarization degree of freedom read: where I n is the n-by-n identity matrix and σ ± are defined as usual σ ± = (σ x ± iσ y )/2 from the 2-by-2 Pauli matrices. The two operators in (14) and (15) commute with each other. The matrix in (12) can be expanded around the K points k → (q x , q y − 4π/(3 √ 3D)), for | q| 1/D at the first order. After a gauge transformation ‡ The first term in (16) is a constant and the second term is the Dirac Hamiltonian of a p x,y -band graphene [33]; both of them are qualitatively unaffected by the pretensioning, i.e. by Ω T = Ω L . More interesting are the third and fourth terms (the ones on the last line) which introduce the effective spin-orbit coupling effects. The effect of the former term is similar to that of a Rashba spin-orbit coupling, as discussed in [34], while the latter gives a trigonal warping effect [35]: both of them are proportional to the ∆ parameter quantifying the difference between Ω L,T .

The hexagonal ring
We conclude this theoretical section by discussing the effect of spin-orbit coupling in a benzene-like geometry consisting of a ring of six pendula arranged at the vertices ‡ We perform this transformation in order to recover the Dirac-like Hamiltonian in the second term of (16). of a regular hexagon. As in the previous sections, the springs are assumed to have a rest length 0 shorter than the distance D between the hanging points. In the present spatially-finite geometry, the pre-tensioning means that the equilibrium positions of the pendula lie on a hexagon of reduced side < D, so that they make a non-zero angle with respect to the vertical direction as shown in the left panel of figure 1 and in the left picture of figure 6.

Equations of motion and eigenmodes
In order to write Newton's equations of motion for the system of pendula, it is useful to separate the motion along the L and T directions and use the frequencies defined in (5). For this purpose, we define unit vectors parallel and orthogonal to the direction of the links, as sketched with coloured arrows in figure 3. The explicit form of the longitudinal vectors is: while that of the transverse vectors is: Expanding the elongation of each spring in this basis, Newton's equations of motion for the i = 1, . . . 6 pendulum take the form: where ψ i = (ψ x i , ψ y i ) and ψ x i , ψ y i are the displacements of the i-th pendulum in the x − y directions. Periodic boundary conditions are applied in the form i + 1 → 1 for i = 6 and i − 1 → 6 for i = 1.
We solve the eigenvalue problem, searching for a solution of the type ψ i (t) = ψ i e iΩt . We can cast the equation in (17) in a matrix form with a state vector either  in the x, y basis Ψ xy = (ψ x 1 , ψ y 1 . . . ψ x 6 , ψ y 6 ) or in the circularly polarized +/− basis Regardless of the basis that is used, the system in (17) is: − from which diagonalization of a 12 × 12 dynamical matrix D gives the frequencies of the eigenmodes. In the general case of Ω T = Ω L , the twelve eigenmodes are grouped into a set of eight different frequencies: 2.3.2. Symmetry classification of the eigenmodes To better understand the properties of these eigenmodes, we can make use of the classification in terms of their total angular momentum first introduced for polaritons in [29]. To see this, we consider the transformationT which leaves the system invariant and which combines a translation T that sends the i-th site to i + 1-th together with a rotation R π/3 of angle π/3. For a plane wave of wavevector Q ∈ (−π, π] around the ring (corresponding to an orbital angular momentum l ≡ Q/(2π/6) ∈ (−3, 3]) and uniform circular polarization S = ±, we have: As the rotation operator in the ± basis acts as R π/3 Ψ S = e iSπ/3 Ψ S , we have that: and so we can identify k = l + S as the total angular momentum.
As this total angular momentum k is conserved in our rotationally symmetric system, we can make use of the following Fourier-like transformation to put the system in (18) in a block diagonal form, with each block being a 2by-2 matrix acting on the sub-space with a given value of k. The result of this diagonalization exactly recovers the formulas in (19) and is graphically illustrated in figure 4 and figure 5. Further details can be found in [29]. Panels (a)-(c) in figure 4 show the frequencies of the system as a function of the total angular momentum k for Ω L /ω 0 = 1 for the three cases Ω T /Ω L = 0, 0.5, 1. In open dots we show the eigenvalues obtained for the discrete integer values k ∈ (−3, 3], while the solid lines are just guides to the eye. The eigenstates, shown as open dots, are labelled according to the analytic expression of their frequencies given in (19): the degenerate states are distinguished by the different value of the total angular momentum, which prevents their mixing by symmetry as long as the system is rotationally invariant.
The weak value of the spin-orbit coupling in the polariton experiments in [29] restricted that investigation to the case where |Ω T − Ω L | Ω L,T . In this limit, only four different eigenfrequencies can be clearly spectrally distinguished, as shown in figure 4(c). In the mechanical system, this regime is instead hard to access as it requires springs with a very small rest length, 0 / → 0. The opposite Ω T /Ω L ≈ 0 regime is instead realized when the spring rest length is equal to the distance between the pendula 0 = = D, and there is no pre-tensioning in the equilibrium configuration. This regime is illustrated in figure 4(a). In polariton systems, this regime is achieved in a lattice of pillars when the two degenerate p-wave orbitals of micropillars are considered instead of the single s-wave one used in [29]. In this p-wave system, the tunneling amplitude of orbitals aligned orthogonal to the link is in fact strongly suppressed with respect to the one of transverse orbitals. While this physics was not specifically addressed in [29], the large difference between the coupling amplitudes is apparent in the flatness of the higher p-wave bands of the honeycomb lattice studied in [36].
As a general remark, we point out another difference between the polariton system and the classical system of pendula. In the polariton system of [29], the Hamiltonian is chirally symmetric and, as a result, the eigenfrequencies are symmetrically located around the bare cavity frequency. Chiral symmetry for our pendula system is instead broken and the coupling induced by springs in systems of pendula always increases the frequency of the oscillation modes. This is due to additional terms in (17) where the elastic force acting on the i-th pendulum depends on the position ψ i of the pendulum itself. In the presence of spin-orbit coupling in a hexagonal ring, this results in there being different diagonal elements in the dynamical matrix (18) for the x-polarized mode and y-polarized mode of a given site, breaking chiral symmetry.
In figure 5 we show the motion of the pendula in each of the twelve eigenmodes. This is obtained from the eigenvectors of the dynamical matrix D, corresponding to the eigenfrequencies of (19), and its form does not depend on the specific value of Ω T /Ω L . The panels in the left part of the figure show the theoretical eigenmodes, while in the right part we show a comparison with the experimental modes as obtained from the data analysis. The motion of the pendula is represented within a period of  oscillation around the equilibrium positions, and the colour gradient represents time, where a darker colour stands for earlier time. For the degenerate eigenfrequencies, the experimentally excited mode is a linear superposition of the two modes, whose weights are determined by the specific excitation procedure. More details on these experimental eigenmodes will be given in the following.

Experimental setup
In figure 6 we show a picture of the experimental setup of six coupled pendula connected with pre-tensioned springs. Each pendulum is realized by a sphere of mass m = (0.596 ± 0.001) Kg and radius R = (2.65 ± 0.05) cm attached to a string of length L = (16.0 ± 0.5) cm, hanging from a vertex of the hexagonal transparent plastic roof. As visible in figure 6, the top plastic roof is pierced with sets of holes that allow us to hang the hexagonal system of pendula at five possible distances D, as reported in table 1 from the outermost to the innermost. Under the assumption that the pendula are made of point-like masses attached to a wire of length L + R, their natural frequency is expected to be: where we have used the geometrical parameters given above. The pendula are coupled through springs of rest length 0 = (7.30±0.01) cm. The springs are attached to a hook that is located at the bottom of the sphere. The distance D between the hanging point of the pendula in the five configurations is always larger than the rest length of the spring and, as a result, the springs are pre-tensioned. In the equilibrium configuration, the pendula move in with respect to their hanging points, as visible in the left part of figure 6. The values of for the five configurations is reported in table 1. We characterize the springs by measuring their elongations when subjected to a known force and extracting the constant of the spring as the slope of the resulting curve. We observe that, for extensions ∆ = ( − 0 ) > 15 mm, the springs behave nonlinearly. However, provided that the extension does not exceed ∆ ≈ 35 mm, the elastic response of the springs can still be locally approximated as linear, with a modified spring constant κ eff . When ∆ > 40 mm, the spring is instead in the plastic regime where the deformations are permanent. For each of the experimental configurations, the modified value of the spring constant κ eff , which takes into account the nonlinear behaviour, is given in table 1. The value of the modified spring constant κ eff is obtained from a local linear fit of the force-extension curve of the springs in the relevant range of extensions. In table 1 we also report the corresponding ratio Ω T /Ω L of the longitudinal and transverse frequencies defined from (5) for the experimental parameters.
To excite the system, we displace by hand one of the pendula from its equilibrium position and suddenly release it, as shown in the supplementary video. The local nature of the initial condition means that the subsequent motion of the pendula is a superposition of all the eigenmodes of the system depicted in figure 5. A different choice in the direction of the initial condition would only affect the relative weight of each mode.
The motion of the pendula is recorded with a standard video camera, positioned above the system. This motion is sampled at a frequency ν s = 25 Hz, for a total time T of about 120 s. The video is then digitally analysed to obtain the displacements of the centre of mass of each pendulum. In order to facilitate the measurement of the position of the pendula, white circles of paper have been rigidly attached to the top of the spheres, as shown in the right panel of figure 6, in order to extract the position of the pendula through the analysis of the data, see the supplementary video. In figure 7 we show a typical result for the experimental motion of the six pendula as obtained for the parameters of the fifth configuration in table 1. We then perform a Fourier transform on each set of data to obtain the Fourier amplitudes for both the x and the y components: |F xi (Ω)| and |F yi (Ω)|. The frequencies span from Ω ∈ [−πν s , πν s ], with a step of ∆Ω = 2π/T .

Results and discussions
In figure 8 we show the frequency spectra as calculated from the total Fourier amplitude, for frequencies in the region of interest. Each spectrum has been normalised in such a way that the integral over the whole frequency range is equal to one. Panels   (19) calculated from the experimental uncertainties in the parameters. The relative height of the peaks depends on the initial condition. A different initial condition will excite another superposition of modes, each of them with a different coefficient and hence with different spectral heights, but the frequencies of the peaks are not affected by the initial condition.
In figure 8, we observe a qualitative overall agreement between the experimental spectra and the theoretical predictions. In particular the low frequency peaks in figure 8(c)-(e) match the theoretical predictions well within the experimental error. On the other hand, appreciable deviations are visible for the peaks around 11 rad/s for all the panels. A convincing explanation for these discrepancies will be given in the following of this section.
Before entering into this discussion, it is useful to look at the Fourier transform of the displacements, which allows us to reconstruct the oscillation amplitude pattern of the eigenmodes from the F i (Ω) evaluated for Ω located at a peak. The experimental eigenmodes are plotted in the right part of figure 5, labelled as (Exp) and ordered, from bottom to top, according to the increasing value of the corresponding eigenfrequency. The oscillation patterns of the experimental eigenmodes are in excellent agreement with the ones of the theory depicted in the left part of figure 5 and discussed in section 2.3, especially regarding the symmetry of the pattern. In particular, we notice that the eigenmodes associated with eigenfrequencies Ω 3 and Ω 4 present an azimuthal symmetry for the pattern of oscillation, while the eigenmodes associated with eigenfrequencies Ω 2 and Ω 5 present a radial symmetry pattern. Such well-defined radial and azimuthal patterns are peculiar of spin-orbit coupled systems [29]. The remaining oscillation patterns are the result of a linear superposition of the doubly degenerate eigenmodes: while the oscillation patterns of peaks Ω As a last point, we wish to shine light on the deviations visible in figure 8 between the theoretical frequencies derived in section 2 and the experimental spectra. We can interpret them as a consequence of the finite size of the sphere making up the pendulum and, more precisely, of its rotation around the hook that connects it to the string. Such an additional degree of freedom is in fact not included in the simple model discussed in section 2 and gives extra oscillation modes at higher frequency. To verify this hypothesis, we numerically simulated the system by solving Euler-Lagrange equations that take into account the effect of a non-zero radius R of the mass and its rotation around the hook. More details on such an approach are given in Appendix A.
As a first consistency check, we have verified that the full numerical simulation well reproduces the eigenmodes (19) in the limit of R → 0. As the rotational symmetry of the system is the same in the extended approach, we expect that the symmetry of the mode oscillation patterns is unchanged: this has also been successfully checked for different values of R.
Finally, when R is taken to be equal to the actual experimental value, the resulting eigenfrequencies are found in excellent agreement with the experimental ones. This is displayed in detail in figure 9, where we plot the frequencies of the eigenmodes as a function of the ratio Ω T /Ω L . Dots are the experimental eigenmodes as obtained from the peaks in the spectra of figure 8 for different hanging positions and therefore different values of the spring pre-tensioning. Solid lines are obtained for the theoretical frequencies in (19) for the simple pendulum model for point masses pendula using experimental parameters. Dashed lines are instead the eigenfrequencies obtained from the full numerical simulation including the non-zero radius of the masses and their rotation around the hook.

Conclusions
In this paper we have given experimental evidence for a tunable spin-orbit coupling in classical mechanics using a system of six coupled pendula arranged in a hexagonal geometry and connected by pre-tensioned springs. The experimental results for the oscillation frequencies and the oscillation patterns are compared with theoretical models: while the qualitative agreement with a simple pendulum approximation is already quite good, it becomes quantitatively excellent once we take into account the finite radius of the masses and their rotation around the hook connecting them to the string. By changing the hanging position of the pendula, we have demonstrated how the strength of the spin-orbit coupling is tunable just by varying the amount of pre-tensioning of the springs.
Future developments will include extending the experimental study of the spinorbit coupling to larger two-dimensional lattices of pendula. In this case one could study the topological Lifshitz transition as theoretically proposed in [32], and observe the creation, motion and annihilation of Dirac cones. Moreover, the simulation of an artificial magnetic field for this simple mechanical system is straightforward to realise experimentally by mounting the system on a rotating table [32,37] or by adding a spatially inhomogeneous strain in a honeycomb geometry [38], so to study the interesting interplay of orbital magnetic effects with the spin-orbit coupling induced by the pre-tensioned springs.
Note: Immediately before the submission of the paper, a related theoretical work appeared on spin-orbit coupling in mechanical graphene [39]. Figure A1. Sketch of the coordinates used for representing the motion of the pendula. On the left, a lateral view of the i−th pendulum. θ i is the angle formed by the string with the vertical direction. The distance along the axis between the centre of mass of the sphere and the center of the hook connecting it to the string is R , which is greater than the radius R. ϕ i is the angle between the axis of the sphere and the vertical direction. On the right, a top view of the i−th pendulum. The suspension point of the pendulum is identified by the coordinates (x 0i , y 0i ). The azimuthal angles α i and β i are also indicated.
We include in our theory the effects of a non-zero radius of the sphere making up the pendula. In particular, we have to consider that the centre of mass of the sphere does not coincide neither with the top hook where the string is attached, nor with the bottom hook where the coupling springs are connected. These three points are instead aligned along a direction which defines the axis of the sphere.
The sphere can also freely rotate around the top hook, meaning that the axis of the sphere makes a non-zero angle with the direction defined by string of the pendulum, as visible in the left part of figure A1. Also visible in figure A1 is the radius R , defined as the distance between the centre of mass of the sphere and the top hook. Such a distance is larger than the radius of the sphere R, because it includes the hook itself. Since the two hooks have the same length, R also indicates the distance between the sphere's centre of mass and the bottom hook where the springs are attached. The experimental value of this distance is R = (2.85 ± 0.05) cm.
The i-th pendulum is then represented by four coordinates: θ i , ϕ i , α i , β i , as schematically shown in figure A1. Two coordinates, θ i and α i , are used for defining the position of the point where the string is attached to the sphere. The other two, ϕ i and β i , define the position of the centre of mass of the sphere. Moreover, θ i measures the angle between the string of the pendulum and the vertical direction, while ϕ i measures the angle between the axis of the sphere and the vertical direction, as visible in the left part of figure A1. α i measures the angle between the projection of the string on the x − y plane and the x−axis, while β i is measured between the projection of the axis of the sphere on the x − y plane and the x−axis. The two azimuthal angles α i and β i are indicated in the right part of figure A1.
We can then express the position of the centre of mass of the sphere with these spherical coordinates: x cm i = x 0i + L sin θ i cos α i + R sin ϕ i cos β i y cm i = y 0i + L sin θ i sin α i + R sin ϕ i sin β i z cm i = z 0i + L cos θ i + R cos ϕ i . It is then straightforward to write the kinetic energy of the i-th pendulum as: where we have considered the rigid body rotation of the sphere around the joint with the string. Thanks to the bottom position of the hook connecting the sphere to the springs we do not need to include the rotation of the sphere around its axis, as this rotation is decoupled from the motion of the other pendula. The potential energy of the i−th pendulum is defined as U pot i = −mgz cm i . To write down the expression for the elastic potential energy, we have instead to consider the position of the bottom of the sphere, where the springs are attached: x s i = x 0i + L sin θ i cos α i + 2R sin ϕ i cos β i y s i = y 0i + L sin θ i sin α i + 2R sin ϕ i sin β i z s i = z 0i + L cos θ i + 2R cos ϕ i .

(A.3)
With these coordinates, the elastic potential energy of the spring that connects the i−th pendulum with its nearest-neighbour j−th pendulum is: We then write down 24 Euler-Lagrange equations starting from the Lagrangian: where the sum over j is done for the two neighbouring pendula of the i−th pendulum.
By numerically solving the Euler-Lagrange equations with the experimental parameters we obtain the dynamics of the pendula system. We then Fourier transform the solutions to obtain numerical frequency spectra. We observe two sets of frequency peaks. The first set is located at low-frequency (6 − 16 rad/s), while the second set is at higher frequency (35 − 50 rad/s). The low-frequency peaks of the numerical spectra are used in figure 9 to draw the numerical dashed lines. As expected, they correspond to modes where the motion of the whole sphere follows that of its upper hook. On the other hand, the high-frequency modes correspond to oscillations where the spheres have a significant rotation around the top hook, so that the motion of the center of mass is somehow out of phase with that of the top hook. In our set-up these latter modes have a lower spectral height and are more damped than the low-frequency modes. Although in figure 8 we focus on the low-frequency part of the spectra, both sets of peaks are however visible in the experimental spectra and are found in good agreement with the numerical results.