Microscopic Aspects of Magnetic Lattice Demagnetizing Factors
Mikael Twengström Laura Bovo Michel J. P. Gingras Steven. T. Bramwell Patrik Henelius



 
Department of Physics, Royal Institute of Technology, SE-106 91 Stockholm, Sweden
  
 London Centre for Nanotechnology and Department of Physics and Astronomy, University College London, 17-19 Gordon Street, London, WC1H OAH, U.K.
  
  Department of Physics and  Astronomy, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada 
  Canadian Institute for Advanced Research, 180 Dundas St. W., Toronto, Ontario, M5G 1Z8, Canada
Perimeter Institute for Theoretical Physics, 31 Caroline St. N., Waterloo, Ontario, N2L 2Y5, Canada

London Centre for Nanotechnology and Department of Physics and Astronomy, University College London, 17-19 Gordon Street, London, WC1H OAH, U.K.
 
Department of Physics, Royal Institute of Technology, SE-106 91 Stockholm, Sweden



  
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 Dy_2Ti_2O_7 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.


Microscopic Aspects of Magnetic Lattice Demagnetizing Factors
    P. Henelius
    Received December 14, 2016; accepted January 23, 2017
=============================================================



§ 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^-α where α≤ d <cit.>. 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 <cit.>.  This raises the question of whether approximating other shapes to corresponding ellipsoids <cit.> 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 <cit.>. Demagnetizing effects are also important in superconductors, while analogues occur, for example, in electric systems <cit.> (depolarizing factor), in the problem of strain fields around inclusions <cit.>, and in the treatment of avalanching systems in confined geometries <cit.>.

In an applied magnetic field H_ ext = B_ ext/μ_0, the thermodynamic energy of an ellipsoid of volume V and magnetic moment m acquires a contribution E_mag=(μ_0/2) N m^2/V, where N is the demagnetizing factor. 
After subtracting E_mag from the total energy, differentiation with respect to the magnetization, M ≡ m/V, defines the internal field ≡+, where H_ d = - NM is the demagnetizing field. The intrinsic magnetic susceptibility = ∂ M/ ∂ is a shape-independent material property derived from the experimentally determined  susceptibility ≡∂ M/ ∂ through

    1/= 1/ - N .
 

The determination of N is a fundamental problem that dates back to the work of Poisson and Maxwell <cit.>. In the 1940s, Osborn <cit.> and Stoner <cit.> tabulated N for general ellipsoids, while more recently,  Aharoni <cit.> treated cuboids in the 
→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 χ_int itself <cit.>, and that ii) many experiments are routinely performed not on ellipsoids but on cuboids <cit.>, it is perhaps remarkable 
that it was only very recently that the χ-dependence of N was calculated for cuboids 
away from the →0 limit <cit.>.

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 <cit.> (through  χ_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  <cit.>. Our interest in this problem was spurred by the recent experimental observation of anomalous demagnetizing effects in the spin ice material Dy_2Ti_2O_7 <cit.>.

One may ask whether small differences
in the estimated N really matter for exposing important physics. The answer is found in Eq. (<ref>).
If ≪ 1, then  is insensitive to the precise value of N. However, in many physical systems that display unusual and interesting magnetic phenomena,  is large, and  becomes a sensitive function of N. 
Examples include the spin ice materials Dy_2Ti_2O_7 and Ho_2Ti_2O_7, which support magnetic monopole excitations <cit.>, and LiHo_1-xY_xF_4 which displays 
ultra-slow relaxation <cit.>. Important demagnetizing effects are manifest when an accurately directed field is required:
for example in experiments on the elusive Kasteleyn transition <cit.>, sub-lattice pinning <cit.> and multiple field-driven transitions <cit.>;  or else for disentangling the in- and out-of-phase frequency response <cit.>. In such cases, quantitative conclusions and accurate tests of theory depend, through χ_ 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 <ref> we discuss how to determine N experimentally. In <ref> we introduce an iterative method for obtaining N, and we consider in <ref> a Monte Carlo calculation of N. Finally, we close the paper with a discussion in <ref>. For details regarding the experimental and numerical procedures we refer the reader to Appendices <ref>-<ref>. The effects of short-range interactions are considered in some detail in Appendix <ref>.







§ 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 Dy_2Ti_2O_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 3̅ m) with no evidence of crystal distortion <cit.>, and can be accurately cut into high-quality single crystal samples of different shape. Since its spin Hamiltonian has been established in great detail <cit.>, it is convenient to adopt Dy_2Ti_2O_7 as a model system for studying the demagnetizing factor.


