Non-monotonic Dependence of Magnetic Ordering on Geometry Regulated by Intra-Unit Interactions in Macroscopic Square Ice

1. Introduction

Geometric frustration refers to constraints that prevent simultaneous minimization of all interactions, leading to a highly degenerate ground state and the emergence of novel phenomena. The square ice model is a representative classical frustrated system and reveals the mathematical essence of frustration in spin ice—the impossibility of graph-coloring in certain cases. In the square-ice model, at each vertex the four nearest-neighbor spins satisfy the “two-in, two-out” ice rule; the ground-state manifold is highly degenerate, providing a suitable theoretical framework to study emergent excitations and phase behavior in frustrated systems. Artificial spin ice systems based on nanoscale magnetic islands have experimentally realized the square ice model’s ground-state manifold and enabled observation of magnetic monopole excitations, Dirac strings, and Coulomb phases. However, due to ultrafast spin dynamics and limitations of nanoscale characterization, direct tracking of spin flips and vertex rearrangements remains challenging, leaving gaps in understanding the mechanisms and evolution of square ice. Macroscopic frustrated systems with limited degrees of freedom have traditionally struggled to reproduce the complete ground-state manifold and its evolution. In response, we designed a macroscopic square ice platform comprising hundreds to thousands of NdFeB magnets confined in vertical cylindrical wells on a PMMA base, each unit containing two magnets of different heights. The platform offers three advantages: (1) room-temperature optical observation enables direct recording; (2) unit gap l and vertex gap d can be independently tuned to study intra-unit vs vertex interactions; (3) vertex types and magnetic-domain evolution can be tracked throughout, providing a visual method to study evolution in frustrated systems. This study focuses on how intra-unit interactions regulate the non-monotonic evolution in the square ice model. By fixing d = 1.1 mm and varying l from 9.7 mm down to 0.7 mm, we observe the evolution of vertex-type proportions and magnetic-domain sizes. We find that as l decreases (increasing intra-unit interaction), the system experiences a non-monotonic evolution. At l ≈ 9d, the system is excited ordered due to abundant excited units and mostly ice-rule vertices; at l ≈ 1.4d and 2.1d, long-range antiferromagnetic order emerges; at l ≈ d, the Type-I vertex fraction plunges and disordered states appear; for l = 0.64d, the system re-enters an antiferromagnetic order. Magnetic structure factor analysis shows qualitative characteristic features that motivate a possible connection to constrained-string physics. This work provides the direct optical imaging of assembled configurations for frustrated system ordering at the macroscopic scale and offers a tunable platform for exploring topological excitations and evolution.

2. Experimental Methods

2.1. Sample Design and Geometry

Macroscopic square ice consists of hundreds to thousands of NdFeB magnets arranged in a square lattice on a PMMA base. Each magnet is confined in a cylindrical well slightly larger than the magnet diameter, permitting only two vertical states: up or down. The well depth is 10 mm to ensure adequate vertical space. Each unit contains two magnets with different heights (long: 6 mm, short: 5 mm), representing a pseudo-spin. There are two possible ground states: (a) left magnet up, right magnet down (pseudo-spin left) and (b) left long magnet down, right short magnet up (pseudo-spin right). All magnets are oriented with the North pole up; the long and short magnets are staggered in the lattice. The minimum horizontal spacing within a unit defines the unit gap l (Figure 1(b)); l values used are 9.7, 2.3, 1.5, 1.1, and 0.7 mm to tune intra-unit magnetic interactions.

The four neighboring units meet at a vertex; the minimum distance between the two adjacent units defines the vertex gap d, fixed at 1.1 mm, enabling analysis of intra-unit vs vertex interactions. The pseudospin direction is determined based on the absolute height of the magnets and the relative vertical displacement between the centers of two adjacent magnets. The difference between the center heights of the first and second magnets is denoted as Z. Taking the trap bottom as zero, the height of the magnet top relative to the trap bottom is called the absolute height. When a magnet exhibits the up state, the threshold (minimum) of the absolute height is 9.5 mm; when it exhibits the down state, the threshold (maximum) is 6.5 mm. In principle, a magnet could have an intermediate state, with its absolute height between 4.5 mm and 5.5 mm. Two up magnets form an up-up unit, with Z ranging from -0.5 mm to 0.5 mm. Its pseudospin is represented by two outward arrows. Two down magnets form a down-down unit, with Z ranging from -0.5 mm to 0.5 mm, and its pseudospin is represented by two inward arrows. The spin magnitude is 1. Magnets in intermediate states are extremely rare in the square lattice; in such cases, and the intermediate-state magnet is represented by a fractionalized spin, marked as two small arrows, with spin magnitude 1/2. At this point, all possible states of a unit have been assigned pseudospins. According to the real parameters, positions and magnetic flux density, unit energy (interaction energy of two magnets in a unit) can be obtained by Comsol for different Z (Figure 1(d)).

