text:
# Dipolar interaction between two-dimensional magnetic particles

## I Introduction

In this paper we are interested in two-dimensional magnetic particles interacting through the long range dipolar forces. These particles have a two-dimensional character in two respects: first, because they are platelet shaped, that is to say their thickness $`t`$ is much smaller than their linear size $`L`$; second, because they are arranged on a two dimensional substrate.

They can be obtained, e.g., growing by Molecular Beam Epitaxy (MBE) a magnetic element on a high symmetry substrate. In this case, growth is driven by surface diffusion, nucleation and aggregation: growth\_mbe  in the submonolayer regime, the magnetic overlayer is made up of an ensemble of atomically thick islands that generally are not uniform in size neither arranged in a regular array. The distribution of islands may be regular if nucleation (and therefore island formation) takes place on a reconstructed surface, ricostruzione  or on a network of dislocations. dislocazioni  Alternatively, particles can be produced via lithographic techniques:litografia  in this case, they are much bigger in size and their distribution is generally uniform.

In the following, small particles obtained by MBE growth in the submonolayer regime will also be called islands and large particles obtained by lithographic techniques will also be called dots. Particle is a generic term for both cases.

Each particle is made up of a large number $`𝒩`$ of spins which interact ferromagnetically through the strong intra-particle exchange interaction. In an island, $`𝒩10^210^4`$, while in a dot $`𝒩`$ may be several orders of magnitude greater.

Sufficiently small particles are expected to be in a single domain state, even if their actual magnetic state may depend on several factors: the shape of the particle, the strength of the anisotropies, the single crystal or polycrystalline character of the particle, and so on. In this paper we are assuming that particles are in a single domain state and have a crystalline structure. Within these hypotheses, the magnetic state of an isolated particle is fixed, first of all, by the balancing between dipolar interaction (which has an easy-plane effect) and possible anisotropies favouring the direction perpendicular to the plane ($`z`$ direction, in the following). In the absence of quartic and higher order anisotropies a canted configuration is impossible and the resulting effect may be easy-axis or easy-plane only.

The effect of dipolar interaction between spins belonging to the same (ultrathin) particle has been studied in a previous paper. stripes  We showed that in-plane shape anisotropy is weak: this means that if the magnetization of the particle is within the film plane, its orientation is expected to be settled by the symmetry of the underlying lattice (through magnetocrystalline anisotropies) rather than by the shape of the particle (through the intra-particle dipolar coupling). Such a feature was indeed experimentally observed in MBE-grown Co on Cu(100) ultrathin particles.Science

In this paper we aim to study the interparticle dipolar interaction and to address the following questions:

i) Assuming the particle to be in a single domain state, how is the dipolar interaction between particles modified by their finite size?

ii) Assuming the particles to be arranged in a regular array, is the dipolar ferromagnetic (FM) state a stable configuration?

iii) If the ground state is ferromagnetic, is it stable at finite temperature and what is the value of the Curie temperature?

The previous questions are addressed in Sections II, III and IV, respectively; in Section V the conclusions are drawn.

## II Effective dipolar interaction between particles

In the following we are going to consider two particles $`I_1,I_2`$ of any shape (see Fig. 1), with linear sizes $`L_1,L_2`$ and thicknesses $`t_1,t_2`$. Each particle is a discrete collection of spins and it is supposed to be in a single domain state: each spin of the two particles is indicated by $`\stackrel{}{s}_1`$ and $`\stackrel{}{s}_2`$, respectively. The effective dipolar interaction between the two particles is evaluated by taking the interaction between a spin $`\stackrel{}{s}_1`$ located in $`\stackrel{}{R}_1`$, a spin $`\stackrel{}{s}_2`$ located in $`\stackrel{}{R}_2`$ and summing up on all them:

$$E_{\text{dip}}=\frac{1}{2}\mathrm{\Omega }\frac{t_1t_2}{c_0^2}\underset{\stackrel{}{R}_1}{}\underset{\stackrel{}{R}_2}{}\left[\frac{\stackrel{}{s}_1\stackrel{}{s}_2}{R_{12}^3}3\frac{(\stackrel{}{s}_1\stackrel{}{R}_{12})(\stackrel{}{s}_2\stackrel{}{R}_{12})}{R_{12}^5}\right]$$
(1)
where $`\mathrm{\Omega }=g^2\mu _B^2`$ ($`g`$ is the gyromagnetic factor and $`\mu _B`$ the Bohr magneton) and $`c_0`$ is the interplane distance in the $`z`$ direction. In the previous expression, we have supposed that the thickness of each particle is much smaller than its linear size, $`tL`$. In this hypothesis, $`E_{\text{dip}}`$ is just linear in the numbers of atomic planes, $`t_1/c_0`$ and $`t_2/c_0`$, and $`\stackrel{}{R}_{12}`$ are two dimensional vectors.

As explained in Appendix A, the two quantities appearing in square brackets in Eq. (1) can be expanded in the ratio $`r/D`$, where $`\stackrel{}{D}`$ is the distance between the centers of mass of the particles and $`\stackrel{}{r}=\stackrel{}{r}_2\stackrel{}{r}_1`$ (see Fig. 1). The effective dipolar interaction between two particles at distance $`\stackrel{}{D}`$ takes the approximate form:

$$E_{\text{dip}}E_{\text{dip}}^{(0)}+E_{\text{dip}}^{(2)}$$
(2)
where $`E_{\text{dip}}^{(0)}`$ is the zero-order coupling and $`E_{\text{dip}}^{(2)}`$ takes into account the finite size of the particles:

$$E_{\text{dip}}^{(0)}=\frac{1}{2}\mathrm{\Omega }\left[\frac{\stackrel{}{𝒮}_1\stackrel{}{𝒮}_2}{D^3}3\frac{(\stackrel{}{𝒮}_1\stackrel{}{D})(\stackrel{}{𝒮}_2\stackrel{}{D})}{D^5}\right]$$
(3)
$`E_{\text{dip}}^{(2)}`$ $`=`$ $`\mathrm{\Omega }{\displaystyle \frac{9_{12}}{4}}{\displaystyle \frac{\stackrel{}{𝒮}_1^{}\stackrel{}{𝒮}_2^{}}{D^5}}`$ (4)
$`+\mathrm{\Omega }{\displaystyle \frac{3_{12}}{4}}\left[{\displaystyle \frac{\stackrel{}{𝒮}_1^{}\stackrel{}{𝒮}_2^{}}{D^5}}5{\displaystyle \frac{(\stackrel{}{𝒮}_1^{}\stackrel{}{D})(\stackrel{}{𝒮}_2^{}\stackrel{}{D})}{D^7}}\right]`$
Each particle behaves as a single spin $`\stackrel{}{𝒮}=𝒩\stackrel{}{s}`$, where $`𝒩`$ is the total number of spins in the particle. If $`L`$ is the linear dimension of a particle and $`t`$ its thickness, denoting by $`a_0`$ the in-plane atomic distance, one has $`𝒩=\gamma (L/a_0)^2(t/c_0)`$, where $`\gamma `$ is a geometric factor, depending on the shape of the particle and the lattice structure. The correction terms included in $`E_{\text{dip}}^{(2)}`$ decay with distance as $`1/D^5`$, whilst the usual dipolar interaction decays as $`1/D^3`$. More precisely, $`E_{\text{dip}}^{(2)}`$ is a factor $`_{12}/D^2`$ smaller than $`E_{\text{dip}}^{(0)}`$, where $`_{12}=\frac{1}{2}(_1+_2)`$ is the semisum of the “moments of inertia” $`_i`$ of the two particles (see Appendix A).

In the continuum approximation, we have the following expressions. For a square particle of side $`L`$, $`=L^2/6`$; for a circular particle of radius $`\rho `$, $`=\rho ^2/2`$; for a triangular (equilateral) particle of side $`L`$, $`=L^2/12`$.

