Microscopic Aspects of Magnetic Lattice Demagnetizing Factors

TL;DR

This paper investigates how microscopic details like spin type, orientation, and shape influence the demagnetizing factor N in magnetic lattices, combining theory, simulations, and experiments for a comprehensive understanding.

Contribution

It reveals that N depends on microscopic spin and shape details, challenging the common ellipsoid approximation in magnetic systems with long-range interactions.

Findings

N varies with spin dimensionality and orientation

Monte Carlo simulations confirm theoretical predictions

Experimental N for Dy2Ti2O7 matches calculations closely

Abstract

The demagnetizing factor N is of both conceptual interest and practical importance. Considering localized magnetic moments on a lattice, we show that for non-ellipsoidal samples, N depends on the spin dimensionality (Ising, XY, or Heisenberg) and orientation, as well as the sample shape and susceptibility. The generality of this result is demonstrated by means of a recursive analytic calculation as well as detailed Monte Carlo simulations of realistic model spin Hamiltonians. As an important check and application, we also make an accurate experimental determination of N for a representative collective paramagnet (i.e. the Dy2Ti2O7 spin ice compound) and show that the temperature dependence of the experimentally determined N agrees closely with our theoretical calculations. Our conclusion is that the well established practice of approximating the true sample shape with "corresponding ellipsoids" for systems with long-range interactions will in many cases overlook important effects stemming from the microscopic aspects of the system under consideration.

Tables

Table 1. Table D.1: RMS error and extrapolated N 𝑁 N for various fitting functions applied to the MC susceptibility calculated with open boundary conditions for χ int = 1.82 subscript 𝜒 int 1.82 \chi_{\text{int}}=1.82 (red squares in Fig. 3 in the main text). The data points for the seven largest system sizes are included in the fit.

function	RMSE (${10}^{-6}$)	N
$a+b/L^1+c/L^2$	4.87	0.3238
$a+b/L^2+c/L^3$	9.63	0.3203
$a+b/L^3+c/L^6$	83.1	0.3177
$a+b/L^3+c/L^4$	23.1	0.3187
$a+b/L^3+c/L^5$	28.9	0.3184
$a+b/L^2+c/L^4$	13.1	0.3198
$a+b/L^1+c/L^3$	4.92	0.3235

Table 2. Table D.2: RMS error and calculated N 𝑁 N for various fitting functions applied to the MC susceptibility calculated with periodic boundary conditions for χ int = 1.82 subscript 𝜒 int 1.82 \chi_{\text{int}}=1.82 (red squares in Fig. 3 in the main text). The data points for the seven largest system sizes are included in the fit.

function	RMSE (${10}^{-6}$)	N
$a+b/L^1+c/L^2$	7.60	0.3235
$a+b/L^2+c/L^3$	3.87	0.3238
$a+b/L^3+c/L^6$	3.68	0.3238
$a+b/L^3+c/L^4$	3.76	0.3238
$a+b/L^3+c/L^5$	3.72	0.3238
$a+b/L^2+c/L^4$	4.02	0.3242
$a+b/L^1+c/L^3$	3.93	0.3238

Equations

$$\frac{1}{{\chi }_{\text{int}}}=\frac{1}{{\chi }_{\text{exp}}}-N.$$

$$B_i^{\parallel }=B_i^{\parallel ,\text{dip}}+B_i^{\parallel ,\text{ext}}+B_i^{\parallel ,\text{self}},$$

$$B_i^{\parallel ,\text{dip}}=\frac{{\mu }_0{\mu }_{\text{B}}}{4\pi }\sum _{j\ne i}\left(\frac{3(\hat{\text{ȷ}}\cdot {\hat{\mathbf{r}}}_{ij})(\hat{\text{ı}}\cdot {\hat{\mathbf{r}}}_{ij})-\hat{\text{ȷ}}\cdot \hat{\text{ı}}}{r_{ij}^3}\right)m_j.$$

$$B_i^{\parallel ,\text{self}}=B^{1,\parallel }+B^{2,\parallel }={\mu }_0{\mu }_B\frac{\mathcal{N}}{V}\left[1-N_0{\cos }^2\theta \right]m_i,$$

$$m_i=\frac{V}{\mathcal{N}}\left(\frac{{\chi }_{\text{loc}}}{{\chi }_{\text{loc}}+1}\right)\frac{B_i^{\parallel }}{{\mu }_0{\mu }_B},$$

$${\chi }_{\text{exp}}={\chi }_{\text{exp}}^{zz}={\left(\frac{\partial M^z}{\partial H_{\text{ext}}^z}\right)}_T=\frac{{\mu }_0{\mu }_B}{V{B}^{\text{ext}}}\sum _{i=1}^{\mathcal{N}}{m}_i\cos \theta ,$$

