Screening and the Pinch Point Paradox in Spin Ice

TL;DR

This paper investigates the nature of pinch point singularities in spin ice, revealing a paradox between simulation/theory and experimental observations, and exploring potential explanations involving screening, disorder, or quantum effects.

Contribution

It provides high-resolution simulations and new analytical theory of pinch points in dipolar spin ice, highlighting a paradox between theoretical predictions and experimental results.

Findings

Simulations and existing theories predict infinitely sharp pinch points.

Experiments show broadened pinch points, indicating screening effects.

The paradox suggests possible roles for quantum fluctuations or fundamental theory corrections.

Abstract

A pinch point singularity in the structure factor characterizes an important class of condensed matter that is a counterpoint to the paradigm of broken symmetry. This class includes water ice, charge ice and spin ice. Of these, dipolar spin ice affords the the pre-eminent model system because it has a well-established Hamiltonian, is simple enough to allow analytical theory and numerical simulation, and is well represented in experiment by Dy$_{2}$Ti$_{2}$O$_{7}$ and Ho$_{2}$Ti$_{2}$O$_{7}$. Nevertheless it is a considerable challenge to resolve the pinch points in simulation or experiment as they represent a very long range correlation. Here we present very high resolution simulations of the polarized neutron scattering structure factor of dipolar spin ice and new analytical theory of the pinch point profiles. We compare these with existing theory and experiment. We find that our simulations are consistent with theories that reveal the pinch points to be infinitely sharp, as a result of unscreened dipolar fields. However, neither simulation nor these theories are consistent with experiments, which instead is quantitatively captured by a theory that allows for screening of the dipolar fields and consequent strong broadening of the pinch points. This striking paradox is not easily resolved: broadening of the pinch points by random disorder seems to have been ruled out by existing theory, while deficiencies in the Hamiltonian description are not relevant. Intriguingly, we are left to consider the role of quantum fluctuations or the possibility of a fundamental correction to either the standard method of simulating dipolar systems, or the theory of polarized neutron scattering. More generally, our results may have relevance far beyond ice systems. For example, spin ice is a model Debye-H\"uckel (magnetic) electrolyte, so our basic observation that the screening length may diverge while...

Equations

$$\mathcal{H}={\mathcal{H}}_{\mathrm{exchange}}+{\mathcal{H}}_{\mathrm{dipolar}}.$$

$$\frac{\text{d}{\sigma }^{\alpha \beta }}{\text{d}\Omega }=A\left({\delta }_{\alpha \beta }-{\hat{Q}}_{\alpha }{\hat{Q}}_{\beta }\right)S^{\alpha \beta }(\mathbf{Q}),$$