Finally, we would like to remark that the isotropic coupling term, proportional to $`(\stackrel{}{𝒮}_1\stackrel{}{𝒮}_2)`$, has different correction terms according to the orientation of the spins: see Eq. (4). In other words, such a term is no more isotropic once the finite size of the particle is taken into account.

In Fig. 2 we compare the exact dipolar coupling $`E_{\text{dip}}`$ (symbols) with the zero order approximation $`E_{\text{dip}}^{(0)}`$ (dashed lines) and with the second order approximation $`(E_{\text{dip}}^{(0)}+E_{\text{dip}}^{(2)}`$) (full lines), in two cases: i) the particles are magnetized perpendicularly to the plane, along the $`\widehat{z}`$ axis (positive coupling energies), and ii) the particles are magnetized in plane along the $`\widehat{x}`$ axis (negative coupling energies). We can see that the second order approximation is fairly good except at very small distances: the smallest allowed distance between centers (without superposing the particles) is $`D_{\text{min}}=\sqrt{2}(L+a_0)=29.7a_0`$, where $`L=20a_0`$ is the side of the square particle. The smallest value of $`D`$ plotted in Fig. 2 is $`D=30a_0`$.

In order to understand which configurations are energetically favoured, let us start by considering just a couple of spins $`\stackrel{}{S}_1,\stackrel{}{S}_2`$, located in plane along the $`\widehat{x}`$ axis. The dipolar coupling can be generally written as

$`E_{12}`$ $``$ $`E_{12}^{(0)}+E_{12}^{(2)}`$
$``$ $`\stackrel{~}{\mathrm{\Omega }}_0(\stackrel{}{S}_1\stackrel{}{S}_23S_1^xS_2^x)`$
$`+\stackrel{~}{\mathrm{\Omega }}_2(2S_1^zS_2^z+\stackrel{}{S}_1\stackrel{}{S}_25S_1^xS_2^x)`$
where the explicit expressions of $`\stackrel{~}{\mathrm{\Omega }}_0`$ and $`\stackrel{~}{\mathrm{\Omega }}_2`$, depending on the distance between spins and on the size of the particles, are irrelevant. We observe that $`\stackrel{}{S}_1,\stackrel{}{S}_2`$ may be either veritable spins $`\stackrel{}{s}_1,\stackrel{}{s}_2`$ or they may represent the effective spins of two particles $`\stackrel{}{𝒮}_1,\stackrel{}{𝒮}_2`$: in the former case one has $`\stackrel{~}{\mathrm{\Omega }}_2=0`$, while in the latter $`\stackrel{~}{\mathrm{\Omega }}_2`$ is the correction due to the finite sizes of the particles.

Both $`E_{12}^{(0)}`$ and $`E_{12}^{(2)}`$ are the sum of competitive interactions. For the sake of definiteness, let us consider $`E_{12}^{(2)}`$. The sum $`2S_1^zS_2^z+\stackrel{}{S}_1\stackrel{}{S}_2`$ favours an antiferromagnetic (AFM) alignment of the spins in the $`\widehat{z}`$ direction, perpendicular to their joining vector. The term $`5S_1^xS_2^x`$ favours a ferromagnetic (FM) alignment along the $`\widehat{x}`$ axis. The energy of the former configuration is $`3\stackrel{~}{\mathrm{\Omega }}_2`$ and the energy of the latter one is $`4\stackrel{~}{\mathrm{\Omega }}_2`$. As proved in Appendix B, where a more detailed discussion is given, the latter configuration is the ground state indeed.

In conclusion, two spins interacting through the dipolar coupling minimize their energy by ordering ferromagnetically along the joining line.

We are now going to discuss the more complex case of a two dimensional lattice of spins. It is well known Maleev  that in the presence of a direct exchange interaction, the system is ferromagnetic and magnetized in the plane, because of an easy-plane effect of $`E_{\text{dip}}`$. If such an exchange interaction is absent, the easy-plane effect survives, but the actual configuration in the plane strongly depends on the lattice structure.Roz91  It is useful to explain the origin of such a dependence on the spin arrangement.

We have seen that spins would like to point along the line joining them: in a two dimensional lattice it is impossible, of course, to fulfill this requirement for all couples of spins. It is possible, however, for a chain of spins: so, we can start by addressing the nature of the coupling between chains.Malo ; Frae99  Let us consider a ferromagnetic chain of spins along the $`\widehat{y}`$ axis of the plane and evaluate the dipolar field $`\stackrel{}{H}_{\text{dip}}`$ generated at a point at distance $`d`$. Spins are oriented in the $`+\widehat{y}`$ direction and $`\stackrel{}{H}_{\text{dip}}=H_{\text{dip}}\stackrel{}{y}`$. In the continuum approximation,