$$N={\left[\frac{{\mu }_0{\mu }_B}{V{B}^{\text{ext}}}\sum _{i=1}^{\mathcal{N}}{m}_i\cos \theta \right]}^{-1}-\frac{1}{{\chi }_{\text{loc}}{\cos }^2\theta }.$$

$$\mathcal{H}=\frac{{\mu }_0{\mu }^2}{4\pi }\sum _{i>j}{\Lambda }_{ij}{\sigma }_i{\sigma }_j,$$

$${\chi }^{zz}=\frac{\partial M^z}{\partial H^z}=\frac{{\mu }_0{\mu }^2}{k_{\text{B}}TV}⟨{\left(\sum _{i=1}^{\mathcal{N}}{\sigma }_i\cos \theta \right)}^2⟩.$$

$${\mathcal{H}}_{\text{exchange}}=\sum _{i>j}{J}_{ij}\hat{\text{ı}}\cdot \hat{\text{ȷ}}{\sigma }_i{\sigma }_j,$$

Full text

I Introduction

Long-range interactions are important in many areas of science, from cosmology, through the gravitational interaction, to biology, through Coulomb’s law. A long-range interaction may be defined in $d$ spatial dimensions by its two-body potential $V(r)$ scaling with distance $r$ as $r^{-\alpha}$ where $\alpha\leq d$ Campa et al. (2009). The paramount problem in such systems is how to integrate $V(r)$ over an extended system. Following Newton and Euler, the analysis of general systems has been largely based on the exact solutions for spheres and ellipsoids Osborn (1945); Stoner (1945); Beleggia et al. (2006a, b); Di Fratta (2016). This raises the question of whether approximating other shapes to “corresponding ellipsoids” Aharoni (1998) just neglects uninteresting details or whether there are crucial properties that are lost in the approximation. The demagnetizing problem in magnetic systems is a natural setting for exploring this question since it is accessible and of intrinsic importance in experiments, and constitutes a paragon for exploring the thermodynamics of long-range interacting systems Campa et al. (2009). Demagnetizing effects are also important in superconductors, while analogues occur, for example, in electric systems Ponomareva et al. (2005) (depolarizing factor), in the problem of strain fields around inclusions Eshelby (1957), and in the treatment of avalanching systems in confined geometries Durin and Zapperi (2000); Csikor et al. (2007); Laurson et al. (2013).

In an applied magnetic field $H_{\rm ext}=B_{\rm ext}/\mu_{0}$, the thermodynamic energy of an ellipsoid of volume $V$ and magnetic moment $m$ acquires a contribution $E_{\textrm{mag}}=(\mu_{0}/2)Nm^{2}/V$, where $N$ is the demagnetizing factor. After subtracting $E_{\textrm{mag}}$ from the total energy, differentiation with respect to the magnetization, $M\equiv m/V$, defines the internal field $H_{\text{int}}\equiv H_{\text{ext}}+H_{\text{d}}$, where $H_{\rm d}=-NM$ is the demagnetizing field. The intrinsic magnetic susceptibility $\chi_{\text{int}}=\partial M/\partial H_{\text{int}}$ is a shape-independent material property derived from the experimentally determined susceptibility $\chi_{\text{exp}}\equiv\partial M/\partial H_{\text{ext}}$ through

[TABLE]

The determination of $N$ is a fundamental problem that dates back to the work of Poisson and Maxwell Chen et al. (1991). In the 1940s, Osborn Osborn (1945) and Stoner Stoner (1945) tabulated $N$ for general ellipsoids, while more recently, Aharoni Aharoni (1998) treated cuboids in the $\chi_{\text{int}}\!\rightarrow\!0$ limit. These highly cited papers bear witness to the importance of accurately computable and easily accessible demagnetizing factors. Given that i) it was realized already in the 1920s that $N$ for a non-ellipsoidal sample is a function not only of the sample shape, but also of $\chi_{\text{int}}$ itself Würschmidt (1923); Stäblein and Schlechtweg (1935), and that ii) many experiments are routinely performed not on ellipsoids but on cuboids Quilliam et al. (2008); Higashinaka and Maeno (2005), it is perhaps remarkable that it was only very recently that the $\chi$-dependence of $N$ was calculated for cuboids away from the $\chi_{\text{int}}\!\rightarrow\!0$ limit Chen et al. (2002, 2005).

