Collective dipole reorganization in magnetostructures

Neodymium spherical magnets are inexpensive objects that demonstrate how dipolar particles self-assemble into various structures ranging from 1D chains to 3D crystals. Assemblies of these magnets are nicknamed magnetostructures and this paper focuses on a variety called magnetotubes, which are some curved square lattices forming cylinders. We experimentally and numerically observe that such magnetotubes can self-buckle, above a critical aspect ratio. In fact, the underlying dipolar ordering of such structures is found to exhibit a collective reorganization, altering the mechanical stability of the entire system. We identify the conditions in which these phenomena occur, and we emphasize that metastable states coexist. This suggests that a wide variety of magnetostructures, including chains and magnetocrystals, may collapse due to the coexistence of multiple ground states and global reorientation of dipoles.


Introduction
The self-assembly of dipolar particles has emerged as a compelling topic over recent decades, largely due to its promising implications in phase separation [1] and the formulation of novel materials like photonics crystals [2]. One remarkable characteristic of magnetic particles, notably magnetic colloids [3], lies in their ability to be remotely manipulated via an externally applied magnetic field. Recent breakthroughs illustrate that metachronal waves [4] and locomotion [5] can be effectively initiated within these magnetic self-assemblies.
Magnetic balls are the macroscopic counterpart of these microscopic particles, and they share the same physical behavior. They can be manipulated to reveal the physical nature of dipole-dipole interactions. By assuming that such magnets are uniformly magnetized, neodymium spheres behave as point-like dipoles [6,7]. The interaction energy U αβ between two point-like dipoles ⃗ µ α and ⃗ µ β is then given by where⃗ r αβ =⃗ r β −⃗ r α is the vector joining particle centers α and β and µ 0 is the vacuum magnetic permeability. Looking carefully at equation (1), one notes that, depending on relative positions and orientations, two dipoles can attract or repel each other. Minimizing U αβ is therefore achieved when both spheres are in contact and their dipoles are aligned which is why such magnets are naturally forming chains [8,9]. Nevertheless, cubic crystals [10][11][12] as well as tubular [13] and planar structures [14] can also be created using such permanent magnets. Such assemblies are called magnetostructures and some of them are illustrated in figures 1(a)-(c): rings and tubes. In order to compare the cohesive energy of structures made with different amounts of beads, we introduce the dimensionless energy u per particle defined by where µ = |⃗ µ|, D is the sphere diameter i.e. the distance between two contacting beads in equation (1), and N is the total amount of beads in the structure. When creating magnetostructures, one realizes rapidly that the underlying dipole arrangement is as significant as the contact network between particles. This has been explained using equations (1) and (2) [9][10][11][12][13][14][15][16][17][18].
Recent experimental and numerical studies highlighted some remarkable features of magnetostructures, feeding the growing interest about this topic. As already recalled, dipolar particles naturally self-assemble into 1D chains. Those latter can be bent in order to create a ring that minimizes cohesive energy u compared to a chain for the same amount N > 5 of beads. Moreover, defects can be created in the perfect ordering of dipoles. These defects act as attractive or repulsive magnetic monopoles [9]. Nonetheless, more complex structures have also been investigated. In a 2D square lattice, dipoles orients in chains that are antiparallel to their neighbors. Conversely, in a triangular lattice, dipoles orients in chains that are parallel to their neighbors. These latter configurations are both the most compact and the more stable ones compared to the square lattice [14,15]. Please note that folding such 2D lattices leads to the formation of magnetotubes as shown in figure 1. A study of 3D magnetocrystals reveals also some unusual physical features: the highest cohesion is obtained for a hybrid structure made for a stack of hexagonal planes while the hexagonal compact (hcp) structure, known as the densest structure, is less cohesive and nearly unstable [16][17][18]. Moreover it has been proved that there exists an infinite number of stables configurations for the orientation of 8 dipoles at the vertices of a cube, similar to Goldstone modes [10,11]. These multipolar magnetic cubes may also interact with extremely short range interactions [12].
However when large and more complex magnetostructures are built by adding particles one by one, the system may suddenly self-buckle or collapse into a disordered assembly, for some reason. In other words, dipoles can not only change their orientation, but also translate and form a new geometric structure. The latter observation was the starting point of our work. The goal of this paper is therefore to study the dipolar arrangement in a growing magnetostructure, emphasizing such instabilities. First, experimental observations and measurements are given. Then, we introduce the numerical methods used to capture these results. Finally, we discuss all results and confront them with other systems.

