© C. Reiss et al., Published by EDP Sciences, 2018

1 Introduction

The physics of metastable states is a classical topic of statistical physics [1]. Many periodic systems can relax towards equilibrium when a large enough disturbance appears, such as supercooled liquids [2], dewetting cells [3], snow avalanches [4] and nucleating bubbles in liquids [5]. Of particular note are instabilities driven by a competition between weight and magnetic forces created by permanent magnets. They have already been studied in the case of 2D vertical [6] and horizontal [7] systems submitted to a vertical acceleration.

If such horizontal systems, made of magnets on a flat surface with their north pole facing up, are submitted to little or no vertical acceleration, a 2D crystal-like array appears. These magnets repel each other and self-organize in hexagonal structures (See Fig. 1).

This state is metastable, and no vertical acceleration is necessary to cause an instability: when the density of magnets is increased one magnet jumps up and attracts all the other magnets, which makes the structure collapse (See Fig. 2).

The goal of the study is to determine the cause of the collapse and the maximum density of magnets before it, which we will call the critical density. We also look into the dependence of this density on the physical parameters of the magnets. To answer these questions, we combine an experimental approach consisting in increasing the density of magnets up to the collapse with the predictions of a dipole model.

Fig. 1 2D hexagonal crystal-like magnet array. In this configuration, there are 58 ferrite magnets (radius 10 mm, height 3 mm) in a ∅25 cm circle.

Fig. 2 Chronophotography of the collapse of the array. One magnet is pushed into the frame through the slot (green circle). The critical density is reached, and a magnet jumps up between T = 0 ms and T = 100 ms (red circle) and attracts all the others afterwards. In this configuration, there are 33 ferrite magnets (R = 15 mm, h = 5 mm) in a ∅[30] cm circle.

2 Method

During this study, we conduct experiments on a single table made of laminated chipboard, a nonmagnetic material with far less friction than wood which enables the arrays to form.

We use only cylindrical magnets to preserve symmetry. 16 types of magnets are used: 13 made of ferrite and 3 of neodymium. For the ferrite (neodymium) magnets, the residual flux density is 0.405 T (1.305 T) and the mass density 4.840 g cm⁻³ (7.427 g cm⁻³). The radii R of our magnets range from 5 to 15 mm and their heights h vary from 3 to 15 mm.

We implement two ways of increasing the density of magnets: raising the number of magnets in a given surface (Method 1 and Fig. 1 for ferrite magnets) or reducing the surface with a given number of magnets inside (Method 2 and Fig. 3 for ferrite and neodymium magnets). For both methods, some magnets would crack on impact during the collapse and broken magnets were discarded.

We also develop a theoretical model based on small movements around the equilibrium. We suppose that having one magnet jump is sufficient to cause the collapse, as this has always been verified experimentally.

Fig. 3 Photography of the lasso used in Method 2 for ferrite magnets (R = 20 mm, h = 3 mm).

2.1 Experimental Method 1

In this method, we use 9 different circular wooden frames with diameters ranging from 10 to 30 cm to confine the magnets (See Fig. 1). The frames are circular to limit the effects of corners on the shape of the array and contain a slot that enables us to easily add magnets.

We push magnets through the slot and into the frame until the structure collapses (See Fig. 2), then count the number of magnets inside the frame.

2.2 Experimental Method 2

This method consists in reducing the support area with a given number of magnets inside. We put 7 magnets in a lasso and reduce its perimeter until the structure collapses (See Fig. 3). We chose to use 7 magnets so they would place themselves in a hexagonal array, as observed in Method 1. We measure the length of the lasso to deduce the distance between two magnets. Considering L_lasso the length of the lasso, d the distance between two magnets and R the radius of a magnet (See Fig. 4): $$d=\frac{L_{\text{Lasso}}-2\pi R}{6}\text{.}$$(1) We can then calculate the area $\mathcal{A}$ occupied by the central magnet, and the density of magnets $\mathcal{D}$: $$\mathcal{A}=\frac{\sqrt{3}}{2}{d}^2=\frac{1}{\mathcal{D}}\text{.}$$(2)

For each type of magnet, we carried out 5 or 10 measurements of the length of string surrounding the magnets at the moment of the collapse. We use the average length to determine the area occupied by the central magnet which we inverse to obtain the critical density using equation (2). We use the standard deviation to determine uncertainty.