The existence of demagnetizing factors for cuboids suggests that their thermodynamics may be formulated in terms of an internal field, with corrections that become dependent on both shape and temperature Chen et al. (2002) (through $\chi_{\text{int}}$). In this work, we have found that, for magnetic lattices, the demagnetizing factor of cuboids depends also on the local spin symmetry and allowed orientations of the magnetic moments. With reference to the question posed at the very beginning, our result illustrates a case where a long-range interaction integrates in a qualitatively different way for a cuboid and an ellipsoid, such that the discrete microscopic nature of the system matters in the former case but not in the latter. We are aware of only a few previous studies where effects of such discreteness have been discussed  Christensen and Mørup (1983); Vedmedenko et al. (2003); Millev et al. (2003); Ponomareva et al. (2005). Our interest in this problem was spurred by the recent experimental observation of anomalous demagnetizing effects in the spin ice material Dy2Ti2O7 Bovo et al. (2013).

One may ask whether small differences in the estimated $N$ really matter for exposing important physics. The answer is found in Eq. (1). If $\chi_{\text{exp}}\ll 1$, then $\chi_{\text{int}}$ is insensitive to the precise value of $N$. However, in many physical systems that display unusual and interesting magnetic phenomena, $\chi_{\text{exp}}$ is large, and $\chi_{\text{int}}$ becomes a sensitive function of $N$. Examples include the spin ice materials Dy2Ti2O7 and Ho2Ti2O7, which support magnetic monopole excitations Castelnovo et al. (2008), and LiHo1-xYxF4 which displays ultra-slow relaxation Biltmo and Henelius (2012). Important demagnetizing effects are manifest when an accurately directed field is required: for example in experiments on the elusive Kasteleyn transition Fennell et al. (2007), sub-lattice pinning Ruff et al. (2005); Higashinaka and Maeno (2005); Sato et al. (2006, 2007) and multiple field-driven transitions Schiffer et al. (1994); or else for disentangling the in- and out-of-phase frequency response Quilliam et al. (2011). In such cases, quantitative conclusions and accurate tests of theory depend, through $\chi_{\rm int}$, on an accurate knowledge of $N$. Our work illustrates how this may be achieved.

The rest of the paper is organized as follows. In Section II we discuss how to determine $N$ experimentally. In Section III we introduce an iterative method for obtaining $N$, and we consider in Section IV a Monte Carlo calculation of $N$. Finally, we close the paper with a discussion in Section V. For details regarding the experimental and numerical procedures we refer the reader to Appendices A-D. The effects of short-range interactions are considered in some detail in Appendix E.

II Experimental determination of $N$

To illustrate the importance of the demagnetizing correction, and to test the theory presented below, we first present the experimental determination of $N$ for a particular case. The localized-moment paramagnet ${\rm Dy_{2}Ti_{2}O_{7}}$ (a spin ice) is well-suited to this purpose as it has a large susceptibility, is crystallographically well-defined (in the cubic space group $Fd\bar{3}m$) with no evidence of crystal distortion Ruminy et al. (2016), and can be accurately cut into high-quality single crystal samples of different shape. Since its spin Hamiltonian has been established in great detail Yavors’kii et al. (2008); Henelius et al. (2016), it is convenient to adopt ${\rm Dy_{2}Ti_{2}O_{7}}$ as a model system for studying the demagnetizing factor.

A sphere of diameter 4 mm and a cube of dimensions $2\times 2\times 2$ mm3, with edges precisely oriented along the cubic crystallographic axes $[100]$, $[010]$ and $[001]$ directions, were commercially hand-cut from different larger crystals of Dy2Ti2O7 provided by D. Prabhakaran Prabhakaran and Boothroyd (2011) (see Ref. Bovo et al. (2013)). The cube was epi-polished on all sides epi . Crystal shape, orientation, and experimental conditions were carefully controlled to minimize measurement errors; see Appendix A. The experimental susceptibilities of both the sphere and the cube ($\chi_{\textrm{exp}}^{\textrm{sphere}}$, $\chi_{\textrm{exp}}^{\textrm{cube}}$) were determined from measurements of the magnetic moment.

Setting the demagnetizing factor of the sphere to $N_{\rm sphere}=1/3$, that of the cube was determined through Eq. (1), i.e., $N_{\rm cube}=1/\chi_{\textrm{exp}}^{\textrm{cube}}-1/\chi_{\textrm{exp}}^{\textrm{sphere}}+N_{\rm sphere}$. In order to match the susceptibility of the cube and sphere in the high-$T$ limit, $\chi_{\textrm{exp}}^{\textrm{cube}}$ was shifted by about 1% ($\chi_{\textrm{exp}}^{\textrm{cube}}\rightarrow\chi_{\textrm{exp}}^{\textrm{cube}}/1.0074$) before calculating $N_{\rm cube}$. Fig. 1 shows how the experimental $N_{\rm cube}$ departs significantly from the 1/3 value when $\chi\gtrsim 1$. This is the main experimental result of our study. The inset of the figure compares the uncorrected susceptibility data and the data derived from assuming $N=1/3$ for both samples. The predicted theoretical continuation of the experimental data below 2 K (dashed curves) is based on a generalized version of the dipolar spin ice model Yavors’kii et al. (2008); Henelius et al. (2016).