2.2. Vertex Types and Classification Criteria

At a vertex, four pseudo-spins (one per unit) only with unit ground states create 16 configurations. Following square ice and artificial spin ice conventions, these configurations are grouped into four classes: Type I (2-in-2-out ice rule; 2 configurations; lowest energy), Type II (2-in-2-out ice rule; 4 configurations; next-lowest energy), Type III (3-in-1-out or 1-in-3-out; 8 configurations; monopole excitations; higher energy), Type IV (4-in or 4-out; 2 configurations; highest energy; rare). I and II comply with ice rule; III and IV violate it. The four vertex types are illustrated with dark magnets down and red magnets up; arrows denote pseudo-spin directions. In a fully random array, the proportions are 12.5%, 25%, 50%, and 12.5% for Type I–IV, respectively.

Figure 2. Vertex types in macroscopic square ice. (a) Schematic of Type-I vertices. Vertex gap d is fixed to be 1.1 mm. (b) Schematic of Type-II vertices, where pseudospins share the same direction for each row and each column. (c) Schematic of Type-III vertices, which are elementary excitations in macroscopic square ice. (d) Schematic of Type-IV vertices, which rarely appear in the system.

Figure 2. Vertex types in macroscopic square ice. (a) Schematic of Type-I vertices. Vertex gap d is fixed to be 1.1 mm. (b) Schematic of Type-II vertices, where pseudospins share the same direction for each row and each column. (c) Schematic of Type-III vertices, which are elementary excitations in macroscopic square ice. (d) Schematic of Type-IV vertices, which rarely appear in the system.

When vertex types include up-up and down-down units, the number of vertex configurations can be expanded to 256, among which there are 32 configurations for Type-I vertices, 64 for Type-II vertices, 128 for Type-III vertices, and 32 for Type-IV vertices. The vertex configuration is mainly determined by the states of the inner four magnets of the vertex, while the outer four magnets of the vertex can be in arbitrary states, thus multiplying the original number of configurations by 16 and generating a large number of new configurations. The correlation U between the two magnets in up-up and down-down units is 1, which does not affect the nearest-neighbor correlation I of the vertex; the nearest-neighbor correlations of Type-I and Type-IV vertices remain -1 and 1, respectively, and those of Type-II and Type-III vertices remain 0. The correlation function L is also unaffected; the L correlations of Type-I, Type-II, Type-III, and Type-IV vertices remain 1, -1, 0, and 1, respectively. Besides affecting the unit correlation function U, up-up and down-down units also change the next-nearest-neighbor correlation L. At this point, the next-nearest-neighbor correlation L cannot be calculated by simple rules and must be analyzed vertex by vertex. In the calculation of the magnetic structure factor, to describe up-up and down-down units and distinguish them from ordinary up-down units, this study uses 11 to represent an up-up unit, 00 to represent a down-down unit, and 01 and 10 to represent different up-down units. The magnet states in each row of horizontal units are represented by binary symbols, with 1 denoting a magnet up and 0 denoting a magnet down. The states of each row of horizontal units are converted into row vectors, and the overall state of the horizontal units is transformed into a matrix. The magnet states in each column of vertical units are also represented by the same binary symbols; then the states of each column of vertical units form a row vector, and the states of all vertical units are transformed into another matrix. These two matrices are then fed into a computer program, which can automatically generate the magnetic structure factor. The magnetic structure factor is computed according to the real positions of magnets.

2.3. Sample Preparation and Assembly

To control variables, all samples use a consistent assembly. In the PMMA base, horizontal unit rows are filled with magnets in a sequence as shown in Figure 3(a); neighboring units have opposite magnet order, and adjacent rows have opposite horizontal order. Vertical units follow a similar scheme. This assembly tends to promote Type-I vertices and enhances near-neighbor interactions, while stronger intra-unit interactions may induce Type-II and Type-III vertices spontaneously. Vertex energy is shown in Figure 3(c).