A sphere of diameter 4 mm and a cube of dimensions 2× 2× 2 mm^3, with edges precisely oriented along the cubic crystallographic axes [100], [010] and [001] directions, were commercially hand-cut from different larger crystals of Dy_2Ti_2O_7 provided by D. Prabhakaran <cit.> (see Ref. <cit.>). The cube was epi-polished on all sides <cit.>. Crystal shape, orientation, and experimental conditions were carefully controlled to minimize measurement errors; see Appendix <ref>.  The experimental susceptibilities of both the sphere and the cube (χ_exp^sphere, χ_exp^cube) were determined from measurements of the magnetic moment. 

Setting the demagnetizing factor of the sphere to N_ sphere = 1/3, that of the cube was determined through Eq. (<ref>), i.e., N_ cube=1/χ_exp^cube -1/χ_exp^sphere + N_ sphere. In order to match the susceptibility of the cube and sphere in the high-T limit, χ_exp^cube was shifted by about 1% (χ_exp^cube→χ_exp^cube/1.0074) before calculating N_ cube. Fig. <ref> shows how the experimental N_ cube departs significantly from the 1/3 value when χ≳ 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 <cit.>.



§ 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  equals  and calculate the induced local magnetization for an assumed . This magnetization generates a demagnetizing field that, in turn, modifies . 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  magnetic moments. As a first case, we focus on Ising moments _i=m_iμ_B, where  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, =_i ·, which is the sum of three contributions:

    = +  +  ,
 
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, 
≡_i^dip·, is given by
the familiar form <cit.>

    =
    μ_0 μ_B/4π∑_j i(3 (·_ij)(·_ij) -·/r_ij^3 )m_j.


Second, we consider an external field in the global  direction, ^ext=B^ext, with = B^extcosθ_i, where cosθ_i≡·, the angle between the direction of the Ising axis  at site i and the direction of ^ext.



Third, is the contribution from the self-field, . In the classic case of a single point dipole <cit.>, a term 2/3μ_0 μ_Bδ() 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 =2/3μ_0 for a uniformly magnetized sphere or cube <cit.>. 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 χ measurement. For a non-ellipsoidal sample, the induced magnetization is in general  non-uniform, except in the χ→0 limit. In this limit,  vanishes and, as a consequence,  and M are uniform. Our goal is therefore to determine the self-field so that the magnetic field has the expected value in the χ→0 limit. We demonstrate the basic idea with two examples. 



We first take a cubic sample with all moments aligned in the global  direction. In this case ,  and  are all aligned with the  direction for which the field equation =μ_0(+) reduces to B^z=μ_0(M^z-N_0M^z)=2/3μ_0M^z, where N_0=1/3 is the χ→0 limit of N for a cube <cit.>. If we consider a simple cubic lattice, it is well known that the lattice sum  vanishes <cit.>. This implies that B^z,self=2/3μ_0M^z must be incorporated to ensure the expected net B^∥ field value.


As a second example, we consider the case of a lattice where all the Ising axes are tilted by the same angle θ_i=θ 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  or  directions. The total , , and  fields are again in the -direction, but what should the B^ field parallel to the magnetic moments be? From =μ_0(+), it follows that , is generated by two terms, which we discuss separately. We begin with the term generated directly by , namely ^1=μ_0, or B^1,z=μ_0 M^z=μ_0M^cosθ, where M^ is the magnetization in the local Ising directions, M^=V^-1∑_i=1^_i·. This equation is satisfied by B^1,=μ_0M^. The second term, B^2,z=μ_0H^z=-μ_0N_0M^z=-μ_0N_0M^cosθ is generated by .
The field along the magnetic moment is thus B^2,=-μ_0N_0M^cos^2θ, and the net self-field becomes

    = B^1,+B^2,= μ_0μ_B/V[1- N_0cos^2θ]m_i,

 which is valid when the dipolar lattice sum,  Eq. (<ref>), vanishes
and  when the average  is along ^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. 

