# Continuous degeneracy of the fcc kagome lattice with magnetic dipolar interactions

###### Abstract

Results are presented on analytic and computational analyses of the spin states associated with a 3D fcc lattice composed of ABC stacked kagome planes of magnetic ions with only long-range dipole-dipole interactions. Extending previous work on the 2D kagome system, where discrete six-fold discrete degeneracy of the ground state was revealed [Holden et al. Phys. Rev. B 91, 224425 (2015)], we show that the 3D lattice exhibits a continuous degeneracy characterized by just two spherical angles involving six sublattice spin vectors. Application of a Heat Bath Monte Carlo algorithm shows that thermal fluctuations reduce this degeneracy at very low temperature in an order-by-disorder process. A magnetic field applied along directions of high symmetry also results in lifting the continuous degeneracy to a subset of states from the original set of ground states. Metropolis Monte Carlo simulation results are also presented on the temperature and system size dependence of the energy, specific heat, and magnetization, providing evidence for a phase transition at T 0.38 (in units of the dipole strength). The results can be relevant to a class of magnetic compounds having the AuCu crystal structure.

## I Introduction

Geometrically frustrated spin systems are typically associated with short-range Heisenberg-like antiferromagnetic exchange interactions.frustration () A more subtle geometry-induced frustration can also occur with only dipole coupling between spin vectors. Features which lead to this frustration include an antiferromagentic (AF)-like contribution to the dipole interaction, the first term in Eq. (1) (where D is the dipole strength), as well as the second term which involves explicit coupling of spin and lattice vectors. The long-range nature of the dipole coupling gives additional aspects (equivalent to a decaying higher-neighbor exchange in the first term) which are usually revealed only through numerical calculations. We adopt the dipole interaction of the form

(1) |

where is the dimensionless vector connecting sites and , and is a unit spin vector. For convenience, all calculated quantities involving the energy have units relative to the dipole strength. In the present work, spin structures without and with an applied magnetic field, as well as at finite temperature, are explored in the case of a 3D fcc lattice composed of ABC stacked kagome planes along cubic directions with classical spin vectors interacting only through the dipole coupling energy Eq. (1). The structure is inspired by magnetic compounds, such as IrMn,tomeno (); kren () which adopt a AuCu crystal structure, as shown in Fig. 1. Unlike the regular fcc lattice, the impact of the inherent AF frustration plays a more important role in the dipolar spin states of the fcc kagome lattice.

In many cases, a consequence of geometrical frustration is a degenercy of spin configurations in the ground state, normally when only nearest-neighbor (NN) antiferromagentic exchange interactions are included. This is particularly relevant for the 2D kagome lattice composed of corner-sharing triangles where the degeneracy is macroscopic, involving all the structures with the basic 120 spin configuration around each triangle.mekata2003 (); harris1992 (); zhito2008 (); schnabel2012 (); chern2013 (); taillef2014 () In the case of only dipole coupling between spins on the 2D kagome lattice, this degeneracy is reduced to six-fold due to the direct spin-lattice coupling.maksy2015 (); holden2015 (); maksy2017 () In general, the effect of thermal fluctuations or an applied magnetic field on frustrated spin systems is to reduce the degeneracy. For the kagome lattice with NN antiferromagnetic exchange, thermal fluctuations select co-planar q=0 spin structures at very low temperatures. With only dipole coupling, the six-fold (in-plane) spin degeneracy is largely unaffected by temperature except on cooling below T 0.43 (in units of the dipole strength ) where the system tends to lock-in to one (or more) of the six spin configurations at T 0.2. At all temperatures below T the 2D kagome dipole system displays a net magnetization.