To quantitatively characterize the degree to which the system satisfies the ice rule, this work defines two derived physical quantities: (1) the ice-rule satisfaction rate R_ice = (N_I + N_II)/N_total, which equals the sum of the proportions of Type-I and Type-II vertices, reflecting the overall proportion of vertices satisfying the ice rule; (2) the monopole excitation density ρ_mono = (N_III + 2N_IV)/N_total, reflecting the proportion of vertices that break the ice rule, where N_I, N_II, N_III, and N_IV denote the numbers of the corresponding vertex types, and N_total is the total number of vertices. Since Type-IV vertices are rare, the monopole excitation density is approximately equal to the proportion of Type-III vertices.

2.4. Experimental Parameter Characterization

To support the consistency of sample preparation, this study systematically measured and statistically analyzed the key geometric and magnetic parameters of the traps and magnets used in the experiments. As shown in Figure 4(a), the key geometric parameter of the trap is the trap diameter D′. The distribution of the trap diameter D′ is shown in Figure 4(b). One hundred traps were randomly selected for measurement; their diameters range from 3.88 to 4.88 mm, with a mean value of 4.37 mm and a standard deviation of 0.15 mm, approximately following a Gaussian distribution. Optical images and schematic diagrams of the two types of magnets are shown in Figure 4(c). Both types of magnets have the same diameter D, whose distribution is shown in Figure 4(d), with a mean value of 3.97 mm and a standard deviation of 0.03 mm. Their heights are denoted as h₁ and h₂. The magnetic flux densities and masses are denoted as B₁, B₂, m₁, and m₂, respectively.

One hundred short magnets and one hundred long magnets were selected, and their magnetic flux densities, masses, and heights were measured individually. The measurement results are shown in Figure 5, Figure 6 and in Table 1 and Table 2. The mean values and standard deviations of each parameter are given in the second and third rows of Table 1 and Table 2, respectively. Gaussian fitting (red curves) shows that all parameters approximately follow Gaussian distributions, indicating good consistency of the magnets. The mean values of the magnetic flux density B₁, mass m₁, and height h₁ of the short magnets are 253.7 mT, 0.47 g, and 4.96 mm, respectively; the corresponding values for the long magnets are 278.2 mT, 0.55 g, and 5.77 mm. The actual heights of the two types of magnets are in good agreement with the design values (5 mm and 6 mm).

For ease of description, this study classifies the system states into three categories based on the characteristic features of the magnetic structure factor: states in which bow-tie diffuse characteristic features occur are called excited ordered states (EO); states in which the magnetic structure factor shows only antiferromagnetic Bragg peaks are called antiferromagnetic ordered states (AFM); states in which antiferromagnetic and ferromagnetic Bragg peaks coexist are called disordered states (DS). It should be noted that the term "state" here is used only to describe the experimentally observed orderliness, and does not correspond to a strict thermodynamic phase.

3. Results

As l decreases, intra-unit interactions become stronger. Figure 7 shows the spatial distribution of vertex types as a function of l. At l = 9.7 mm, nearly all vertices are Type-I, forming a large connected Type-I domain. Because most units are in excited configurations, the system is identified as an excited ordered state. When l is reduced to 2.3 mm and 1.5 mm, Type-I domains remain dominant, but Type-II and Type-III domains emerge in some regions, yielding an antiferromagnetic-like ordered state. At l = 1.1 mm, the Type-I vertex fraction drops to 36.36%, while Type-II and Type-III rise significantly, placing the system in a disordered state with a high monopole density. With a further reduction of l to 0.7 mm, Type-I domains re-grow and Type-II / Type-III domains shrink, partially restoring antiferromagnetic order. This “order–disorder–order” evolution is referred to as non-monotonic behavior.

It is found that the ratio between the unit gap l and the vertex gap d determines the evolution of system order. When l gradually decreases from far greater than d to less than d, the system shows a non-monotonic geometric dependence: when l ≈ 9d, Type-I vertices dominate, although a large fraction of units remain in excited states, resulting in an excited ordered state; l ≈ 1.4d and 2.1d, Type-I vertices dominate and most units are in their ground states, exhibiting long-range antiferromagnetic order; when l ≈ d, the proportion of Type-I vertices drops sharply to 36.36%, while the proportions of Type-II and Type-III (monopole excitation) vertices rise significantly, and the system enters a highly disordered state; when l = 0.64d, the proportion of Type-I vertices rebounds to 78.79%, and antiferromagnetic order is restored. The above evolution originates from interaction competition.

3.1. Vertex Proportions and Evolution of Order