<ref> 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 . Using =χ (linear media), we get =μ_0(+/χ)=μ_0χ+1/χ, leading to

    m_i=V/(/+1)/μ_0μ_B,

where  is the local susceptibility in the  direction, M^= H^. 

We can now proceed to iterate the expressions for  in Eq. (<ref>) and m_i in Eq. (<ref>) until convergence, and then calculate N from Eq. (<ref>), where  is given by 

    =^zz= (∂ M^z/∂ H^z_ext)_T =μ_0μ_B/V B^ext∑_i=1^m_icosθ,

where we are still considering site independent tilt angles, θ_i=θ. The intrinsic susceptibility, , expresses the relation between ^ext and induced  under Ewald, or tin foil, boundary conditions <cit.>, which eliminate demagnetizing fields and correspond to the N = 0 limit. As a result, both  and  are responses to an internal field.  While  measures the response in the direction of ,  measures the response along the local Ising axis . With · =cosθ, H_ext^zcosθ induces a magnetization M^= H_ext^zcosθ.
This magnetization, in turn, has a component M^z=M^cosθ= H_ext^zcos^2θ along , and therefore =cos^2θ.

To sum up, once the converged  and m_i distributions have been determined, N is calculated using Eq. (<ref>),

    N =[μ_0μ_B/V B^ext∑_i=1^m_icosθ]^-1 - 1/cos^2θ.





In Fig. <ref>, we include N calculated for the pyrochlore lattice using the iterative method (red line), and the main theoretical result is shown in Fig. <ref>, where N is displayed as a function of  for cubic samples of the simple cubic (sc) and body centered cubic (bcc) lattices with the Ising direction parallel to ^ext (cosθ=1). Results (not shown) for a tetragonal lattice, relevant to LiHoF_4 <cit.>, are found to be identical to the sc case. We also display results for a bcc lattice with spins pointing in the [101] and [1̅ 0 1 ] directions (cosθ=1/√(2)), and a pyrochlore lattice (cosθ=1/sqrt3) built from the conventional cubic unit cell <cit.>. Finally, we include results for the dipolar model with spherically symmetric Heisenberg spins on an sc lattice. 



§ DETERMINATION OF N VIA MONTE CARLO SIMULATIONS


With the iterative method, we are able to reach relatively large system sizes of 𝒪(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.  <ref>). For a single data point, the MC approach requires 𝒪(10^5) core hours <cit.> to reach the necessary precision for 𝒪(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 χ_int,  so that  the MC susceptibility calculated using Ewald boundary conditions, χ_int^MC, matches the susceptibility from the iterative calculation. We have chosen to adjust the MC temperature, T, in order to tune  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 <ref> -  <ref>.



For definitiveness, we use the magnetostatic dipolar Hamiltonian

    ℋ=μ_0μ^2/4π∑_i > jΛ_ijσ_iσ_j,

where σ_i=± 1, μ is the magnetic moment and Λ_ij=[ (·)-3 (·_ij) (·_ij)]/r_ij^3, and χ^zz, in zero field, is determined according to

    χ^zz=∂ M^z/∂ H^z=μ_0μ^2/k_TV⟨(∑_i=1^σ_icosθ)^2⟩.

Using Ewald boundary conditions, we obtain , while open boundary conditions yield , with N obtained from Eq. (<ref>).


Results for the MC method are shown in  Fig. <ref>.
All MC and iterative results have been extrapolated to infinite system size and, in Fig. <ref>, 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 <ref>. 



§ DISCUSSION


The key results of this study are threefold. First, we find quantitative agreement between two theoretical methods – iterative and MC (Fig. <ref>) – and experiment (Fig. <ref>), demonstrating that our methods are sound. Second, the explicit T dependence of N for a cuboid has been verified for a real material (Fig. <ref>). Finally, N is found to depend on the symmetry and direction of the moments (Fig. <ref>). The sc, bcc and LiHoF_4 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  causes a more rapid decrease of N with increasing . The pyrochlore lattice with tilt angle cosθ=1/√(3) yields a smaller N than the bcc lattice with cosθ=1/√(2) for >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. <cit.>, and we conjecture that models with isotropic χ will generally follow this behavior <cit.>. 

Exchange interactions, even when known in detail (e.g., for Dy_2Ti_2O_7 <cit.>),  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 <cit.>, and inclusion of short-range interactions does not alter the thermodynamic limit results for N; see Appendix <ref>. 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. <ref>), 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  is required, then the corrections to N(χ→0) identified here may be dramatic for χ≳ N. For example, in  the  case of Dy_2Ti_2O_7,  T χ_ (T)
features a peak, which is easily shifted outside the experimental temperature window by application of the ordinary  χ=0 demagnetizing correction (see Fig. <ref> and Ref.  <cit.>).
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 <cit.> and LiHoF_4 <cit.> – 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 <cit.>.

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_mag = (μ_0/2) V N(T) M^2 <cit.> where M is thermodynamically conjugate to H_ int. By going beyond Maxwell's continuum theory, we show that N depends not only on sample shape and χ, 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.


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. 




 