Method 1 is closer to the problem that we wish to solve, i.e. determining the critical density of an array of magnets; but Method 2 is much quicker to implement experimentally.

Fig. 4 Geometry of the array and of the lasso used in Method 2; in blue the area occupied by the central magnet.

2.3 Theoretical model: two perfectly flat magnets

In this section, we consider magnets as magnetic dipoles [8] on a flat surface (See Fig. 5). In the past, the dynamics of close-range (d ≈ 3R) interacting magnets of varying shape has been predicted accurately using the dipole model [9,10,11]. Considering the magnetic moment ${\displaystyle \overrightarrow{\mu }}$ and magnetic moment norm μ of a single magnet, the vacuum permeability μ₀ and the distance between two magnets d, the magnetic energy of magnet (2) placed in the field ${\displaystyle \overleftarrow{B}}$ generated by magnet (1) is: $$E_{magnetic}=-{\displaystyle \overrightarrow{\mathit{\mu }}}\cdot {\displaystyle \overrightarrow{B}}=\frac{{\mu }_0{\mu }^2}{4\pi d^3}\text{.}$$(3)

Equation (3) shows that to minimize their interaction energy, the magnets must be as far away as possible and repel each other. Furthermore, as the energy only depends on one horizontal parameter there can be no vertical force, no jump and no collapse. In this case, the maximum density will be for a compact hexagonal packing of magnets inside our surface, which leads to a mathematical optimization problem [12]. But this is insufficient to explain the phenomenon we observe. This is why we introduce irregularities in the next section.

Fig. 5 Situation used to calculate the energy of interaction between two flat magnets.

2.4 Theoretical model: two magnets with irregularities

We will now suppose that there are irregularities: the magnetic dipoles that represent our magnets can be at an altitude z ≪ d with a moment tilted with and angle θ ≪ 1 from the vertical (See Fig. 6). They can be askew in this way because of irregularities on the surface or inside the magnets. In the small irregularities region, the dipolar magnetic field created by (1) can be approximated with a second order Taylor expansion [8]: $${\displaystyle \overrightarrow{B}}\approx \frac{{\mu }_0\mu }{4\pi d^3}\left(1-\frac{3}{2}\frac{z^2}{d^2}\right)\left(3\frac{z}{d}{\displaystyle \overrightarrow{u_x}}+\left(-1+3\frac{z^2}{d^2}\right){\displaystyle \overrightarrow{u_z}}\right)\text{.}$$(4) We then compute a second-order Taylor expansion for the magnetic moment of magnet (2): $${\displaystyle \overrightarrow{\mu }}=\mu \left(\left(1-\frac{1}{2}{\theta }^2\right){\displaystyle \overrightarrow{u_z}}+\theta \cos \varphi {\displaystyle \overrightarrow{u_x}}+\theta \sin \varphi {\displaystyle \overrightarrow{u_y}}\right)\text{.}$$(5)This enables us to calculate the magnetic energy of our system: $$E_{magnetic}=-{\displaystyle \overrightarrow{\mathit{\mu }}}\cdot {\displaystyle \overrightarrow{B}}=\frac{{\mu }_0{\mu }^2}{4\pi }\frac{1}{d^3}\left(1-\frac{9}{2}{\left(\frac{z}{d}\right)}^2-\frac{1}{2}{\theta }^2-3\theta \frac{z}{d}\cos \varphi \right)\text{.}$$(6)For z = 0 and θ = 0, we obtain the same formula as equation (3).

Fig. 6 Positions of both magnets and physical meaning of the quantities used in the model: θ, φ and z.

2.5 Theoretical model: one tilted magnet surrounded by six flat magnets