III Determination of $N$ via an iterative method

In this section we introduce an iterative method to calculate the on-site field distribution inside a linear magnetic material placed in a uniform magnetic field. In the iterative algorithm we first assume that $H_{\text{int}}$ equals $H_{\text{ext}}$ and calculate the induced local magnetization for an assumed $\chi_{\text{int}}$. This magnetization generates a demagnetizing field that, in turn, modifies $H_{\text{int}}$. The resulting field-magnetization equations are iterated until convergence. With the converged field and magnetization distributions in hand, one then computes $N$.

To proceed, we consider a sample of volume $V$ with $\mathcal{N}$ magnetic moments. As a first case, we focus on Ising moments $\bm{m}_{i}=m_{i}\mu_{B}{{\bm{\hat{\textnormal{\i}}}}}$, where ${{\bm{\hat{\textnormal{\i}}}}}$ is the unit vector in the local Ising direction at site $i$, and $m_{i}$ is dimensionless. We first determine the component of the local field along the Ising moment at site $i$, $B_{i}^{\shortparallel}=\bm{B}_{i}\cdot{{\bm{\hat{\textnormal{\i}}}}}$, which is the sum of three contributions:

[TABLE]

which we now discuss one by one.

First, the dipolar field at site $i$ produced by all the other point magnetic dipoles within the sample, $B_{i}^{\shortparallel,\text{dip}}\equiv\bm{B}_{i}^{\text{dip}}\cdot{{\bm{\hat{\textnormal{\i}}}}}$, is given by the familiar form Griffiths (1999)

[TABLE]

Second, we consider an external field in the global $\bm{\hat{z}}$ direction, $\bm{B}^{\text{ext}}=B^{\text{ext}}\bm{\hat{z}}$, with $B_{i}^{\shortparallel,\text{ext}}=B^{\text{ext}}\cos\theta_{i}$, where $\cos\theta_{i}\equiv\bm{\hat{z}}\cdot{{\bm{\hat{\textnormal{\i}}}}}$, the angle between the direction of the Ising axis at site $i$ and the direction of $\bm{B}^{\text{ext}}$.

Third, is the contribution from the self-field, $B_{i}^{\shortparallel,\text{self}}$. In the classic case of a single point dipole Griffiths (1999, 1982), a term $\frac{2}{3}\mu_{0}\mu_{\textrm{B}}\delta(\bm{r})$ must be added to ensure that the average magnetic field in a sphere containing the dipole gives the correct macroscopic field. Similarly, we add a self-field to ensure that the internal magnetic field in a uniformly magnetized sample has the expected value, for example $\bm{B}\!=\!\frac{2}{3}\mu_{0}\bm{M}$ for a uniformly magnetized sphere or cube Chen et al. (2005). Note that one should, in general, treat the limit of a uniformly magnetized non-ellipsoidal sample with some care. In this work, we are primarily concerned with paramagnetic samples in the linear response regime, where a weak magnetic field induces a magnetization proportional to it, as in a typical $\chi$ measurement. For a non-ellipsoidal sample, the induced magnetization is in general non-uniform, except in the $\chi\!\rightarrow\!0$ limit. In this limit, $H_{\text{d}}$ vanishes and, as a consequence, $H_{\text{int}}$ and $M$ are uniform. Our goal is therefore to determine the self-field so that the magnetic field has the expected value in the $\chi\!\rightarrow\!0$ limit. We demonstrate the basic idea with two examples.

We first take a cubic sample with all moments aligned in the global $\bm{\hat{z}}$ direction. In this case $\bm{B}$, $\bm{M}$ and $\bm{H}$ are all aligned with the $\bm{\hat{z}}$ direction for which the field equation $\bm{B}\!=\!\mu_{0}(\bm{M}+\bm{H})$ reduces to $B^{z}\!=\!\mu_{0}(M^{z}-N_{0}M^{z})\!=\!\frac{2}{3}\mu_{0}M^{z}$, where $N_{0}\!=\!\frac{1}{3}$ is the $\chi\!\rightarrow\!0$ limit of $N$ for a cube Chen et al. (2005). If we consider a simple cubic lattice, it is well known that the lattice sum vanishes Jackson (1975). This implies that $B^{z,\text{self}}\!=\!\frac{2}{3}\mu_{0}M^{z}$ must be incorporated to ensure the expected net $B^{\parallel}$ field value.