Experiments
In our experiment, we use neodymium (NeFeB) magnetic sphere of diameter D = 5 mm and a remanent field B r ≃ 1.2 T. The magnetic momentum can be retrieve by the relation where V is the volume of the sphere. The coating of the beads is made of an alloy of nickel and copper. Those values are also used in the simulations described in next section.

Building magnetotubes
Dipoles in a ring orient head to tail forming a loop with all dipoles tangent to the ring curvature. By stacking two rings to form a tube of square lattice, the respective loops are in opposite directions, i.e. one loop of dipoles is turning clockwise and the other one is turning anti-clockwise. By stacking m rings made of l beads, it is possible to form a magnetotube with a curved square lattice as illustrated on figures 1(b) and (c). Each magnetotube is therefore determined by the numbers ℓ and m, and the total amount of beads composing the tube is simply N = ℓm and we adopt the notation (ℓ × m). We consider rings as the building blocks of our tubes.
Parameter values from ℓ = 5 to ℓ = 14 and m = 1 to m = 60 have been considered. Therefore, the aspect ratio of the magnetotube can be either ℓ/m < 1 (long tube) or ℓ/m > 1 (short tube). In order to avoid the collapse of long tube due to gravity, they are lying onto a horizontal rigid plate. By stacking the rings on the extremities of a magnetotube, one may consider that all dipoles are tangent to the ring curvature, even for long magnetotubes. We will see below that the situation is more complex than expected.

Acoustic measurements
A long tube is initially assembled and placed horizontally. One extremity is pinched and rolled in order to change the orientations of the beads. The tube is not crushed, the geometry barely changes, but the orientations do. This is often enough to provoke transition. One can hear a creaking sound before the elastic properties changes such that the entire structure becomes less rigid. After the sound production, the tube becomes more flexible and can easily be squeezed into a parallelepiped. This is illustrated in the first video of supplementary materials. The origin of emitted sound is due to the solid friction between rotating beads. If beads are rotating, one understands that their dipoles are reorienting. A microphone is placed close to the tube before disturbing the dipoles at one extremity. In figure 2, typical resulting sound patterns are presented for two different tube lengths and for ℓ = 6. Plot presents the sound amplitude in arbitrary units as a function of time, assuming that the disturbance is done at t = 0. For larger tubes, transition is harder to provoke and we sometimes need to roll all beads. This was the case for ℓ = 7 in the second video.
From figure 3, one can see that the duration of the whole creaking sounds is increasing linearly with the tube length m, meaning that the sequence of dipole rotations is propagating from the disturbed extremity to the other one at constant speed. As previously mentioned, larger tubes are harder to disturb and are therefore not illustrated. The acoustic measurements prove that a global reorientation of the dipoles, i.e. sphere rotations, takes place when the magnetotube is disturbed. While long tubes (m > ℓ) show bead rotation patterns and a drastic modification of their mechanical properties, short tubes (m < ℓ) do not experience such an event. The relevant question is therefore to find the physical origin of this global reorientation of dipoles and why it happens only for long tubes. In the next section, simulations will allow to look at the dipolar orientations during such an event.  . Three pictures of the same (7 × 10) magnetotube that has self-buckled. One can clearly see the triangular lattice formed by the two parallel chains along the tube axis. One also sees that contacts are lost between beads located at the extremities.

Self-buckling
The previous experimental results were obtained when long tubes are built. In general, this can be achieved when ℓ is an even number. For an odd ℓ number of beads composing successive rings, the above dipole reorientation provokes the collapse of the entire structure and the magnetotube is distorted. Three pictures of a (7 × 10) magnetotube are shown in figure 4, illustrating this self-buckling. After such an event, the local structure seen along the lateral curved surface is a triangular lattice instead of a square lattice. Rows of beads modified their relative positions such that the entire structure is irreversibly changed. This phenomenon is also illustrated in supplementary materials.
The fundamental question is therefore why this collapse appears only for odd ℓ values. In order to tackle this second major observation, we will in the next section perform simulations to reveal the dipolar organization during the above phenomena.