Figure 8(a) shows the states for five sample sets: Type-I, Type-II, and Type-III vertex proportions vs l. When l decreases from 9.7 mm to 0.7 mm, a pronounced non-monotonic evolution is observed. At l ≈ 9d, Type-I vertices dominate with almost all vertices obeying ice rules, but many units are excited (inelastic to gravity), producing an excited ordered state. When l ≈ 1.5 mm (≈ 1.4d), Type-I remains dominant and most units are in ground states, yielding long-range antiferromagnetic order. At l ≈ d (1.1 mm), the Type-I vertex fraction decreases sharply to 36.36%, while II and III rise to 24.52% and 39.12%, respectively, indicating a disordered state with a high monopole density. When l decreases to 0.7 mm, the Type-I rises to 78.79% and antiferromagnetic order partly recovers. The evolution is attributed to interaction competition.

The states and correlation data (Figure 8) show that the nearest-neighbor correlation of Type-I vertices and next-nearest-neighbor correlation L vary with l in a manner that tracks vertex-type proportions, with a notable minimum near l = d, indicating maximal competition. Figure 8(a) shows the states of macroscopic square ice, presenting the variations of the proportions of Type-I, Type-II, and Type-III vertices with l for five groups of samples. When l decreases from 9.7 mm to 0.7 mm, the system shows a significant non-monotonic evolution. When l = 9.7 mm (about 9d), the proportion of Type-I vertices is 97.22%, while the proportions of Type-II and Type-III vertices are both close to zero. At this point, almost all vertices satisfy the ice rule; however, it is noteworthy that a large number of the two magnets within each unit are simultaneously in the up or down excited states.

This can be attributed to the intra-unit magnetic interaction and vertex interaction. For the first row of horizontal units, all magnets are down because the intra-unit magnetic interaction is too weak to overcome the gravitational weight of the magnets. The next row of vertical units will be in the ground state because the upper half of the magnets is stabilized by vertex interactions to be up while the lower half is governed by gravity, resulting in a downward alignment. For the remaining horizontal units, the vast majority of macroscopic magnets are supported by vertex interactions and thus align upward, with only a small number at the boundaries deviating from this state. For the remaining vertical units, the magnets are simultaneously influenced by gravity and the downward repulsive vertex interaction exerted by the horizontal-unit magnets; consequently, except for a few boundary magnets, all are aligned downward. The system is highly ordered with many up-up units and down-down units, so this state is called an “excited ordered state (EO)”.

As l decreases to 1.5 mm (about 1.4d), the proportion of Type-I vertices drops to 84.35%, while the proportions of Type-II and Type-III vertices rise slowly. At this point, the intra-unit interaction has become strong enough to stabilize most units in their ground states, and the system enters ground states at both the vertex and unit levels, thereby establishing long-range antiferromagnetic order. When l is further reduced to 1.1 mm (i.e., l = d), the vertex proportions change dramatically: The Type-I vertex fraction decreases sharply to 36.36%, while the proportions of Type-II and Type-III vertices jump to 24.52% and 39.12%, respectively. At this point, the proportions of both Type-I and Type-II vertices are below 50%, the ice-rule satisfaction rate (I + II) is only 55.62%, and the magnetic monopole density is as high as 39.12%, indicating that the system has entered a highly disordered state with a high monopole density. This state arises because the intra-unit interaction and the inter-vertex magnetic interaction become comparable in strength; the friction between the first magnet of a unit and the trap wall is large enough to overcome gravity and keep it in the up state, while the second magnet is down, causing the pseudospin of the unit to reverse. This process is called unit flipping. One unit flip can turn a Type-I vertex into a Type-III vertex, while two unit flips turn it into a Type-II vertex; therefore, the proportion of Type-III vertices is usually higher than that of Type-II vertices, which is consistent with the experimental results. When l decreases from 1.5 mm to 1.1 mm, macroscopic square ice evolves from an antiferromagnetic state to a disordered state; when l decreases from 1.1 mm to 0.7 mm, it returns from the disordered state to the antiferromagnetic state [71].

When l decreases to 0.7 mm (0.64d), the proportion of Type-I vertices rebounds to 78.79%. The proportion of Type-III vertices drops to 15.15%, and the ice-rule satisfaction rate recovers to 73.37%, with antiferromagnetic order partially restored. From an energy perspective, Type-I vertices have the lowest interaction energy, so when only one type of interaction is present, the system preferentially forms Type-I vertices. However, the system involves interaction competition. Therefore, as l decreases, the system undergoes the non-monotonic evolution. For l takes the values of 9.7, 2.3, 1.5, 1.1 and 0.7 mm, the numbers of samples are all three. Repeated experiment shows that the standard deviations of vertex proportions are all below 13%. With the number of vertices fixed to be 121, non-monotonic evolution occurs in representative samples. For the proportion of Type-I vertices, the largest standard deviation 12.83% occurs at l = 1.1 mm and the least standard deviation 3.91% occurs at l = 2.3 mm (Table 3). For the proportion of Type-II vertices, the largest standard deviation 9.29% also occurs at l = 1.1 mm, but the least standard deviation 0.95% occurs at l = 0.7 mm (Table 4). For the proportion of Type-III vertices, the largest standard deviation 9.02% occurs at l = 0.7 mm and the least 1.52% occurs at l = 9.7 mm (Table 5).