As a second example, we consider the case of a lattice where all the Ising axes are tilted by the same angle $\theta_{i}=\theta$ with respect to the $z$-axis, with half the spins tilted to the right and half to the left so that there is no net magnetization in the $\bm{\hat{x}}$ or $\bm{\hat{y}}$ directions. The total $\bm{B}$, $\bm{M}$, and $\bm{H}$ fields are again in the $\bm{\hat{z}}$-direction, but what should the $B^{\shortparallel}$ field parallel to the magnetic moments be? From $\bm{B}\!=\!\mu_{0}(\bm{M}+\bm{H})$, it follows that $\bm{B}$, is generated by two terms, which we discuss separately. We begin with the term generated directly by $\bm{M}$, namely $\bm{B}^{1}\!=\!\mu_{0}\bm{M}$, or $B^{1,z}\!=\!\mu_{0}M^{z}\!=\!\mu_{0}M^{\shortparallel}\cos\theta$, where $M^{\shortparallel}$ is the magnetization in the local Ising directions, $M^{\shortparallel}\!=\!V^{-1}\sum_{i=1}^{\mathcal{N}}\bm{m}_{i}\cdot{\bm{\hat{\textnormal{\i}}}}$. This equation is satisfied by $B^{1,\shortparallel}\!=\!\mu_{0}M^{\shortparallel}$. The second term, $B^{2,z}\!=\!\mu_{0}H^{z}\!=\!-\mu_{0}N_{0}M^{z}\!=\!-\mu_{0}N_{0}M^{\shortparallel}\cos\theta$ is generated by $H_{\text{d}}$. The field along the magnetic moment is thus $B^{2,\shortparallel}\!=\!-\mu_{0}N_{0}M^{\shortparallel}\cos^{2}\theta$, and the net self-field becomes

[TABLE]

which is valid when the dipolar lattice sum, Eq. (3), vanishes and when the average $\bm{M}$ is along $\bm{B}^{\text{ext}}$. For the case in which the lattice sum does not vanish, it must be subtracted from the self-field in order to ensure the expected net field value.

Eqs. 2, 3 and 4 give the local field in terms of the set of local magnetizations, $\{m_{i}\}$. With the local fields available we next consider the reverse relation that yields the $\{m_{i}\}$ induced by $B_{i}^{\shortparallel}$. Using $\bm{M}\!=\!\chi\bm{H}$ (linear media), we get $\bm{B}\!=\!\mu_{0}(\bm{M}+\bm{M}/\chi)\!=\!\mu_{0}\frac{\chi+1}{\chi}\bm{M}$, leading to

[TABLE]

where $\chi_{\text{loc}}$ is the local susceptibility in the ${{\bm{\hat{\textnormal{\i}}}}}$ direction, $M^{\shortparallel}=\chi_{\text{loc}}H^{\shortparallel}$.

We can now proceed to iterate the expressions for $B_{i}^{\shortparallel}$ in Eq. (2) and $m_{i}$ in Eq. (5) until convergence, and then calculate $N$ from Eq. (1), where $\chi_{\text{exp}}$ is given by

[TABLE]

where we are still considering site independent tilt angles, $\theta_{i}=\theta$. The intrinsic susceptibility, $\chi_{\text{int}}$, expresses the relation between $\bm{B}^{\text{ext}}$ and induced $\bm{M}$ under “Ewald”, or “tin foil”, boundary conditions Ewald (1921), which eliminate demagnetizing fields and correspond to the $N\!=\!0$ limit. As a result, both $\chi_{\text{int}}$ and $\chi_{\text{loc}}$ are responses to an internal field. While $\chi_{\text{int}}$ measures the response in the direction of $H_{\text{ext}}$, $\chi_{\text{loc}}$ measures the response along the local Ising axis ${{\bm{\hat{\textnormal{\i}}}}}$. With $\bm{\hat{z}}\cdot{\bm{\hat{\textnormal{\i}}}}=\cos\theta$, $H_{\text{ext}}^{z}\cos\theta$ induces a magnetization $M^{\shortparallel}=\chi_{\text{loc}}H_{\text{ext}}^{z}\cos\theta$. This magnetization, in turn, has a component $M^{z}=M^{\shortparallel}\cos\theta=\chi_{\text{loc}}H_{\text{ext}}^{z}\cos^{2}\theta$ along $\bm{\hat{z}}$, and therefore $\chi_{\text{int}}=\chi_{\text{loc}}\cos^{2}\theta$.

To sum up, once the converged $B_{i}^{\shortparallel}$ and $m_{i}$ distributions have been determined, $N$ is calculated using Eq. (1),

[TABLE]

In Fig. 1, we include $N$ calculated for the pyrochlore lattice using the iterative method (red line), and the main theoretical result is shown in Fig. 2, where $N$ is displayed as a function of $\chi_{\text{int}}$ for cubic samples of the simple cubic (sc) and body centered cubic (bcc) lattices with the Ising direction parallel to $\bm{B}^{\text{ext}}$ $(\cos\theta=1)$. Results (not shown) for a tetragonal lattice, relevant to LiHoF4 Gingras and Henelius (2011), are found to be identical to the sc case. We also display results for a bcc lattice with spins pointing in the $[101]$ and $[\bar{1}01]$ directions $(\cos\theta=1/\sqrt{2})$, and a pyrochlore lattice $(\cos\theta=1/sqrt3)$ built from the conventional cubic unit cell Melko and Gingras (2004). Finally, we include results for the dipolar model with spherically symmetric Heisenberg spins on an sc lattice.