Numerical approaches
In order to capture the above observation, we use two different numerical approaches to simulate the magnetotubes. They are described below.

Simulated annealing (SA)
Stochastic optimization algorithms such as annealing, can be used to find the ground states of equation (2), as proposed in [9]. However, this approach requires long computation times due to the high number of degrees of freedom. Since each dipole α has five degrees of freedom, being the positions (x α , y α , z α ) and angular orientations (θ α , φ α ), minimizing the cohesive energy u, given by equation (2), requires to explore a complex configuration space with 5N degrees of freedom. In order to speed up our simulations, we consider that the positions of the dipoles are fixed, thus reducing configuration space to 2N degrees of freedom.
The structure is first initialized in a specific geometry where dipoles are oriented randomly. The dimensionless energy u as defined by equation (2) is then computed taking all dipoles into account. It is stored in a variable u min . This corresponds to the minimum value found so far. In addition, all dipole orientations corresponding to this energy are stored in an array. Then, we change the dipole's orientations in a random direction with some amplitude ∆θ and ∆φ between 0 and A. This latter is the maximum amplitude and is chosen large enough at the beginning of the simulation (typically ∆θ = ∆φ = 16 • for our simulations). The energy u is computed and compared to u min . If u min < u, the newly computed orientations are rejected, and we start back from the orientations corresponding to u min . Otherwise the newly computed energy u is saved as u min and dipoles orientations are stored replacing the old ones. This procedure is iterated thousands times before dividing ∆θ and ∆φ by a factor of 2. The simulation is stopped when ∆θ = ∆φ < 1/2. With this method, a local minimum of u is finally reached. In order to get the global minimum of energy, several simulations are performed in parallel using Open MPI over hundreds processors.

Discrete element method (DEM)
Of course, minimizing energy u does not provide information about the mechanical stability of the dipole arrangement. For exploring these characteristics, DEM simulations were realized. Our home made algorithm [19,20] considers contact and magnetic forces. Conversely to SA, the system has 5N degrees of freedom, so each bead can freely move in space. Moreover, DEM allows us to follow the dynamics of the system whereas the SA does not. Starting from any initial configuration, the system evolves towards mechanical equilibrium by moving and reorienting the magnetic spheres. At each time step of about 10 −6 s, the system computes all forces acting on beads, and update their positions and velocities solving Newton's equations of motion using Leapfrog integration scheme.
Contacts are modeled taking tangential and normal forces. Two beads are in contact if there is a tiny overlap between them. A repulsive normal force is modeled by a linear spring and a dashpot for dissipation with the relation where k N is the spring stiffness, δ is the overlap of both beads, γ N is the viscous damping and v N is the relative normal velocity [21]. Of course if the considered beads are not in contact, F N = 0. The tangential force is also modeled by a spring dashpot model, but sliding regime is added. Indeed, a criterion called Coulomb's law that couple normal and tangential forces is implemented as follow where µ s is the static friction coefficient and is chosen between 0.6 and 0.9. If tangential force satisfies equation (5), then it is a static regime. Otherwise it is a sliding regime and the tangential force is bounded by µ s F N . In summary, DEM provides the dynamics of a structure from specific initial conditions. While SA is able to find the minimum of u, DEM provides information about its stability against fluctuations. They are complementary methods.

Dipole reorientation in growing magnetotubes
Consider an arbitrary magnetotube (ℓ × m). The magnetic energy of each ring made of ℓ dipoles, is found to be minimized when each dipole is tangent to the circle, as calculated exactly in a previous study [9]. Consequently, one expects a circular orientation, even when m → ∞. Let us test this hypothesis by looking at simulations results.  However, above the critical size m c , the dipoles align with the tube's axis and a global reorganization of all dipoles, as described in [22], takes place. In figure 5, the SA curve diverges from the circular case as m increases. Furthermore, the energy u does not reach the case featuring only dipoles oriented along the tube's axis, since boundary effects are seen at the extremities of the magnetotube, as illustrated in figures 5(d) and (e). Nevertheless, extremely long tubes (1/m → 0) tend toward the same energy u for annealing and fixed dipoles along tube's axis since boundary effects become negligible in regard to the tube's length. Open diamond and dot lines meet for 1/m → 0 around u = −2.55 on figure 5(a). This asymptotic value is slightly larger than the cohesive energy for a square lattice [16].
Our numerical results confirm therefore the existence of a global reorientation of all dipoles in magnetotubes. Size and parameter effects will be investigated in the next subsections.