Spatial dimensionality of the lattice also plays an important role on both the ground state degeneracy and impact of thermal fluctuations. Previous studies of the 3D fcc kagome lattice with NN antiferromagnetic exchange J (involving 8 neighboring sites), for both XY and Heisenberg models, showed that it also exhibits 120 spin configurations associated with q=0 order as in the 2D case (but with reduced degeneracy), and a phase transition to long-range order at T = 0.760J and T = 0.476J, respectively.hemmati (); leblanc1 (); leblanc2 () In the case of the triangular lattice with only dipole interactions, the ground state is ferromagnetic in 2D as well as for the 3D hexagonal structure.mckeehan (); tomita (); johnston () The impact of higher dimensionality on ground state dipolar spin structures in going from the 2D square lattice to 3D cubic structures is more complicated. The dipolar 2D square lattice exhibits AF orderdebell (); johnston () whereas the regular fcc lattice shows ferromagentic order.bouchard (); johnston () In contrast, and of particular relevance to the present work, the simple cubic dipolar lattice is characterized by a four sublattice spin configuration with continuous degeneracy involving two angles,belobrov () also seen in cubic clusters.schonke (); kure () Monte Carlo simulations on the cubic system suggest a discontinuous phase transition at T 0.56.romano ()

In this work we examine ground-state spin structures of the fcc kagome dipolar lattice using an Effective Field Methodwalker () (EFM) with lattice sums performed using Ewald techniques.ewald () The results reveal degenerate ground-state configurations characterized by six sublattice spins, composed of three spins per adjacent (111) kagome plane, with antiferromagnetic alignment between planes, yielding a zero net magnetic moment (in contrast with the 2D case). The six spin vectors are characterized by just two angles. This continuous degeneracy is shown to be removed by a order-by-disorder process using a Heat Bath algorithm. Similar degeneracy reduction is also shown to be achieved with magnetic field cycling. Finally, Metropolis Monte Carlo (MC) simulations are used to demonstrate a phase transition to long-range magnetic order at T 0.38 with a lock-in transition at T 0.3.

The remainder of this work is organized as follows. We describe the model used in Sec. II as well as the characterization of the ground state. The selection of particular ground states by thermal fluctuations at very low temeprature is described in Sec. III. The impact of an applied magnetic field is described in IV. MC simulation results on various thermodynamic quatities are presented in Sec. V, with a summary and conclusions given in Sec. VI.

## Ii Ground State Characterization

EFM simulations were performed on three-dimensional lattice consisting of ABC stacked 2D kagome planes (along cubic axes), with each plane occupied by () unit spin vectorshemmati (); leblanc1 (); leblanc2 () at lattice sizes 6, 12, and 18 using periodic boundary conditions. Several thousand runs were executed with random initial spin configurations to determine the lowest energy state, where each run involved several hundred sweeps through the lattice. Analysis of the results reveals a many-fold degenerate ground state allowing configurations with a six sublattice spin structure, having three sublattice spins per kagome layer, alternating in sign along the axes (i.e., adjacent layers are AF aligned). The basis spin vectors themselves can also be viewed as occupying the sites of three simple cubic lattices imbedded in the fcc kagome structure, as shown in Fig. 2.

By inspection of numerous resulting ground-state spin vectors, it was determined that every configuration is found to be characterized by the following set of equations:

(2) |

These equations act as elementary building blocks for the components of the six ground state spin vectors. There are only two parameters that characterize the resulting ground states, and , which are the polar and azimuthal angles of spin “1” (), where the polar axis lies along the positive z axis [001] of a Cartesian reference frame. Only six spins are required to fully characterize the ground states. The six ground state spin vectors themselves may be constructed as follows:

(3) |

Here, , , and are the three spins on a given triangle in one (111) kagome plane and , , and are spins on a triangle above or below on an adjacent (111) kagome plane. In all figures, appears in red, appears in blue, and appears in green.

We note that the relations between components of , , and can be obtained by directly calculating the system’s dipolar energy with a finite radial cutoff as a function of the components of the three independent spins, and assuming the six-sublattice structure given in Eq. 3. Consider in Fig. 2 as a central spin. It resides on a vertex in the red simple cubic lattice. The dipolar energy is zero for this sublattice (this is straight forward to show for a given spin and its six nearest neighbours). The central spin also resides in the middle of a face of the blue simple cubic lattice of alternating and spins, and similarly for the green simple cubic lattice of alternating and spins. On the level of nearest neighbours, is in the middle of two orthogonal squares, one blue and one green. Because neighbouring spins alternate signs, the AF-like contribution to the energy is zero (the first term in the sum in Eq. 8). This is in contrast to the 2D kagome dipolar system where the AF-like contribution from nearest neighbours is not zero, and hence the sublattice structure of the ground state emerges only after considering next-nearest neighbours.holden2015 () Thus, only the lattice-coupling term of Eq. 8 contributes to the ground state energy in the 3D kagome system. In terms of independent spin components of , , and the per particle dipolar energy of the six-spin system is,

