From confinement to deconfinement of magnetic monopoles in artificial rectangular spin ices

In recent years, a great deal of effort has been dedicated to the study of materials with frustrated interactions in an attempt to find and understand new states of matter [1–9]. In particular, for ferromagnetic materials, two geometric structures are often investigated: the three-dimensional (3D) pyrochlore spin ice [2–5, 10] and the two-dimensional (2D) analogue thereof, the artificial kagome spin ice [8, 9, 11–14]. These systems share the common phenomenon of a highly degenerate classical ground state. Besides, they also have excitations that are similar to Dirac magnetic monopoles [5, 8, 15]. The artificial 2D square spin ice is another very well studied structure, where monopole excitations also appear [6, 16–20]; however, the square ice does not exhibit a degenerate ground state, since the two topologies that obey the ice rule have different energies [6]. The main reason for this difference is the fact that, unlike the case of a tetrahedron in the natural 3D spin ice, the six bonds between the four islands belonging to a vertex are not all equivalent. As a first consequence, the square array provides a different type of magnetic monopole-like excitation: a kind of Nambu pair of monopole–antimonopole [16, 19, 21], in which the opposite charges are effectively interacting by means of the usual Coulombic law plus a linear confining potential, the latter being related to a stringlike excitation binding the monopoles [7, 16]. Therefore, these monopoles are confined by a string similar to quark confinement in quantum chromodynamics [21].

In a recent paper Möller and Moessner [22] showed that the problem related to the two inequivalent bond energies in the square ice could be remedied by introducing a height displacement h between magnetic islands pointing in the horizontal and vertical directions. In this modified lattice, the ground state changes its configuration at a critical value of h given by h = h₁ ≈ 0.444b (b is the lattice spacing) [16, 22]. Indeed, the string tension connecting the opposite monopoles decreases as h increases from zero, almost vanishing at h₁. At this particular value of h, the ground state should become degenerate and the monopoles would become deconfined, similar to those found in the natural 3D spin ice compounds [5]. However, there is a price to be paid. The modified array is also 3D and, in principle, it is considerably more difficult to be realized artificially. To our knowledge, such a system has not been built yet.

On the other hand, there is a simpler way of dealing with this difficulty without transforming the array into a 3D one. Indeed, one could mimic techniques used for natural systems in which the lattice spacing can even be tuned continuously by applying pressure along a particular direction [23]. Here, for a 2D array, one does not need to apply any pressure to deform the lattice; it is sufficient to fabricate samples with an intentional change in either the horizontal or the vertical lattice spacing on the original square spin ice ('compressing' or 'stretching' the lattice, see figures 1(a)–(d)). Really, such a deformation can tune the ratios of the interactions between neighboring elements resulting in a different magnetic ordering of the system as shown by Li et al [24]. Then, the ground state of this system is modified if the ratio of vertical to horizontal lattice spacing is $\sqrt {3}$ (or equivalently $1/\sqrt {3}$), which is equivalent to placing the centers of the nanoislands in a triangular lattice. This critical value leads to a degenerate ground state, suggesting a residual entropy at absolute zero temperature similar to what happens to the natural and artificial 3D spin ice materials. In this paper we study the ground state, elementary excitations and thermodynamics of rectangular lattices which are characterized as a function of the assigned lattice anisotropy.

Zoom In Zoom Out Reset image size

At each site (x_i,y_i) of the rectangular array (figure 1(d)), two spin variables are defined: $\skew3\vec{S}_1 = \pm (1,0)$ located at $\vec {r}_1 = (a x_i + a/2, b y_i)$ and $\skew3\vec{S}_2 = \pm (0,1)$ located at $\vec {r}_2 = (a x_i,b y_i + b/2)$. The parameters a and b denote the horizontal and vertical lattice spacings, respectively (here, b will be the standard lattice spacing, while a is varied to give different rectangular arrays). When the parameters are equal (a = b) we then recover the square lattice [6, 7, 17, 19]. Therefore, in a lattice of area A = l²b² one gets 2 × l²b/a spins. In this work we study the ratio γ = a/b from 0.4 to 4.0 and we have fixed b = 1. Representing the spins by $\skew3\vec{S}_i$, then the system is described by the following Hamiltonian:

$$\begin{equation} H_{SI} = D \sum_{i \neq j}\left[\frac{\vec{S}_i \cdot \vec{S}_j}{r_{ij}^3}- \frac{3 (\vec{S}_i \cdot \vec{r}_{ij})(\vec{S}_j \cdot \vec{r}_{ij})}{r_{ij}^5}\right], \end{equation}\noindent \tag{ 1 }$$

where D = μ₀μ²/4π is the coupling constant of the dipolar interactions.

Each rectangular lattice vertex has 16 possible magnetic moment configurations that can be characterized by five distinct topologies denoted by T_ν (ν = 0,1,2,3,4, see figure 1(c)). We remark that the square lattice admits only four different topologies [6, 7, 17]. Indeed, for the particular case in which γ = 1, topologies T₂ and T₃ have the same energy; on the other hand, for any γ ≠ 1, T₂ and T₃ have different energies.

Another important difference between the cases γ ≠ 1 and the square lattice (γ = 1) is the emergence of residual magnetic charges even for topology T₀, which obeys the ice rule where two spins point inward and two spins point outward (2-in, 2-out) on a vertex. Although the dumbbell picture proposed in [5] cannot be simply transposed to the artificial square spin ice (since it does not describe the system quantitatively), it can help us to understand (qualitatively) some differences between the square and rectangular spin ices. The dumbbell picture [5] is obtained when each spin in the lattice is replaced by a pair of magnetic charges placed on the adjacent vertices. Due to the lattice anisotropy, the horizontal spins should fix opposite magnetic charges (at adjacent vertices) given by q_h = ± μ/a. On the other hand, the vertical spins would give adjacent opposite charges given by q_v = ± μ/b. Therefore, at a vertex (x_i,y_i) described by topology T₀, we can associate a residual charge with modulus Q = |2μ/a − 2μ/b| (see figures 1(a) and (c)). Of course, the residual charges present in topology T₀ vanish for the square lattice in which a = b. Note that topology T₁, which also obeys the ice rule, does not exhibit a residual charge (figure 1(b)); it has, however, a residual magnetic moment (here, referred to as 'SPIN'), not present in topology T₀. There is therefore an interesting asymmetry for the two topologies that obey the ice rule: on the vertices of T₀, one can find residual magnetic charges but not 'SPINS', while on the vertices of T₁ one finds residual 'SPINS' but the neutrality of magnetic charges is preserved. As a/b is increased, the energy of the magnetic charge Q (for T₀) increases while the energy of the effective 'SPIN' (for T₁) decreases. Such properties of the vertices obeying the ice rule (charge for T₀ or 'SPIN' for T₁) may play important roles in the ground state of the rectangular array. There may eventually exist a balance of the values of charge in T₀ and 'SPIN' in T₁ in such a way that the energies of these two topologies become equal for a particular a/b, causing a degenerate ground state. This possibility is studied in the next section. In figure 1(c) the residual vertex charge Q (in the dumbbell picture) is shown as a function of the lattice parameters a and b. An interesting point is that topologies T₂ and T₃ are degenerate in the square lattice (γ = 1) and can be grouped in a single topology. However, as long as γ ≠ 1 they have different energies and charges. For instance, considering γ > 1 (a > b), vertices on topology T₂ have more charge than vertices on topology T₃. This difference has important consequences on charge interactions as will be shown soon. Of course, as stated above, the dumbbell picture cannot be used to quantitatively describe artificial spin ices such that we do not expect that the charge dependences on lattice spacings a and b are indeed those shown in figure 1(c).

It is well known that, on a square lattice (γ = 1), the topology T₀ is more energetically favorable than the others. Consequently, the ground state of an artificial square spin ice has a configuration with all vertices obeying the ice rule with topology T₀ as illustrated in figure 1(d) state-0. In this case the ground state is only two-fold degenerate. It is also valid for an appreciable range of values (γ ≠ 1) of the rectangular ice. Nevertheless, the energies of topologies T₀ and T₁ approach each other when we set the dipoles on a rectangular geometry such that γ → 0.556 or γ → 1.797 (see figure 2(a)). Indeed, the analytical expression for the energy difference (ΔE) between a vertex on topology T₀ and a vertex on topology T₁ as a function of γ is