3.2. Supporting Evidence from Correlation Functions, Unit Proportions, and Domain Sizes

The changes in correlation functions further support the above evolution trend. Figure 8(b) shows the variation of the correlation strengths of nearest-neighbor I and next-nearest-neighbor L with l. The inset shows the positional relationship among nearest-neighbor I, next-nearest-neighbor L, and the gray magnet: all three are located at the same vertex; the two I magnets are respectively at the upper-left and lower-left positions of the gray magnet, and the L magnet is collinear with the unit containing the gray magnet. When l decreases from 2.3 mm to 1.1 mm and then to 0.7 mm, the nearest-neighbor I correlation changes from -0.98 to -0.47 and then recovers to -0.82; the next-nearest-neighbor L correlation changes from 0.92 to 0.26 and then rebounds to 0.77. The extrema of the correlation functions occur at l = d, in full agreement with the changes in vertex proportions. Overall, as l decreases, the correlation functions exhibit an evolution trend similar to that of the vertex proportions, reflecting changes in the system order, with weak correlations occurring exactly at l = d.

When the two magnets have equal lengths, the magnets with flush tops have horizontal interactions, which tend to form excited states. To minimize the occurrence of excited states, this work adopts magnets of different lengths. When the magnet tops are flush, the interactions along the vertical direction between magnets can reduce excited states. The up and down states of the magnets correspond to the north and south poles of ferromagnetic nanoislands, respectively. To avoid a significant difference in magnetic flux density at the two ends of a pseudospin, the length difference between the long and short magnets is only 1 mm. The interaction energy of a unit depends only on the relative vertical displacement between the centers of the two magnets; the vertical displacement, mass, and magnetic flux density between the long and short magnets do not break the energy degeneracy of the two pseudospin directions. Odd numbers represent short magnets, marked in blue, and even numbers represent long magnets, marked in red. In horizontal units of odd-numbered rows, the magnet order follows a 1243 cycle; in even-numbered rows, the corresponding two units show a 2134 cycle. For vertical units, from the first column to the last column, the magnet order of every two units in odd-numbered columns follows a 2134 cycle, while that in even-numbered columns follows a 1243 cycle. For horizontal units in adjacent rows, the magnet orders are opposite. For vertical units in adjacent columns, the magnet orders are also opposite. At the boundaries, the magnet orders of the first row of horizontal units and the first column of vertical units are given, while those of the last row of horizontal units and the last column of vertical units should be determined according to the parity of the row and column numbers.

As shown in Figure 8(c), when l is 9.7 mm, the proportion of up-down units is only 15.11%, while the proportions of up-up and down-down units are 39.97% and 44.92%, respectively. When l decreases to 2.3 mm, the proportion of up-down units reaches 93.91%, and the proportions of up-up and down-down units drop to 2.24% and 3.85%, respectively. As the intra-unit interaction strengthens, the proportions of up-up and down-down units decrease significantly.

The magnetic domain size is defined as the maximum number of vertices within a single domain. The changes in the sizes of Type-I, Type-II, and Type-III domains at l = 9.7 mm, 2.3 mm, 1.5 mm, 1.1 mm, and 0.7 mm are detailed in Table 6. The maximum size of Type-I domains gradually decreases from 118 vertices at l = 9.7 mm to 30 vertices at l = 1.1 mm, and then rebounds to 68 vertices at l = 0.7 mm. The size of Type-III domains increases monotonically from 0 to 17 vertices, and does not decrease significantly when l = 0.64d, indicating that monopole excitations are more likely to nucleate and aggregate under strong interaction conditions, but remain smaller than Type-I domains (Table 7).