$`H_{\text{dip}}=`$ $``$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }𝒮{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑y\left({\displaystyle \frac{1}{r^3}}3{\displaystyle \frac{y^2}{r^5}}\right)`$
$``$ $`\frac{3}{4}\mathrm{\Omega }𝒮{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑y\left({\displaystyle \frac{1}{r^5}}5{\displaystyle \frac{y^2}{r^7}}\right)`$
where $`r=\sqrt{y^2+d^2}`$.

Both integrals have the form $`𝑑y(\frac{1}{r^{3+n}}(3+n)\frac{y^2}{r^{5+n}})`$. It is sufficient to integrate by parts the term $`1/r^{3+n}`$ to prove that the integral vanishes for any $`n`$ and for any $`d`$. Therefore, in the continuum approximation $`H_{\text{dip}}=0`$.

An exact calculation on a discrete lattice gives a finite value for $`H_{\text{dip}}`$, but its sign depends on the actual lattice structure. Therefore, the reason for the sensitivity of the ground state on the spin arrangement is clear. Dipolar interaction favours the formation of spin chains magnetized ferromagnetically along the chain: these chains are very weakly coupled, in the continuum approximation being even uncoupled. The sign of the coupling and therefore the nature of the ground state do depend on the lattice structure.

In the next Section we are going to analyze the nature of the ground state and to study the spin-wave spectrum with respect to a ferromagnetic alignment.

## III Ground states and spin-wave spectra in two-dimensional lattices

In the following we consider an orthorombic lattice (see Fig. 3a) with primitive vectors

$$\stackrel{}{u}_1=\frac{1}{2}(a,b),\stackrel{}{u}_2=\frac{1}{2}(a,b)$$
(5)
where $`a=2D_0\mathrm{sin}(\varphi /2)`$, $`b=2D_0\mathrm{cos}(\varphi /2)`$, and $`D_0`$ is the rhombus side. For $`\varphi =\frac{\pi }{3}`$ and $`\varphi =\frac{\pi }{2}`$, the orthorombic lattice reduces to a triangular and to a square lattice, respectively. The angle $`\varphi `$ can be supposed smaller than $`\frac{\pi }{2}`$, because if $`\varphi >\frac{\pi }{2}`$ there is just an interchange between $`a`$ and $`b`$ (see Fig. 3a).

It is known Roz91 ; Frae99  that the ground state is ferromagnetic for $`\varphi <\varphi _c80^{}`$ and antiferromagnetic for larger values of the angle $`\varphi `$. It is noteworthy that we find $`\varphi _c`$ to be almost unaffected when $`E_{\text{dip}}^{(2)}`$ is added to $`E_{\text{dip}}^{(0)}`$: as a first approximation, the finite size of the particle has therefore negligible effects on the ground state. We can say as well that the finite size of the particles (i.e. the term $`E_{\text{dip}}^{(2)})`$ reinforces the FM character of the ground state for a triangular lattice, $`\varphi =\frac{\pi }{3}`$, and the AFM one for a square lattice, $`\varphi =\frac{\pi }{2}`$. This statement will be substantiated by the analysis of the spin-wave spectra at the end of this Section.

Let us now calculate the frequency of spin-wave excitations with respect to a ferromagnetic configuration with magnetization directed in plane. We remind that, depending on the value of the rhombic angle $`\varphi `$, the ferromagnetic ground state is different:Roz91  for $`\varphi \pi /3`$, the dipoles are oriented along the short diagonal of the rhombus (i.e. along the $`\widehat{x}`$ axis in Fig. 3a) while for $`\pi /3\varphi \varphi _c`$, along the long diagonal ($`\widehat{y}`$ axis).

The effective dipolar interaction between particles of spins whose centers of mass are located on the sites $`\stackrel{}{D}_i`$ of an orthorombic lattice is $`_{\text{dip}}=_{\text{dip}}^{(0)}+_{\text{dip}}^{(2)}`$; it is obtained from Eqs. (3,4) summing over all sites

$`_{\text{dip}}^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }{\displaystyle \underset{\stackrel{}{D}_i}{}}{\displaystyle \underset{\stackrel{}{D}_j}{}}\left[{\displaystyle \frac{\stackrel{}{𝒮}_i\stackrel{}{𝒮}_j}{D_{ij}^3}}3{\displaystyle \frac{(\stackrel{}{𝒮}_i\stackrel{}{D}_{ij})(\stackrel{}{𝒮}_j\stackrel{}{D}_{ij})}{D_{ij}^5}}\right]`$
$`_{\text{dip}}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }{\displaystyle \underset{\stackrel{}{D}_i}{}}{\displaystyle \underset{\stackrel{}{D}_j}{}}\{3_{ij}\left[{\displaystyle \frac{\stackrel{}{𝒮}_{i}^{}{}_{}{}^{}\stackrel{}{𝒮}_{j}^{}{}_{}{}^{}}{D_{ij}^5}}\right]`$ (6)
$`+{\displaystyle \frac{3_{ij}}{2}}[{\displaystyle \frac{\stackrel{}{𝒮}_i\stackrel{}{𝒮}_j}{D_{ij}^5}}5{\displaystyle \frac{(\stackrel{}{𝒮}_i\stackrel{}{D}_{ij})(\stackrel{}{𝒮}_j\stackrel{}{D}_{ij})}{D_{ij}^7}}]\}`$
where $`\stackrel{}{D}_{ij}=\stackrel{}{D}_j\stackrel{}{D}_i`$ and $`_{ij}=\frac{1}{2}(_i+_j)`$ is the semisum of the “moments of inertia” of two particles whose centers of mass are located in $`\stackrel{}{D}_i`$ and $`\stackrel{}{D}_j`$, respectively. In the following we are assuming to have an array of identical particles, so that $`_i`$.

Taking the $`\widehat{y}`$ direction as quantization axis, we perform the Holstein-Primakoff transformation from spin to boson operators

$$𝒮_j^x=i\sqrt{\frac{𝒮}{2}}(a_j^{}a_j),𝒮_j^y=𝒮a_j^{}a_j,𝒮_j^z=\sqrt{\frac{𝒮}{2}}(a_j+a_j^{})$$

Next, exploiting the translational invariance in the film plane, we introduce the Fourier transform

$$a_j=\sqrt{\frac{1}{N_{}}}\underset{\stackrel{}{k}}{}a_\stackrel{}{k}\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}},a_j^{}=\sqrt{\frac{1}{N_{}}}\underset{\stackrel{}{k}}{}a_\stackrel{}{k}^{}\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}},$$

where $`\stackrel{}{k}=(k_x,k_y)`$ is the two-dimensional in-plane wavevector ranging over the first Brillouin zone, generated by the primitive vectors (see Fig. 3b)

$$\stackrel{}{v}_1=2\pi (\frac{1}{a},\frac{1}{b}),\stackrel{}{v}_2=2\pi (\frac{1}{a},\frac{1}{b})$$
(7)
and $`N_{}`$ is the total number of spins in the two-dimensional lattice. The spin-wave Hamiltonian takes the form

$$_{\text{dip}}=\underset{\stackrel{}{k}}{}\left[A_\stackrel{}{k}a_\stackrel{}{k}^{}a_\stackrel{}{k}+\frac{1}{2}B_\stackrel{}{k}(a_\stackrel{}{k}a_\stackrel{}{k}+a_\stackrel{}{k}^{}a_\stackrel{}{k}^{})\right]$$
(8)
The coefficients $`A_\stackrel{}{k}=A_\stackrel{}{k}^{(0)}+A_\stackrel{}{k}^{(2)}`$ and $`B_\stackrel{}{k}=B_\stackrel{}{k}^{(0)}+B_\stackrel{}{k}^{(2)}`$ can be expressed (see Appendix C) through the dipolar sums ($`\alpha ,\beta =x,y,z`$)

$$𝒟_{\alpha \beta }^{(n)}(\stackrel{}{k})=\underset{\stackrel{}{D}_j}{}\frac{1}{D_{ij}^{3+n}}\left[1(3+n)\frac{D_{ij}^\alpha D_{ij}^\beta }{D_{ij}^2}\right]\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}}$$
(9)
where for $`n=0`$ and $`n=2`$ one has respectively the zero and second order expressions in the ratio $`L/D_0`$ between the linear dimension of the particle and the interparticle distance. The spin-wave energy is

$$ϵ_\stackrel{}{k}=\left[\left(A_\stackrel{}{k}B_\stackrel{}{k}\right)\left(A_\stackrel{}{k}+B_\stackrel{}{k}\right)\right]^{1/2}.$$
(10)
In the continuum limit $`\stackrel{}{k}0`$, one finds the approximate analytic expression

$`ϵ_\stackrel{}{k}`$ $``$ $`\mathrm{\Omega }𝒮[(\mathrm{\Delta }_1(\varphi ){\displaystyle \frac{4\pi k}{ab}}+O(k^2))`$ (11)
$`\times (\mathrm{\Delta }_0(\varphi )+{\displaystyle \frac{4\pi k}{ab}}\mathrm{cos}^2\theta _k+O(k^2))]^{1/2}`$
where $`\mathrm{\Delta }_0`$ and $`\mathrm{\Delta }_1`$, defined in Appendix C, depend on the angle $`\varphi `$. In particular, $`\mathrm{\Delta }_0`$ vanishes for the highly symmetric cases $`\varphi =\pi /3`$ (triangular lattice) and $`\varphi =\pi /2`$ (square lattice). The angle $`\theta _k`$ defines the orientation of the in-plane wavevector $`\stackrel{}{k}`$ with respect to the $`\widehat{x}`$ axis.

It is worth stressing that $`ϵ_\stackrel{}{k}`$ has the same structure both at zero and second order, because $`H_{\text{dip}}^{(2)}`$ does not contribute to the linear term in $`k`$. In particular, we have that the linear term in the quantity $`(A_\stackrel{}{k}+B_\stackrel{}{k})`$ vanishes when $`\stackrel{}{k}`$ is oriented along the magnetization ($`\theta _k=\pi /2`$). Its sign is therefore decided by the quadratic term, which is different for different lattices. It is positive for the triangular lattice (see Fig. 4) and it may be either positive (Fig. 5b) or negative (Fig. 5a) for the square lattice, depending on the orientation of the magnetization.

In Figs. 4,5 we show the dispersion relation of the spin waves, as obtained from Eq. (10), in the case of a triangular lattice ($`\varphi =\pi /3`$) and of a square one ($`\varphi =\pi /2`$), for different orientations $`\theta _k`$.

In the remaining part of this Section we are making general comments on these results.

Dipolar interaction, since it couples the spins to the underlying lattice, is not rotationally invariant, contrary to, e.g., the Heisenberg interaction ($`=J\stackrel{}{s}_1\stackrel{}{s}_2`$). Therefore, the energy of a generic spin configuration depends on how the spins are oriented with respect to the lattice. However, lattices with four-fold ($`\varphi =\pi /2`$) or six-fold ($`\varphi =\pi /3`$) symmetries are special, in the sense that the energy of certain configurations (the ferromagnetic one, for example) is rotationally invariant. This means that the dipolar ferromagnetic ground state of a triangular lattice has a continuous degeneracy.

The coupling between spins and crystal lattice manifests itself in two ways (we are considering the highly symmetric cases, $`\mathrm{\Delta }_0(\varphi )=0`$). First, the spin-wave energy $`ϵ_\stackrel{}{k}`$ depends on the orientation of the wavevector in the continuum limit $`\stackrel{}{k}0`$ as well. According to Eq. (11), $`ϵ_\stackrel{}{k}|\mathrm{cos}\theta _k|k^{1/2}`$, showing that the $`\theta _k`$-dependence of $`ϵ_\stackrel{}{k}`$ is maintained in the limit $`\stackrel{}{k}0`$.

Second, even if the energy of the FM state is degenerate with respect to its orientation in the plane, the spin-wave energy is not (see Fig. 5). This is true for the triangular lattice as well (not shown), but it has most striking consequences for the square lattice. In this case, the FM configuration is known to have a higher energy than the structure where spins are ferromagnetically coupled along lines and antiferromagnetically coupled between neighbouring lines.Roz91  The FM state, if oriented along an axis of the square lattice (the typical configuration considered in the literature) is not even locally stable. This is clearly shown in Fig. 5a where the square of the spin-wave energy is plotted, in the case of $`\stackrel{}{k}`$ and $`\stackrel{}{M}`$ parallel to the a side of the square lattice. The square of the energy is negative, signaling that the corresponding FM state is unstable. However, if the magnetization is oriented along the diagonal, the FM state is no more unstable, as shown in Fig. 5b (we show just the case of $`\stackrel{}{k}`$ parallel to $`\stackrel{}{M}`$, but no instability appears for any value of $`\theta _k`$).

Figures 4 and 5 display the effect of $`E_{\text{dip}}^{(2)}`$ on the spin-wave spectrum: the energy increases for any value and orientation of $`\stackrel{}{k}`$. This result corroborates our statement that the finite size of the particles reinforces ferromagnetism in the triangular lattice (Fig. 4). For the square lattice we have showed that ferromagnetism is locally stable if $`\stackrel{}{M}`$ is oriented along a diagonal of the square lattice (Fig. 5b) and locally unstable if $`\stackrel{}{M}`$ is parallel to a side of the square lattice (Fig. 5a): $`E_{\text{dip}}^{(2)}`$ reinforces the stability in the former case and the instability in the latter one.

## IV The Curie temperature

The Curie transition temperature of an orthorombic lattice of particles with magnetization directed in plane along the $`\widehat{y}`$ axis can be estimated in the framework of spin-wave theory. The relative deviation of the magnetization from the saturation value takes the formMaleev

$`{\displaystyle \frac{\delta 𝒮}{𝒮}}`$ $`=`$ $`{\displaystyle \frac{1}{N_{}𝒮}}{\displaystyle \underset{\stackrel{}{k}}{}}a_\stackrel{}{k}^{}a_\stackrel{}{k}={\displaystyle \frac{𝒱_2}{(2\pi )^2𝒮}}`$ (12)
$`\times `$ $`{\displaystyle d^2\stackrel{}{k}\left[\frac{A_\stackrel{}{k}}{ϵ_\stackrel{}{k}}\frac{1}{\mathrm{e}^{ϵ_\stackrel{}{k}/T}1}\frac{A_\stackrel{}{k}ϵ_\stackrel{}{k}}{2ϵ_\stackrel{}{k}}\right]}`$
where $`𝒱_2`$ is the volume of the two-dimensional unit cell and the integration is over the first Brillouin zone. The first term on the r.h.s. of Eq. (12) gives the temperature dependence of the magnetization, while the second term represents the zero-point spin deviation, which can be safely neglected. A rough estimate of the Curie temperature $`T_C`$ is obtained by imposing that $`\delta 𝒮(T_C)/𝒮1`$.

Since the only energy scale in the problem is given by the dipolar interaction, $`T_C`$ is expected to be of order $`W_{\text{eff}}=\mathrm{\Omega }𝒮^2/D_0^3(g\mu _Bs)^2(L/a_0)^4(t/c_0)^2/D_0^3`$, the effective dipolar interaction between different particles.

The convergence of the integral (12) in $`\stackrel{}{k}=0`$ is easily proved.Malo  For generic orthorombic lattices, it is guaranteed by the gap in the dispersion curve, $`ϵ_00`$. For square and triangular lattices, the spin-wave energy vanishes for $`k0`$ and one can expand the exponential on the denominator of Eq.(12), because the possible divergence is infrared-like. Thus, the integral giving the temperature dependent spin deviation is found to converge, provided that the positive nota  $`O(k^2)`$ terms in the spin-wave energy are taken into account. Eq. (12) can be rewritten as follows:

$$\frac{\delta 𝒮}{𝒮}\frac{T}{W_{\text{eff}}}_0^{2\pi }𝑑\theta _0^{q_M}\frac{qdq}{c_1q\mathrm{cos}^2\theta +c_2q^2}$$
(13)
where $`q=kD_0`$, $`\theta =\theta _k`$, $`c_1`$ and $`c_2`$ are positive constants, and $`q_M1`$. The double integral can be easily evaluated, giving

$$T_CW_{\text{eff}}\frac{c_2}{2\pi }\left[\mathrm{ln}\left(\frac{q_M+\sqrt{q_M^2+(c_1/c_2)^2}}{c_1/c_2}\right)\right]^1.$$
(14)
In the limit $`c_1=0`$, $`ϵ_\stackrel{}{k}=c_2q^2`$ reproduces the dispersion curve of the Heisenberg ferromagnet and the Curie temperature $`T_C`$ vanishes, in agreement with the Mermin-Wagner theorem.MW

Finally, it is interesting to discuss the case where a uniaxial single-ion anisotropy $`\lambda `$, favouring the in-plane $`\widehat{y}`$ axis, is present in the system. The effect of the finite size of the particle is straightforward, in this case:

$$E_{\text{ani}}=\lambda 𝒩\underset{\stackrel{}{D}_i}{}(s_i^y)^2=\stackrel{~}{\lambda }\underset{\stackrel{}{D}_i}{}(𝒮_i^y)^2$$
(15)
where $`\stackrel{~}{\lambda }=\lambda /𝒩`$ and, as usual, $`\stackrel{}{𝒮}=𝒩\stackrel{}{s}`$.

The effect of $`E_{\text{ani}}`$ on the spin-wave energy is simply to add the quantity $`2\lambda `$ to $`A_\stackrel{}{k}`$; for large particles, this constant factor dominates the dipolar terms in Eq. (11). Because of that, spin-wave approximation is no more suitable and Eq. (12) can not be used to evaluate $`T_C`$: in fact, it would give a Curie temperature $`T_C2𝒩\lambda `$, which diverges when $`\lambda \mathrm{}`$.

On physical grounds, we expect that $`T_C`$ is an increasing function of $`\lambda `$, but in the limit of strong anisotropy, $`T_C`$ is always of order $`W_{\text{eff}}`$. An analogy can be done with the three dimensional Heisenberg model in the presence of an easy-axis anisotropy $`\lambda `$: $`T_C`$ increases with $`\lambda `$, but $`T_C(0)`$ and $`T_C(\mathrm{})`$ are of the same order of magnitude and are both of order $`J`$, the exchange coupling constant. In our case, $`W_{\text{eff}}`$ replaces $`J`$ and the dimension of the system is two instead of three (the two-dimensional Heisenberg model is not ordered at finite temperature in the absence of anisotropy).

In order to corroborate our argument, $`T_C`$ will be calculated in the mean field approximation, starting from the Hamiltonian

$$^{MF}=\underset{\stackrel{}{D}_i}{}\left[\xi 𝒮_i^y\stackrel{~}{\lambda }\left(𝒮_i^y\right)^2\right]$$
(16)
where $`\xi =\frac{1}{2}\mathrm{\Omega }𝒮^y𝒟_{yy}^{(0)}(0)`$. The magnetization is given by

$$𝒮^y=\frac{\underset{M=𝒮}{\overset{𝒮}{}}M\mathrm{e}^{(\xi M\stackrel{~}{\lambda }M^2)/T}}{_{M=𝒮}^𝒮\mathrm{e}^{(\xi M\stackrel{~}{\lambda }M^2)/T}}$$
(17)
where $`M`$, $`M^2`$ denote the eigenstates of $`𝒮^y`$, $`(𝒮^y)^2`$ respectively. As $`TT_C`$ one has $`𝒮^y0`$ so that $`\xi 0`$ and the exponential can be expanded. Thus $`𝒮^y\frac{\xi }{T_C}R`$ where $`R`$ is defined as

$$R=\frac{\underset{M=𝒮}{\overset{𝒮}{}}M^2\mathrm{e}^{\frac{\stackrel{~}{\lambda }M^2}{T_C}}}{_{M=𝒮}^𝒮\mathrm{e}^{\frac{\stackrel{~}{\lambda }M^2}{T_C}}}=(𝒮^y)^2_{T_C}.$$
(18)
For $`\stackrel{~}{\lambda }=0`$ one has $`R=\frac{1}{3}𝒮(𝒮+1)`$, so that $`T_C(\stackrel{~}{\lambda }=0)=\frac{1}{6}\mathrm{\Omega }[𝒟_{yy}^{(0)}(0)]𝒮(𝒮+1)W_{\text{eff}}`$.

For $`\stackrel{~}{\lambda }0`$ the mean field critical temperature is

$$\frac{T_C(\stackrel{~}{\lambda })}{T_C(\stackrel{~}{\lambda }=0)}=\frac{3}{𝒮(𝒮+1)}(𝒮^y)^2_{T_C},$$
(19)
where the mean on the r.h.s. must be calulated numerically. Now we observe that in the limit $`\stackrel{~}{\lambda }\mathrm{}`$, one has $`(𝒮^y)^2_{T_C}𝒮^2`$, so that the ratio $`T_C(\mathrm{})/T_C(0)`$ tends to the finite value $`3𝒮/(𝒮+1)`$. We conclude that for any value of $`\stackrel{~}{\lambda }`$, the Curie transition temperature remains of the order of the effective dipolar interaction $`W_{\text{eff}}`$, with a prefactor changing by a factor three as $`\stackrel{~}{\lambda }`$ increases from 0 to $`\mathrm{}`$.

## V Discussion

Let us get back to the three questions formulated at the end of the Introduction.

i) What is the effective dipolar interaction between single domain particles?

In the hypothesis that spins are strongly coupled ferromagnetically inside each particle, it is straightforward to define an effective dipolar coupling between (microscopic) spins $`\stackrel{}{s}`$ of two particles at distance $`D`$: $`W_{\text{eff}}=\mathrm{\Omega }(s𝒩)^2/D^3`$, where $`\mathrm{\Omega }=(g\mu _B)^2`$ and $`𝒩`$ is the number of spins in each particle. If $`L`$ and $`t`$ are respectively its linear size and thickness, $`W_{\text{eff}}(g\mu _Bs)^2(L/a_0)^4(t/c_0)^2/D^3`$. If a single-ion anisotropy $`\lambda `$ is present, its effective value is just nota\_lambda  $`\lambda \text{eff}=\lambda 𝒩\lambda (L/a_0)^2(t/c_0)`$.

The full dipolar interaction between particles can be expanded in (even) powers of $`(L/D)^m`$. The first correction term ($`m=2`$) gives an interaction decaying as $`1/D^5`$ with the distance between particles. A couple of remarks are in order here. First, the ‘purely’ spin term in the dipolar interaction (i.e., the term not coupled to the lattice) is no more rotationally invariant: the $`z`$-components are more strongly coupled than in-plane components. Second, $`E_{\text{dip}}^{(2)}`$ preserves two important features: it is minimized when spins are aligned ferromagnetically along the joining line, and the dipolar field generated by a ‘continuum’ line of spins aligned along the line in a point outside the line, vanishes.

ii) What is the dipolar ground state of an ordered array of magnetic particles?

We have considered the class of orthorombic lattices, which comprises the triangular and the square lattices. It is known that in the case of a lattice of spins, the ground state for the six-fold symmetry is ferromagnetic and for the four-fold symmetry has zero net magnetization. These results are not modified when single-domain particles replace single spins and $`E_{\text{dip}}^{(2)}`$ is considered in addition to $`E_{\text{dip}}^{(0)}`$. The effect of $`E_{\text{dip}}^{(2)}`$ is to reinforce, in some sense, the effect of $`E_{\text{dip}}^{(0)}`$. In particular, therefore, a triangular lattice of two dimensional particles interacting through the dipolar interaction has a ferromagnetic ground state.

iii) What is the finite temperature behaviour of a in-plane ferromagnetically ordered array of particles?

In two dimensional systems, long range order at finite temperature is not certain: however, the long range dipolar forces are known to grant it. In this respect, $`E_{\text{dip}}^{(2)}`$ is of minor importance, because it decays as $`1/D^5`$ and consequently the dipolar sums in the $`\stackrel{}{k}`$-space do not contribute to the terms linear in $`\stackrel{}{k}`$, but to the quadratic terms only.

In the absence of anisotropies, it is elementary that the Curie temperature $`T_C`$ is of order of the effective dipolar coupling, $`T_CW_{\text{eff}}`$, because it is the only energy scale in the problem. If in-plane easy-axis anisotropies are present, we have shown that $`T_C`$ is expected to increase, but not to change in order of magnitude: according to mean field theory, $`T_C`$ increases by a factor three passing from the ‘weak’ anisotropy regime into the ‘strong’ anisotropy regime.

Denoting by $`w=(g\mu _Bs)^2/a_0^3`$ the dipolar coupling between microscopic spins on a two-dimensional lattice with atomic distance $`a_0`$, the effective dipolar interaction between particles of linear dimension $`L`$ and thickness $`t`$ at distance $`D_0`$ can be rewritten as $`W_{\text{eff}}=w(L/a_0)^4(t/c_0)^2/(D_0/a_0)^3`$. Thus, the Curie temperature $`T_CW_{\text{eff}}`$ of an ensemble of magnetic particles may be significantly larger than the Curie temperature of a two-dimensional lattice of microscopic spins, which is of order $`w`$. However, $`L`$ cannot be made larger than $`D_0`$, because $`D_0`$ scales with $`L`$: at the best, therefore, $`T_Cw(L/a_0)(t/c_0)^2`$.

Even if we have no definite evaluation of the numerical prefactor appearing in the previous estimate for $`T_C`$, we suggest that regular arrays of two-dimensional particles magnetized in plane might sustain long range order at experimentally relevant temperatures (see also Ref. JP, ). However, this collective phenomenon may be masked by the superparamagnetic behaviour of the single particle, appearing below the blocking temperature $`T_B\lambda \text{eff}\lambda (L/a_0)^2(t/c_0)`$.

The condition $`T_B>T_C`$ is equivalent to $`\lambda /w>\eta L^2t/D_0^3`$, where $`\eta `$ is an unknown numerical factor. If such a condition is satisfied, the Curie phase transition at $`T=T_C`$ is not visible because thermodinamic equilibrium cannot be attained below $`T_B`$. So, if we decrease temperature from the high-$`T`$ region each particle becomes superparamagnetic when dipolar forces are still unable to induce a long range order in the system. When $`T=T_C`$ dipolar interaction comes into play, but the magnetization of each particle is frozen and the system is unable to attain equilibrium. In the opposite case, $`T_C>T_B`$, a phase transition at $`T_C`$ should be visible.

Recently, the system Co/Cu(001) has drawn the attention because it has been suggested Co\_JMMM ; Co\_PRB  that for $`t<1.8`$ML (ML=monolayer) this system displays a dipolar induced ferromagnetic order. $`T_C(t)`$ is seen Co\_JMMM  to be finite and increase from $`T_C`$(1ML)=25K to $`T_C`$(1.8ML)=200K, where it has a sudden jump, attributed to the percolation in the second layer.

The growth morphology of the system Co/Cu(001) is complicated by the alloying Co\_SS ; Co\_PRB  between the two elements, which mainly takes place in the first layer (75$`\%`$ in Co, 25$`\%`$ in Cu). However, with such a high percentage of Cobalt, it is hard to suppose that the first layer is made up of an ensemble of disconnected Co-islands which interact only through long range forces (in Ref. Co\_PRB,  the percolation threshold is theoretically estimated to be of order 60$`\%`$). More likely, an infinite cluster of Cobalt does exist in the first layer and therefore the magnetic behaviour of the system (and the value of $`T_C`$) follows from the combining effect Co\_PRB  of the direct exchange interaction between spins and long range forces.

Finally, we would like to mention an additional difficulty in the interpretation of experimental data concerning an array of magnetic islands: the random character of deposition and diffusion gives rise to a non-uniform distribution of sizes and positions.JP  This fact, along with frustration due to dipolar interaction, makes difficult even the determination of the ground state, because the system has a glassy behaviour with a lot of metastable states.

## Appendix A Multipolar expansion of interparticle interaction

We start from Eq. (1), assuming single monolayers.nota\_spessore  Summations on $`\stackrel{}{R}_{1,2}`$ are replaced by sums on $`\stackrel{}{r}_{1,2}`$ (see Fig. 1 for notations):

$$E_{\text{dip}}=\frac{1}{2}\mathrm{\Omega }\underset{\stackrel{}{r}_1}{}\underset{\stackrel{}{r}_2}{}\left[\frac{\stackrel{}{s}_1\stackrel{}{s}_2}{R_{12}^3}3\frac{(\stackrel{}{s}_1\stackrel{}{R}_{12})(\stackrel{}{s}_2\stackrel{}{R}_{12})}{R_{12}^5}\right]$$

We make the expansions:

$$\frac{1}{R_{12}^3}\frac{1}{D^3}\left[13\frac{\stackrel{}{r}\stackrel{}{D}}{D^2}\frac{3}{2}\frac{r^2}{D^2}+\frac{15}{2}\frac{(\stackrel{}{r}\stackrel{}{D})^2}{D^4}\right]$$
(20)
$`{\displaystyle \frac{(\stackrel{}{s}_1\stackrel{}{R}_{12})(\stackrel{}{s}_2\stackrel{}{R}_{12})}{R_{12}^5}}`$ $`{\displaystyle \frac{1}{D^5}}\{(\stackrel{}{s}_1\stackrel{}{D})(\stackrel{}{s}_2\stackrel{}{D})[15{\displaystyle \frac{\stackrel{}{r}\stackrel{}{D}}{D^2}}{\displaystyle \frac{5}{2}}{\displaystyle \frac{r^2}{D^2}}+{\displaystyle \frac{35}{2}}{\displaystyle \frac{(\stackrel{}{r}\stackrel{}{D})^2}{D^4}}]+(\stackrel{}{s}_1\stackrel{}{r})(\stackrel{}{s}_2\stackrel{}{r})`$ (21)
$`+`$ $`[(\stackrel{}{s}_1\stackrel{}{D})(\stackrel{}{s}_2\stackrel{}{r})+(\stackrel{}{s}_2\stackrel{}{D})(\stackrel{}{s}_1\stackrel{}{r})][15{\displaystyle \frac{\stackrel{}{r}\stackrel{}{D}}{D^2}}]\}`$
The following expressions are easily calculated:

$`{\displaystyle \underset{\stackrel{}{r}_1}{}}1=𝒩_1{\displaystyle \underset{\stackrel{}{r}_2}{}}1=𝒩_2`$
$`{\displaystyle \underset{\stackrel{}{r}_1}{}}{\displaystyle \underset{\stackrel{}{r}_2}{}}\stackrel{}{r}=𝒩_1𝒩_2(\stackrel{}{r}_2\stackrel{}{r}_1)`$
$`{\displaystyle \underset{\stackrel{}{r}_1}{}}{\displaystyle \underset{\stackrel{}{r}_2}{}}r^2=𝒩_1𝒩_2(r_1^2+r_2^22\stackrel{}{r}_1\stackrel{}{r}_2)`$
$`{\displaystyle \underset{\stackrel{}{r}_1}{}}{\displaystyle \underset{\stackrel{}{r}_2}{}}(\stackrel{}{V}_1\stackrel{}{r})(\stackrel{}{V}_2\stackrel{}{r})=𝒩_1𝒩_2[(\stackrel{}{V}_1\stackrel{}{r}_1)(\stackrel{}{V}_2\stackrel{}{r}_1)`$
$`+(\stackrel{}{V}_1\stackrel{}{r}_2)(\stackrel{}{V}_2\stackrel{}{r}_2)(\stackrel{}{V}_1\stackrel{}{r}_2)(\stackrel{}{V}_2\stackrel{}{r}_1)`$
$`(\stackrel{}{V}_1\stackrel{}{r}_1)(\stackrel{}{V}_2\stackrel{}{r}_2)]`$
In the previous expressions, $`𝒩_{1,2}`$ are the number of spins in the particles $`I_{1,2}`$ and $`\stackrel{}{V}_{1,2}`$ are generic vectors.

It is always possible to choose the origin of the reference system for a given particle in its center of mass, so that $`\stackrel{}{r}_i=0`$ and the “moment of inertia” $`_i`$ of a particle is $`_i=r_i^2`$. Consequently, we have the following results

$`(\stackrel{}{D}\stackrel{}{r}_i)^2`$ $`=`$ $`\frac{1}{2}D^2_i`$ (22)
$`(\stackrel{}{V}_1\stackrel{}{r}_i)(\stackrel{}{V}_2\stackrel{}{r}_i)`$ $`=`$ $`\frac{1}{2}_i\stackrel{}{V}_1^{}\stackrel{}{V}_2^{}`$ (23)
where $`\stackrel{}{V}^{}`$ is the in-plane component of the generic vector $`\stackrel{}{V}`$.

If we define $`_{12}=\frac{1}{2}(_1+_2)`$ and $`\stackrel{}{𝒮}_i=𝒩_i\stackrel{}{s}_i`$, we obtain the following expression for the effective dipolar interaction between two particles at distance $`\stackrel{}{D}`$:

$`E_{\text{dip}}=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }{\displaystyle \frac{\stackrel{}{𝒮}_1^{}\stackrel{}{𝒮}_2^{}}{D^3}}\left(1+{\displaystyle \frac{9}{2}}{\displaystyle \frac{_{12}}{D^2}}\right)`$ (24)
$`+`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }{\displaystyle \frac{\stackrel{}{𝒮}_1^{}\stackrel{}{𝒮}_2^{}}{D^3}}\left(1+{\displaystyle \frac{3}{2}}{\displaystyle \frac{_{12}}{D^2}}\right)`$ (25)
$``$ $`{\displaystyle \frac{3}{2}}\mathrm{\Omega }{\displaystyle \frac{(\stackrel{}{𝒮}_1^{}\stackrel{}{D})(\stackrel{}{𝒮}_2^{}\stackrel{}{D})}{D^5}}\left(1+{\displaystyle \frac{5}{2}}{\displaystyle \frac{_{12}}{D^2}}\right)`$ (26)
## Appendix B Minimization of the dipolar coupling between two spins