IV Determination of $N$ via Monte Carlo simulations

With the iterative method, we are able to reach relatively large system sizes of $\mathcal{O}(10^{6})$ spins. To verify that this method, which is mean-field like and does not include fluctuations in the $m_{i}$’s, gives the same result as a full statistical calculation for a given spin Hamiltonian, we have also calculated $N$ using Monte Carlo (MC) simulations for several representative cases (see Fig.  2). For a single data point, the MC approach requires $\mathcal{O}(10^{5})$ core hours pro to reach the necessary precision for $\mathcal{O}(10^{4})$ moments. Since the iterative formulation contains an internal susceptibility, but no explicit temperature, $T$, it is necessary to tune either the MC $T$, or the iterative method $\chi_{\text{int}}$, so that the MC susceptibility calculated using Ewald boundary conditions, $\chi_{\text{int}}^{\text{MC}}$, matches the susceptibility from the iterative calculation. We have chosen to adjust the MC temperature, $T$, in order to tune $\chi_{\text{int}}$ to the desired value. In other words, and to emphasize, we do not compare a temperature-dependent mean-field theory calculation with a MC calculation at the same nominal temperature, a calculation which would not generally yield the same $N$ in the thermodynamic limit. For details concerning the numerical methods, we refer the reader to Appendices B - D.

For definitiveness, we use the magnetostatic dipolar Hamiltonian

[TABLE]

where $\sigma_{i}=\pm 1$, $\mu$ is the magnetic moment and $\Lambda_{ij}=\left[({\bm{\hat{\textnormal{\i}}}}\cdot{\bm{\hat{\textnormal{\j}}}})-3({\bm{\hat{\textnormal{\i}}}}\cdot\bm{\hat{r}}_{ij})({\bm{\hat{\textnormal{\j}}}}\cdot\bm{\hat{r}}_{ij})\right]/r_{ij}^{3},$ and $\chi^{zz}$, in zero field, is determined according to

[TABLE]

Using Ewald boundary conditions, we obtain $\chi_{\text{int}}^{\text{MC}}$, while open boundary conditions yield $\chi_{\text{exp}}^{\text{MC}}$, with $N$ obtained from Eq. (1).

Results for the MC method are shown in Fig. 2. All MC and iterative results have been extrapolated to infinite system size and, in Fig. 3, we compare the system-size dependence of the iterative and MC methods. Results for open boundary conditions are extrapolated using the form $a+b/L+c/L^{2}$ where the leading $1/L$ term represents a surface to volume ratio effect, while we use $a+b/L^{3}+c/L^{6}$ for periodic boundary conditions, with the leading $1/L^{3}$-term representing the inverse volume of the system. These functions yield the best fit to the data, but we find that the extrapolated value of $N$ is rather insensitive to the precise fitting function, see Appendix D.

V Discussion

The key results of this study are threefold. First, we find quantitative agreement between two theoretical methods – iterative and MC (Fig. 2) – and experiment (Fig. 1), demonstrating that our methods are sound. Second, the explicit $T$ dependence of $N$ for a cuboid has been verified for a real material (Fig. 1). Finally, $N$ is found to depend on the symmetry and direction of the moments (Fig. 2). The sc, bcc and LiHoF4 lattices with collinear Ising spins yield the same $N$, indicating that $N$ is not directly sensitive to the lattice. However, turning the local Ising axes away from $B_{\text{ext}}$ causes a more rapid decrease of $N$ with increasing $\chi_{\text{int}}$. The pyrochlore lattice with tilt angle $\cos\theta=1/\sqrt{3}$ yields a smaller $N$ than the bcc lattice with $\cos\theta=1/\sqrt{2}$ for $\chi_{\text{int}}>0$. The spin ice pyrochlore lattice and the dipolar model with Heisenberg spins yield the same result as the continuum method of Chen et al. Chen et al. (2005), and we conjecture that models with isotropic $\chi$ will generally follow this behavior Chen et al. (2005).

Exchange interactions, even when known in detail (e.g., for Dy2Ti2O7 Yavors’kii et al. (2008); Henelius et al. (2016)), have not been included in our theoretical models. This is because demagnetizing fields arise solely from the long-range dipolar interactions. The thermodynamic limit for short-range models is well-defined Levy (1968); Ruelle (1963); Campa et al. (2014), and inclusion of short-range interactions does not alter the thermodynamic limit results for $N$; see Appendix E. Thermal fluctuations also appear irrelevant in this limit. For ellipsoids, $N$ is calculated from averaged macroscopic fields that do not include thermal fluctuations and, similarly, our mean-field like iterative method captures the essential demagnetizing effects also for cuboids. However, in the non-universal approach to the thermodynamic limit (Fig. 3), there is an expected and significant finite-size difference between the iterative and the MC methods.