$$\begin{eqnarray*} &&\Delta E=4\left(1+\frac{1}{\gamma^3}-\frac{24\gamma}{(1+\gamma^2)^{5/2}}\right), \end{eqnarray*}$$

which gives the above numerical results. Then, there are two special ratios γ that may degenerate the ground state of the system, since the topologies satisfying the ice rule have the same energy for γ = 0.556 and 1.797. However, the values of γ in which the system is degenerate may not be those obtained in a single vertex analysis, since long range interactions are present. To get the correct values we have evaluated the energy (see equation (1)) for the four different ordered configurations of spins as shown in figure 1(d) using different lattice sizes. A typical curve, obtained for a lattice with 800 spins, is shown in figure 2(b). From a finite size scaling analysis we see that the values are $\gamma _{{\mathrm { c}}1} =1.732\,87(5)\approx \sqrt {3}$ and $\gamma _{{\mathrm { c}}2}=0.577\,07(2)\approx 1/\sqrt {3}$. Moreover, for γ_c1 < γ < γ_c2, the ground state is composed by topology T₀ (possessing residual charges alternating the signs at adjacent vertices) and for γ < γ_c1 or γ > γ_c2, the ground state is dominated by topology T₁ (possessing residual 'SPINS' located at the vertices). There are then two equivalent ground state transitions on a rectangular lattice. State-0 (see figure 1(d)) is the ground state of a square array (and also of a rectangular array for γ_c1 < γ < γ_c2) with all vertices described by topology T₀; state-1 is the ground state of a rectangular array with γ < γ_c1 or γ > γ_c2 with all vertices described by topology T₁. There are other configurations that also obey the ice rule with topology T₁ but they are not ground states since they have higher energies than that of state-1. Note also that state-1 has null total magnetization while state-2 is magnetized along the 'vertical' direction and state-3 is magnetized along the diagonal. Therefore, the ground state of the rectangular lattice is either, a state with residual charge at each vertex but with null net charge or a state with residual spins at each vertex but with null net magnetization.

Zoom In Zoom Out Reset image size