Now that we have the interaction energy for two magnets, we place ourselves in a situation closer to the one just before the jump. We will suppose that a tilted magnet is surrounded by 6 perfectly flat magnets in a hexagonal array, all placed at the same distance d. This cancels the cosφ term and adds up the other terms in equation (6). We can therefore derive from the energy a vertical force and a torque exerted by the surrounding magnets on the middle of the central magnet (Point A on Fig. 6): $$F_{\text{total},A}=-\frac{\partial }{\partial z}(E_{\text{total}})=6\frac{{\mu }_0{\mu }^2}{4\pi }\frac{1}{d^3}\frac{9z}{d^2}\text{,}$$(7) $$T_{\text{total,A}}\mathrm{}=-\frac{\partial }{\partial \theta }\mathrm{}(E_{\text{total}})=6\frac{{\mu }_0{\mu }^2}{4\pi }\frac{1}{d^3}\theta \text{.}$$(8) This enables us to compute the total torque on Point B, at the basis of the central magnet: $$T_{\text{total},B}=T_{\text{total},A}+R\cdot F_{\text{total},A}=6\frac{{\mu }_0{\mu }^2}{4\pi }\frac{1}{d^3}\left(\theta +R\frac{9z}{d^2}\right)\text{.}$$(9)The equilibrium is a result of a competition between this magnetic torque and weight. Considering the mass m and the radius R of the magnet, the central magnet will jump when: $$T_{\text{total},B}=6\frac{{\mu }_0{\mu }^2}{4\pi }\frac{1}{d^3}\left(\theta +9\frac{Rz}{d^2}\right)>Rmg\text{.}$$(10)A jump cannot be caused by the force calculated in equation (7): since θ > 0, the total torque in equation (10) always becomes greater than weight before force.

If a magnet's displacement is sufficient for it to start rotating, z and θ will increase, and so will the magnetic force and the magnet's angular velocity. This makes the equilibrium metastable: if a magnet undergoes a tiny movement it will stay in place, but if it undergoes a big enough movement the structure will collapse.

We separate two distinct cases whether θ or $9\frac{Rz}{d^2}$ is dominant. We assume that θ and $\frac{z}{R}$ will be of the same order of magnitude: $\theta \approx \frac{z}{R}$. This enables us to determine whether the torque or the force applied in A is mainly responsible for the jump.

If d_crit is the critical distance between two magnets, force causes the jump if: $$9\frac{Rz}{d_{crit}^2}>\theta \approx \frac{z}{R}\text{,}$$(11) Which amounts to: $$d_{crit}<3R\text{.}$$(12)Force causes the jump at close range, and torque at long range.

We can use equation (10) in order to determine an approximate critical distance at the moment of the jump if force is dominant: $$d_{\text{force}}\approx {\left(\frac{27}{2\pi }\frac{{\mu }_0{\mu }^2}{mg}z\right)}^{\frac{1}{5}}$$(13) Considering magnets of height h, mass density ρ, and residual flux density M, and knowing that μ = MπR²h/μ₀, we can determine a relation between the critical density $\mathcal{D}$ and the other parameters thanks to equation (2): $${\mathcal{D}}_{\text{force}}=\alpha z^{-\frac{2}{5}},\alpha =\frac{2}{\sqrt{3}}{\left(\frac{27}{2\pi }\frac{M^2\pi R^{2}h}{{\mu }_0\rho g}\right)}^{-\frac{2}{5}}\text{.}$$(14)

In the same way, if torque causes the jump: $${\mathcal{D}}_{\text{torque}}=\beta {\theta }^{-\frac{2}{3}},\beta =\frac{2}{\sqrt{3}}{\left(\frac{3}{2\pi }\frac{M^{2}R\pi h}{{\mu }_0\rho g}\right)}^{-\frac{2}{3}}\text{,}$$(15) where α and β only depend on the properties of the magnets.

3 Results

3.1 The critical density is independent from the area

For Method 1, we plot the number of magnets at the collapse according to the surface of the frame and compute an affine regression. The relation is linear, with a systematically positive intercept that is very small compared with the typical number of magnets in our system. We use the slope as critical density (See Fig. 7) and the standard deviation to determine the uncertainty.

Fig. 7 Experimental results for Method 1 using three ferrite magnets (Magnet 1: h = 3 mm, R = 5 mm; Magnet 2: h = 3 mm, R = 10 mm, Magnet 3: h = 5 mm, R = 15 mm). The slopes of the fits, which we will use as the critical density, are respectively 0.148 cm−2, 0.118 cm−2 and 0.059 cm−2.

3.2 Force or torque causes the jump depending on the critical density

For Method 2, we plot the measured critical density of magnets according to α and β using the data sheets of our magnets (See Fig. 8a and b).