§ 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 <cit.>. 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 μ_0H_0 = 0.0025, 0.005 and 0.01 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 K in zero field, applying the weak magnetic field, measuring the susceptibility while warming up to 350 K, cooling to 1.8 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 <cit.>, 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.




§ 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 <cit.> is used in order to achieve good strong scaling <cit.> when run on many processors, a necessity in order to reach 𝒪(10^6) spins used in this study. The number of iterative steps required to reach convergence increases with increasing susceptibility but is 𝒪(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× 10^6 spins) and 1024 cores <cit.> takes around 6 hours and requires roughly 400 communications when the intrinsic susceptibility, ∼ 10.



§ MONTE CARLO METHOD

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



§ 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 (χ_int^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 <ref>, 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 χ_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 <ref>, 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 <ref>. 

In Table <ref>, 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 <ref>. 









§ 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 Dy_2Ti_2O_7  and Ho_2Ti_2O_7  dipolar spin ice materials <cit.>. The  Hamiltonian for this model consists of the dipolar term defined in <ref> and an exchange term of the form

    ℋ_=∑_i>j J_ij· σ_iσ_j,

where the strength of the dipolar interaction is given by D=μ_0μ^2/4π r_nn^3k_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. <ref> 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. <cit.>) 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. <ref>, while the difference of the inverse susceptibilities (demagnetizing factor N) is independent of boundary conditions, as shown in Fig. <ref>. Hence, we expect no entangling between the dipolar part and the exchange part for the determination of N when both are present.



55
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Campa et al.(2009)Campa,
  Dauxois, and Ruffo]campa09
  author author A. Campa, author T. Dauxois,   and author S. Ruffo,  @noop  journal journal Phys. Rep.  volume 480, pages 57 (year
  2009)NoStop
[Osborn(1945)]osborn45
  author author J. A. Osborn, @noop  journal journal Phys.
  Rev. volume 67, pages 351 (year 1945)NoStop
[Stoner(1945)]stoner45
  author author E. C. Stoner, @noop  journal journal
  Philos. Mag. volume 36, pages 803
  (year 1945)NoStop
[Beleggia et al.(2006a)Beleggia, De Graef, and Millev]beleggia
  author author M. Beleggia, author M. De Graef,
   and author Y. T. Millev,  @noop  journal journal Philos. Mag.  volume 86, pages 2451 (year
  2006a)NoStop
[Beleggia et al.(2006b)Beleggia, De Graef, and Millev]graef
  author author M. Beleggia, author M. De Graef,
   and author Y. T. Millev,  @noop  journal journal J. Phys. D  volume 39, pages 891 (year
  2006b)NoStop
[Di Fratta(2016)]difratta16
  author author G. Di Fratta, @noop  journal journal
  Proc. R. Soc. A volume 472, pages
  0197 (year 2016)NoStop
[Aharoni(1998)]aharoni98
  author author A. Aharoni, @noop  journal journal J.
  Appl. Phys. volume 83, pages 3432
  (year 1998)NoStop
[Ponomareva et al.(2005)Ponomareva, Naumov, Kornev, Fu, and Bellaiche]ponomareva05
  author author I. Ponomareva, author I. I. Naumov, author I. Kornev,
  author H. Fu,  and author L. Bellaiche, @noop 
  journal journal Phys. Rev. B volume 72, pages 140102 (year
  2005)NoStop