Let us consider two unitary spins $`\stackrel{}{S}_1,\stackrel{}{S}_2`$ located at a distance $`R_{12}`$ along the in-plane $`\widehat{x}`$ axis, taken as the polar axis, while the $`\widehat{z}`$ axis is perpendicular to the plane. Their orientations are defined by the polar and azimuthal angles $`\theta _i,\phi _i`$.

The dipolar coupling can be generally written as

$`E_{12}`$ $``$ $`E_{12}^{(0)}+E_{12}^{(2)}`$
$``$ $`\stackrel{~}{\mathrm{\Omega }}_0(\stackrel{}{S}_1\stackrel{}{S}_23S_1^xS_2^x)`$
$`+\stackrel{~}{\mathrm{\Omega }}_2(2S_1^zS_2^z+\stackrel{}{S}_1\stackrel{}{S}_25S_1^xS_2^x)`$
where $`\stackrel{~}{\mathrm{\Omega }}_0=\frac{1}{2}\mathrm{\Omega }S^2/R_{12}^3`$ and $`\stackrel{~}{\mathrm{\Omega }}_2=\frac{3}{4}_{12}\mathrm{\Omega }S^2/R_{12}^5`$ (see Eqs. 3,4).

Both terms $`E_{12}^{(0)}`$ and $`E_{12}^{(2)}`$ are minimized by the same configuration. First, let us treat the zero-order term. We have to minimize the function