where the nearest neighbour distance is unity and is a factor that starts at a value of 2 when only considering nearest neighbours and decreases with increasing the radial cut-off distance (to 1.7030969 when the cutoff is 22). A minimization of this energy with respect to the spin components subject to normalization of the spin vectors (through the method of Lagrange multipliers, for example) yields the expressions for , , and in Eqs. 2-3, and a minimum value of the quantity in brackets of . This ground state energy of compares well with our Ewald result of -2.55458.

Any pair chosen from the plane in Fig. 3 will give rise to a valid ground state of the same energy with the exception of those pairs of that lie within the bound region of the graph. Within the bound region of the graph, in Eq. 2 has an imaginary part and thus spin states characterized by those angle pairs () are not allowed and do not appear in the EFM results. At each node of each bounded area, in Eq. 2 becomes undefined, but the limit as approaches for , , is well defined.

Note that the pattern in Fig. 3 repeats in with a period of due to crystal symmetry considerations. Nodes occur at where is an integer. The plane is simplified further due to reflections about the axes in the middle of each bounded region. The three spins at node states lie in a plane. For example, for and ,

(5) |

are in the plane.

Analysis of the relations Eq. 2 reveals that the region of allowed pairs satisfy,

(6) |

Thus, the entirety of the ground state is characterized by Eqs. 2-3 and the parameters such that they satisfy Eq. 6.

By generating the spin vectors described by Eqs. 2-3, spin configurations from what we will call planar and non-planar states can be obtained. An example of a non-planar state is shown in Fig. 4.

To gain insight into what choices of pairs give rise to particular states, consider the contour plot in Fig. 5. It illustrates the volume of a parallelepiped formed by the , , and spin vectors corresponding to a particular chosen from Fig. 3. It is from this illustration that one can obtain an understanding of what choices of result in planar states or non-planar states. Planar states are those having zero volume of the parallelepiped which appear as the darkest regions in Fig. 5 as nearly circular rings with a radius of roughly . Note that the nodes indicated on Fig. 3 and Fig. 5 represent planar states. In terms of the previous polar coordinates planar states can be shown to satisfy trigonometric relations, and, e.g., for the locus of planar states is given by,

(7) |

outlining a section of the darkest rings in Fig. 5.

## Iii Order-by-Disorder at Low Temperature

In this section results are presented on the evolution of ground states under the influence of thermal fluctuations. A low temperature Monte Carlo Heat Bath method,heatbath () that uses the local field determined through the EFM algorithm, was employed for this purpose. Beginning with the particular ground state shown in Fig. 8, the temperature was increased from T=0 to T=0.008 with increments of T=0.00001. After each increase in temperature, the spin configurations were run through the T=0 EFM to obtain the nearest corresponding ground state. Each of the states obtained following each zero temperature EFM run were recorded and are displayed in Fig. 9.

The initial state is clearly a non-planar state. As temperature increases, the spin configurations gradually transition from the non-planar region to the planar region in Fig. 9, illustrated also by the snapshot of the state at shown in Fig. 10. This order-by-disorder process appears to select the planar states, near the “node” of the bounded region, as the minimum free energy configuration from all possible states after exposing the spin lattice to thermal fluctuations.

Increasing the temperature results in the spin configuration moving towards a planar node state regardless of its initial state. Temperature has the effect of lifting the degeneracy of the system from one that includes all possible (valid) points in the ground-state manifold to those that are planar (or near planar).

## Iv Applied Magnetic Field

In this section the effect of an applied magnetic field at zero temperature on the spin configurations is presented. Our first set of results stem from seven EFM simulations performed on lattices of size , all starting with the same non-planar ground state, with the Zeeman term

