Order in the ground state of a simple cubic dipole lattice in an external field
aa r X i v : . [ c ond - m a t . o t h e r] S e p Order in the ground state of a simple cubic dipole lattice in anexternal field ∗ S. Ashhab, M. Carignano, and M. E. Madjet
Qatar Environment and Energy Research Institute,Hamad Bin Khalifa University, Qatar Foundation, Doha, Qatar (Dated: September 4, 2019)
Abstract
Motivated by the presence of a lattice of rotating molecular dipoles in the high temperaturephase of methylammonium lead iodide, we investigate the ground state of a simple cubic latticeof dipoles interacting with each other via the dipole-dipole interaction and with an external fieldvia the Zeeman interaction. In the absence of an external field, the ground state is infinitelydegenerate, and all the configurations in the ground state manifold are periodic along the threelattice axes with period 2. We numerically determine the ground state of a 1000-dipole latticeinteracting with an external field, and we analyze the polarization, dipole orientation statistics andcorrelations in this state. These calculations show that for some special directions of the externalfield the two-site periodicity in the dipole configurations is preserved, while in the general case thisperiodicity is lost and complex dipole configurations form in the presence of the external field. ∗ This is the pre-peer reviewed version of an article that will appear in the International Journal of QuantumChemistry. The manuscript was initially submitted under the title “Ground state structure of a simplecubic dipole lattice in an external field”. This article may be used for non-commercial purposes inaccordance with Wiley Terms and Conditions for Use of Self-Archived Versions. . INTRODUCTION Spin models are ubiquitous in theoretical physics in the study of a variety of physicalsystems, and they can be used to understand a wide range of phenomena such as orderedstates of matter, domain boundaries and phase transitions [1]. Recently, they have alsoserved as natural models for the study of entanglement and related properties in the studyof condensed matter systems from the point of view of quantum information [2, 3]. In thefield of light harvesting devices, the recently emerging hybrid organic-inorganic lead-halideperovskites have opened a novel application for spin models because the organic componentoften has a net dipole moment. The most representative perovskite material in this areais CH NH PbI , which in its high temperature phase contains a simple cubic lattice ofdipoles from the polar molecule CH NH that occupies the A site in the perovskite crystalstructure [4–6]. Spin models have been used to explore the possible origin of hysteresisobserved in solar cells [7, 8], to understand the dynamics of the cation rotations [9] and alsoto explain enhancement of carriers conductivity at the interface between domains [10]. Morerecently, spin-models have been used to explore the long-range order in perovskites [11] andthe structure of the interface between domains with different ordering patterns across theinterface [6].Spin models apply more generally to other systems in nature, some examples being theordering of the electric dipoles of water molecules in ice [12], the dipolar ordering of smallmolecules trapped in regular cage structures [13–15]. Moreover, the application of spin mod-els to dielectric and magnetic properties of materials has a long history [17–23]. Artificiallyfabricated magnetic superstructures can also be described by similar models [16]. One par-ticularly interesting result, found by Luttinger and Tisza [19], is that in the absence of anexternal field a simple cubic lattice of dipoles has an infinitely degenerate ground state withno net polarization. There have been several studies on the polarization properties of theground state of a dipole lattice under the influence of an external field [19–22, 24–28]. Thesestudies have generally considered cases where the external field points in some special direc-tions that allow the derivation of analytical results. In the general case of an arbitrary fielddirection, one cannot perform fully analytical derivations to our knowledge, and numericalmethods are necessary. Here we perform such numerical calculations on a dipole lattice in anexternal field, analyzing several different field directions. By evaluating dipole orientation2tatistics and correlation functions, we show that for a field direction outside any of themain lattice planes (i.e. xy, yz and xz) of a simple cubic lattice the simple structure of theground state is lost, in addition to the lifting of the degeneracy that exists at zero externalfield.We emphasize from the outset that an accurate model of CH NH PbI and similar ma-terials would include the inorganic lattice and its interactions with the organic moleculardipoles. However, our purpose in this work is to gain an understanding and insight intothe physical mechanisms related to the dipole-dipole interactions in the system. The resultsthat we obtain in this work can then be used to contribute to our understanding of the rolethat dipolar interactions could play in such a complex material. II. GROUND STATE IN THE ABSENCE OF AN EXTERNAL FIELD
We consider a simple cubic lattice of dipoles interacting via the dipole-dipole interaction.For definiteness, we shall use the language of electric dipoles, although the results are generaland apply e.g. to magnetic dipoles. The total dipole-dipole interaction energy of the latticeis given by U int = 18 πǫ X i,j ~p i · ~p j − ~p i · ˆ r ij ) ( ~p j · ˆ r ij ) | r ij | , (1)where ǫ is the permitivity, ~p i is the dipole vector of dipole i , ˆ r ij is the unit vector pointingfrom the location of dipole i to that of dipole j , and r ij is the distance between dipoles i and j . In order to transform the Hamiltonian into a more universal form, we define the newvariables ~p i = p ˆ p i r ij = r ˜ r ij (2) U = p πǫr , where p is the magnitude of the dipole moments (which is assumed to be fixed and equal forall the dipoles), and r is the lattice parameter. The vector ˆ p i is now the unit vector in thedirection of dipole i , and ˜ r ij is the renormalized distance between dipoles i and j in units ofthe lattice parameter. We emphasize that the problem under study is a classical one, andthe hat symbols should not be confused as indicating quantum mechanical operators. With3hese definitions the dipole-dipole total interaction energy is given by U int = U X i,j ˆ p i · ˆ p j − p i · ˆ r ij ) ( ˆ p j · ˆ r ij ) | ˜ r ij | . (3)We keep the factor of 2 in the coefficient on the right-hand side to indicate that this factoris used to avoid double-counting the interaction energy of a single dipole pair.The ground state of the above-described dipole lattice is infinitely degenerate [19]. If wechoose one dipole in the lattice and designate it as the origin of our system of coordinates,this dipole can point in any direction on the unit sphere with Cartesian components p x , p y and p z . The orientations of all other dipoles are then given by the following simple rule:the dipole at location ( x, y, z ) [each of these being an integer] has Cartesian components( − y + z p x , ( − x + z p y and ( − x + y p z . This dipole configuration has periodicity 2 along eachof the crystal axes (i.e. x, y and z). The periodicity in any of these three directions becomes1 in the special case where the dipoles are oriented along that direction. x yz FIG. 1: Schematic diagram of the LT ground state. In the absence of an external field, one has atwo dimensional manifold of degenerate ground states. The dipole at the origin can be chosen topoint in any direction, which then determines the directions of all the dipoles in the lattice. The LTground state has two-site translation symmetry in the x, y and z directions. The net polarizationis zero in the absence of an external field.
Figure 1 illustrates one example of the manifold of degenerate Luttinger-Tisza (LT)4round states, all of which possess the same energy ( − . N U for an interaction cutoffof 3 r , with N being the total number of dipoles) and no net polarization. III. POLARIZATION AND DIPOLE CONFIGURATIONS IN THE PRESENCEOF AN EXTERNAL FIELD
We now calculate and analyze the ground state of the dipole lattice under the influenceof an externally applied field.The energy arising from the interaction between the dipole lattice and the external fieldis given by U ext = − ~E · X i ˆ p i , (4)where the vector ~E is the external field, including any factors that convert the field into theappropriate energy units. The total energy is now given by U total = U int + U ext .Before analyzing the case of a finite external field, we mention that in the absence of theexternal field we can take the energy per dipole, i.e. − . U , and infer from it the valueof the internal field induced by all neighbouring dipoles at the location of any given dipole.Since all the dipoles in the lattice are equivalent to each other in any of the LT configurations,all of the dipoles will have the same energies and experience the same absolute value of thelocal field. Using the formula U int = − ~E i, int · ˆ p i × N/
2, where the unspecified index i indicatesthat we can choose any dipole in the lattice, we find that the internal field is 5 . U pointingin the same direction as the dipole (which is true for every dipole in the lattice). Note thatwe again use the factor of 2 here as in Eq. (3), because the interaction energy of a pair iscalculated as the energy of only one of the two dipoles in the field induced by the otherdipole.In our calculations, we generally use a 10 × ×
10 lattice of dipoles. In each calculation wesearch for the dipole configuration that minimizes the energy for a given field direction andstrength. We do so by initializing the dipole configuration in a state of randomly orienteddipoles and allowing the dipoles to relax to lower energy states: in each iteration of thecalculation we evaluate the total field at each dipole location and then rotate each dipoles inthe direction that maximizes the reduction in its energy. The rotation angle for each dipole isproportional to the torque felt by the dipole multiplied by an overall step size. This step sizeis reduced whenever the above procedure stops reducing the total energy, and the calculation5s stopped when the step size reaches 10 − , at which point we expect that we have reacheda good approximation of the ground state. We find that, with some exception explainedbelow, we generally obtain the same dipole configuration for all randomly generated initialstates. After obtaining a given optimized dipole configuration, we analyze its polarizationand statistical properties.
0 0.2 0.4 0.6 0.8 1 E ext / E P FIG. 2: The polarization P as a function of external field strength for different field directions:along one of the main crystal axes [i.e. x, y or z] (red squares), in one of the main planes andmaking an angle π/ θ, φ ) = ( π/ , π/
6) (magenta diamonds). The filled symbols areobtained using a 10 × ×
10 lattice, while the open symbols are obtained using a 2 × × The polarization, defined as P = | P i ˆ p i | /N , as a function of external field strength isshown in Fig. 2. When the external field lies in one of the three planes xy, xz and yz(to which we refer as the main planes below), the polarization has a linear dependence onthe field from E = 0 to E = E , where E is the value of the externally applied field atwhich the system reaches a state of full polarization along the direction of the applied field.6hen we take different field directions, i.e. directions where the external field has finitecomponents along x, y and z, we find that the polarization does not deviate much from thelinear dependence, as also shown in Fig. 2. We note here that although this result agreeswith intuitive expectations for dielectric and paramagnetic materials, to our knowledge ithas not been shown to be true for a LT dipole lattice before. For all field directions, thepolarization is parallel to the direction of the external field, except for small deviations atintermediate values of the external field.It is interesting that the value of the external field E at which the system becomes fullypolarized is given by 5 . U , which is exactly the value of the internal field induced byneighbouring dipoles at zero field, even though these two fields occur in different situations.In other words, the dipole-induced field is given by 5 . U only when the external fieldstrength is zero; the dipole-induced field then decreases and reaches the value zero as theexternal field is increased to the value 5 . U . Based on this observation, one might wonderwhether the total field at any dipole location remains constant between the two above limits.We have confirmed numerically that this is not the case. The local field decreases from5 . U to a smaller value as the external field strength is increased from zero, and thenabove a certain value of the external field the local field reverses its tendency and startsincreasing to come back to the value 5 . U when the external field reaches the value5 . U (meaning that the dipole-induced field vanishes at that point, as discussed in moredetail in the Appendix).In the following three subsections, we analyze the dipole configurations in some moredetail. A. Weak external field
Let us take the infinitely degenerate ground states in the case of zero field discussed inSec. II, i.e. the LT state, and consider the effect of an infinitesimally weak field.We start with the case where the externally applied field is parallel to one of the threelattice axes, and for definiteness we take this axis to be the z axis. In physical systemsthat have a degenerate ground state, a weak external field generally plays the role of aperturbation that lifts the degeneracy. We therefore ask: what dipole configuration doesthis external field favour? If we had a single dipole, which can point in any direction in7he absence of an external field, the dipole would align with the external field. Somewhatcounter-intuitively, in the case of the dipole lattice under study, the external field pointingin the z direction favours dipole configurations in which the dipoles point perpendicular tothe z axis [19, 24]. This result can be understood by considering that any given dipole feelsan internal field that is induced by the neighbouring dipoles, and this field is parallel to thedirection of the dipole under consideration. The weak external field acts as a small pertur-bation that slightly modifies the field felt by each dipole. In order for the infinitesimallyweak external field to make the largest possible change in the total energy, it is favourableto have the external field perpendicular to the internal field. In this case, each dipole rotatesby a small amount towards the direction of the external field, leading to an energy reductionrelative to the unperturbed ground state. This situation can now be contrasted with thatof a configuration in which the dipole orientations are all parallel to the z axis and we add aweak external field along the z axis. In this case, since the internal field points along the zaxis, the external field slightly modifies the value of the field at each dipole location, but itdoes not change the direction of the field. As a result, each dipole will still point along thetotal field at its location, making each dipole feel that it is in its lowest-energy orientation,and no energy reduction can be obtained by rotating any dipole. Note that although theapplication of a weak field pointing along the z axis breaks the two-dimensional continuoussymmetry of LT configurations, the rotation symmetry about the z axis is preserved, andone therefore still has a one-dimensional manifold of ground states, which is confirmed byour numerical simulations.A somewhat similar result is obtained when the external field lies in one of the threemain planes. In the limit of weak external field, it is energetically favourable for the dipolesto point perpendicular to the external field. However, no infinite degeneracy survives inthis case. For example, if the external field lies in the xy plane (but not along x or y), thebroken-symmetry configuration will have the dipoles pointing in the z direction. Choosinga configuration where the dipoles have a finite xy component and following the rules ofconstructing the LT dipole configuration would result in a situation where some dipoles arenot perpendicular to the external field, which is not as energetically favourable as having allthe dipoles perpendicular to the direction of the external field.The situation becomes more complicated when the external field has finite componentsin all three directions. One can see that this case is trickier than the above two by noting8hat if we take any given dipole and set it to any direction perpendicular to the externalfield, the LT rules determining the orientations of the other dipoles in the lattice dictatethat some dipoles will point in the same general direction as the external field (even ifnot exactly parallel to it). In our simulations, we find that the dipoles point along one ofthe three crystal axes, specifically the one that makes the largest angle with the externalfield. For example if we take an external field direction defined by the spherical coordinates( θ, φ ) = ( π/ , π/ B. Strong external field
An infinitely strong external field would obviously dominate over the dipole-dipole inter-action and polarize all the dipoles parallel to the external field. We therefore take as thestrong-field state the configuration where the dipoles are ordered in a fully ferroelectric typeof state, with the orientation vectors ( p x , p y , p z ) being uniform for all the dipoles. In thisstate, as shown in the Appendix, the internal field induced by all neighbouring dipoles atthe location of a central dipole vanishes, independently of the orientation direction of thedipoles. This result means that the internal field does not favour any rotation of any dipoleaway from the external field direction and confirms that for sufficiently strong external fieldsall the dipoles will align with the external field.Note that the vanishing of the internal field in the ferroelectric state means that theferroelectric configuration with axis parallel to the external field must be a dynamicallystable state for any value of the field. The reason is that in this case every dipole is parallelto the total field at the dipole’s location, which is a dynamically stable state. However, forfields smaller than a certain value (which turns out to be 5 . U as discussed above) otherconfigurations become energetically more favourable than the ferroelectric configuration.This situation can also lead to a hysteresis effect, as discussed e.g. in Ref. [25]. C. Intermediate field values
In the case where the external field lies in one of the main planes, the dipoles start offpointing perpendicular to the external field in the limit of very weak field. As the field9trength increases, all the dipoles rotate by the same angle to become partially aligned withthe external field. The angle is determined by the linear increase in polarization as thefield is increased from zero to 5 . U . For field values above 5 . U , the dipoles are allaligned with the external field. For all values of the external field strength, the two-sitetranslation symmetry is clearly preserved in all three directions [19]. In other words, itsuffices to perform calculations on a 2 × × × ×
10, give different results from simulations on a 2 × × × ×
10 configurations therefore do not have the two-site translation symmetry of theLT configurations. Otherwise the results calculated using the 10 × ×
10 lattice and thosecalculated using a 2 × × α i =ˆ p i · ~E ext / | E ext | . When the field lies in one of the main planes, all the dipoles make thesame angle with the field, and the histogram therefore contains only one peak that movesfrom α = 0 to α = 1 as the field strength is increased from E ext = 0 to E ext = E . Whenthe field lies outside the main planes, the histogram starts off with two peaks at E ext = 0,as would be expected for a LT configuration in which the dipoles point along one of themain axes. When the field strength is increased and the two-site translation symmetry isbroken, the histogram in general contains a large number of α values, which reflects the factthat one no longer has at most eight dipole orientations that are repeated throughout thelattice. The fact that the histogram turns into a broad distribution at intermediate values ofthe external field is also interesting because it demonstrates that the dipole configurationsbecome complex, with no simple obvious pattern, even in the ground state. It is expectedthat such complex patterns also arise for larger lattices, such that the ground state wouldlook disordered at intermediate field values.Figure 4 shows the dipole-dipole correlation function C ijk = h ˆ p abc · ˆ p a + i,b + j,c + k i , (5)where the average is taken over all the dipoles in the lattice. Interestingly, for the cases shown10 α P r ob ( α ) (a) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (b) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (c) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (d) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (e) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (f) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (g) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (h) 00.20.40.60.81 -1 -0.5 0 0.5 1 α P r ob ( α ) (i) FIG. 3: Histograms showing the distribution of values for the individual dipole components α alongthe direction of the external field. From left to right the field direction is varied: ( θ, φ ) = ( π/ , π/ θ, φ ) = ( π/ , π/
6) (right column). Fromtop to bottom the field strength is varied: E ext = 0 . E (top row), 0 . E (middle row) and 0 . E bottom row. in the figure, the correlation function remains at its maximum value C ijk = 1 (meaningthat the two-site translation symmetry is preserved) up to finite fraction of E (specifically E ext = 0 . E and E ext = 0 . E for the two cases plotted in Fig. 4). This result indicatesthat a LT-like state persists for a finite range of external field values. At somewhat highervalues of the field, the correlation function drops below one. The initial drop is ratherquick, such that the dipole configuration changes suddenly from a perfectly ordered stateto a highly disordered state. As the external field becomes stronger, the dipoles becomeincreasingly aligned with the external field until they are fully polarized along the directionof the external field, which occurs when the external field E ext reaches E . At this point C ijk E ext / E C ij k (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 E ext / E C ij k (b) FIG. 4: The correlation function C ijk for different values of ( i, j, k ): (2,0,0) [filled red squares],(0,2,0) [filled red circles], (0,0,2) [filled red triangles], (4,0,0) [open green squares], (0,4,0) [opengreen circles] and (0,0,4) [open green triangles]. Two field directions are used in the calculations:the direction (1,1,1) [Panel a] and ( θ, φ ) = ( π/ , π/
6) [Panel b]. In the case of the direction (1,1,1),there is a rotation symmetry in which the three main axes are equivalent to each other. Thissymmetry is clearly not present in the dipole configuration, since for example the red symbols donot all coincide with each other. returns to its maximum value C ijk = 1, because all the dipoles are aligned with each otherby virtue of their alignment with the external field. The dipole-dipole correlations in thiscase have nothing to do with internal interactions, since a strong field would lead to the samecorrelations regardless of any details related to the internal interactions. Note that thesesomewhat artificial correlations can be suppressed in the theoretical analysis by defining thecorrelation function as h ˆ p abc · ˆ p a + i,b + j,c + k i − P . However, since our purpose of analyzing thecorrelation function is to investigate the periodicity in the dipole configuration, we use thedefinition given in Eq. (5). 12 V. CONCLUSION
We have analyzed the behaviour of the ground state of a simple cubic dipole latticewith dipole-dipole internal interactions and interacting with an external field. This systemresembles the high temperature phase of methylammonium lead iodide, although it neglectsthe interaction of the dipolar cations with the inorganic lattice. In the special case wherethe external field lies in one of the main lattice planes, our results agree with the knownresult that the dipole configurations remain rather simple and give an exactly linear relationbetween the polarization and the external field. When the external field does not lie inany of the main lattice planes, we find a number of interesting results. The ground statebecomes a complex dipole configuration with no simple regular structure at intermediatevalues of the external field, even though the polarization deviates only slightly from thelinear dependence on the strength of the external field. This result is unchanged when weuse somewhat different sample sizes, up to 20 × ×
20 in some test simulations, suggestingthat these results are not consequences of finite size effects. We also find that the perfectLT ordering persists for relatively weak external fields until we reach a certain value of theexternal field at which the system suddenly turns to a disordered state, before returning toan ordered state at strong fields. Another interesting result that we have found is the factthat the value of the external field at which the system becomes fully polarized coincideswith the value of the internally induced field in the absence of any external field. Ourresults help elucidate the microscopic physics of dielectric and magnetic materials in externalfields, including hybrid organic-inorganic perovskites, and can find applications in recentlyemerging artificial magnetic materials and superstructures.
Acknowledgment
We would like to thank S. Rashkeev for useful discussions. This work was made possibleby NPRP grant ppendix: Internal field in the ferromagnetic state
Here we show that the neighbouring-dipole-induced field vanishes in the ferromagneticstate, independently of the dipole orientation.In the ferromagnetic state, all dipoles have the same components ( p x , p y , p z ). Becauseof the cubic symmetry, for every dipole at location ( x, y, z ) there are dipoles at locationsgiven by the permutations [e.g. ( z, y, x )] as well as mirror image locations [e.g. ( x, − y, z )].Taking the dipole at the origin and the dipole at ( x, y, z ), the first term in the dipole-dipoleinteraction energy (Eq. 1) has the factor ~p i · ~p j = p x + p y + p z = p . (A1)For the second term in Eq. (1), taking the dipoles at the origin and at ( x, y, z ), we find( ~p i · ˆ r ij ) ( ~p j · ˆ r ij ) = ( xp x + yp y + zp z ) r ij = 1 r ij (cid:16) x p x + y p y + z p z + 2 xyp x p y + 2 xzp x p z + 2 yzp y p z (cid:17) (A2)The cross terms, i.e. the last three terms inside the brackets, can be ignored because for everynonzero coordinate, e.g. x , there will be a partner term with the coordinate − x giving thesame term with the opposite sign. Each one of the first three terms has the square of a dipoleCartesian component multiplied by the square of one Cartesian component of the relativeposition vector. When combined with similar terms for permutations of ( x, y, z ), we obtainproducts such as y p x and z p x , such that the sum over all the permutations results in thevalue (cid:16) p x + p y + p z (cid:17) multiplied by the total number of permutations of ( x, y, z ) and dividedby 3. When we take into consideration the factor of 3 in the second term in Eq. (1), we findthat the two terms in the sum exactly cancel each other, independently of the orientationsof the dipoles. We therefore find that the total energy and hence the local field at any dipolemust vanish regardless of the orientation of the dipoles. [1] See e.g. R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, Lon-don, 1982); H. E. Stanley,
Introduction to Phase Transitions and Critical Phenomena (OxfordUniversity Press, New York, 1987); S. Sachdev,
Quantum Phase Transitions (Cambridge Uni-versity Press, Cambridge, 2011).
2] L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. Mod. Phys. , 517 (2008).[3] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. , 277 (2010).[4] T. Baikie,Y. Fang, J. M. Kadro, M. Schreyer, F. Wei, S. G. Mhaisalkar, M. Graetzeld, and T.J. White, J. Mater. Chem. A , 5628 (2013).[5] C. C. Stoumpos, C. D. Malliakas, and M. G. Kanatzidis, Inorg. Chem. , 9019 (2013).[6] S. Ashhab, M. Carignano, and M. E. Madjet, J. Appl. Phys. , 163103 (2019).[7] J. M. Frost, K. T. Butler, and A. Walsh, APL Mater. , 081506 (2014).[8] J. M. Frost, K. T. Butler, F. Brivio, C. H. Hendon, M. van Schilfgaarde, and A. Walsh, NanoLett. , 2584 (2014).[9] A. M. A. Leguy, J. M. Frost, A. P. McMahon, V. Garcia Sakai, W. Kockelmann, C. Law, X.Li, F. Foglia, A. Walsh, B. C. O‘Regan, J. Nelson, J. T. Cabral, and P. R. F. Barnes, NatureCommun. , 7124 (2015).[10] S. N. Rashkeev, F. El-Mellouhi, S. Kais, and F. H. Alharbi, Sci. Rep. , 11467 (2015).[11] J. Lahnsteiner, R. Jinnouchi, and M. Bokdam, arXiv:1905.12540 (2019).[12] S. T. Bramwell, Nature , 212 (1999).[13] J. Cioslowski and A. Nanayakkara, Phys. Rev. Lett. , 2871 (1992).[14] S. Aoyagi, N. Hoshino, T. Akutagawa, Y. Sado, R. Kitaura, H. Shinohara, K. Sugimoto, R.Zhang, and Y. Murata, Chem. Commun. , 524 (2014).[15] B. P. Gorshunov et al ., Nature Commun. , 12842 (2016).[16] R. Skomski, J. Phys.: Condens. Matter , R841 (2003); D. J. Sellmyer and R. Skomski(Eds.), Advanced Magnetic Nanostructures (Springer, New York, 2006).[17] K. Honda and J. Okubo, Phys. Rev. , 705 (1917).[18] L. S. Ornstein and F. Zernike, KNAW Proceedings , 911 (1919).[19] J. M. Luttinger and L. Tisza, Phys. Rev. , 954 (1946); ibid. , 257 (1947).[20] M. Lax, J. Chem. Phys. , 1351 (1952).[21] T. H. Berlin and J. S. Thomsen, J. Chem. Phys. , 1368 (1952).[22] B. R. A. Nijboer and F. W. De Wette, Physica , 422 (1958).[23] A. Horvat, L. Pourovskii, M. Aichhorn, and J. Mravlje, Phys. Rev. B , 205115 (2017).[24] V. V. Kukhtin and O. V. Shramko, Phys. Lett. A , 271 (1988).[25] A. M. Abu-Labdeh, A. B. MacIsaac, J. P. Whitehead, K. DeBell, and M. G. Cottam, Phys.Rev. B , 094412 (2006).
26] A. Yu. Galkin and B. A. Ivanov, JETP Letters , 383 (2006).[27] P. V. Bondarenko, A. Yu. Galkin, B. A. Ivanov, and C. E. Zaspel, Phys. Rev. B , 224415(2010).[28] D. C. Johnston, Phys. Rev. B , 014421 (2016)., 014421 (2016).