An energy competition model can qualitatively explain this phenomenon. Vertex proportions are related not only to the interaction energy but also to the assembly method. Because magnets interact with each other during the process, which is similar to the mutual collisions of microscopic molecules, the effect of the assembly method can be attributed to an energy barrier. The energy barrier is related to the number of unit flips. Under a given assembly method, the fewer the unit flips required for a vertex, the lower its energy barrier. Type-II vertices require more unit flips than Type-III vertices, so the energy barrier of Type-II vertices is higher than that of Type-III vertices, and the proportion of Type-III vertices is generally higher than that of Type-II vertices. As shown in Table 8 and Table 9, Type-II domains are generally smaller than Type-III domains. Since Type-III vertices can be regarded as elementary excitations, macroscopic square ice provides new macroscopic experimental evidence for studying elementary excitations in frustrated systems.

3.3. Magnetic Structure Factor

In addition to vertex proportions, correlation functions, and domain sizes, the magnetic structure factor can also effectively characterize the magnetic order of the system. The image of the magnetic structure factor can reveal information about magnetic order, evolution, and quasiparticle excitations; therefore, its calculation is of great significance. For l equals 2.3 mm, 1.5 mm, 1.1 mm and 0.7 mm, the proportion of the ground state of units is over 80%. Most units can be regarded as an Ising spin according to real positions of magnets. Following the conventional calculation method for the magnetic structure factor in artificial spin ice, the magnetic structure factor I(q) of macroscopic square ice is defined as follows.

Definition 1.

Let S_i be the pseudospin of the i-th unit, r_i be the position of the i-th unit, q be the reciprocal lattice vector, and N be the total number of pseudospins. Then the magnetic structure factor I(q) is expressed as

$$I(\mathit{q})=(1/N){\sum }_{i,j }{e}^{i\mathit{q}\cdot (\mathit{r}i- \mathit{r}j)} ⟨\mathit{S}i\cdot \mathit{S}j⟩$$

(1)

In Eq. (1), S_i denotes the spin of the i-th site. r_i is the position of site i, and q is the scattering vector. ⟨···⟩ denotes the thermal average. When q = 0, I (0) can be expressed as

I (0) = (1/N)∑_i,j⟨S_i·S_j⟩ = (1/N) ⟨(∑_iS_i)²⟩

(2)

The calculation method of the magnetic structure factor is given by Eqs. (1) and (2), and the results are shown in Figure 9, reflecting the evolution of system order from the reciprocal-space perspective.

When l = 9.7 mm, the magnetic structure factor cannot be calculated by Eqs. (1) and (2). This suggests that the evolution from this excited ordered state to the antiferromagnetic state may be related to string-like excitations. When l decreases to 2.3 mm or even 1.5 mm (Figure 9(a)), the magnetic structure factor shows sharp antiferromagnetic Bragg peaks at the Brillouin zone boundaries such as (π, π), consistent with the real-space antiferromagnetic order.

When l = 1.1 mm (Figure 9(b)), the intensity of the antiferromagnetic Bragg peaks is significantly weakened, and ferromagnetic Bragg peaks also appear at the zone center. The two types of Bragg peaks coexist, both with relatively low intensity, indicating that the system exhibits a disordered state. Notably, this state shows two types of Bragg peaks, which is distinctly different from the excited ordered state at l = 9.7 mm in terms of reciprocal-space characteristics. This indicates that the disordered state at l = 1.1 mm shows coexistence of multiple ordered components.

For l = 0.64d, Type-I vertex fraction recovers to 78.79%, partial AFM ordering reemerges. This evolution results from interaction competition. By using a uniform assembly method, we demonstrate that interaction competition governs order in macroscopic frustrated systems. The macroscopic platform enables direct optical imaging of assembled configurations, room-temperature visualization of frustrated behavior and monopole-like excitations, and it can be extended to field control, disorder effects, and finite-temperature studies. Future work could incorporate Monte Carlo simulations and analogies to artificial spin ice analyses to deepen understanding of monopole generation and dynamics, and to explore broader classes of evolution and topological excitations in frustrated systems.

When l = 0.7 mm (Figure 9(c)), the magnetic structure factor again shows a single sharp antiferromagnetic Bragg peak, consistent with the restoration of antiferromagnetic order. In the figures, r.l.u. denotes reciprocal lattice units, and a.u. denotes arbitrary units. From a physical perspective, the above non-monotonic geometric dependence can be explained by an energy competition model, namely the interaction competition. When l is large, the intra-unit interaction are weak and units tend to remain in excited states, and the system behavior is dominated by inter-vertex interactions. As l decreases, the intra-unit interaction becomes strong enough to stabilize the unit ground states, and the system establishes antiferromagnetic order. When l further approaches d, the two types of interactions become comparable in strength, and the competition reaches its maximum; at this point, the friction is relatively large, and the pseudospins of the units are prone to flipping, thereby disrupting the long-range order. When l is smaller than d, the intra-unit interaction becomes absolutely dominant, and instead drives the system to rearrange into larger antiferromagnetic domains. The above non-monotonic evolution corresponds to three cases [71]: the first corresponds to the shrinking of strings, from an excited ordered state to an antiferromagnetic state; the second is from the antiferromagnetic state to a disordered state; and the third is from the disordered state back to the antiferromagnetic state.