[Eshelby(1957)]eshelby57
  author author J. D. Eshelby, @noop  journal journal Proc.
  R. Soc. A volume 241, pages 376
  (year 1957)NoStop
[Durin and Zapperi(2000)]durin00
  author author G. Durin and author S. Zapperi, @noop  journal journal J.
  Appl. Phys. volume 87, pages 7031
  (year 2000)NoStop
[Csikor et al.(2007)Csikor,
  Motz, Weygand, Zaiser, and  Zapperi]csikor07
  author author F. F. Csikor, author C. Motz,
  author D. Weygand, author M. Zaiser,  and author
  S. Zapperi, @noop  journal journal Science volume
  318, pages 251 (year 2007)NoStop
[Laurson et al.(2013)Laurson, Illa, Santucci, Tallakstad, Måløy, and Alava]alava13
  author author L. Laurson, author X. Illa,
  author S. Santucci, author K.-T. Tallakstad, author K. J. Måløy,  and  author M. J. Alava, @noop
   journal journal Nat. Comm. volume 4, pages 2927 (year
  2013)NoStop
[Chen et al.(1991)Chen,
  Brug, and Goldfarb]chen91
  author author D.-X. Chen, author J. A. Brug,   and author R. B. Goldfarb,  @noop  journal journal IEEE Trans.
  Magn. volume 27, pages 3601 (year 1991)NoStop
[Würschmidt(1923)]wuhr23
  author author J. Würschmidt, @noop  journal journal
  Z. Phys. volume 12, pages 128
  (year 1923)NoStop
[Stäblein and Schlechtweg(1935)]stab35
  author author F. Stäblein and author H. Schlechtweg, @noop  journal journal
  Z. Phys. volume 95, pages 630
  (year 1935)NoStop
[Quilliam et al.(2008)Quilliam, Meng, Mugford, and  Kycia]quil08
  author author J. A. Quilliam, author S. Meng,
  author C. G. A. Mugford,  and  author J. B. Kycia, @noop
   journal journal Phys. Rev. Lett. volume 101, pages 187204 (year
  2008)NoStop
[Higashinaka and Maeno(2005)]Higa_112
  author author R. Higashinaka and author Y. Maeno, @noop  journal journal Phys.
  Rev. Lett. volume 95, pages 237208
  (year 2005)NoStop
[Chen et al.(2002)Chen,
  Pardo, and Sanchez]chen02
  author author D.-X. Chen, author E. Pardo,  and  author A. Sanchez, @noop
   journal journal IEEE Trans. Magn. volume 38, pages 1742 (year
  2002)NoStop
[Chen et al.(2005)Chen,
  Pardo, and Sanchez]chen05
  author author D.-X. Chen, author E. Pardo,  and  author A. Sanchez, @noop
   journal journal IEEE Trans. Magn. volume 41, pages 2077 (year
  2005)NoStop
[Christensen and Mørup(1983)]morup83
  author author P. H. Christensen and author S. Mørup, @noop  journal journal J.
  Magn. Magn. Mater. volume 35, pages
  130 (year 1983)NoStop
[Vedmedenko et al.(2003)Vedmedenko, Oepen, and Kirschner]vedmed02
  author author E. Y. Vedmedenko, author H. P. Oepen,  and author J. Kirschner, @noop  journal journal J.
  Magn. Magn. Mater. volume 256, pages
  237 (year 2003)NoStop
[Millev et al.(2003)Millev,
  Vedmedenko, and Oepen]millev03
  author author Y. T. Millev, author E. Vedmedenko,
   and author H. P. Oepen,  @noop  journal journal J. Phys. D  volume 36, pages 2945 (year
  2003)NoStop
[Bovo et al.(2013)Bovo,
  Jaubert, Holdsworth, and Bramwell]bovo13
  author author L. Bovo, author L. D. C. Jaubert, author P. C. W. Holdsworth,  and author S. T. Bramwell, @noop  journal journal J.
  Phys.: Condens. Matter volume 25, pages
  386002 (year 2013)NoStop
[Castelnovo et al.(2008)Castelnovo, Moessner, and Sondhi]Castel08
  author author C. Castelnovo, author R. Moessner,  and author S. L. Sondhi, @noop  journal journal
  Nature volume 451, pages 42
  (year 2008)NoStop