$$\mathcal{H}=J_1\sum _{⟨i,j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+D{a}^3\sum _{i>j}\frac{{\mathbf{S}}_i\cdot {\mathbf{S}}_j-3\left({\hat{\mathbf{r}}}_{ij}\cdot {\mathbf{S}}_i\right)\left({\hat{\mathbf{r}}}_{ij}\cdot {\mathbf{S}}_j\right)}{r_{ij}^3}.$$

$$\mathbf{M}=g{\mu }_{\mathrm{B}}{v}^{-1}\sum _{i=1}^4{\mathbf{s}}_i,$$

$${\mathbf{s}}_1$$

$${\mathbf{s}}_2$$

$${\mathbf{s}}_3$$

$${\mathbf{s}}_4$$

$${\mathbf{r}}_1$$

$${\mathbf{r}}_2$$

$${\mathbf{r}}_3$$

$${\mathbf{r}}_4$$

$$\frac{\text{d}{\sigma }_{hh2}^{\mathrm{NSF}}}{\text{d}\Omega }$$

$$\frac{\text{d}{\sigma }_{hh2}^{\mathrm{SF}}}{\text{d}\Omega }$$

$$G[\mathbf{M}]={\mu }_0\int \left(\frac{{\mathbf{M}}^2}{2\chi }+\frac{{\xi }^2(\nabla \cdot \mathbf{M}{)}^2}{2\chi }-\mathbf{H}\cdot \mathbf{M}\right){\text{d}}^{3}r,$$

$$S^{\mathrm{T}}(q)=\chi T/3C{S}^{\mathrm{L}}(q)=\frac{\chi T/3C}{1+{\xi }^2{q}^2}.$$

$$\gamma (T)=\chi T/C=\frac{2(1+e^{2J_{\mathrm{eff}}/T})}{2+e^{2J_{\mathrm{eff}}/T}+e^{-6J_{\mathrm{eff}}/T}},$$

$$n=\frac{4e^{-2J_{\mathrm{eff}}/T}}{3+4e^{-2J_{\mathrm{eff}}/T}},$$

$${\xi }^2=\frac{\chi kT{v}_0}{{\mu }_0{Q}^{2}n}.$$

$$U_{\mathrm{mag}}=-\frac{{\mu }_0}{2}\int \mathbf{M}\cdot \tilde{\mathbf{H}}{\text{d}}^{3}r,$$

$$\mathbf{M}-{\xi }^2\nabla \left(\nabla \cdot \mathbf{M}\right)=\chi (\mathbf{H}+\tilde{\mathbf{H}}),$$

$$\tilde{\mathbf{H}}(\mathbf{r})=-\frac{{\mu }_0}{4\pi }\nabla \int \frac{\rho ({\mathbf{r}}^{\prime })}{|{\mathbf{r}-\mathbf{r}}^{\prime }|}{\text{d}}^3{r}^{\prime }.$$

$$S^{\mathrm{T}}=\gamma /3S^{\mathrm{L}}=\frac{\gamma /3}{(1+\chi )+q^2{\xi }^2}.$$

$$\nabla \cdot \nabla \varphi (\mathbf{r})+{\kappa }^2\varphi (\mathbf{r})=0,$$

$$-\nabla \left(\nabla \cdot \tilde{\mathbf{H}}(\mathbf{r})\right)-{\kappa }^2\tilde{\mathbf{H}}(\mathbf{r})=\left[\nabla \left(\nabla \cdot \mathbf{M}\right)+{\kappa }^2{N}_{\mathrm{d}}\mathbf{M}(\mathbf{r})\right]\delta (\mathbf{r}),$$

$$-\nabla \left(\nabla \cdot \tilde{\mathbf{H}}(\mathbf{r})\right)-{\kappa }^2\tilde{\mathbf{H}}(\mathbf{r})=$$

$$\int \left[\nabla \left(\nabla \cdot \mathbf{M}({\mathbf{r}}^{\prime })\right)+{\kappa }^{2}Nd\mathbf{M}({\mathbf{r}}^{\prime })\right]\delta ({\mathbf{r}-\mathbf{r}}^{\prime }){\text{d}}^3{r}^{\prime }.$$

$$\tilde{\mathbf{H}}(\mathbf{q})=-\left(\mathbf{q}\left(\mathbf{q}\cdot \mathbf{M}(\mathbf{q})\right)+{\kappa }^2{N}_{\mathrm{d}}\mathbf{M}(\mathbf{q})\right)/({\kappa }^2+q^2),$$

$$S^{\mathrm{L}}=\frac{(\gamma /3)\left({\kappa }^2+q^2\right)}{{\kappa }^2(\chi N_{\mathrm{d}}+1)+q^2(1+\chi +{\kappa }^2{\xi }^2)+q^4{\xi }^2}$$

$$S^{\mathrm{T}}=\frac{(\gamma /3)\left({\kappa }^2+q^2\right)}{\chi {\kappa }^2{N}_{\mathrm{d}}+({\kappa }^2+q^2)}.$$

Full text

I Introduction

The decomposition of the scattering amplitude into the product of form factor and structure factor is a key device in the theory of scattering from condensed matter. While the form factor represents the scattering amplitude of an object, the structure factor represents the correlation between objects. Arrays of dissimilar objects may exhibit similar structure factors because their correlations are similar. Hence the classification of structure factors plays a key role in classifying the types of physical state that may occur in condensed matter systems and beyond.

The various classes of structure factor $S^{\alpha\beta}({\bf q})$ may be elucidated by considering possible singularities in the function. The pinch point Fig. 1 a is a type of singularity that is expected, and may be observed, in at least two types of system: firstly, dipolar systems, such as ferromagnets Krivoglaz (1964); Als-Nielsen et al. (1976), and secondly, ice-rule systems, including hydrogen bonded ferroelectrics Paul et al. (1970); Havlin et al. (1982); Youngblood and Axe (1978); Youngblood et al. (1980), water ice Nield and Whitworth (1995); Wehinger et al. (2014), spin ice Moessner and Sondhi (2003); Fennell et al. (2007, 2009); Chang et al. (2010); Sen et al. (2013); Ryzhkin et al. (2013), ionic ice Anderson (1956); Fennell et al. (2019), artificial spin ice Perrin et al. (2016); Östman et al. (2018); Farhan et al. (2019), quantum spin ice Benton et al. (2012); Fennell et al. (2012); Sibille et al. (2018) and antiferromagnetic spin liquids Zinkin et al. (1997); Conlon and Chalker (2009); Ballou et al. (1996); Canals and Lacroix (2000); Henley (2005). Such systems have the notable feature that they enter a highly correlated low-temperature state without any symmetry breaking. As such, they provide a valuable contrast to the usual Landau paradigm of condensed matter, where formation of low temperature states is invariably accompanied by symmetry breaking.

The pinch point structure factor is the principal characteristic and diagnostic of such states (in some, but not all, of the above examples, Pauling entropy Pauling (1935); Ramirez et al. (1999); Giblin et al. (2018) is another key characteristic). However, despite the paradigm nature and widespread relevance of the phenomenon, there have been very few direct confrontations of theory and experiment for the pinch point structure factors. The reason seems to be that any approach to the singularity is difficult to achieve in both numerical simulation and experiment: the correlation that gives rise to the pinch point is very subtle and long ranged, meaning that the amplitude of $S^{\alpha\beta}({\bf q})$ is everywhere very small and high resolution is required to distinguish its fine features as ${\bf q\rightarrow 0}$. In addition, well-defined model systems are very scarce. Spin ice, whose pinch points arise from the combined effects of dipolar and ice-rule correlations Bramwell and Gingras (2001), stands out as the exception to this particular rule: it is a nearly ideal model system which lends itself well to advanced experiment and numerical modelling.

The idea of spin ice (see Fig. 2) was invented Harris et al. (1997); Bramwell and Harris (1998) largely to describe the unusual correlations and absence of long range order observed in neutron scattering experiments on ${\rm Ho_{2}Ti_{2}O_{7}}$. Subsequently, a very deep and consistent understanding of the many interesting properties of spin ice has been developed in terms of the microscopic dipolar spin ice (DSI) model den Hertog and Gingras (2000), with spin Hamiltonian consisting of exchange and dipolar terms:

[TABLE]

A truncation of the dipole-dipole interaction to near neighbor in Eq. 1 gives a near neighbor ferromagnetic coupling and the original near–neighbor spin ice (NNSI) model of Refs. Harris et al. (1997); Bramwell and Harris (1998), which remains a valuable point of reference. Today, DSI (along with its extensions Yavors’kii et al. (2008)) is fully accepted as the basic microscopic description of spin ice, one that leads eventually to some surprising properties, including its residual (Pauling) entropy Ramirez et al. (1999), emergent electromagnetism Ryzhkin (2005); Castelnovo et al. (2008), fractionalization Castelnovo et al. (2008) and fragmentation Brooks-Bartlett et al. (2014); Petit et al. (2016). Yet, despite this progress, the neutron scattering structure factor that first inspired the discovery of spin ice, is not yet fully understood, as we discuss in this paper. In particular, we find that the pinch point profile, that summarizes the long ranged correlations in spin ice, presents a clear contradiction between theory and experiment.

Motivations for studying the pinch points in classical spin ices like Ho2Ti2O7 and Dy2Ti2O7 include not only their above-mentioned status as model representatives of a large class of condensed matter, but also their more recently-established status as model Coulomb fluids Castelnovo et al. (2008, 2011); Kaiser et al. (2018, 2015); Bovo et al. (2013a); Paulsen et al. (2016). Here they afford perhaps unique examples of a model lattice Coulomb fluid (of magnetic monopoles) that quantitatively obeys Debye-Hückel theory in both the dilute and concentrated regimes Kaiser et al. (2018). Neutron scattering affords a rare opportunity to directly image the field correlations and screening in such a model experimental Coulomb fluid. To do this may shed some light on foundational issues of the theory of Coulomb fluids, for example the highly topical problem of ‘underscreening’ in concentrated ionic liquids, where the experimentally measured screening length appears to diverge, while the Debye length remains finite Lee et al. (2015, 2017).

Structure factor of spin ice – A complete characterization of the static correlations in spin ice can only be achieved using polarized neutron scattering Fennell et al. (2009). The reason is that the dipole-dipole interaction breaks the rotational invariance of the correlation function, such that the structure factor tensor $S^{\alpha\beta}({\bf q})$ is characterized by two numbers, corresponding to the longitudinal (to ${\bf q}$) and transverse eigenvalues, $S^{\rm L}$ and $S^{\rm T}$ respectively. Polarized neutrons are required to separate these numbers, which in unpolarized neutron scattering appear only in combination.

Here we define ${\bf Q}={\bf q}+{\bf G}$, where ${\bf Q}$ is the scattering vector and ${\bf G}$ is a reciprocal lattice vector (Fig. 1c). The neutron scattering differential cross section measures the transverse projection of the structure factor on the scattering vector Marshall and Lowde (1968):

[TABLE]

(we set $A=1$ henceforth). Within an unimportant caveat discussed below, the structure factor is the same in each Brillouin zone (hence characterized as $S^{\alpha\beta}({\bf q})$), but its projection may vary from zone to zone, as seen in Fig. 1c.

The origin of the pinch points may be visualized by representing the inverse $S^{\alpha\beta}({\bf q})$ tensor as an ellipsoid of dimensions $\sqrt{S^{\rm L}({\bf q})}\times\left(\sqrt{S^{\rm T}({\bf q})}\right)^{2}$, oriented with principal (L) axis parallel to ${\bf q}$ Bramwell (2017); Nye (1985). In a dipolar ferromagnet (say) the ellipsoid is eccentric (i.e. $S^{\rm L}\neq S^{\rm T}$) only at small $q$ Als-Nielsen et al. (1976), but in spin ice, the conspiracy of ice rules and very strong dipolar forces impose a non-spherical structure factor tensor over nearly the entire Brillouin zone. Bramwell (2017). Considering the in plane projection of this ellipsoid on ${\bf Q}$ (Fig. 1c), the pattern shown in Fig. 1b is easily inferred. If the ellipsoid remains eccentric as $q\rightarrow 0^{+}$ then the pinch points are infinitely sharp; if it evolves to a sphere in that limit they are broadened (Fig. 1c). Precisely at $q=0$, in either case, the structure factor ellipsoid is a sphere, in accord with the crystal symmetry.

We focus on the case of ${\bf G}=002$ in the $1\bar{1}0$ scattering plane used in the experiments Fennell et al. (2009); Chang et al. (2010). If we inspect wavevectors ${\bf Q}=hh(l=2)$ within the $002$ zone, then $h\bar{h}0$ and $hh0$ define two principal axes of the tensor (Fig. 1d). Non-spin flip (NSF) scattering then measures the out of plane, transverse (to ${\bf q}$) eigenvalue, while spin flip (SF) scattering measures a mixture of longitudinal and transverse eigenvalues in plane: with the projection perpendicular to ${\bf Q}$, this mixture creates the pinch points in the SF channel (Fig. 1b), while the NSF scattering, as a pure eigenvalue, remains a periodic function. However, it turns out (see below) that a SF scan along $hh2$ picks out only the longitudinal eigenvalue up to about $h=0.5$, so this particular scan ‘across’ the pinch point isolates the two eigenvalues needed to characterize the correlations in spin ice: $S^{\rm L}$ in the SF channel and $S^{\rm T}$ in the NSF channel. Thus, except very near the zone boundary, the entire correlation function of spin ice is contained in the two channels (SF,NSF) of the single line scan across the $002$ zone centre in the $1\bar{1}0$ scattering plane.

The transverse and longitudinal structure factors represent the correlations of the Helmholtz decomposed components of the magnetization Bramwell (2017). In the low temperature limit the inverse structure factor tensor becomes a disc ($S^{\rm L}=0$), reflecting a divergence free state, the so-called Coulomb phase Henley (2005). In this special case there is a singularity in the correlation function as one tensor eigenvalue, $S^{\rm L}$, disappears. The more general case of an oblate spheroid that is invariant with ${\bf q}$ has been termed the ‘harmonic phase’ Bramwell (2017) as it is implied by Onsager’s theory Onsager (1936) of dipolar fluids, which contain only harmonic (divergence and curl free) fields. In this context, we ask the question, what happens at finite temperature? As $q\rightarrow 0^{+}$, does the structure factor tensor of spin ice becomes a sphere (hence broadened pinch points) or does it become an oblate spheroid (hence infinitely sharp pinch points)? These are hard questions to answer as the long range correlations that give rise to the pinch points make meaningful simulations and experiments quite challenging. Nevertheless, we discuss these questions with respect to new, high resolution, numerical simulations of NNSI and DSI, as well as with respect to existing experimental data Fennell et al. (2009); Chang et al. (2010) and analytical theories Sen et al. (2013); Ryzhkin et al. (2013).

For completeness we note the subtle caveat mentioned above – that happily proves to be unimportant. This arises from the glide symmetry that emerges in the Fd$\bar{3}$m space group of spin ice. With respect to the face centred cubic reciprocal lattice $\{{\bf G}\}$ the function $S^{\alpha\beta}{(\bf Q})$ is modulated, such that unit cell in reciprocal space is doubled. This means that every second zone is inequivalent: for example the function is different in the zones defined by ${\bf G}=002,004$. We discuss this point in detail elsewhere Twengström et al. (2019) where we refer to zones like $002$ as ‘false zones’ and zones like $004$ as ‘true’ zones. The fact that there are two types of zone is the only reason that experiment can clearly resolve the pinch points in spin ice: at the true zone centres the nuclear Bragg peak has its full intensity and this will always obscure the subtle magnetic diffuse scattering, even in a polarized experiment. In contrast, at a false zone centre (e.g. $002$), the nuclear intensity is zero, making the magnetic diffuse scattering accessible to precise measurement. However this advantage potentially comes at a price: at a false centre $S^{\alpha\beta}({\bf Q})$ is not necessarily a measure of magnetization fluctuations, the quantity of most immediate interest. Yet there is a fortunate crystallographic coincidence that renders this point irrelevant to the present study: scans across the pinch point at $002$ (surprisingly) do yield the magnetization fluctuations as required (see Section III).

Plan of the paper – The paper is organized as follows. In section II we describe details of our numerical simulations. In Section III we present an analytical theory of the correlations that allows for different scenarios – particularly as regards screening of the dipolar fields – and we compare these with our simulations and other analytical theory Sen et al. (2013); Ryzhkin et al. (2013). We find that this theory accounts very accurately for the simulated pinch point profiles in the absence of field screening. In section IV we compare these results with experiment Fennell et al. (2009); Chang et al. (2010) and expose a very clear disagreement between experiment and simulation. However, we find that we can capture the experimental data by allowing screening of dipolar fields in the theory. This defines the pinch point paradox. In the final section, V, we summarize and discuss our main results, including possible resolutions of the paradox.

II Numerical Simulations

In the spin ice materials ${\rm R_{2}Ti_{2}O_{7}}$, rare earth ions R = Ho or Dy are located on the points of a cubic pyrochlore lattice of corner-sharing tetrahedra Subramanian et al. (1983) (see Fig. 2). The trigonal crystal field enforces a doublet ground state for each ion Rosenkranz et al. (2000); Bertin et al. (2012); Jana et al. (2002); Rau and Gingras (2015a) and establishes a local Ising-like Rau and Gingras (2015a) confinement with effective two-state spins ${\bf S}_{i}$ pointing between the centers of each pair of adjacent tetrahedra. The associated magnetic moments are very large, $\mu\approx 10~{}\mu_{\rm B}$, with the consequence that dipole-dipole interactions are particularly strong in these materials. For modelling purposes we consider classical spins of unit length. The dipolar spin ice model (DSM) combines the long-ranged dipolar interaction with short-ranged exchange terms den Hertog and Gingras (2000):

[TABLE]

Here $D$ is the dipolar interaction constant, $r_{ij}$ the distance between spin $i$ and $j$ and $J_{1}$ the nearest-neighbor exchange interaction Bramwell and Gingras (2001). Considering only nearest-neighbor interactions, the model reduces to the nearest-neighbor spin ice model (NNSI), which gives an accurate approximation to DSI down to about 0.6 K Melko and Gingras (2004). The NNSI model has effective ferromagnetic exchange parameter $J_{\rm eff}$, which takes the values 1.1 K for ${\rm Dy_{2}Ti_{2}O_{7}}$ and 1.9 K for ${\rm Ho_{2}Ti_{2}O_{7}}$ Bramwell and Gingras (2001). The NNSI model maps exactly Bramwell and Harris (1998) to Pauling’s model of (cubic) water ice, hence exhibiting a degenerate ground state with strong correlations, but no long range order (see Fig. 2). By introducing the dipolar interaction, the degeneracy of the NNSI model is weakly broken and an ordering transition at very low temperature is induced. However, this is not relevant at the temperatures ($T>1$ K) considered in this paper, where DSI behaves qualitatively like NNSI. In addition to the nearest-neighbor exchange interactions $J_{1}$, the generalized spin ice model (g–DSM) contains second and third nearest-neighbor interactions $J_{2},J_{3a}$ and $J_{3b}$. A set of parameter values were previously determined for ${\rm Dy_{2}Ti_{2}O_{7}}$ ($J_{1}=3.41$ K, $J_{2}=-0.14$ K, $J_{3a}=J_{3b}=0.025$ K) which models unpolarized neutron scattering and bulk thermodynamic properties at a quantitative level Yavors’kii et al. (2008).

We previously reported that demagnetizing effects are a pure outcome of dipolar interactions in highly correlated systems Twengström et al. (2017), i.e. the exchange interactions do not alter the shape-dependent physics at the zone centre. Nevertheless the diffuse scattering elsewhere in the Brillouin zone may depend upon exchange, as discussed subsequently. To model the neutron structure factor in more detail we reconstruct the Fourier transform of the spin-spin correlation function in $\mathbf{Q}$-space. Following the use of a parallel Monte Carlo code Giblin et al. (2018) that exploits the symmetry of the dipolar interactions Ewald (1921), we use periodic Ewald boundary conditions Ewald (1921) and a loop algorithm Melko et al. (2001) to speed up equilibration when needed at low temperature.

The surface dependence is treated via spherical boundary conditions by considering the addition of a microscopic term to the Ewald sum Twengström et al. (2019); Twengström (2018). The choice of this boundary term is the foundation of the conditional convergence of the dipolar sum and exploits a way to manage the demagnetizing effects.

High resolution is key to this study. The allowed $\mathbf{Q}$-points were determined by considering only the set of points in which the discrete Fourier transform is defined, i.e. the inverse cubical system size, $1/L$. For example, the high resolution of $L=16$ which we reach corresponds to $65536$ particles and a resolution of $1/L=0.0625$ in units of $2\pi$. Considering that the dipolar systems under study occupy an ordered lattice, parallel computing is essential, given that a system of $65536$ particles corresponds to almost 4.3 billion interactions per Monte Carlo step. Note that we do not employ any interpolation schemes such as the Nyquist–Shannon sampling theorem as this would defeat our purpose of establishing the actual line shapes and their comparison with experiment.

III Analytical theory of $S^{\alpha\beta}({\bf Q})$

In this section we derive closed-form analytical expressions for the eigenvalues of the structure factor tensor $S^{\alpha\beta}({\bf Q})$ and compare them with our simulations. For dipolar spin ice, the scattering functions we derive turn out to be consistent with those of Refs. Sen et al. (2013); Ryzhkin et al. (2013) but in addition, we have quantified all parameters and controlled boundary effects. Further, our derivation makes a connection with Ref. Bramwell (2017) where the structure factor tensor of spin ice was derived from a very different viewpoint: that of the Onsager dipolar fluid Onsager (1936).

The theory presented below assumes a continuum magnetization ${\bf M}({\bf r})$. This approximation can only really be valid with respect to magnetization defined on a Bravais lattice, such that a primitive, unit cell can be defined, which, on repetition just fills all space. The pyrochlore lattice (Fig. 2) is a face centred cubic ($F$) Bravais lattice with a four point tetrahedral basis. Hence ${\bf M}$ should be interpreted as the magnetization averaged over the four $\langle 111\rangle$ spins of the spin ice basis:

[TABLE]

where $g$ is the g-factor, $v$ is the primitive unit cell volume and the spins

[TABLE]

are located at

[TABLE]

respectively in the conventional cubic unit cell. It is straightforward to demonstrate that with respect to the axes $x=[1\bar{1}0]$, $y=[110]$, $z=[001]$, the structure factors probed in the NSF and SF channels, $S^{xx}$ and $S^{yy}$, measure the fluctuations $\langle M^{x}({\bf q})M^{x}(-{\bf q})\rangle$ and $\langle M^{y}({\bf q})M^{y}(-{\bf q})\rangle$. Because $002$ is a ‘false’ zone centre (see above) this is not true of the third $(zz)$ component. The difference arises because $zz$ structure factor samples fluctuations of four spins per tetrahedron in the combination $\sqrt{1/3}(s_{1}^{z}+s_{2}^{z}-s_{3}^{z}-s_{4}^{z})$, which is not the same as $M^{z}$, while the $xx,yy$ structure factors sample two spin combinations e.g. $\sqrt{2/3}(s_{1}^{x}+s_{2}^{x})$ and $\sqrt{2/3}(s_{3}^{y}+s_{4}^{y})$, which are the same as the corresponding projections of $M^{x},M^{y}$ (this is because two of the spins in each case have zero component on the given axes $x,y$). Below, we focus on a ‘cut’ through the ${\bf G}=002$ zone centre such that ${\bf q}=hh0$ in reduced units and the observed differential cross sections are

[TABLE]

where $\theta$ is the angle between ${\bf Q}={\bf G}+{\bf q}$ and ${\bf G}=002$. Fortunately, the phase factors of the $zz$ component act to suppress the second term in Eq. 7 such that it may be neglected up to $h\approx 0.5$. The net result is that the SF and NSF channels measure the longitudinal (L) and transverse (T) eigenvalues of the magnetization structure factor to an excellent approximation. We can therefore conveniently compare $\textup{d}\sigma^{\rm NSF}_{hh2}/\textup{d}\Omega$ with $S^{\rm T}_{hh0}$ and $\textup{d}\sigma^{\rm SF}_{hh2}/\textup{d}\Omega$ with $S^{\rm L}_{hh0}\cos^{2}\theta=S^{\rm L}_{hh0}\times 4/(4+2h^{2})$ in this range.

Near Neighbor Spin Ice – The scattering function of near neighbor spin ice (NNSI) may be calculated from the free energy functional Bramwell (2012):

[TABLE]

which comes from an expansion of the free energy in powers of magnetization ${\bf M}$. Its first term in the susceptibility $\chi$ and its third term in the field ${\bf H}$ are common to all classical spin systems, while the second term in $(\nabla\cdot{\bf M})^{2}$ expresses the free energy cost of divergence and is the lowest order gradient term for ice systems Youngblood and Axe (1978); Youngblood et al. (1980). The Euler-Lagrange equation ${\bf M}-\xi^{2}\nabla\left(\nabla\cdot{\bf M}\right)=\chi{\bf H},$ may be solved in Fourier components to give the wave vector dependent susceptibility $\chi({\bf q})$ and, by the classical fluctuation dissipation theorem, the eigenvalues of the scattering tensor:

[TABLE]

Here T/L indicate transverse/longitudinal to the wavevector ${\bf q}$. These depend on the susceptibility $\chi$ and ‘diffusion’ length $\xi$ Ryzhkin et al. (2013).

The susceptibility may be expressed as $\chi=\gamma(T)C/T$, with $C=\mu_{0}m^{2}/(3k_{\textup{B}}v_{0}/2)$, where $m$ is the rare earth moment and $v_{0}=v/2$ is the volume per tetrahedron. The moment is related to the monopole charge $Q$ by $Q=2m/a$ where $a$ is the diamond lattice constant Castelnovo et al. (2008). The factor $\gamma(T)$ has been calculated by Jaubert et al. Jaubert et al. (2013) in a Husimi tree approximation:

[TABLE]

where $J_{\rm eff}$ is the effective exchange defined above Bramwell and Gingras (2001).

The diffusion length $\xi(T)$ depends on the densities of single and double charge monopoles (which have no Coulomb interaction in NNSI). We make a low temperature approximation by setting the chemical potential of the double charge monopoles infinite. We write the entropy Ryzhkin (2005); Kaiser et al. (2018) in terms of the densities of positive and negative monopoles, pick out a component of the free energy $G=-\mu{n_{+}}-\mu{n_{-}}-TS$ (with $\mu=2J_{\rm eff}$) and then substitute $n_{+}\rightarrow(n+\delta_{n})/2$ and $n_{-}\rightarrow(n-\delta_{n})/2$ where $n$ is the equilibrium monopole density and $\delta_{n}$ a charge fluctuation. A Taylor expansion of the free energy in powers of $\delta_{n}$, followed by setting the linear term to zero gives the equilibrium density Ryzhkin (2005):

[TABLE]

and a second order fluctuation term $G^{\prime\prime}=\mu_{0}kT(\delta_{n})^{2}/2nv_{0}$. The latter may be related to Eq. 8 using $(\delta_{n})^{2}=(\nabla\cdot{\bf M})^{2}/(Q/v_{0})^{2}$ to give Bramwell (2012):

[TABLE]

Eqs. 9, 10, 11 and 12 combine to give an analytic calculation of the full scattering function tensor. This is a low temperature approximation: for ${\rm Ho_{2}Ti_{2}O_{7}}$ parameters ($2J_{\rm eff}=-\mu=3.8$ K), it should be valid at $T=2$ K and below, where double charge monopoles are very scarce.

In Fig. 3 we compare the calculated scattering functions of NNSI with the numerically simulated one: the parameter-free comparison is excellent, showing that the continuum theory is a nearly exact description of NNSI over the wave vectors considered. This gives us great confidence that the theory can be adapted to describe dipolar spin ice (DSI) as well.

Inclusion of dipolar interactions – Given the excellent comparison of theory and simulation for NNSI, we seek to describe DSI by a natural extension of the NNSI theory. In particular, we simply add a term to Eq. 8 to describe the magnetostatic energy:

[TABLE]

where $\tilde{\bf H}({\bf r})$ represents the local internal far field arising from the dipoles at points ${\bf r}^{\prime}$ in the sample. Treating $\tilde{\bf H}$ as an independent variable, the Euler-Lagrange equation becomes:

[TABLE]

which can be solved for the structure factor as above if $\tilde{\bf H}$ can be expressed as a function of ${\bf M}$.

It is well known (see, for example, Ref. Morrish (1980)) that the field at a point external to a system of dipoles (e.g. at ${\bf r}$ where the dipoles are at ${\bf r}^{\prime}$) may be expressed as a sum over the fields arising from surface magnetic charge density $\sigma={\bf M}\cdot\hat{\bf n}$ (where $\hat{\bf n}$ is the unit vector normal to the surface) and volume magnetic charge density $\rho=-\nabla\cdot{\bf M}$. Without yet being explicit about surface effects, we write

[TABLE]

By Helmholtz’ theorem, the dipolar field is then equal to the irrotational (longitudinal) component of the magnetization, $\tilde{\bf H}=-{\bf M}^{\rm L}$. Writing ${\bf M}={\bf M}^{\rm L}+{\bf M}^{\rm T}$ thus renders the dipolar correction entirely local and we may proceed as above to find the structure factors:

[TABLE]

These are the same as Eq. 9 but with the the longitudinal structure factor suppressed by a factor $1/(1+\chi)$ at the zone centre. Also, the characteristic length becomes the Debye length $\xi/\sqrt{1+\chi}$ rather than the diffusion length $\xi$ Ryzhkin et al. (2013). In addition (not shown in Eq. 16) there is a delta function at the zone centre: this is elucidated for the case of different boundary conditions below.

It has been shown Bramwell (2017) that treating the dipolar field in this way is equivalent to Onsager’s cavity construction Onsager (1936) for dipolar fields, where the magnetic charge resides on the surface of a spherical cavity cut in a continuous polarizable medium. Alternatively, we can see that the dipolar integral is equivalent to a discrete summation over the effective magnetic monopoles of spin ice. These two different physical pictures give the same result only if the fields are unscreened Bramwell (2017). Hence we refer to this as the ‘unscreened’ model going forwards.

Following Eq. 7 and the discussion thereafter, we may once again compare $\textup{d}\sigma^{\rm NSF}_{hh2}/\textup{d}\Omega$ with $S^{\rm T}_{hh0}$ and $S^{\rm SF}_{hh2}$ with $S^{\rm L}_{hh0}\cos^{2}\theta$ up to $h=0.5$. Fig. 4 compares the DSI simulation with the unscreened dipolar calculation at $T=2$ K. We see that the simulated SF scattering function is immediately well described by the unscreened model, the NSF slightly less well. Since $\chi$ is large (e.g. $\chi\approx 6$ at $T=2$ K Bovo et al. (2013b)), the Lorentzian peak becomes very broad and flat compared to that of NNSI.

There are two small modifications to the unscreened model that are worth considering. First, we notice that the NSF scattering is better described by $S^{\rm T}\sim 2/3-S^{\rm L}$ (see Fig. 4). As discussed in Ref. Bramwell (2017), this correction is both plausible, because it satisfies the total moment sum rule, and justifiable, because the longitudinal fields were altered independently of the transverse. Second, the scattering function Eq. 21 apparently depends on $2J_{\rm eff}$ rather than the magnetic monopole chemical potential, which exceeds $2J_{\rm eff}$ in magnitude owing to the Coulomb interaction between monopoles. This difference increases the equilibrium defect density $n$, altering the specific heat Kaiser et al. (2018). As Eq. 21 depends on $n$ through $\xi$ (Eqs. 11, 12 and 21), it might be appropriate to replace $n$ by the true equilibrium density Kaiser et al. (2018). A fit to the ensuing function is tested in Fig. 4 but in fact it makes the description of the DSI simulation slightly poorer.

Boundary conditions and screening – To account for different boundary conditions, we add a ‘self field’ Griffiths (1982) to the dipolar field which generates the correct average when the field is summed over a ellipsoidal volume: in this way the demagnetizing factor $N_{\rm d}$ ($=1/3,0$ for spherical and Ewald boundaries respectively) naturally enters the equations, obviating the need to treat surface charge explicitly. This factor is conveniently introduced by treating the dipolar integral as a sum over thermally generated monopoles and allowing these to screen each other as in a Debye-Hückel gas Kaiser et al. (2018). In this approximation, the dipolar field is not simply equal to minus ${\bf M}^{\rm L}$, but in addition has a longer-range harmonic component. At the very least, such a component is expected from the demagnetizing field that arises from the surface charge (surface monopoles). It was suggested in Ref. Bramwell (2017) that it could also arise from thermally generated monopoles at mesoscopic distances Bramwell (2017), but we will examine the validity of this proposal subsequently.

We first consider the Debye-Hückel equation for the potential $\phi$:

[TABLE]

which applies outside a region of size $\sim a$ at the origin. By taking the gradient of this equation and introducing a delta function to extend it to the origin we find for the field at all ${\bf r}$:

[TABLE]

where $-N_{\rm d}{\bf M}({\bf r})\delta({\bf r})$ is the self-field.

The Debye-Hückel solution $\kappa\rightarrow\kappa_{\rm DH}$ then gives the Greens function solution for the more general equation where charge density is widely spread:

[TABLE]

In Fourier components the field becomes

[TABLE]

and the solution of the Euler-Lagrange equation is then:

[TABLE]

[TABLE]

Setting $\kappa\rightarrow 0$ in Eq. 21 results in the unscreened scattering functions (Eq. 16). In Fig. 4, to compute the scattering functions, we left $\kappa$ small but finite, which broadens the delta function slightly to render it visible to the eye.

It might seem surprising that we have to remove the screening ($\kappa\rightarrow 0$) in order to describe the numerical data: after all, it is well established from analysis of the specific heat, that monopoles are screened with a finite $\kappa$ Kaiser et al. (2018). In fact, it was shown in Ref. Ryzhkin et al. (2013) that the monopole potential $\phi$ consists of both a screened part and an unscreened part (see Fig. 5). The energy and hence specific heat is determined by the short range, screened, part of the potential whereas the pinch point is determined by the long ranged, unscreened part. Our analysis using the Poisson-Boltzmann equation above forces the potential into a perfectly exponential decay. Hence Debye-Hückel theory can only describe the pinch point if the screening length $\kappa^{-1}$ is set to be infinite, recovering the Poisson equation. In contrast, to calculate the specific heat Kaiser et al. (2018), it incurs negligible error to neglect the power law tail of the potential, and equate $\kappa$ with the Debye length of the monopole fluid.

Eq. 21 suggests further possibilities to test. Varying $\kappa$ and $\xi$ generates a family of curves, as illustrated in Fig. 6. We see that finite $\kappa$ broadens the central delta function into a Lorentzian in the longitudinal channel, plus a inverse Lorentzian ‘dip’ in the transverse channel. We consider the model in which $\kappa$ retains its identity as the inverse Debye length for single and double-charge magnetic monopoles. This may be calculated as in Ref. Kaiser et al. (2018), which has been shown to give a very accurate calculation of specific heat in both simulation and experiment. Additionally, setting $\xi=0$ results in an entirely flat background to this central Lorentzian. This flat background is the scattering function of a dipolar fluid in the Onsager picture - the ‘harmonic phase’ of Ref. Bramwell (2017). Hence the right hand curve of Fig. 6 is the scattering function of dipoles screened by magnetic monopoles, the ‘screened dipolar fluid’ Bramwell (2017).

Infinitely sharp pinch points in DSI simulation– We have already established (Fig. 4) that the unscreened dipolar model describes our simulation well, but can we actually rule out the screened model (finite $\kappa$)? To test this, we examine the simulation at a much higher temperature (10 K) where $\kappa$ should be sufficiently large to make the central Lorentzian (broadened delta function) easily visible within the resolution of our simulation. The result, shown in Fig. 7, does indeed rule out any pinch point broadening on the expected scale. (Note that, to obtain the fit to the unscreened model, we adjusted $\xi$ slightly to account for double charge monopoles). We also note that our simulation largely rules out the screened dipolar fluid model of Ref. Bramwell (2017). The ‘harmonic phase’ (flat background), is, however, present to a certain approximation, given that the ‘diffusion Lorentzian’ is very broad and flat at higher temperatures.

The physical origin of the infinitely sharp pinch points in dipolar spin ice simulations may now be traced to the unscreened part of the monopole potential calculated in Ref. Ryzhkin et al. (2013) and already mentioned. This feature was also noted as a dipolar tail in the correlation function in Ref. Sen et al. (2013). Referring to Fig. 5, the long-ranged unscreened part of the potential gives rise to a ‘gap’ between longitudinal and transverse eigenvalues of the scattering tensor at $q\rightarrow 0^{+}$, and hence the zone centre delta function that we observe. However, in view of the experimental result (see below), it was suggested in Ref. Ryzhkin et al. (2013) that the gap is unrealistic. This is an issue that we consider further below.

IV The pinch point paradox

In the preceding Section we have established a parameterized form of the pinch point profile that accurately and quantitatively describes high resolution simulations of the dipolar spin model. This is valuable for making a bridge to experiment as through parameter adjustment it allows for other physical possibilities and physical insight. As already advertised, it is in this confrontation of theory and experiment that the pinch point paradox arises.

The pinch point lineshape of polarized neutron scattering in ${\rm Ho_{2}Ti_{2}O_{7}}$ spin ice has been independently studied by Fennell et al. Fennell et al. (2009) and by Chang et al. Chang et al. (2010). Both sets of authors found that it can be described by a central Lorentzian plus a temperature-dependent flat (i.e. $q$-independent) component at all temperatures. An explanation of the flat component was put forward by one of us in Ref. Bramwell (2017) in terms of the screened dipolar fluid model described above, i.e. $\xi=0$, finite $\kappa$ as calculated in Ref. Kaiser et al. (2018). This quantitatively produced the observed temperature dependence of the flat component as well as producing a Lorentzian central peak (note that in subsequent figures, the ‘flat’ background is not completely flat because of the $\cos^{2}\theta$ factor discussed above). However, the screened model of Ref. Bramwell (2017) is already ruled out as a model for the DSI simulations.

The essence of the pinch point paradox is shown in Fig. 8 where we compare the SF scattering at $T=2$ K measured by Chang et al. Chang et al. (2010) with the DSI simulation, and the screened and unscreened theories (the latter for $\xi=0$, the dipolar fluid case Bramwell (2017)). The experimental data were scaled such that the $q=0$ point is coincident with the theory: i.e. we assume that theory and experiment are consistent as regards the bulk susceptibility (this is an approximation, but a reasonably accurate one). We see very clearly that the experiment rules out the unscreened DSI but is remarkably well described by the screened dipolar fluid. Indeed, the screened line shape can be transformed Bramwell (2017) into independent flat and Lorentzian components, where the Lorentzian component has the same width parameter as NNSI. This is in perfect agreement with the analysis of Chang et al. who fitted this data using Monte Carlo simulations of NNSI plus an arbitrary flat background: our theory eliminates the adjustable parameter used by these authors. They also drew a central Lorentzian through the points of their simulations of DSI – a reasonable extrapolation in view of knowledge available at the time, but one that we have now shown to be incorrect through higher resolution simulations. Because DSI has a delta function rather than a Lorentzian central peak, the discrepancy between simulation and experiment for spin ice is a very large one – yet both simulation and experiment should a-priori be considered correct, hence the paradox.

In the case of the 1.7 K SF, NSF and total neutron scattering data of Fennell et al. Fennell et al. (2009), the DSI again fails badly, while the screened dipolar fluid is again close to the experimental data (see Fig. 9 – the fit has been improved slightly by adjusting the demagnetizing factor to  0.2). Notably, the screened dipolar fluid model captures the striking ‘zone center’ dip in the SF scattering as well as the small ‘zone center bump’ in the total scattering. In general, we found that the screened model outperforms the unscreened one at all temperatures, in every case describing the flat (harmonic phase) background well, though at higher temperatures overestimating the Lorentzian width by a significant factor (e.g. $\sim 5$). Despite this, we can certainly conclude that, in the temperature range where we are most sure of the theory, the screened dipolar fluid model is fully consistent with experiment, while the DSI simulation and unscreened models are ruled out as descriptions of experiment.

V Discussion and Conclusion

Our study exposes two significant results: (i) the pinch point profile in the dipolar spin ice simulation is infinitely sharp, and (ii) experiment disagrees with simulation in that the profiles are broad, rather than sharp. Yet both are described quantitatively by different limits of the same theory.

To elaborate on (i), we recall that the subtle long-ranged correlations in spin ice are entirely revealed by neutron polarization analysis of the pinch point profile, as explained in the introduction (section I). In the near neighbor spin ice model the pinch points are broad, but addition of the long range dipole-dipole interaction causes them to be infinitely sharp with a delta function singularity at $q=0$ Sen et al. (2013); Ryzhkin et al. (2013). The structure factor (inverse) tensor remains an oblate spheroid (see Fig. 1 c,d), except precisely at $q=0$, where it becomes a sphere, in keeping with the cubic space symmetry Nye (1985).

As dipolar spin ice represents a model Coulomb gas, this raises an interesting and topical point. The Debye length for magnetic monopole interactions $(=\xi/\sqrt{1+\chi})$ is always finite, and any long-ranged correction to the exponential screening is essentially negligible as regards monopole correlations or specific heat. However the field-screening length ($=1/\kappa$) is, at the same time, divergent. Thus dipolar spin ice – a polarizable Coulomb gas that is relevant even as a model of water ice Ryzhkin (1984) – presents a very different picture as regards charge correlations and field correlations: the screening length for fields diverges while the Debye length remains finite. This is not a contradiction as it depends on how a particular experimental measure picks out the details of what is a rather subtle and long-ranged correlation function Ryzhkin et al. (2013) (see Fig. 5). In this context, we may make a connection with an interesting recent discussion of the experimental observations of ‘underscreening’ in dense ionic liquids Lee et al. (2015, 2017). We see that dipolar spin ice affords a tangible example where the screening length diverges, while the Debye length remain constant. This suggests that the observations of Ref. Lee et al. (2015, 2017) may be consistent with the dense ionic liquid having a stable polarization.

To elaborate on (ii) we have constructed flexible theories of the line shape that accurately and quantitatively account for both our numerical results and the experimental results: but only, paradoxically, under different limits. Specifically, in experiment magnetic monopoles appear to screen the dipolar fields, leading to broad pinch points, while in simulation they do not screen, leading to sharp pinch points.

In seeking a resolution of the pinch point paradox, we believe that a first obvious possibility – pinch point broadening by structural defects and disorder in the samples – may immediately be discarded. Although real samples do contain defects and disorder (see Refs. Sala et al. (2014); Revell et al. (2012); Ghasemi et al. (2018) for a discussion of this and how to prepare defect free samples), theory appears to rule this out as a cause of pinch point broadening. Thus, the role of defects and disorder has been considered theoretically in Ref. Sen et al. (2013). The case of non-magnetic dilution was shown to lead to no broadening of the dipolar pinch points. Since our results agree with this theory in other respects, it therefore seems highly unlikely that a small concentration of defects and disorder in the experimental samples can give the very strong pinch point broadening that is actually observed.

We are also confident that a second obvious possibility – broadening by extra terms in the spin Hamiltonian – may also be ruled out. We have checked this by simulating the refined $g-$DSI model (see Section II) which has been shown to give a most accurate description of ${\rm Dy_{2}Ti_{2}O_{7}}$ Yavors’kii et al. (2008). We find differences with DSI only near zone boundaries (see Fig. 7). This is to be expected: in this dipolar paramagnet, only long-ranged (i.e dipolar) interactions should affect the correlation function at small $q$. We also tested the ‘dumbbell model’ Castelnovo et al. (2008) (as implemented in Ref. Bovo et al. (2018)), but found this to be qualitatively the same as DSI.

Similarly, we have considered the possibility of various experimental corrections, but none of these give a convincing explanation of the paradox. For example, incomplete beam polarization, corrections to the static approximation, or other sources of background would not lead to the systematics that we observe. The fact that neutron data taken under different conditions on different instruments and sources are in excellent agreement further points against experimental corrections being complicit. It is also crucial to emphasize, as discussed above, that Chang et al.’s analysis of their data (Fig. 8) is quantitatively consistent with ours: this is a crucial proof that there is no inadvertent bias in our analysis (for example a scaling of the experimental data onto one theoretical curve rather than the other in Fig. 8).

These considerations leave us to face the likelihood of more fundamental causes of the pinch point paradox, including quantum fluctuations, corrections to neutron scattering theory and corrections to dipolar simulation methods.

**Quantum fluctuations – ** We have treated spin ice as a classical spin system. A careful appraisal of the problem Rau and Gingras (2015b) showed that the experimental systems considered here should approximate classical dipolar spin ice very closely. The only correction would be at the monopole sites where there is experimental evidence of a second order Zeeman effect, indicating the superposition of ‘spin up’ and ‘spin down’ states associated with magnetic monopoles Bovo et al. (2013a); Paulsen et al. (2019). However, it is hard to see this having a significant effect on the correlation function or structure factor measured in the static approximation of neutron scattering.

**Corrections to neutron scattering theory – ** The theory of neutron polarization analysis is well established following the pioneering work of Halpern and Holstein Halpern and Holstein (1941) and Moon et al. Moon et al. (1969). It is customary to map the orbital contributions to the magnetic moment on to effective spin operators, but the experimental systems considered here have largely orbital, rather than spin magnetism. We suggest it would be worth closely examining the polarization analysis theory for spin ice. For example, it is certain that there would be some degree of interaction of the neutron ‘within’ the orbital dipole, i.e. with the self-field, and we speculate that this might have some effect on the structure factor. Similarly, we might imagine some rather complex depolarization effects in spin ice, but these would have to go well beyond the standard picture of depolarization to explain our observations.

**Corrections to dipolar simulation method – ** Clearly, the method used agrees with experiment at large $q$ and at $q=0$, but does not do so at small but finite $q$. If we take the analytical theory to be a guide, we can expect the pinch point broadening to occur, regardless of the boundary conditions, so that it would even be present for periodic (Ewald) boundaries. Although it seems most unlikely that a tried and tested method like Ewald summation would have some basic flaw, it is also unlikely that it has ever been subject to such a stringent experimental test, so this question is worth examining in much greater detail. In this context, it is worth emphasizing that the simulation does agree with experiment in so far as the total scattering or trace of correlation function: the discrepancy is at the level of the Helmholtz decomposition or polarization analysis. As above, it would be worth reconsidering the role of ‘self field’ (the ‘internal structure of the dipole’) as an obvious point of difference between theory and experiment. It is not present in simulation, the demagnetizing effect being implemented via a demagnetizing tensor, but it remains a very real and potentially important Griffiths (1982) experimental property.

**Conclusion – ** In conclusion, pinch points typify an important and widespread correlation in condensed matter. As we have stressed, they do not necessarily imply a Coulomb phase or even a singular structure factor, but they remain a key diagnostic of long ranged, or critical correlations in diverse systems. Topics of current theoretical interest include the coexistence of pinch points with ordered states Banks and Bramwell (2012); Brooks-Bartlett et al. (2014); Benton (2016), extensions of the pinch point concept Yan et al. (2018); Benton et al. (2016a), pinch points as a diagnostic of quantum fluctuations Benton et al. (2012, 2016b) and pinch points in Coulomb systems McClarty et al. (2014); Slobinsky et al. (2019); Gray (2019). In general a sharp pinch point indicates that the structure factor tensor remains eccentric as $q\rightarrow 0$, typically as as a result of unscreened power law correlations. These pinch points will broaden if the long range correlations are exponentially damped, as occurs with the screening of long ranged interactions. Using very high resolution dipolar simulations and by comparing to existing experiment, we have exposed an unusual paradox in dipolar spin ice: the dipolar fields are not screened in simulation but do appear to be screened in experiment.

From a purely theoretical perspective, the pinch point paradox may be summarized as follows. In description of experiment, the Poisson-Boltzmann equation for magnetic monopoles appears where the Poisson equation is expected: it is as if the magnetic monopoles of spin ice behave much more like ‘real charges’ (Dirac monopoles) than expected Castelnovo et al. (2008).

More generally, given the special status of spin ice as a model ice rule system, a model dipolar system and a model Coulomb fluid, it is important to find the resolution of this paradox. When considered alongside recent results on ionic liquids Lee et al. (2017) it seems that the screening of long ranged interactions in condensed matter presents many aspects that remain to be understood.

Acknowledgements.

It is a pleasure to thank T. Fennell, P. C. W. Holdsworth and C. Gray for many useful discussions and related collaborations and P. C. W. Holdsworth for critical comments on the manuscript. 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). We gratefully acknowledge the NVIDIA Corporation for the donation of GPU resources. M.T. is supported by Stiftelsen Olle Engkvist Byggmästare (187-0013) with support from Magnus Bergvalls Stiftelse (2018-02701). The authors declare no competing financial interests.