(8) |

added to the dipole energy, Eq. (1). The magnetic field was increased in steps of up to and then larger steps of until the field reached and the spins became nearly saturated. Subsequently, the magnetic field for each lattice was decreased in steps of until magnitude was zero (field cycling). Following every change in magnetic field, the spins were subjected to the EFM method using 3000 iterations. The magnetic field directions used for the simulations were along directions of high crystal symmetry; [001], [010], [100], [011], [101], [110], and [111].

Illustrative results for the field along principal cube axes are shown in Fig. 11. The magnetization as a function of magnetic field indicates several distinct phases. For the example state in Fig. 11 with [001], the magnetization lingered close to zero at low magnetic field.

At some critical value (H 0.012), a sudden change in magnetization occurred, signalling the formation of a planar state. A similar jump to a planar state occurred for all the field directions. A linear change in magnetization was observed after this critical field magnitude and at the magnetization reaches near saturation. (Note there is a slight positive slope in the plateau regime.) Snapshots of the six sublattices spin vectors orientations with increasing field are shown in Fig. 12. For the field in the [001] direction, there is a pronounced hysteresis near the transition to the saturation plateau (), which is largely absent for H along [110] or [111]. These MH loops can also be dependent upon which of the crystallographically equivalent high symmetry directions are chosen. For example, three different MH curves can occur for a given starting ground state with H increasing along [100], [010], and [001]. However, all of the crystallographically equivalent directions show identical MH curves upon decreasing the field from saturation. These behaviours indicate that the energetically equivalent ground states experience a complex energy landscape as a function of applied field.

These simulations reveal how a sufficiently large external magnetic field causes a spin configuration to transition from a non-planar state to a planar state. When decreasing the magnetic field magnitude, the spin vectors orient themselves away from the magnetic field direction and form a planar state at zero magnetic field. The general response of the spins to the magnetic field were consistent with the exception of one point; spin configurations that are already planar at zero field do not experience an abrupt change into a planar state, but rather a linear change similar to non-planar states after transitioning to a planar state.

Our second set of results show that all possible ground states can be transformed into planar states through field cycling. Twenty-six different initial ground states spanning and were prepared, as indicated in Fig. 13. Each of these states was field cycled with a field along the [111] direction.

Twenty-four of the initial non-planar states resulted in the same final state upon field cycling, while only two resulted in a different final state, a node state, as indicated in Fig. 13. Both of these are planar states. Thus, it is clear that by cycling an external magnetic field, the degeneracy of the ground state is reduced to include only states that are planar.

## V Monte Carlo Simulations

### v.1 Energy and Specific Heat

As in our previous MC simulations on the 2D kagome lattice,holden2015 () three main types of temperature runs were performed: isolated temperature, cooling and heating. All simulations were done at zero magnetic field on lattices with 6, 12, and 18. Isolated temperature data are collected by a simulation at every initialized with a random configuration. These initial configurations have no dependence on the configurations of the neighbouring temperatures. Furthermore, each temperature may be run simultaneously on separate processors, greatly decreasing the real time length of simulations. Heating and cooling runs, however, are completed at consecutive using the final spin configuration of the previous as the initial configuration. In the case of isolated temeprature runs, at each temperature, at least 910 MC steps (MCS) were used with the initial MCS discarded for thermalization giving about MCS for averaging. A single MCS involves, on average, an attempt to flip each spin in the system. For cooling and heating runs, these numbers were reduced by a factor of ten. Preliminary MC results on the 3D kagome system using a small lattice ( = 8) with a limited number of MCS were previously reportedholden2015 (); shanethesis () and indicated anomalous behaviour at a temperature near 0.35.

Figure 14 shows results for the specific heat, given by,

(9) |

for = 6, 12, and 18, where the sum is over the MC generated configurations.

In each set there is a peak in the range which becomes larger, narrower, and appears at higher , with increasing system size. The largest and most well defined peak occurs in the case near . At low temperatures, the results suggest , consistent with earlier simulation results.holden2015 ()