We optimize the threshold value d_t of the ratio $\frac{d_{\text{crit}}}{R}$ to have the best fit on the two figures. We choose d_t to maximize the average of the Pearson correlation coefficients between α (β) and experimental points for force-driven (torque-driven) magnets with $\frac{d_{\text{dcri}}}{R}<d_t$ $\left(\frac{d_{\text{crit}}}{R}\mathrm{}>d_t\right)$. This yields d_t = 3.6 which is close to the prediction of 3 given by equation (12). We then separate the magnets into two groups: force and torque driven. In Figure 8d, the torque-driven jumps are far from the fit on force. In Figure 8b, force-driven magnets are closer to the torque-driven magnets. This means that torque still has a role to play at close range, while force does not at long range.

3.3 The critical density can be computed thanks to the parameters of the magnets and of the irregularities

Using the results of the previous section, we compute a theoretical total critical density by solving in equation (10) by injecting the values for z and θ given by the fits on Figures 8a and b and using a numerical solver. This gives us a theoretical density that takes into account force and torque, for which the results are displayed on Figure 8c. The Pearson coefficient of correlation for this fit is 0.92.

Fig. 8 Experimental results and comparison with theory (a) Measured density as a function of α defined in (14); triangles are for force driven jumps $\left(\frac{d_{\text{crit}}}{R}<d_t=3.6\right)$, circles for torque driven jumps $\left(\frac{d_{\text{crit}}}{R}>d_t\right)$, and the color scale gives the distance between two magnets over their radius at the time of the jump (dcrit/R); the fit yields z = 15 mm. (b) Measured density as a function of β defined in (15); the conventions used are the same as in (a); the fit yields θ = 9∘. (c) Measured density as a function of theoretical density calculated thanks to equation (10); the values used for z and θ are computed from the fits in (a) and (b); the conventions used are the same as in (a).; (d) Comparison between results from method 1 and method 2; data gathered from method 1 are more scattered.

3.4 The lasso method yields slightly higher critical densities than the circles

We then plot the experimental results for Method 1 and 2 on the same graph (See Fig. 8d). We can see that the measured density is smaller for Method 1 than for Method 2. Furthermore, it appears that the results for Method 1 are more dispersed than the results for Method 2.

4 Discussion

Torque has a greater role to play at close range than force at long range because the distance between two magnets cannot drop under 2R, as they cannot overlap. In practice, the smallest measured $\frac{d_{\text{crit}}}{R}$ was 2.5, which is close to $\frac{d_{\text{crit}}}{R}=3.6$, so torque still plays a role. We can also see in Figure 8a that the greater $\frac{d_{\text{crit}}}{R}$, the further the points are from the fit.

The experimental separation we find for torque or force causing the jump is for $\frac{d_{\text{crit}}}{R}=3.6$ and our theoretical separation is for $\frac{d_{\text{crit}}}{R}=3$. This difference is quite small, as the approximation $\theta \approx \frac{z}{R}$ is an equality only for specific topologies. That being said, it is normal that our experimental separation is greater than our theoretical separation: due to irregularities in the crystal-like structure, for the same average distance between magnets, some will be closer in practice than in theory. This will lead to force causing the jump.

The density for Method 1 being lower than for Method 2 can be explained in part by imperfections in the crystalline matrix that were not present on method 2 (where we forced the structure to be hexagonal). It can also be explained by the fact that for method 1, only one out of all magnets needs to be askew to cause a collapse.

It is surprising to find values of irregularities as important as those that we have here: z = 1.5 mm and θ = 9^∘. There are three possible causes for these irregularities: a bumpy surface, geometrically imperfect magnets and an inhomogeneous magnetization. As we were able to have z and θ independent of the magnets, we believe that a bumpy surface is the most probable cause of irregularities. The surface we used is visibly wavy, so these values are physically possible. But the different causes could also add up.

It is possible to understand why magnets form a cylinder when one jumps by comparing the interaction energy between two magnets in an array and in a collapsed cylinder (See Fig. 2). This energy is typically $E=-{\displaystyle \overleftarrow{\mu }}\cdot {\displaystyle \overleftarrow{B}}\approx \pm \frac{{\mu }_0{\mu }^2}{4\pi d^3}$ [8]. The energy of a cylinder is negative while the energy of an array is positive. Furthermore, the distance between two magnets in an array is of the order of magnitude 10 times that in a collapsed cylinder, so the norm of the total energy of a cylinder is typically 10³ that of an array. Magnets will strive to go into a collapsed state and join the cylinder once it has started to form, so one magnet jumping is sufficient for a collapse.