$$E_{12}^{(0)}/\stackrel{~}{\mathrm{\Omega }}_0=\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{cos}(\phi _1\phi _2)2\mathrm{cos}\theta _1\mathrm{cos}\theta _2$$
(27)
By taking the derivatives with respect to $`\phi _{1,2}`$ we find that $`(\phi _1\phi _2)=0,\pi `$ or that one $`\theta _i`$ at least must vanish. If, e.g., $`\theta _1=0`$, it is straightforward to derive that $`\theta _2=0`$ as well. If both $`\theta _{1,2}`$ are not vanishing, taking the derivatives with respect to them implies $`\mathrm{cos}\theta _1=\mathrm{cos}\theta _2=0`$, i.e. $`\theta _{1,2}=\frac{\pi }{2}`$. In simple words, it is sufficient to consider two kinds of configurations: i) the configuration where both spins are perpendicular to the joining vector and they are parallel or antiparallel; ii) the configuration where both spins are aligned along the joining vector. The former configuration, with antiparallel spins, corresponds to the minimization of $`(\stackrel{}{S}_1\stackrel{}{S}_2)`$ and its energy is $`\stackrel{~}{\mathrm{\Omega }}_0`$. The latter configuration, with parallel spins, corresponds to the minimization of $`(3S_1^xS_2^x)`$ and its energy is $`\stackrel{~}{\mathrm{\Omega }}_03\stackrel{~}{\mathrm{\Omega }}_0=2\stackrel{~}{\mathrm{\Omega }}_0`$. We conclude that $`E_{12}^{(0)}`$ is minimized by the ferromagnetic configuration with spins aligned along their joining vector.