Figure 15 shows the energy at = 12 performed at isolated temepratures as well as cooling and heating runs. There appears to be little difference in the simulation types. This observation is important to the results shown below for the sublattice magnetization, which do exhibit a strong dependence on simulation type.

### v.2 Magnetization and Susceptibility

The total ferromagnetic magnetization is often used in classical spin systems as an order parameter, especially for those associated with dipolar coupling. In the present case, the ground state exhibits a zero net moment due to the alternating spin vectors along [111] axes. For the 2D kagome system, was shown to increase with decreasing temperature with a maximum value of for a pure, single domain, ground state.holden2015 () Figure 16 shows MC results for decreasing with decreasing temperature, tending towards zero at low . This behaviour holds true for all lattice sizes, and diminishes rapidly with increasing at all . Figure 17 shows the magnetization as a function of the inverse square root of number of spins at selected temperatures, showing a positive, linear scaling that grows with temperature. The magnetization approaches zero at all temperatures as lattice size increases while the slope of the fit decreases with temperature, as expected for an antiferromagnetic system. The ground state magnetization is expected to have a slope of zero (i.e. perfect antiferromagnetic order for all lattice sizes), which corresponds with this behaviour.

For the present system, as in the 2D case,holden2015 () the sublattice magnetization serves as an order parameter associated with the phase transition to long range order. The value of this quantity in each individual sublattice of size is calculated as

(10) |

where represents the subset of all spins belonging to any given sublattice, and is selected such that adjacent spins, , carry an opposing sign (i.e. for lattice position ). The spins of each cubic sublattice are expected to approach perfect ordering () at zero temperature. By taking a thermal average of each sublattice magnetization the total sublattice magnetization is calculated as

(11) |

Figure 18 shows the results for the three different simulation types. It can be seen that in cooling simulations, some domains of “frozen-in” spins of a different orientation occur, as in the 2D case.holden2015 () These domains reduce the total sublattice magnetization. In both cases of cooling and isolated temperature runs, for , for some temperature , the system becomes non-ergodic and is no longer able to sample states effectively.

The frozen-in states occur when it is no longer energetically favorable for single spins, even on the border of the domains, to change orientation. Thermal averaging results for the three individual sublattice magnetization are displayed in Fig. 19 for cooling and heating simulations. Two of the cooled system sublattices do not uniformly approach unity. A depiction of frozen-in states in the sublattice is presented in the cross section of Fig. 20 with a portion of differently-oriented spins highlighted. These frozen-in states occur in domains which are symmetrically equivalent on the lattice, spanning the entire lattice in one dimension (i.e. individual plane), and appear any number of times. That is, these domains are related to each other by the crystal symmetry. When compared to a single-domain state, the energy differences at are on the order of . Such a small energy difference leads to domain occurrence being commonplace in both isolated, cooling, and, on rare occasion, heating simulations.

The ferromagnetic susceptibility, , of the lattice is defined as the variance in ferromagnetic magnetization and is given by

(12) |

Figure 21 shows from isolated temperature simulations with varying lattice sizes. Figure 22 shows for isolated temperature, and cooling and heating simulations at only. The majority of simulations produce a narrow, noisy peak at followed by a weak, broad peak around . The first, narrow peak occurs at the temperature at which the lattices become locked into their configurations, as seen in the total sublattice magnetization, is related the the non-ergodicity of the system below this temperature. These per site values of the susceptibility do not appear to scale with respect to the system size, which is consistent with a lack of a ferromagnetic transition in the system.

The obvious exceptions to this generalization are the isolated temperature runs where finite-size effects are strongest. This may facilitate the sampling of more states at lower , i.e., increase ergodicity. The heating simulation results in Fig. 22 at also show signifcant noise at lower , as with the heating results in Fig. 18, where spins remain locked into a ground state until the transition begins at . This behaviour is witnessed in all simulations, including the cooling simulations in which domains are not realized. Heating simulation results, in which the spin lattice quickly develops domains, align closely with the cooling simulation susceptibility suggesting the domains serve to lock spins into an orientation and reduce the susceptibility of the system.

Similarly, the antiferromagnetic susceptibility, , can be defined in terms of the variance in the total sublattice magnetization defined in Eq. 11. This is given as