At the theoretical level, this evolution driven by the interaction competition, accompanied by partial breaking of the ice rule and significant changes in the monopole excitation proportion, shares several similarities with the physical picture of strings in spin ice systems [70,72,73,74,75]. In this work, the evolution from the excited ordered state at l = 9d to the antiferromagnetic ordered state at l = 2.1d, as well as the restoration from the strongly disordered state at l = d to the antiferromagnetic ordered state at l = 0.64d, exhibit several characteristic features—including the mutual suppression and enhancement of diffuse scattering and Bragg peaks, partial relaxation of topological constraints (ice rule), and significant changes in the monopole excitation concentration.

The magnetic structure factor I(q) is calculated from the perpendicular components of unit pseudo-spins, following standard artificial spin ice methodology. At l = 9.7 mm, I(q) cannot be calculated by the above definition. At l = 2.3 mm and 1.5 mm, sharp antiferromagnetic Bragg peaks appear near (π, π) in the Brillouin zone, consistent with emergent long-range AFM order. At l = 1.1 mm, AFM Bragg peaks weaken, and FM Bragg peaks indicate a disordered state with coexisting AFM and FM components. At l = 0.7 mm, a single sharp AFM Bragg peak re-emerges, consistent with AFM order restoration.

The observed evolution aligns with a phenomenological energy-competition model: when intra-unit interactions are weak (large l), units are excited and inter-vertex interactions govern behavior, favoring excited ordered ice-rule states; as l decreases, intra-unit interactions stabilize unit ground states and promote AFM order; near l ≈ d, competing interactions plus boundary/friction effects lead to unit flips and disruption of long-range order; when l = 0.64d, intra-unit interactions dominate and the system again orders antiferromagnetically.

The observed sequence shares qualitative similarities with constrained strings. However, due to lacking definitive statistics on winding numbers, energy-entropy competition quantification, and precise critical behavior, we do not claim a strict topological invariant. The magnetic structure factor exhibits diffuse patterns and evolving Bragg intensities, providing a real-space visualization of topological and energetic competition in a macroscopic frustrated system [78,79].

4. Discussion

The observed non-monotonic evolution can be understood from the competition between intra-unit and inter-vertex interactions. This work currently lacks the following key pieces of evidence: (1) systematic statistics of Dirac string winding numbers; (2) quantitative analysis of the energy-entropy competition; and (3) precise determination of critical behavior and scaling laws. Moreover, the features in the magnetic structure factors are not as pronounced as those observed in artificial spin ice. The significance of this study lies in providing clear experimental evidence and a tunable research platform for exploring string-related physical phenomena through real-space visualization.

We constructed a macroscopic square ice platform with a standardized assembly method and systematically measured vertex-type proportions, correlation functions, magnetic-domain sizes, and magnetic structure factors as a function of l. The geometric and magnetic parameter distributions (trap diameters, magnet diameters, magnetic flux densities, masses, and heights) were approximately Gaussian, indicating good sample uniformity. The key finding is that local intra-unit interaction strength provides an effective route to control global magnetic ordering in macroscopic square ice, leading to a non-monotonic evolution. As l decreases from values well above d to below d, the system displays non-monotonic order–disorder–order behavior:

At l ≈ 9d, Type-I vertices dominate, although a large fraction of units remains in excited states, leading to an excited ordered state. At l = 1.5 mm and 2.3 mm (≈ 1.4d and 2.1d), Type-I vertices dominate and units are mainly in ground states, yielding long-range AFM order. At l ≈ d, Type-I vertex fraction drops sharply to 36.36%, while Type-II/Type-III fractions rise, producing a disordered state.

5. Conclusions

In conclusion, we have developed a macroscopic square ice platform that enables direct visualization and controllable manipulation of magnetic configurations. By tuning the unit gap l, we systematically investigated how the strength of intra-unit interactions influences the collective ordering behavior of the frustrated system.

We demonstrate that the magnetic ordering does not vary monotonically with geometry. Instead, the system exhibits a non-monotonic order–disorder–order evolution as the intra-unit interaction strength changes. This behavior originates from the competition between intra-unit interactions and inter-vertex interactions, where local interaction regulation modifies the balance between ordered and disordered configurations.

