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

We study a frustrated two-dimensional array of dipoles forming an artificial rectangular spin ice with horizontal and vertical lattice parameters given by $a$ and $b$ respectively. We show that the ice regime could be stabilized by appropriate choices for the ratio $\gamma \equiv a/b$. Our results show that for $\gamma \approx \sqrt{3}$, i.e., when the center of the islands form a triangular lattice, the ground state becomes degenerate. Therefore, while the magnetic charges (monopoles) are excitations connected by an energetic string for general rectangular lattices (including the particular case of a square lattice), they are practically free to move for a special rectangular lattice with $\gamma \approx \sqrt{3}$. Besides that, our results show that for $\gamma>\sqrt{3}$ the system is highly anisotropic in such a way that, even for this range out of the ice regime, the string tension almost vanishes along a particular direction of the array. We also discuss the ground state transition and some thermodynamic properties of the system.


Introduction
In the recent years, a great deal of efforts has been dedicated to the study of materials with frustrated interactions in an attempt to find and understand new states of matter [1,2,3,4,5,6,7,8,9]. In particular, for ferromagnetic materials, two geometric structures are often investigated: the three-dimensional (3d) pyrochlore spin ice [2,3,4,5,10] and the two-dimensional (2d) analog thereof, the artificial kagome spin ice [11,8,9,12,13,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 [15,5,8]. The artificial two-dimensional square spin ice is another structure very well studied, where monopoles excitations also appear [6,16,17,18,19,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 [21,16,19], 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 1 ≈ 0.444b (b is the lattice spacing) [22,16]. Indeed, the string tension connecting the opposite monopoles decreases as h increases from zero, almost vanishing at h 1 . 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 three-dimensional and, in principle, it is considerable more difficult to be realized artificially. To our knowledge, such a system was not built yet.
On the other hand, there is a simpler way of dealing with this difficulty without transforming the array in a three-dimensional 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 two-dimensional 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 Figs.1(a)-(d)). Really, such a deformation can tune the ratios of the interactions between neighboring elements resulting in different magnetic ordering of the system as shown by Li at al. [24]. Then, the ground state of this system is modified if the ratio of show the 16 possible magnetic moment configurations on a vertex and the five distinct topologies. In (b) e (c) we illustrate the T 0 and T 1 topologies and the residual charge due the dipole fracionalization. In (d) we show the configuration of a state denoted as state-0; the configuration of the dipoles looks like the ground state of the usual square spin ice. One can also see other configurations denoted as state-1, state-2 and state-3. The state-1 represents the ground state for a/b > √ 3.
vertical to horizontal lattice spacing is √ 3 (or even 1/ √ 3). 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.

Model and Methods
At each site (x i , y i ) of the rectangular array ( Fig. 1(d)), two spin variables are defined: S 1 = ±(1, 0) located at r 1 = (ax i + a/2, by i ) and S 2 = ±(0, 1) located at r 2 = (ax i , by 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 equals (a = b) we then recover the square lattice [6,7,19,17]. Therefore, in a lattice of area A = l 2 b 2 one gets 2 × l 2 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 S i , then the system is described by the following Hamiltonian where D = µ 0 µ 2 /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 Fig. 1(a)). We remark that the square lattice admits only four different topologies [6,7,17]. Indeed, for the particular case in which γ = 1, topologies T 2 and T 3 have the same energy; on the other hand, for any γ = 1, T 2 and T 3 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 0 , 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 Ref. [5] can not 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 0 , we can associate a residual charge with modulus Q =| 2µ/a − 2µ/b | (see Fig.1(b)). Of course, the residual charges present in topology T 0 vanish for the square lattice in which a = b. Note that topology T 1 , which also obeys the ice rule, does not exhibit a residual charge ( Fig.1(c)); it has, however, a residual magnetic moment (here, referred to as "SPIN"), not present in topology T 0 . There is therefore an interesting asymmetry for the two topologies that obey the ice rule: on the vertices of T 0 , one can find residual magnetic charges but not "SPINS" while on the vertices of T 1 one finds residual "SPINS" but the neutrality of magnetic charges is preserved. As a/b is increased, the magnetic charge Q (for T 0 ) increases while the effective "SPIN" (for T 1 ) decreases. Such properties of the vertices obeying the ice rule (charge for T 0 or "SPIN" for T 1 ) may play important roles in the ground state of the rectangular array. It may eventually exist a balance of the values of charge in T 0 and "SPIN" in T 1 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. Table  (1) shows the residual vertex charge Q (in the dumbbell picture) as a function of the lattice parameters a and b, for all topologies presented in Fig.1a. An interesting point is that topologies T 2 and T 3 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 2 has more Table 1. Correspondence between topologies and charges in the dumbbell picture. In this work b is kept fixed while a is varied in such a way that the charges of topologies T 0 , T 2 and T 4 varies as long as γ = a/b changes.

T opology
Charges charge than vertices on topology T 3 . This difference has important consequences on charges interactions as will be shown soon. Of course, as stated above, the dumbbell picture can not be used to quantitatively describe artificial spin ices such that we do not expect that the charges dependence on lattice spacings a and b are indeed those shown on Table 1.
Here, we perform standard Monte Carlo techniques to obtain the thermodynamics 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 [25,19], 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 2 /γ single spin flips and 1 loop move. Usually, 10 4 MCS were shown to be sufficient to reach equilibrium configuration and we have used 10 4 − 10 6 MCS to get thermodynamic averages.

Ground states
It is well known that, on a square lattice (γ = 1), the topology T 0 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 0 as illustrated in Fig. 1(d) state-0. In this case the ground state is only twofold degenerate. It is also valid for an appreciable range of values (γ = 1) of the rectangular ice. Nevertheless, the energy of topologies T 0 and T 1 approaches each other when we set the dipoles on a rectangular geometry such that γ → 0.559(1) or γ → 1.79(1) (see Fig. 2(a)). 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 (γ c1 , γ c2 ). Moreover, for γ c1 < γ < γ c2 , the ground state is composed by topology T 0 (possessing residual charges alternating the signs at adjacent vertices) and for γ < γ c1 or γ > γ c2 , the ground state is dominated by topology T 1 (possessing residual "SPINS" located at the vertices). There are then two equivalent ground state transitions on a rectangular lattice. Figure 2(b) shows the energies of some special configurations of spins obeying the ice rule (referred to as states 0, 1, 2 and 3) as a function of γ. The configurations of these four states can be seen in Fig. 1(d). State-0 is the ground state of a square array (and also of a rectangular array for γ c1 < γ < γ c2 ) with all vertices described by topology T 0 ; State-1 is the ground state of a rectangular array with γ < γ c1 or γ > γ c2 with all vertices described by topology T 1 . There are other configurations that obey the ice rule with topology T 1 but they are not ground states since they have higher energies than that of the 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. A calculation considering the dipolar interaction between all dipoles of the lattice shows that the values of γ are slightly modified to γ c1 ≈ 0.576 (5) and γ c2 ≈ 1.735(5) (see Fig 2( 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 √ 2b  between the two adjacent spins pointing along perpendicular directions. However, if the lattice is deformed toward the rectangular shape, there are two relations concerning the parameters a and b which give the same distances for the above mentioned situations, i.e., a = √ 3b ≈ 1.732b or a = b/ √ 3 ≈ 0.577. Choosing the parameters a and b in such a way that the spins around a vertex become equidistant should lead to interaction energies approximately equal. The values found by this simple argument are very near the numerical results. 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. For this reason we shall consider in this study 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 the Fig.  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 favorable than residual "SPINS" found on topology T 1 .
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 Fig. 1(d). For this range, the residual "SPINS" are more energetically favorable than the residual charges of topology T 0 . The balance between local residual charges (found in topology T 0 ) and local residual "SPINS" (found in topology T 1 ) must be eventually established (from the energetic point of view) when γ = γ c2 , making these two topologies to become degenerate.

Excitations and monopole deconfinement
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 in 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 4 in the Fig. 1(a)) and the lower energetic ones, in which the spins are in the (3 − in, 1 − out) or (3 − out, 1 − in) states (topologies T 2 and T 3 shown in Fig. 1(a)). It is important to remark that topologies T 2 and T 3 have different charges and energies for general γ = 1 (see Table 1). 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 Ref. [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 2 and excitation 1h creates two adjacent vertices in topology T 3 ). 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 Fig. 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 4P h e 4P v 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 on Table 1. However, as said before, the dumbbell picture does not quantitatively describe the system and thus, the dependence on γ is not really that of Table 1. Indeed, we could not find a simple law relating the energetic cost with the magnetic charge of the excitations (which will be considered latter).
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 monopoles interaction should be highly dependent of 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 2 and T 3 ). However, the quantity that measures the interaction between monopoles involves the product between these charges, i.e., the quantity q of Refs. [7,16] should be somewhat like q = q 1 q 2 where q 1 and q 2 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 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-axes of the array, and q is the product between monopoles 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 Refs. [7,16], two different shapes (2L and 4P ) as the building blocks of the string, such that X = √ 2R 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 an horizontal spin, a vertical one must be flipped and so on. In this context, using the same string shapes of Refs. [7,16] (2L and 4P ), we always have at the end points of a string, a type T 2 vertex and a type T 3 vertex for excitations above state-0. So, in this case, the constant q is the product of different types of monopoles (q = q T 2 q T 3 ). We remark that in the separation process, the moving monopole changes the topology (from T 2 to T 3 and vice-versa) in such a way that during this procedure, there will be interactions between the same kind of vertices (T 2 with T 2 for example) alternating with interactions between different kinds of vertices; however, we only keep in the expression of the potential used to fit Eq. 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 4P h and 4P v strings.
Besides that, the string in the rectangular lattice is composed by T 1 vertices, which has neutral charge. However, since the ground state (with the configuration of state-0) is a charge's 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 2 and T 3 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 1 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 2 or T 3 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 1 vertices, while, along the region it bends (making a corner) the topology of the string changes to T 0 and thus the monopole changes from topology T 2 to T 3 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 excitations) as shown in Fig. 5 picture (3) of Ref. [16]. Figure 4 shows the dependence of the parameters q and κ on γ. The parameter q for 1 < γ < γ c2 is q = q T 2 q T 3 since, as discussed above, the monopoles generated by the string shapes 4P h, 4P v and 2L are described by topologies T 2 and T 3 . This quantity has a clear dependence on γ; 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 1 , 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 2 and thus one may expect a 1/γ 2 dependence, as it really is. For string 1v both monopoles are in topology T 3 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 Fig.4).

Thermodynamics
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 Ref. [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 B T ≈ 7.2D as discussed at Ref. [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 Fig. 5(b)). The fact that the string tension also decreases as γ increases, practically vanishing at γ = √ 3, corroborates the argument of Refs. [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 B T /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 Figs. 4 and 5, their dependence on γ is very similar.

Conclusions
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 = √ 3 and a/b = 1/ √ 3 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 = √ 3 or a/b = 1/ √ 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 two-dimensional artificial ices not due to thermal effects, but due to the lattice construction itself. Moreover, we have investigated some aspects of the thermodynamic 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 in the search for artificial spin ice systems where monopoles excitations can be controlled. Indeed, the dependence of the excitation's energies and consequently of the constants q(φ, γ), κ(φ, γ) and c(φ, γ) on γ (see Figs. 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 Refs. [27,28,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 [22,11,16] and now, as more one example, in the artificial rectangular ice with special ratio a/b = √ 3.