(13) |

Figure 23 shows results for isolated temperature simulations at = 6, 12, and 18. A single, sharp peak is apparent only for , corresponding with the inflection point in the value of , and peak in the specific heat, producing the same estimated critical temperature . Also note that for = 18, the results show reduced noise relative to the smaller lattice simulations below a temperature around 0.3, also consistent with the above discussion of lock-in states.

## Vi Summary and Conclusions

This work explores an example of the interplay between the geometrical frustration of a 3D kagome lattice involving the antiferromagnetic and anisotropic character of the dipole interaction, as well as the long-range nature of dipole coupling. The ABC stacking of kagome planes gives rise to a truly 3D example of a kagome lattice, with eight NN sites, four in-plane and two connecting each adjacent plane. In contrast with the 2D case where discrete degeneracy of the three-sublattice ground state spin structures reflects the six-fold hexagonal anisotropy,maksy2015 (); holden2015 (); maksy2017 () and the regular fcc lattice showing ferromagnetic order, the 3D fcc dipolar kagome lattice is shown here to exhibit a continuous degeneracy involving two parameters (the polar angles and ), similar to the simple cubic dipolar lattice.belobrov () Certain regions in the - plane are excluded from the domain of allowed ground states characterized by six sublattice spins, involving three spins around a triangle in the (111) plane, and three spins around the adjacent plane triangle pointing in opposite directions. The total magnetization is thus zero, unlike the 2D case.

This continuous degeneracy is removed by thermal fluctuations at low temperatures in an order-by-disorder process, and through magnetic field cycling, with both processes yielding planar states for which all six sublattice spins lie in a single plane. Thermal fluctuations result in the system approaching any one of a set of discrete “node states” (, and odd multiples thereof). Field cycling appears to pick out a set of discrete set of states along the locus of planar states, including node states.

The evolution of degenerate ground state spin configuration with an applied magentic field is dependent on the field direction, with M vs H curves showing features for H along [100] and [110] directions, but only a continuous rotation of spins in the case of H along [111] up to saturation. The detailed evolution of a state with increasing field is dependent on which of the ground states is used, but always returning to a planar state as the field decreases to zero.

MC simulations show evidence for long-range order at and confirm the absence of a finite magnetization. In addition, the MC results suggest that for , the system can lock-in to a mix of the degenerate ground state configurations and exhibit non-ergodic theromodynamics.

The realization of such a purely dipolar magnetic system appears limited to artificial nanostructuresartificial () but the relevance of the present work my be related to a number of atomic kagome-based structures which also exhibit exchange interactions.dun (); hayashida () Among AB compounds where magnetic ions form the fcc kagome structure, IrMn with strong AF exchange is the most well known for its technologically important role in spin valves. As mentioned in Ref. [holden2015, ], dipole coupling may be expected to be important at the interface with an adjacent thin-film ferromagnet when IrMn is used as an exchange bias material. The present work may also be relevant to the magnetism in a sister compound, GeMn, which exhibits exchange-driven ferromagnetism where dipole coupling can also be important.gemn ()

## Vii Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada, and the Compute Canada facilities of the Atlantic Computational Excellence network (ACEnet) and the Western Canada Research Grid (WestGrid).

## References