What are the experimental implications of our results? If an accurate measurement of $\chi_{\text{int}}$ is required, then the corrections to $N(\chi\!\rightarrow\!0)$ identified here may be dramatic for $\chi\gtrsim N$. For example, in the case of Dy2Ti2O7, $T\chi_{\textup{int}}(T)$ features a peak, which is easily shifted outside the experimental temperature window by application of the ordinary $\chi\!=\!0$ demagnetizing correction (see Fig. 1 and Ref.  Bovo et al. (2013)). More generally, while the demagnetizing correction is readily controlled for needles or ellipsoids, it is not always easy to prepare real samples with these ideal shapes. This is particularly true of non-metallic and often brittle samples – e.g., spin ice Gingras (2011) and LiHoF4 Gingras and Henelius (2011) – which have become of significant interest in recent years. Therefore, insofar as cuboidal samples are often the most practical to prepare and control, the best approach may be to use them alongside the theoretical corrections identified in this work. Our methods are general and valid for localized-moment magnets independently of details like interaction range and spin dimensionality, and the iterative method can be generalized to non-cuboids. The iterative method could also prove useful for calculating demagnetizing effects in aggregate systems, such as biomedically relevant dispersions of magnetic nanoparticles Sánchez et al. (2017).

In conclusion, considering the demagnetizing problem as a paradigm for the study of long-range interactions, our results confirm that $N$ may be defined for cuboids such that their free energy includes a term $F_{\textrm{mag}}=(\mu_{0}/2)VN(T)M^{2}$ Rhodes and Rowlands (1954) where $M$ is thermodynamically conjugate to $H_{\rm int}$. By going beyond Maxwell’s continuum theory, we show that $N$ depends not only on sample shape and $\chi$, but also on microscopic factors: the spin dimensionality and local spin anisotropy. Given that microscopic details affect even such a fundamental and well-studied macroscopic property as $N$, it is interesting to ask how they could affect the thermodynamics of more general long-range interacting systems.

Acknowledgements.

We thank D. Prabhakaran for providing crystals from which the samples were cut, and Tom Fennell and Jeffrey Rau for useful discussions. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Center for High Performance Computing (PDC) at the Royal Institute of Technology (KTH). M.T. and P.H. are supported by the Swedish Research Council (2013-03968), M.T. is grateful for funding from Stiftelsen Olle Engkvist Byggmästare (2014/807), and L.B. is supported by The Leverhulme Trust through the Early Career Fellowship program (ECF2014-284). The work at the University of Waterloo was supported by the Canada Research Chair program (M.J.P.G., Tier 1). This research was supported in part by the Perimeter Institute for Theoretical Physics. Research at the Perimeter Institute is supported by the Government of Canada through Innovation, Science, and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation, and Science.

Appendix A Susceptibility measurement

The magnetic susceptibility was measured using a Quantum Design SQUID magnetometer and the crystals were positioned in a cylindrical plastic tube to ensure a uniform magnetic environment. Measurements were performed in the RSO (reciprocating sample option) operating mode to achieve better sensitivity by eliminating low frequency noise. The position of the sample was carefully optimized to minimize misalignment with respect to the applied magnetic field. In particular, the sphere was measured at different positions and orientations in order to confirm the isotropic response and to fully reproduce the results of Bovo et al. (2013). Similarly, the cube, with edges cut along $[001]$, $[010]$ and $[001]$, was measured with the field aligned along all three orientations giving equivalent results, as would be expected.

Different measurements were made on each sample and orientation: low-field susceptibility (at $\mu_{0}H_{0}=0.0025$, $0.005$ and $0.01$ ${\rm T}$) in field-cooled (FC) versus zero-field-cooled (ZFC) protocol. In addition, magnetic field sweeps at fixed temperature were performed in order to evaluate the susceptibility accurately and confirm the low-field linear response approximation. The FC versus ZFC susceptibility measurements involved cooling the sample to $1.8$ ${\rm K}$ in zero field, applying the weak magnetic field, measuring the susceptibility while warming up to $350$ ${\rm K}$, cooling to $1.8$ ${\rm K}$ again and finally re-measuring the susceptibility while warming. Before switching the magnetic field off, field scans with small steps were performed in order to estimate the absolute susceptibilities. As expected, and previously reported Bramwell et al. (2000), no difference between field-cooled and zero-field-cooled magnetization was observed in this temperature range. The magnetization of each sample was averaged over all six repetitions (three fields, two scans each) to minimize the influence of noise.