The minimization of $`E_{12}^{(2)}`$ proceeds along the same lines, with the minor difference that configurations with both spins perpendicular to the joining line are no more degenerate with respect to a global rotation around the $`\widehat{x}`$ axis: the lowest energy one corresponds to antiparallel spins along the $`\widehat{z}`$ direction and its energy is $`3\stackrel{~}{\mathrm{\Omega }}_2`$. On the other hand, the ferromagnetic configuration with both spins parallel to the $`\widehat{x}`$ axis has the energy $`4\stackrel{~}{\mathrm{\Omega }}_2`$, so the conclusion is unchanged.

## Appendix C Dipolar sums for the orthorombic lattice

For spins ferromagnetically oriented along the $`y`$ axis, the coefficients of the spin-wave Hamiltonian in Eq. (8) take the form

$`{\displaystyle \frac{A_\stackrel{}{k}^{(0)}}{\mathrm{\Omega }𝒮}}`$ $`=`$ $`\left[{\displaystyle \frac{1}{2}}𝒟_{zz}^{(0)}(\stackrel{}{k})+{\displaystyle \frac{1}{2}}𝒟_{xx}^{(0)}(\stackrel{}{k})𝒟_{yy}^{(0)}(0)\right]`$
$`{\displaystyle \frac{B_\stackrel{}{k}^{(0)}}{\mathrm{\Omega }𝒮}}`$ $`=`$ $`\left[{\displaystyle \frac{1}{2}}𝒟_{zz}^{(0)}(\stackrel{}{k}){\displaystyle \frac{1}{2}}𝒟_{xx}^{(0)}(\stackrel{}{k})+i𝒟_{zx}^{(0)}(\stackrel{}{k})\right]`$ (28)
and