- (1) Frustrated Spin Systems, edited by H.T. Diep (World Scientific, Singapore, 2004); Introduction to Frustrated Magnetism, edited by C. Lacroix, P. Mendels, and F. Mila, Springer Series in Solid-State Sciences, Vol. 164 (Springer, Berlin, 2011).
- (2) Krén, G. Kádár, L. Pál, J. Sólyom, and P. Szabó, Phys. Lett. 20, 331 (1966); E. Krén, G. Kádár, L. Pál, J. Sólyom, P. Szabó, and T. Tarnóczi, Phys. Rev. 171, 574 (1968).
- (3) I. Tomeno, H. N. Fuke, H. Iwasaki, M. Sahashi, and Y. Tsunoda, J. Appl. Phys. 86, 3853 (1999).
- (4) M. Mekata, Physics Today 56, 2, 12 (2003).
- (5) A. B. Harris, C. Kallin, and A. J. Berlinsky, Phys. Rev. B. 45, 2899 (1992.)
- (6) M. E. Zhitomirsky, Phys. Rev. B 78, 094423 (2008).
- (7) S. Schnabel, and D. P. Landau Phys. Rev. B 86, 014413 (2012).
- (8) G-W. Chern, and R. Moessner, Phys. Rev. Lett. 110, 077201 (2013).
- (9) M. Taillefumier, J. Robert, C. L. Henley, R. Moessner, and B. Canals Phys. Rev. B 90, 064419 (2014).
- (10) M. Maksymenko, V. R. Chandra, R. Moessner, Phys. Rev. B 91, 184407 (2015).
- (11) M. S. Holden, M. L. Plumer, I. Saika-Voivod, and B. W. Southern, Phys. Rev. B 91, 224425 (2015).
- (12) M. Maksymenko, R. Moessner, K. Shtengel, Phys. Rev. B 96, 134411(2017).
- (13) V. Hemmati, M. L. Plumer, J. P. Whitehead, B. W. Southern Phys. Rev. B 86, 104419 (2012)
- (14) M. D. LeBlanc, M. L. Plumer, and J.P. Whitehead, Phys. Rev. B 88, 094406 (2013).
- (15) M.D. LeBlanc, B.W. Southern, M.L. Plumer and J.P. Whitehead, Phys. Rev. B 90, 14440 (2014).
- (16) L. W. McKeehan, Phys. Rev. 43, 1025 (1933).
- (17) Y. Tomita, J. Phys. Soc. Jpn. 78, 114004 (2009).
- (18) D. C. Johnston, Phys. Rev. B 93, 14421 (2016).
- (19) K. De’Bell, A. B. MacIsaac, I. N. Booth, and J. P. Whitehead, Phys. Rev. B 55, 15108 (1997).
- (20) J. P. Bouchard and P. G. Zérah, Phys. Rev. B 47, 9095 (1993).
- (21) P. I. Belöbrov, R. S. Gekht, and V. A. Ignatchenko, Sov. Phys. JETP 57, 636 (1983).
- (22) J. Schonke, T. M. Schneider, and I. Rehberg, Phys. Rev. B 91, 20410(R) (2015).
- (23) M. Kure, M. Beleggia, and C. Frandsen, J. Appl. Phys. 122, 133902 (2017).
- (24) S. Romano, Il Nuovo Cimento 7, 717 (1986).
- (25) S. Holden, Frustration of Dipolar Spins on the Kagome Lattice. Honours Thesis, Department of Physics and Physical Oceanography, Memorial University of Newfoundland (2015).
- (26) L. R. Walker and R. E. Walstedt, Phys. Rev. B 22, 3816 (1980).
- (27) M. Born and H. Kun, Dynamical Theory of Crystal Lattices (Oxford University Press, Oxford, 1985); M. Enjalran1 and M. J. P. Gingras, Phys. Rev. B 70, 174426 (2004).
- (28) Y. Miyatake, M. Yamamoto, J. J. Kim, M. Toyonaga, and O. Nagai, J. Phys. C 19, 2539 (1986).
- (29) I. Gilbert, C. Nisoli, and P. Schiffer, Physics Today 69, 54 (2016); L. Sun et al., Phys. Rev. B 96, 144409 (2017); P. Gypens, J. Leliaert, and B. Van Waeyenberge, Phys. Rev. Appl. 9, 034004 (2018).
- (30) Z. L. Dun et al., Phys. Rev. Lett. 116, 157201 (2016).
- (31) S. Hayashida et al., Phys. Rev. B 97, 054411 (2018).
- (32) H. Takizawa, T. Yamashita, K. Uheda, and T. Endo, J. Phys.: Condes. Matter 14, 11147 (2002); Y-F. Kong, S-Y. Wang, M. Xu, G. yin, and L-Y. Chen, J. Korean Phys. Soc. 49, 2188 (2006); E. Arras, D. Caliste, T. Deutsch, F. Lancon, and P. Pochet, Phys. Rev. B 83, 174103 (2011).