Metastability
SA simulations predict a reorientation of dipoles if the tube is long enough. Nonetheless, experiments described in section 2 shows that long tubes with circular orientation can be assembled, transition only occurring if we disturbed the system. This mismatching between simulations and observations is the evidence of metastability in such a system. It also evidence that SA lacks of an important ingredient to fit the observations. To figure out what happen in those systems, we increase the size of a short tube by adding rings one by one in a DEM simulation which takes friction into account. We observe the transition for a critical size m larger for DEM than for SA.
Moreover, by removing rings one by one, one can reach a critical height for which the dipoles will orientate along loops again. This transition is not observed at the same aspect ratio whether one increases or decreases the tube's length. This hysteretic transition is illustrated for a (10 × m) tube on figure 6. Energy is plotted as a function of 1/m, showing transition as we increase or decrease the height. This hysteresis is a clear signature of metastability, emphasizing the discontinuous nature of the dipolar reorientation.

Self-buckling instability
Additional SA and DEM simulations have been done to study the dipolar instability for other sizes of magnetotubes. Figure 7 shows similar numerical results obtained by SA for a various number ℓ of magnets for the rings constituting the magnetotube. The results are plotted as a function of ℓ/m. The previous results (ℓ = 10) shown in figure 5 are still plotted for comparison using the same color and symbol. For all magnetotubes, the global reorientation of the dipoles is observed when the aspect ratio, ℓ/m, is in between 1 and 1.5. For small tubes (large aspect ratio), on the right of figure 7, at a given value of m, the magnetic energy u decreases as ℓ increases. However, once the reorientation of the dipoles is obtained, for long tubes (smaller aspect ratios), the different curves cross each other. In particular, the energy u of odd ℓ curves are systematically above the ℓ − 1 and ℓ + 1 ones. This is due to the frustration between dipole chains. Indeed, an odd number of beads per rings implies that two neighboring chains along the tube's axis present a parallel alignment, which is not favorable in a square lattice. This results in an increase of the energy u. That can be seen in figure 7(b) where two parallel chains are obtained for ℓ = 7 in a SA simulation. The magnetic energy is therefore much higher than for a pattern with a perfect antiparallel ordering. In fact, this frustrated state is unstable and DEM simulations provide some insights of the time evolution of such collapsing magnetotubes, as seen in figures 7(c)-(e). In addition to a global reorientation of the dipoles, a geometric restructuring of the tube occurs, and parallel chains takes a triangular lattice configuration which minimize the energy u in this case: the tube has self-buckled. On figure 7(a), one observe that when ℓ increases, the transition is less and less pronounced, meaning that the energy difference between two states becomes less important.

Conclusion
In summary, our experimental and numerical work evidences collective dipole reorientation in curved square lattices and in particular in magnetotubes with increasing aspect ratios. We propose numerical approaches providing explanation for this phenomenon driven by energy minimization. On the one hand, SA allows us to find the orientations of dipoles that minimize magnetic cohesive energy. DEM on the other hand, allows us to study the stability and the dynamics of such a system. Numerical simulations show that friction seems to delay the reorganization, providing metastable states. Moreover, magnetic frustration is seen to induce self-buckling of the magnetostructure after dipole reorientation.
Even if reorientation of dipoles has been evidenced in tubes, it should also exists in other kind of structures. One can consider for example a square lattice plane or a parallelepiped. Changing the aspect ratio of these systems should trigger the transition. Minimization of the magnetic cohesive energy u is the main ingredient of the transition and is valid for all structures. Investigating dipole reorientation in other magnetostructures is left for future researches.
All these numerical results and experimental observations emphasize that the mechanical properties of the magnetostructures take root in the underlying dipole arrangement as well as the contact network between magnets. Moreover, magnetic particles can be found in many topics of physics such as magnetic colloids [3,23], micromachines [4,5], self-assembly in microfabrication methods [24][25][26], one understands that our findings may find future applications.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.