$`{\displaystyle \frac{A_\stackrel{}{k}^{(2)}}{\mathrm{\Omega }𝒮\frac{3}{2}}}`$ $`=`$ $`\left[{\displaystyle \frac{1}{2}}𝒟_{zz}^{(2)}(\stackrel{}{k})+{\displaystyle \frac{1}{2}}𝒟_{xx}^{(2)}(\stackrel{}{k})𝒟_{yy}^{(2)}(0)+^{(2)}(\stackrel{}{k})\right]`$
$`{\displaystyle \frac{B_\stackrel{}{k}^{(2)}}{\mathrm{\Omega }𝒮\frac{3}{2}}}`$ $`=`$ $`\left[{\displaystyle \frac{1}{2}}𝒟_{zz}^{(2)}(\stackrel{}{k}){\displaystyle \frac{1}{2}}𝒟_{xx}^{(2)}(\stackrel{}{k})+i𝒟_{zx}^{(2)}(\stackrel{}{k})+^{(2)}(\stackrel{}{k})\right]`$
where $`𝒩`$ is the number of spins in each particle and the dipolar sums are defined as ($`\alpha ,\beta =x,y,z`$)

$`𝒟_{\alpha \beta }^{(n)}(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{D}_j}{}}{\displaystyle \frac{1}{D_{ij}^{3+n}}}\left[1(3+n){\displaystyle \frac{D_{ij}^\alpha D_{ij}^\beta }{D_{ij}^2}}\right]\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}}`$
$`^{(n)}(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{D}_j}{}}{\displaystyle \frac{\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}}}{D_{ij}^{3+n}}}`$ (30)
We observe that in the ultraflat particle limit, one has $`^{(n)}(\stackrel{}{k})D_{zz}^{(n)}(\stackrel{}{k})`$ and $`D_{zx}^{(n)}(\stackrel{}{k})0`$ and the $`𝒟_{\alpha \alpha }`$’s can be approximately expressed

$`𝒟_{xx}^{(n)}(\stackrel{}{k})`$ $``$ $`Y^{(n)}(\stackrel{}{k})(2+n)X^{(n)}(\stackrel{}{k})`$
$`𝒟_{yy}^{(n)}(\stackrel{}{k})`$ $``$ $`X(n)(\stackrel{}{k})(2+n)Y^{(n)}(\stackrel{}{k})`$
$`𝒟_{zz}^{(n)}(\stackrel{}{k})`$ $``$ $`X^{(n)}(\stackrel{}{k})+Y^{(n)}(\stackrel{}{k})`$ (31)
in terms of the dipolar sums