A simple argument can lead to the above approximate quantities: by replacing the net magnetic moment of the islands by a point-like dipole at their centers, one can easily see that, for any vertex of a square array, the distance b between the two adjacent spins pointing along the same direction is smaller than the distance $\sqrt {2}b$ between the two adjacent spins pointing along perpendicular directions. However, if the lattice is deformed toward the rectangular shape, the lines connecting the spins pointing along perpendicular directions form a rhombus with a short diagonal of length b and a long diagonal of length a. For the special ratios $\gamma =\sqrt {3}$ or $1/\sqrt {3}$, two back to back equilateral triangles are formed, in such a way that the distance between spins pointing along different directions (the length of the sides of the rhombus) becomes equal to the distance between one of the two pairs of spins pointing along the same direction (the length of the rhombus' short diagonal). Then, for this special case, the energies of the different spin configurations become degenerate. The square spin ice was thus transposed to a triangular lattice.

As we have seen, the ground state configuration on rectangular lattices has a dependence on the parameter γ. However, due to rotational symmetry of the Hamiltonian, the process of compressing or stretching the lattice is equivalent to a deformation. Indeed, the parameters γ_c1 and γ_c2 are just the reverse of each other, since it is enough to rotate the sample by an angle π/2 to pass from γ to 1/γ. For this reason we shall consider from now on only the stretching process, i.e. γ ⩾ 1. Thus, for 1 < γ < γ_c2, the ground state of a rectangular array assumes the same configuration of the square spin ice ground state (state-0), as illustrated in figure 1(d). In particular, for this range, the ground state exhibits an interesting ordering of alternating positive and negative charges Q on the vertices. In addition, such residual charges are energetically more favorable than residual 'SPINS' found on topology T₁.

On the other hand, for γ > γ_c2, the lattice anisotropy forces the system to assume a modified ground state (state-1) which is also illustrated in figure 1(d). For this range, the residual 'SPINS' are more energetically favorable than the residual charges of topology T₀. The balance between local residual charges (found in topology T₀) and local residual 'SPINS' (found in topology T₁) must be eventually established (from the energetic point of view) when γ = γ_c2, making these two topologies become degenerate.

The inversion of a single spin (for example, between adjacent vertices i and i + 1) violates the ice rule generating an excited state, which implies an excess of opposite magnetic charges placed on i and i + 1. It is a pair of defects similar to monopole quasi-particles, positioned at two adjacent vertices. It is useful here to distinguish the different types of monopole defects: the most energetic ones, in which the spins are in the (4-in) or (4-out) states (T₄ in figure 1(c)) and the lower energetic ones, in which the spins are in the (3-in,1-out) or (3-out,1-in) states (topologies T₂ and T₃ shown in figure 1(c)). It is important to remark that topologies T₂ and T₃ have different charges and energies for general γ ≠ 1 (see figure 1(c)). Therefore, the system may be filled by several elementary excitations that can be classified by the number of flipped moments and a mnemonic character for shape as introduced in [17]. However, here, there is a small difference that introduces a new factor in the classification: separations of the magnetic charges along the horizontal (h) and vertical (v) lines of the array have different energies. As examples, the simplest excitations are symbolized by 1v and 1h, representing two opposite charges separated by one lattice spacing along the vertical and horizontal lines respectively (excitation 1v creates two adjacent vertices in topology T₂ and excitation 1h creates two adjacent vertices in topology T₃). The separation process of the monopoles generates energetic one-dimensional strings of dipoles that can be seen as lines which pass by adjacent vertices that obey the local ice rule. In figure 3 we present the symbol of some pairs of monopoles and their energy cost as a function of γ. The behavior of the energy cost as a function of γ is not trivial, but, in general, for γ < γ_c2 we observe that the energy cost of the defect decreases as the lattice is stretched. Note that we did not show the energies of the excitations 4Ph and 4Pv for γ > γ_c2 since these excitations generate extra charges on the ground state. One might naively think that the energetic cost of these excitations should follow the dependence on lattice spacings shown in figure 1(c). However, as said before, the dumbbell picture does not quantitatively describe the system and thus, the dependence on γ is not really that of figure 1(c). Indeed, we could not find a simple law relating the energetic cost with the magnetic charge of the excitations (which will be considered later).

Zoom In Zoom Out Reset image size

In recent works [7, 16] considering the square array, we have proposed that the most general expression for the total cost of a monopole–antimonopole pair separated by a distance R is the sum of the usual Coulombic term roughly equal to q/R, and a term roughly equal to κX resulting from the string joining the monopoles (here, X is the string length; there is still a constant term associated with the pair creation energy [7, 16]). The parameters q and κ have a small dependence on the direction in which the monopoles are separated in the crystal plane [16]. This anisotropy is therefore manifested in both the Coulomb and linear terms of the potential. Although the monopole interaction should be highly dependent on the stretching of the lattice, one should expect that for the rectangular lattice, in general, things must have some similarities with the square lattice. The main difference is that the square array has only one kind of unit-charged monopole while the rectangular array has two (topologies T₂ and T₃). However, the quantity that measures the interaction between monopoles involves the product between these charges, i.e. the quantity q of [7, 16] should be somewhat like q = q₁q₂ where q₁ and q₂ are the charges of each monopole (the remaining term, associated with the string tension, is supposed to be independent of the charges itself). So, as long as we can always deal with the same kind of charges (considering for example only the interactions between vertices of the same kind) or even between vertices of different types, we expect that the same phenomenology of the square spin ice applies to the artificial rectangular spin ice. Indeed, the function

$$\begin{equation} V(R) = \frac{q(\phi,\gamma)}{R} + \kappa(\phi,\gamma)X + c(\phi,\gamma) \end{equation} \tag{ 2 }$$

is expected to describe the interactions between the monopoles in both cases (here, ϕ is the angle that the line joining the monopole-like defects makes with the x-axis of the array and q is the product between monopole charges).

Before showing the results, a detailed discussion on how they were obtained would be useful. First we would like to remark that the potential V (R) is simply the energy of each configuration of spins minus the ground state energy. Each point of this curve is obtained by evaluating the system's energy for a given configuration for which the distance between the monopoles is R. One very important point of our calculations is that we must establish a connection between the string length (X) and the distance between the monopoles (R), otherwise we would have two variables in our model. This is why we have chosen, in [7, 16], two different shapes (2L and 4P) as the building blocks of the string, such that $X=\sqrt {2}R$ for the 2L strings and X = 2R for 4P strings. Note also that we have to keep the neutrality in the region between the two opposite monopoles (in the context of two-in/two-out, while flipping spins in sequence); then, after flipping a horizontal spin, a vertical one must be flipped and so on. In this context, using the same string shapes of [7, 16] (2L and 4P), we always have at the end points of a string, a type T₂ vertex and a type T₃ vertex for excitations above state-0. So, in this case, the constant q is the product of different types of monopoles (q = q_T2q_T3). We remark that in the separation process, the moving monopole changes the topology (from T₂ to T₃ and vice versa) in such a way that during this procedure, there will be interactions between the same kind of vertices (T₂ with T₂ for example) alternating with interactions between different kinds of vertices; however, we only keep in the expression of the potential used to fit equation (2), the interactions between different kinds of vertices (since this is the only way of establishing a precise relationship between the string length X and the distance between the monopoles R). One extra ingredient about the rectangular lattice is the fact that the separation of monopoles along the x or y direction is different, such that we have introduced the notation 4Ph and 4Pv strings. Besides that, the string in the rectangular lattice is composed by T₁ vertices, which has a neutral charge. However, since the ground state (with the configuration of state-0) is a charge crystal, the string itself must have a net magnetic charge. Indeed, this net charge of the string is responsible for keeping the system's neutrality, once T₂ and T₃ vertices have different charges, in such a way that the composite object, string plus monopoles, is neutral. For excitations above state-1, things are less complicated since the monopoles can be separated along a straight line by a sequence of type T₁ vertices of opposite 'SPIN', which are neutral in the ground state. In this way, only one kind of monopole can be considered (interaction between either T₂ or T₃ vertices) and the string length X is equal to the distance between monopoles. Of course, it is also possible to chose other shapes for the string and separate the charges along other directions. In this case, straight portions of the string are composed by T₁ vertices, while, along the region it bends (making a corner) the topology of the string changes to T₀ and thus the monopole changes from topology T₂ to T₃ or vice versa. Indeed, we have also considered strings where the building block is a 4S excitation (formed by joining an 1h with two consecutive 1v and another 1h excitation) as shown in figure 5 picture (3) of [16].

Figure 4 shows the dependence of the parameters q and κ on γ. The parameter q for 1 < γ < γ_c2 is q = q_T2(q)_T3 since, as discussed above, the monopoles generated by the string shapes 4Ph, 4Pv and 2L are described by topologies T₂ and T₃. This quantity has a clear dependence on γ (the expected dependence from the dumbbell picture is 1/γ as one can see in figure 1(c)); however, for each string shape, a different dependence is observed. This feature is not completely understood at the moment. The string tension is almost the same for these three string shapes in this region and diminishes as γ increases, showing only a small anisotropy that increases with increasing γ. For γ > γ_c2, the string is a sequence of adjacent vertices with topology T₁, with its local 'SPINS' in different directions as discussed above. In this case, the dumbbell picture describes well the charge's dependence on γ. Indeed, for 1h and 4S strings, both monopoles are in topology T₂ and thus one may expect a 1/γ² dependence, as it really is. For string 1v both monopoles are in topology T₃ and no dependence on γ should be expected as shown. The string tension is constant in this region for strings 1h and 1v and it increases with increasing γ for 4S string, such that the system is anisotropic. Then, for a rectangular array with γ ⩾ γ_c2, the monopoles become deconfined for separation along the vertical direction while they still remain a bit confined for separation along the horizontal direction (there is a small but finite string tension along this direction; see figure 4).

Zoom In Zoom Out Reset image size

Here, we perform standard Monte Carlo techniques to obtain the thermodynamic averages of the system defined by Hamiltonian (1). Periodic boundary conditions were implemented and our Monte Carlo procedure includes a combination of single spin flips and loop moves [19, 25], where all spins contained in a closed random loop are flipped according to the Metropolis prescription. In our scheme, one Monte Carlo step (MCS) consists of 2 × l²/γ single spin flips and one loop move. Usually, 10⁴ MCS were shown to be sufficient to reach equilibrium configuration and we have used 10⁴–10⁶ MCS to get thermodynamic averages. All results shown here were obtained for l = 20.

Figure 5(a) shows the specific heat as a function of the temperature and γ. The sharp peak in the specific heat suggest that the system undergoes a phase transition as shown in [19] for the square lattice (γ = 1). Indeed, the height of these peaks increases logarithmically with the system size [19]. The transition seems to be characterized by the appearance of several excitations (Nambu monopoles and string loops). Particularly, the monopoles may become free at k_BT ≈ 7.2D as discussed in [19]. On the other hand, we can see that the lattice anisotropy changes the temperature of the transition and the peak in the specific heat has different positions and intensities, depending on the parameter γ. The critical temperature decreases as γ increases from 1 to γ_c (see figure 5(b)). The fact that the string tension also decreases as γ increases, practically vanishing at $\gamma =\sqrt {3}$, corroborates the argument of [7, 19] that the phase transition argued for the square arrays is associate with the string configurational entropy: for X sufficiently large, the number of configurations of strings connecting two opposite monopoles would be well approximated by the random walk result p^X/b (for a 2D square lattice, p = 3). Then, using a very simple estimate, the string configurational entropy (k_B ln[p^X/b]) is proportional to X [19], and the string free energy can be approximated by F = [κ − (ln 3)k_BT/b]X, which leads to an estimate of the critical temperature as T_c ∼ bκ/ln(3). Therefore, T_c∝κ and, as one can observe in figures 4 and 5, their dependence on γ is very similar.

Zoom In Zoom Out Reset image size

In summary we have studied a possible artificial spin ice array in two dimensions with rectangular geometry. We have carried out the investigation of non-thermally activated order–disorder transitions on these possible artificial spin ice structures. By varying the ratio between vertical and horizontal lattice spacings (a/b) we have obtained several properties of the system, which may depend on the competition between residual charges and residual 'SPINS' in the ground states. Such rectangular arrangement also contains elementary excitations which are similar to Nambu monopoles [21]: opposite charges confined by energetic strings like quarks in quantum chromodynamics. However, the special values $a/b=\sqrt {3}$ and $1/\sqrt {3}$ (which correspond to placing the islands in a triangular lattice) approach the array to the ice regime, increasing the frustration of the system which allows greater mobility for the monopoles. Then, in contrast to the square lattice [7, 16], monopoles may become free to move, even for very low temperatures, in a rectangular array with $a/b=\sqrt {3}$ or $1/\sqrt {3}$, leading to the possibility of monopole controls, opening an avenue for new technologies. So, here we have shown that there is a possibility of the existence of deconfined poles in 2D artificial ices not due to thermal effects, but due to the lattice construction itself. Moreover, we have investigated some aspects of the thermodynamics of this system, which corroborate the idea of a phase transition induced by the string configurational entropy [19]. In addition, we have shown that this particular change in the system's geometry can be interpreted in terms of residual magnetic charges at each vertex. Indeed, this picture is similar to that observed for a square spin ice where one of the islands is deformed, being smaller or greater than the others [26]. This may indicate that the properties of an imperfect system may be modeled by the presence of residual charges on the vertices. Besides, we remark that changes in the system's geometry may be a simple and effective way to search for artificial spin ice systems where monopole excitations can be controlled. Indeed, the dependence of the excitation's energies and consequently of the constants q(ϕ,γ), κ(ϕ,γ) and c(ϕ,γ) on γ (see figures 3 and 4) may change the system's dynamics for different aspect ratios. For instance, the strength of the system's anisotropy depends on γ and this fact may facilitate the dynamic process that drives the system to its ground state in the presence of external magnetic fields (for more details concerning the system dynamics in the presence of external fields, see [27–29]). Finally, we would like to remark that an important fact to stabilize the ice regime is the equidistance between spins around the vertex centers. Indeed, one can also see it in other systems such as the tetrahedron of natural spin ices, in the artificial kagome and modified square spin ices [11, 16, 22] and now, as a further example, in the artificial rectangular ice with special ratio $a/b=\sqrt {3}$.

The authors thank CNPq, FAPEMIG and CAPES for financial support.