Appendix B Iterative method

The iterative method was implemented using a form of “trivial parallelization”, in which the local field at all sites is calculated in parallel for a given magnetic moment distribution. An MPI allgather call mpi is used in order to achieve good strong scaling Kunkel et al. (2016) when run on many processors, a necessity in order to reach $\mathcal{O}(10^{6})$ spins used in this study. The number of iterative steps required to reach convergence increases with increasing susceptibility but is $\mathcal{O}(10^{2})$ regardless of the number of spins. Therefore, internode communication is not a bottleneck even though we gather and broadcast a vector equal to the length of the number of spins at every iterative step. A typical run for the largest system sizes ($2\times 10^{6}$ spins) and 1024 cores pro takes around 6 hours and requires roughly 400 communications when the intrinsic susceptibility, $\chi_{\text{int}}\sim 10$.

Appendix C Monte Carlo method

The Monte Carlo (MC) method used in this study is mostly based on the Metropolis-Hastings single-spin flip algorithm Hastings (1970) applied to Ising spins. The exception is a loop algorithm Melko et al. (2001), which we applied to the dipolar spin ice Hamiltonian in addition to the single-spin flip algorithm.

Appendix D Extrapolation to infinite system size

The approach to the infinite system size limit of the demagnetizing factor $N$ in the iterative and MC methods is illustrated in Fig. 3 in the main text. Fig. 3 was generated by selecting three MC temperatures ($16$ K, $10$ K and $3.5$ K), and calculating the susceptibilities extrapolated to infinite system size for these temperatures ($\chi_{\text{int}}^{\text{MC}}=1.00,1.82$ and $7.53$). The iterative method calculations were performed with these susceptibility values for all system sizes, and the MC $T$ was also kept the same for all system sizes.

The functional forms used for the extrapolation also deserve further comments. For the open boundary case, the leading term is of the form $1/L$, the surface to volume ratio. This is numerically confirmed in Table D.1, where the first column gives the fitting function, the second the root-mean-square error (RMSE), and the third column the extrapolated value of $N$. The MC susceptibility is calculated with open boundary conditions for $\chi_{\text{int}}=1.82$ (red squares in Fig. 3 in the main text). The smallest RMSE is found in the first and last row of Table D.1, both with a leading $1/L$ term. All data for open boundary conditions in this study have been transformed using the form $a+b/L^{1}+c/L^{2}$, marked in red in Table D.1.

In Table D.2, the corresponding data for periodic boundary conditions are shown, and we note that the RMSE and extrapolated $N$ are not very sensitive to the precise form of the extrapolation function, but the minimum RMSE is found for the function $a+b/L^{3}+c/L^{6}$, which represents an expansion in inverse volume of the surface-free system. All data for periodic boundary conditions in this study have been transformed using the form $a+b/L^{3}+c/L^{6}$, marked in red in Table D.2.

Appendix E Short-range exchange interactions

As discussed in the main text, models with short-range interactions have a well defined shape-independent thermodynamic limit and adding exchange interactions to the dipolar Hamiltonian does not alter the demagnetizing factor. We illustrate this explicitly by a numerical MC simulation of the so-called dipolar spin ice model, which has been found to reproduce a number of properties of the Dy2Ti2O7 and Ho2Ti2O7 dipolar spin ice materials den Hertog and Gingras (2000); Yavors’kii et al. (2008); Henelius et al. (2016). The Hamiltonian for this model consists of the dipolar term defined in Eq. 8 and an exchange term of the form

[TABLE]

where the strength of the dipolar interaction is given by $D=\mu_{0}\mu^{2}/4\pi r_{nn}^{3}k_{B}$ with $r_{nn}$ being the nearest-neighbor distance and $k_{B}$ the Boltzmann’s constant. The matrix $J_{ij}$ is the exchange interaction strength between particle $i$ and $j$. Here we consider first ($J_{1}$), second ($J_{2}$) and third-nearest-neighbor ($J_{3}$) exchange interactions.

In Fig. E.1 we show the demagnetizing factor for this model with parameters ($D=1.3224$ K, $J_{1}=3.41$ K, $J_{2}=-0.14$ K, and $J_{3}=0.025$ K, see Ref. Yavors’kii et al. (2008)) and the same model with no exchange interaction ($D=1.3224$ K, $J_{1}=J_{2}=J_{3}=0$ K). We expect the infinite system size susceptibility to be dependent on boundary conditions, as shown in Fig. E.2, while the difference of the inverse susceptibilities (demagnetizing factor $N$) is independent of boundary conditions, as shown in Fig. E.1. Hence, we expect no entangling between the dipolar part and the exchange part for the determination of $N$ when both are present.