$`X^{(n)}(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{D}_j}{}}{\displaystyle \frac{(D_{ij}^x)^2}{D_{ij}^{5+n}}}\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}}`$
$`Y^{(n)}(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{D}_j}{}}{\displaystyle \frac{(D_{ij}^y)^2}{D_{ij}^{5+n}}}\mathrm{e}^{i\stackrel{}{k}\stackrel{}{D}_{ij}}.`$ (32)
The latter sums can be numerically calculated in a very efficient way following a method, developed some years ago by Benson and Mills,BM  similar to Ewald’s one for the evaluation of lattice sums. For $`n=0`$ one has

$$X^{(0)}(k_x,k_y;a,b)=U^{(0)}(k_x,k_y)+V^{(0)}(k_x,k_y)$$
(33)
where $`U`$ is a summation over the sites of the orthorombic lattice which are located at integer multiples of $`a`$ and $`b`$

$`U^{(0)}(k_x,k_y)`$ $`=`$ $`{\displaystyle \frac{16}{3}}{\displaystyle \frac{1}{b}}{\displaystyle \underset{l=1}{\overset{+\mathrm{}}{}}}\mathrm{cos}(k_xal)`$
$`{\displaystyle \underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}}({\displaystyle \frac{\pi m}{b}}`$ $`+`$ $`{\displaystyle \frac{k_y}{2}})^2K_2[2al|{\displaystyle \frac{\pi m}{b}}+{\displaystyle \frac{k_y}{2}}|]`$
while $`V`$ refers to semi-integer multiples

$`V^{(0)}(k_x,k_y)`$ $`=`$ $`{\displaystyle \frac{16}{3}}{\displaystyle \frac{1}{b}}{\displaystyle \underset{l=0}{\overset{+\mathrm{}}{}}}\mathrm{cos}\left[k_xa(l+{\displaystyle \frac{1}{2}})\right]`$
$`{\displaystyle \underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}}(1)^m({\displaystyle \frac{\pi m}{b}}`$ $`+`$ $`{\displaystyle \frac{k_y}{2}})^2K_2[2a(l+{\displaystyle \frac{1}{2}})|{\displaystyle \frac{\pi m}{b}}+{\displaystyle \frac{k_y}{2}}|]`$
For $`n=2`$ one has

$$X^{(2)}(k_x,k_y;a,b)=U^{(2)}(k_x,k_y)+V^{(2)}(k_x,k_y)$$
(34)
where

$`U^{(2)}(k_x,k_y)`$ $`=`$ $`{\displaystyle \frac{32}{15}}{\displaystyle \frac{1}{ab}}{\displaystyle \underset{l=1}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{\mathrm{cos}(k_xal)}{l}}`$
$`{\displaystyle \underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}}({\displaystyle \frac{\pi m}{b}}`$ $`+`$ $`{\displaystyle \frac{k_y}{2}})^3K_3[2al|{\displaystyle \frac{\pi m}{b}}+{\displaystyle \frac{k_y}{2}}|]`$
and

$`V^{(2)}(k_x,k_y)`$ $`=`$ $`{\displaystyle \frac{32}{15}}{\displaystyle \frac{1}{ab}}{\displaystyle \underset{l=0}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{\mathrm{cos}\left[k_xa(l+\frac{1}{2})\right]}{(l+\frac{1}{2})}}`$
$`{\displaystyle \underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}}(1)^m({\displaystyle \frac{\pi m}{b}}`$ $`+`$ $`{\displaystyle \frac{k_y}{2}})^3K_3[2a(l+{\displaystyle \frac{1}{2}})|{\displaystyle \frac{\pi m}{b}}+{\displaystyle \frac{k_y}{2}}|]`$
Here above, $`K_2(z)`$ and $`K_3(z)`$ are the modified Bessel functions of second and third order, respectively. By symmetry reasons, one has

$$Y^{(n)}(k_x,k_y;a,b)=X^{(n)}(k_y,k_x;b,a)$$
(35)
Finally, in the continuum limit $`\stackrel{}{k}0`$, one can obtain analytical expressions for the dipolar sums. Denoting by $`k`$ the modulus of the two-dimensional wavevector $`\stackrel{}{k}=(k_x,k_y)`$ and by $`\theta _k`$ the angle that $`\stackrel{}{k}`$ forms with the $`\widehat{x}`$ axis, we obtain

$`X^{(0)}(k)`$ $``$ $`X^{(0)}(0){\displaystyle \frac{2\pi }{ab}}k\left[1+{\displaystyle \frac{1}{3}}\mathrm{cos}(2\theta _k)\right]`$
$`Y^{(0)}(k)`$ $``$ $`Y^{(0)}(0){\displaystyle \frac{2\pi }{ab}}k\left[1{\displaystyle \frac{1}{3}}\mathrm{cos}(2\theta _k)\right]`$ (36)
and

$`X^{(2)}(k)`$ $``$ $`X^{(2)}(0)+O(k^2)`$
$`Y^{(2)}(k)`$ $``$ $`Y^{(2)}(0)+O(k^2)`$ (37)
where for the triangular lattice one has $`X^{(0)}(0)=5.5170879/D_0^3`$, $`X^{(2)}(0)=3.3809493/D_0^3`$ and for the square lattice $`X^{(0)}(0)=4.5168109/D_0^3`$, $`X^{(2)}(0)=2.5451291/D_0^3`$. From the previous expressions one obtains approximate expansions for the quantities in Eqs. (10,12):

$`{\displaystyle \frac{A_\stackrel{}{k}B_\stackrel{}{k}}{\mathrm{\Omega }𝒮}}`$ $``$ $`\mathrm{\Delta }_0+{\displaystyle \frac{4\pi k}{ab}}\mathrm{cos}^2\theta _k+O(k^2)`$
$`{\displaystyle \frac{A_\stackrel{}{k}+B_\stackrel{}{k}}{\mathrm{\Omega }𝒮}}`$ $``$ $`\mathrm{\Delta }_1{\displaystyle \frac{4\pi k}{ab}}+O(k^2)`$
$`{\displaystyle \frac{A_\stackrel{}{k}}{\mathrm{\Omega }𝒮}}`$ $``$ $`\mathrm{\Delta }_2{\displaystyle \frac{2\pi k}{ab}}\mathrm{sin}^2\theta _k+O(k^2)`$ (38)
where

$`\mathrm{\Delta }_0`$ $`=`$ $`3\left[Y^{(0)}(0)X^{(0)}(0)\right]+{\displaystyle \frac{15}{2}}\left[Y^{(2)}(0)X^{(2)}(0)\right]`$
$`\mathrm{\Delta }_1`$ $`=`$ $`3Y^{(0)}(0)+{\displaystyle \frac{3}{2}}\left[7Y^{(2)}(0)+2X^{(2)}(0)\right]`$
$`\mathrm{\Delta }_2`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left[2Y^{(0)}(0)X^{(0)}(0)\right]+{\displaystyle \frac{9}{4}}\left[4Y^{(2)}(0)X^{(2)}(0)\right]`$
Hence we observe that for a generic orthorombic lattice one has $`X^{(n)}(0)Y^{(n)}(0)`$, so that the dispersion curve $`ϵ_\stackrel{}{k}=\left[(A_\stackrel{}{k}B_\stackrel{}{k})(A_\stackrel{}{k}+B_\stackrel{}{k})\right]^{1/2}`$ presents a gap for $`\stackrel{}{k}0`$. For the special cases of the triangular and the square lattice one has $`X^{(n)}(0)=Y^{(n)}(0)`$ by symmetry reasons; this implies that $`\mathrm{\Delta }_0=0`$, so that a Goldstone mode is present in the dispersion curve: $`ϵ_\stackrel{}{k}0`$ for $`\stackrel{}{k}0`$.

###### Acknowledgements.

We thank Angelo Rettori for his critical reading of the manuscript.