[Biltmo and Henelius(2012)]bilt12
  author author A. Biltmo and author P. Henelius, @noop  journal journal Nat.
  Comm. volume 3, pages 857 (year 2012)NoStop
[Fennell et al.(2007)Fennell, Bramwell, McMorrow, Manuel, and Wildes]fennell07
  author author T. Fennell, author S. T. Bramwell, author D. F. McMorrow, author P. Manuel,   and author A. R. Wildes,  @noop  journal journal Nat. Phys.  volume 3, pages 566 (year
  2007)NoStop
[Ruff et al.(2005)Ruff,
  Melko, and Gingras]Ruff2005
  author author J. P. C.  Ruff, author R. G. Melko,  and author M. J. P. Gingras, @noop  journal journal Phys.
  Rev. Lett. volume 95, pages 097202
  (year 2005)NoStop
[Sato et al.(2006)Sato,
  Matsuhira, Tayama, Hiroi,
  Takagi, and Sakakibara]Sato1
  author author H. Sato, author K. Matsuhira,
  author T. Tayama, author Z. Hiroi, author
  S. Takagi,  and author
  T. Sakakibara, @noop  journal journal Journal of Physics: Condensed Matter  volume 18, pages L297 (year
  2006)NoStop
[Sato et al.(2007)Sato,
  Matsuhira, Sakakibara, Tayama, Hiroi, and Takagi]Sato2
  author author H. Sato, author K. Matsuhira,
  author T. Sakakibara, author T. Tayama, author
  Z. Hiroi,  and author
  S. Takagi, @noop  journal journal Journal of Physics: Condensed Matter  volume 19, pages 145272 (year 2007)NoStop
[Schiffer et al.(1994)Schiffer, Ramirez, Huse, and  Valentino]Schiffer1994
  author author P. Schiffer, author A. P. Ramirez, author D. A. Huse,   and author A. J. Valentino,  @noop  journal journal Phys. Rev. Lett.  volume 73, pages 2500 (year
  1994)NoStop
[Quilliam et al.(2011)Quilliam, Yaraskavitch, Dabkowska,
  Gaulin, and Kycia]Quilliam2011
  author author J. A. Quilliam, author L. R. Yaraskavitch, author H. A. Dabkowska, author B. D. Gaulin,  and author J. B. Kycia, @noop  journal journal Phys.
  Rev. B volume 83, pages 094424
  (year 2011)NoStop
[Ruminy et al.(2016)Ruminy,
  Groitl, Keller, , and Fennell]ruminy16
  author author M. Ruminy, author F. Groitl,
  author T. Keller, ,  and  author T. Fennell, @noop
   journal journal Phys. Rev. B volume 94, pages 174406 (year
  2016)NoStop
[Yavors'kii et al.(2008)Yavors'kii, Fennell, Gingras, and  Bramwell]Yavo08
  author author T. Yavors'kii, author T. Fennell,
  author M. J. P. Gingras,  and  author S. T. Bramwell,  @noop  journal journal Phys. Rev. Lett.  volume 101, pages 037204 (year 2008)NoStop