Our results reveal that intra-unit interactions provide an effective degree of freedom for controlling emergent magnetic order in macroscopic square ice. Beyond simply adjusting geometric parameters, this work establishes a route for manipulating collective states in frustrated systems through local interaction engineering.

The macroscopic square ice platform developed here offers a versatile approach for exploring the relationship between microscopic interaction competition and macroscopic emergent phenomena, providing new insights into the design and control of artificial spin ice systems and square ice model.

This work successfully achieves room-temperature, direct optical imaging of assembled configurations, and direct observation of the non-monotonic evolution in a frustrated system, addressing the challenge of directly visualizing evolution processes in microscopic spin ice systems. This platform can be further extended to research directions such as magnetic-field tuning, disorder effects, and finite-temperature behavior. In the future, combined with Monte Carlo simulations [43,80,81,82,83] and drawing on the analysis methods of artificial spin ice [22], it can be used to deeply explore the generation and interaction dynamics of monopole excitations [33], providing more experimental evidence for the study of other types of evolution and topological excitations in frustrated systems [9,10,23].

Experimental fabrication: Commercial cylindrical permanent magnets were purchased from Zhanggong District Zhicheng Electronic Products Business Department. The grade of NdFeB magnets is N52, and the diameter of magnets is 4 mm. The height of magnets is 5 mm or 6 mm. These magnets are the highest-grade standard magnets, offering maximum magnetic energy of 52 MGOe but limited to low-temperature environments (< 80℃). The square lattice of traps in the PMMA plates was fabricated by continuous CO₂ laser cutting machine (λ = 10.6 μm). Vertical empty wells were etched into a 20 cm × 20 cm PMMA plate to form a square lattice, with the well depth equal to the plate thickness (10 mm). Following the well cutting, nine small holes (diameter smaller than the magnets) were additionally cut at scattered locations in the array for later fixing. Magnets were manually one by one inserted into predefined traps in the PMMA plate. The positional uncertainty of magnets was approximately ±0.05 mm.

Energy calculation: The magnetic dipole model can be used to approximately calculate the interaction energy of the magnets. If a magnet can be considered as a dipole, the interaction energy of two magnets can be calculated by

$$E={\displaystyle \frac{{\mu }_0}{4\mathsf{\pi }{\left|\mathit{r}\right|}^3}}[{\mathit{m}}_1·{\mathit{m}}_2– {\displaystyle \frac{3({\mathit{m}}_1·\mathit{r})({\mathit{m}}_2·\mathit{r})}{{\left|\mathit{r}\right|}^2}}]$$

(3)

where m₁ and m₂ are the moments of two magnets, and r is the distance between the centers of these two magnets. The magnetic moments of the two magnets are parallel and point upward.

Magnetic state imaging: The magnetic configurations were obtained by identifying the orientation of each magnet in optical images. Images were acquired by using an optical cameral. The device (Honor 80, model ANN-AN00) is equipped with a triple-rear camera system. The main camera features a 160-megapixel sensor (1/1.56-inch optical format, f/1.8 aperture) capable of capturing images with a maximum resolution of 10,944 × 14,592 pixels. The secondary rear cameras comprise an 8-megapixel ultra-wide-angle lens (f/2.2) and a 2-megapixel macro lens (f/2.4). For the front, there is a 32-megapixel selfie camera (f/2.4). The main sensor has a native pixel size of 0.56 µm and supports 16-in-1 pixel binning, which yields an effective pixel size of 2.24 µm to enhance low-light performance. For video, the rear camera supports recording up to 4K resolution (3840 × 2160 pixels). The polarization direction of a magnet can be identified by means of the WT-60 NS magnetic pole indicator.

Error estimation: As shown in Extended Figure 1(a), at a positional deviation of -0.1 mm, the vertex energies of the three vertex types slightly increase; conversely, at +0.1 mm in Extended Figure 1(b), they slightly decrease. The energy changes observed in the two panels are nearly identical, with very close numerical values. Taken together, these energy data indicate that the non-monotonic dependence in vertex propositions does not arise from local positional deviations.

Extended Figure 1. Vertex energies for three vertex types at ±0.1 mm positional deviation. (a) Vertex energies for -0.1 mm positional deviation. (b) Vertex energies for +0.1 mm positional deviation.

Extended Figure 1. Vertex energies for three vertex types at ±0.1 mm positional deviation. (a) Vertex energies for -0.1 mm positional deviation. (b) Vertex energies for +0.1 mm positional deviation.