We obtain a good agreement between theory and experiment by using the dipole model at close range (d ≈ 3R). We were not the first to correctly describe magnet dynamics in this manner [9,10,11].

5 Dead ends

5.1 Dependence on the surface irregularities

We succeeded in understanding that the key parameters of the problem, the ones that cause the jump, are the surface and magnet irregularities. They appear in our equations as the z and θ parameters.

We tried to vary surface irregularities in a controlled way by 3D printing bumpy surfaces: we were able to see that the bigger the irregularities, the smaller the critical density. But we did not obtain any quantitative results.

5.2 Influence of friction

Another parameter we tried to investigate is friction between magnets and the surface. Friction prevents magnets from forming a hexagonal array by blocking their movements. Consequently, for the same global density some magnets will have closer neighbors in the presence of friction than in the absence of friction, and the system will have a bigger chance of collapsing. So in the presence of friction, the critical density should diminish.

This is what we observed by conducting experiments on tables with different friction coefficients. But we were not able to determine the average change in the pattern of the array due to friction, which is the cause of the change in critical density.

6 Conclusion

We proposed in this paper two experimental setups to determine the critical density of magnets in a 2D hexagonal array. We also developed a theoretical model based on the dipole model to explain the physics of this jump, and in what way this state is metastable. We found this jump to be caused by irregularities, either of the surface or of the magnets. Furthermore, depending on the ratio of the distance between two magnets and the radius of a magnet, the jump will be caused by a force or a torque applied on the central magnet by the others. We were able to confront theory and experiment and confirm the relevance of our formulas, in particular the split between jumps due to force and torque and the impact of the different parameters of the magnets, namely their heights, radii, magnetizations and mass densities.

Further research on the subject could involve quantifying the impact of irregularities and pattern changes in the array due to friction.

Acknowledgments

This study was carried out for the participation to the 2017 French Physicists' Tournament and International Physicists' Tournament. We would like to thank the organization of both tournaments. We would also like to thank Simon Merminod for his advice, Jacques Lebret, Laurent Bergeon, Paul Lilin, Declan McCavana, Alistair Rowe and Christophe Clanet for their help and of course our team leader Guilhem Gallot.

References

All Figures

	Fig. 1 2D hexagonal crystal-like magnet array. In this configuration, there are 58 ferrite magnets (radius 10 mm, height 3 mm) in a ∅25 cm circle.
In the text

	Fig. 2 Chronophotography of the collapse of the array. One magnet is pushed into the frame through the slot (green circle). The critical density is reached, and a magnet jumps up between T = 0 ms and T = 100 ms (red circle) and attracts all the others afterwards. In this configuration, there are 33 ferrite magnets (R = 15 mm, h = 5 mm) in a ∅[30] cm circle.
In the text

	Fig. 3 Photography of the lasso used in Method 2 for ferrite magnets (R = 20 mm, h = 3 mm).
In the text

	Fig. 4 Geometry of the array and of the lasso used in Method 2; in blue the area occupied by the central magnet.
In the text

	Fig. 5 Situation used to calculate the energy of interaction between two flat magnets.
In the text

	Fig. 6 Positions of both magnets and physical meaning of the quantities used in the model: θ, φ and z.
In the text

	Fig. 7 Experimental results for Method 1 using three ferrite magnets (Magnet 1: h = 3 mm, R = 5 mm; Magnet 2: h = 3 mm, R = 10 mm, Magnet 3: h = 5 mm, R = 15 mm). The slopes of the fits, which we will use as the critical density, are respectively 0.148 cm−2, 0.118 cm−2 and 0.059 cm−2.
In the text

Fig. 8 Experimental results and comparison with theory (a) Measured density as a function of α defined in (14); triangles are for force driven jumps $\left(\frac{d_{\text{crit}}}{R}<d_t=3.6\right)$, circles for torque driven jumps $\left(\frac{d_{\text{crit}}}{R}>d_t\right)$, and the color scale gives the distance between two magnets over their radius at the time of the jump (dcrit/R); the fit yields z = 15 mm. (b) Measured density as a function of β defined in (15); the conventions used are the same as in (a); the fit yields θ = 9∘. (c) Measured density as a function of theoretical density calculated thanks to equation (10); the values used for z and θ are computed from the fits in (a) and (b); the conventions used are the same as in (a).; (d) Comparison between results from method 1 and method 2; data gathered from method 1 are more scattered.

In the text