[Henelius et al.(2016)Henelius, Lin, Enjalran, Hao, Rau, Altosaar, Flicker, Yavors'kii, and Gingras]hene16
  author author P. Henelius, author T. Lin,
  author M. Enjalran, author Z. Hao, author
  J. G. Rau, author J. Altosaar, author F. Flicker, author T. Yavors'kii,  and author M. J. P. Gingras, @noop  journal
  journal Phys. Rev. B volume 93,  pages 024402 (year 2016)NoStop
[Prabhakaran and Boothroyd(2011)]Prabhak
  author author D. Prabhakaran and author A. Boothroyd, @noop  journal journal
  Crys. Growth volume 318, pages 1053
  (year 2011)NoStop
[epi()]epi
  @noop  note SurfaceNet GmbH,
  http://www.surfacenet.deNoStop
[Griffiths(1999)]griffiths
  author author D. J. Griffiths, @noop  title Introduction to
  Electrodynamics, edition 3rd ed. (publisher
  Prentice-Hall, year 1999)NoStop
[Griffiths(1982)]GriffPaper82
  author author D. J. Griffiths, @noop  journal journal
  American Journal of Physics volume 50, pages 698 (year 1982)NoStop
[Jackson(1975)]jack
  author author J. D. Jackson, @noop  title Classical
  Electrodynamics, edition 2nd ed. (publisher
  John Wiley & Sons, year 1975)NoStop
[Ewald(1921)]ewal21
  author author P. P. Ewald, @noop  journal journal Ann.
  Phys. volume 369, pages 253 
  (year 1921)NoStop
[Gingras and Henelius(2011)]gingras11
  author author M. J. P.  Gingras and author P. Henelius, @noop  journal
  journal J. Phys.: Conf. Ser. volume
  320, pages 012001 (year 2011)NoStop
[Melko and Gingras(2004)]melko04
  author author R. G. Melko and author M. J. P. Gingras, @noop  journal journal J.
  Phys.: Condens. Matter volume 16, pages
  R1277 (year 2004)NoStop
[pro()]process
  @noop  note Intel Xeon E5-2698v3 processors with Cray
  Aries interconnect between nodesNoStop
[Levy(1968)]levy68
  author author P. Levy, @noop  journal journal Phys.
  Rev. volume 170, pages 595 (year 1968)NoStop
[Ruelle(1963)]ruelle63
  author author D. Ruelle, @noop  journal journal Helv.
  Phys. Acta volume 36, pages 183
  (year 1963)NoStop
[Campa et al.(2014)Campa,
  Dauxios, Fanelli, and Ruffo]campa14
  author author A. Campa, author T. Dauxios,
  author D. Fanelli,  and  author S. Ruffo, @noop 
  title Physics of long-range interacting systems,  edition 1st ed. (publisher Oxford University
  Press, year 2014)NoStop
[Gingras(2011)]gingras11a
  author author M. J. P.  Gingras, title Introduction to
  frustrated magnetism,  (publisher Springer (Ed. Lacroix,
  Mendels & Mila), year 2011) Chap. chapter Spin
  IceNoStop
[Sánchez et al.(2017)Sánchez, Mendoza Zélis, Arciniegas,
  Pasquevich, and Fernández van
  Raap]sanchez17
  author author F. H. Sánchez, author P. Mendoza Zélis, author M. L. Arciniegas, author G. A. Pasquevich,  and author M. B. Fernández van Raap, @noop  journal journal Phys. Rev. B volume 95, pages 134421 (year 2017)NoStop
[Rhodes and Rowlands(1954)]rhodes54
  author author P. Rhodes and author G. Rowlands, @noop  journal journal
  Proc. Leeds Phil. Lit. Soc. volume 6, pages 191 (year 1954)NoStop
[Bramwell et al.(2000)Bramwell, Field, Harris, and  Parkin]Bramwell00
  author author S. T. Bramwell, author M. N. Field,
  author M. J. Harris,  and  author I. P. Parkin, http://stacks.iop.org/0953-8984/12/i=4/a=308 journal
  journal Journal of Physics: Condensed Matter volume 12, pages 483 (year
  2000)NoStop
[mpi()]mpi
  @noop  note Open MPI v2.0.0 documentation.
  https://www.open-mpi.org/doc/v2.0/NoStop
[Kunkel et al.(2016)Kunkel,
  Balaji, and Dongarra]scaling
  author author J. Kunkel, author P. Balaji,   and author J. Dongarra, https://books.google.se/books?id=kRJkDAAAQBAJ title
  High performance computing: 31st international conference,  (year 2016), note ISC High Performance 2016, Frankfurt,
  Germany, June 19-23, 2016, ProceedingsNoStop
[Hastings(1970)]Hast
  author author W. K. Hastings, @noop  journal journal
  Biometrika volume 57, pages 97
  (year 1970)NoStop
[Melko et al.(2001)Melko,
  den Hertog, and Gingras]loopalgo
  author author R. G. Melko, author B. C. den
  Hertog,  and author M. J. P. Gingras, 10.1103/PhysRevLett.87.067203 journal journal Phys. Rev. Lett. volume 87, pages 067203 (year
  2001)NoStop
[den Hertog and Gingras(2000)]Hertog2000
  author author B. C. den Hertog and author M. J. P. Gingras, @noop  journal journal Phys.
  Rev. Lett. volume 84, pages 3430
  (year 2000)NoStop

