Determining the anisotropy and exchange parameters of polycrystalline spin-1 magnets

Although low-dimensional $S = 1$ antiferromagnets remain of great interest, difficulty in obtaining high-quality single crystals of the newest materials hinders experimental research in this area. Polycrystalline samples are more readily produced, but there are inherent problems in extracting the magnetic properties of anisotropic systems from powder data. Following a discussion of the effect of powder-averaging on various measurement techniques, we present a methodology to overcome this issue using thermodynamic measurements. In particular we focus on whether it is possible to characterise the magnetic properties of polycrystalline, anisotropic samples using readily available laboratory equipment. We test the efficacy of our method using the magnets [Ni(H$_{2}$O)$_{2}$(3,5-lutidine)$_{4}$](BF$_{4}$)$_{2}$ and Ni(H$_{2}$O)$_{2}$(acetate)$_{2}$(4-picoline)$_{2}$, which have negligible exchange interactions, as well as the antiferromagnet [Ni(H$_{2}$O)$_{2}$(pyrazine)$_{2}$](BF$_{4}$)$_{2}$, and show that we are able to extract the anisotropy parameters in each case. The results obtained from the thermodynamic measurements are checked against electron-spin resonance and neutron diffraction. We also present a density functional method, which incorporates spin-orbit coupling to estimate the size of the anisotropy in [Ni(H$_{2}$O)$_{2}$(pyrazine)$_{2}$](BF$_{4}$)$_{2}$.

are checked against electron-spin resonance and neutron diffraction. We also present a density functional method, which incorporates spin-orbit coupling to estimate the size of the anisotropy in [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 .

Introduction
The investigation of low-dimensional quantum magnets is a key thrust of condensed matter physics. Of particular relevance here are quasi-one-dimensional and quasi-twodimensional systems based on S = 1 magnetic moments, which inspire a great deal of contemporary theoretical attention (e.g. [1][2][3][4][5][6][7][8]) and are predicted to display vibrant phase diagrams arising from competing interactions and their interplay with singleion anisotropy. These diagrams encompass quantum critical points [1,2], nematic and supersolid states [5,9,10], as well as topologically interesting gapped and quantum paramagnetic phases [11][12][13][14]. By contrast, because of the difficulty in making real examples of these systems, experimental work in this area moves more slowly and several predictions remain untested. While recent advances have been made with moleculebased magnets [15][16][17][18][19][20][21][22][23], difficulties in obtaining high-quality single crystals of the newest materials continue to hinder progress.
Detailed thermodynamic studies of single-crystal samples can be used to find anisotropy parameters and the size of the primary magnetic interactions, and reveal the ground state of the system. Unfortunately, crystals of sufficient size for such measurements are often hard to come by, particularly in the case of the newest materials, which are frequently synthesized initially as powders. Optimising the synthesis procedures for clean crystallization of a particular material typically requires considerable time and effort. It is therefore advantageous to be able to characterise the basic properties of a powdered anisotropic magnetic material using simple, readily accessible measurement techniques in order to identify the compounds that merit the additional work required for crystal growth. However, the complication of powder-averaging leads to difficulties interpreting the results of bulk thermodynamic measurements. This issue is made worse if the magnitude of the anisotropy is on a similar energy scale to the strength of exchange interactions in the compound [18,22,24].
Here we discuss the interpretation of experiments on polycrystalline samples of S = 1 magnets under the influence of uniaxial and rhombic single-ion anisotropy with energy D and E, respectively, and with the possibility of nearest-neighbour Heisenberg exchange J acting between spins. We describe the effect of powder averaging and discuss to what extent the parameters in the Hamiltonian can be extracted from data, focusing particularly on bulk thermodynamic measurements of susceptibility, magnetization and heat capacity that can, in principle, be performed using commonly available laboratory apparatus without the need to access equipment at a large user facility. We will start by showing that for materials with anisotropic spins, but negligible exchange interactions, a good estimation of the parameters can be readily extracted from powder data for both easy-plane and easy-axis systems, provided low enough temperatures and high enough magnetic fields can be achieved. We next consider the more challenging situation in which antiferromagnetic exchange interactions are finite and similar in energy to the single-ion anisotropy. We test the reliability of our methods by comparing the findings derived from thermodynamic probes with additional facility measurements of neutron diffraction and high-frequency electron-spin resonance (ESR). Finally we describe an approach, using density functional theory and spin-orbit coupling, to provide reliable estimates of single-ion anisotropy.
We will apply the analysis to three new materials in which Ni(II) ions are separated from one another by molecular ligands. These are (1) [Ni(H 2 O) 2 (3,5-lutidine) 4 ](BF 4 ) 2 and (2) Ni(H 2 O) 2 (acetate) 2 (4-picoline) 2 , both of which are found to be dominated by single-ion anisotropy with no evidence of a significant role for exchange interactions at the temperatures measured, and (3) the antiferromagnet [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 . System (1) was designed to have NiN 4 O 2 octahedra similar to that of (3), but without extended interaction pathways, such that the effect of the local environment on the anisotropy could be elucidated in the absence of exchange, and we will discuss to what extent this approach has been successful.
While a subset of the methods outlined have been preliminarily tested in studies by some of the same authors [22,23], we combine the full methodology here for the first time. This work follows from a related investigation of how to extract exchange parameters in low-dimensional S = 1/2 antiferromagnets [25].
The results presented here are of relevance not only to the study of low-dimensional S = 1 magnets, but also to the growing field of Ni(II) single-ion magnets [26][27][28][29][30].

Systems with negligible exchange
In the absence of exchange interactions, the Hamiltonian that governs the magnetic properties in applied magnetic field iŝ where we apply the constraint (discussed below) Here z is defined by the local axial direction, g is the g-tensor = diag (g x , g y , g z ) and S = (Ŝ x ,Ŝ y ,Ŝ z ) are the S = 1 spin operators. A negative D corresponds to easyaxis anisotropy and positive D is easy-plane anisotropy. The Hamiltonian is readily solved and the energy eigenvalues for two values of E are displayed in Fig. 1 in both the easy-axis and easy-plane cases (assuming g to be isotropic). In the absence of E anisotropy, the easy-axis system [ Fig. 1(a)] has a doubly degenerate ground state that splits with applied field, and has no ground-state level crossing for any field direction. The degeneracy of x and y energy levels is lifted in the presence of a non-zero E [ Fig. 1(b)], and a ground-state energy level crossing appears for the field applied parallel to x. In contrast, for easy-plane anisotropy [ Fig. 1(c) and (d)] a crossover from a nonmagnetic to magnetic ground state occurs even if E = 0, but only for the magnetic field parallel to z. Features arising from these crossovers will be observable in the polycrystalline magnetization data.
Values of E outside the constraint (Eq. 2) can be encompassed by an exchange of coordinate axes, without changing the system properties [31]. For example, it can be shown (see appendix A) that a permutation of the coordinate axes leads to equivalent Hamiltonians, whose eigenvalues differ only by a constant energy shift, but which will have different values of D and E. This means that for experiments on polycrystalline samples, where the correlation between the z-axis of Eq. 1 and the crystallographic directions is lost, fitting of powder data (such as magnetic susceptibility or heat capacity, as described below) without constraining 0 < 3E < |D| can yield two apparently conflicting sets of anisotropy parameters, one with a negative axial parameter and one positive. However, only one of these sets will fulfil the constraint. The parameters can be interconverted via the relations given in appendix A.

Effect of powder averaging on magnetometry and heat capacity measurements
In a polycrystalline measurement of an anisotropic magnetic material the mixing of different crystal directions with respect to the applied field leads to a blurring or loss of information as compared to single-crystal experiments. However, at sufficiently low temperatures and high magnetic fields, features visible in the results of thermodynamic measurements can still yield quantitative data. Here we discuss the problem of measuring polycrystals and by looking at the results of powder-averaged simulations of isolated S = 1 systems we show how to draw conclusions about the magnetic parameters.

Magnetic susceptibility
It has been suggested previously that it is not possible to distinguish between easy-axis and easy-plane isolated S = 1 magnets using polycrystalline measurements of magnetic susceptibility alone [32]. Revisiting this subject, we find that it is indeed possible to distinguish these two cases at sufficiently low temperatures and also extract reasonable estimates of both D and E from fitting polycrystalline data. It is true, however, that at high-temperatures all anisotropy information is lost.
Magnetometry measurements performed in the high-temperature, paramagnetic limit show a linear dependence of the inverse susceptibility on temperature. Extrapolating this linear dependence to obtain a temperature-axis intercept (the Weiss temperature) in isotropic magnetically-interacting systems can be used to obtain an estimate of the size of the exchange energy via the familiar Curie-Weiss law. In exchange-free, anisotropic systems the direction-dependent Weiss temperature reveals information regarding the crystal-field parameters [33] and in particular, estimates for D and E could be deduced from single-crystal measurements. For example, if E = 0 and g is isotropic, then solving the Hamiltonian above in the high-temperature region yields the Weiss temperatures Θ xy ≈ D/6 and Θ z ≈ −D/3 for the field applied perpendicular and parallel to the axial direction, respectively (see appendix B). However, in a polycrystalline experiment these values will be averaged such that the measured Weiss temperature will approach zero, and it is necessary to make measurements at lower temperatures to characterize the anisotropy.
The eigenvalues of Eq. 1 are used to construct a partition function, from which can be found the form of the low-field molar susceptibilities for fields applied along the three principal axes [23]: where β = 1/k B T . The expressions are in agreement with those previously published for the case with E → 0 [32]. For these smoothly varying functions it is possible to obtain a  Figure 2. Temperature dependence of magnetic susceptibility for different field directions calculated using Equations 3 for |D| = 10 K, E = 1 K and isotropic g = 2.
The simple polycrystalline average χ av = 1 3 (χ x + χ y + χ z ) is also displayed. (a) and (b) show χ(T ) and χT (T ), respectively, for easy-axis anisotropy. (c) and (d) show the same quantities for easy-plane anisotropy. The insets show the extrapolated values of the high-temperature inverse susceptibility, indicating the zero value of the Weiss temperature for the averaged data. reasonable approximation to the results of a polycrystalline measurement from a simple average, χ av = 1 3 (χ x + χ y + χ z ). (This is in contrast to the magnetization simulations presented below, which display step-like features for certain field directions and so require averaging over more angles to reproduce the experimental data.) The success of the simple average used here will be inspected in more detail later in comparison with experimental data. The functions above are plotted together with χ av in Fig. 2(a) and 2(c) for the easy-axis and easy-plane cases, respectively, with |D| = 10 K and E = 1 K. For these values there is a clear distinction between the easy-axis and easy-plane data, with the susceptibilities in the easy-plane case reaching a saturated value as temperature is reduced, in contrast to the easy-axis case for which the susceptibility continues to rise down to much lower temperatures. The insets to these figures show the values of the inverse susceptibilities extrapolated to where they cross the T -axis, highlighting the polycrystalline averaging to zero of the Weiss temperatures. The same effect can also be seen in a plot of χT vs T [ Fig. 2(b) and 2(d)] in which, while the single-crystal data either strongly increase or decrease on cooling depending on the direction of magnetic field and the sign of D, the polycrystalline value remains roughly constant down to temperatures of the order of the largest term in the Hamiltonian. In principle, it would be possible to obtain an indication of the size of D from the temperature at which χ av T departs from its high-temperature value. However, more reliable estimates can be obtained by direct fitting of experimental data to χ av as described below. In both cases the feature can only be observed at thermal energies low compared to h c = gµ B µ 0 H c . The data in (b) and (d) are offset for clarity.

Magnetization
At sufficiently low temperatures the magnetization M i for fields applied parallel to i = x, y and z will be dominated by the ground state energy level crossings seen in Fig. 1. As such, for easy-axis anisotropy a step will be observed in M x given by which is zero for E = 0, while for easy-plane anisotropy there will be a step in M z at a critical field given by The abrupt changes in M at the critical fields mean that the simple average used above for susceptibility will not reliably reproduce the result of a measurement. Instead we perform a simulation of polycrystalline M (H) using a full angular average over many possible field directions [34], the results of which are shown in Fig. 3(a) and (c) for easyaxis and easy-plane cases respectively. The easy-axis (easy-plane) system will have a sharp increase in M x (M z ) at H = H c and hence a peak in the differential susceptibility. For a powder, this feature is reduced to a small bump, which is hard to discern in M , but is readily observed in dM/dH or d 2 M/dH 2 as seen in the figures. Thermal occupation of excited states obscures the crossing of the ground state and the strength of the features diminishes as temperature is raised. By simulating the differential susceptibility at different temperatures [ Fig. 3(b) and (d)], we find that the peak indicating the level crossing at H c can be observed only if the temperature is lowered below approximately 0.1 × gµ B µ 0 H c /k B in both cases.

Heat capacity
The zero-field magnetic heat capacity C mag (in units of N A k B ) resulting from solving Equation 1 is found to be and is in agreement with the expression previously published for the case with E → 0 [31]. As temperature is reduced the function reproduces anomalies in the heat capacity resulting from the zero-field splittings shown in Fig. 1, and can be used in combination with lattice heat capacity models to fit measured single or polycrystalline data as shown by example later. The evolution of C mag (T ) under applied field for a polycrystalline sample is distinct for the easy-axis and easy-plane cases. Simulations of C mag (T ) for both situations are obtained via a full polycrystalline average at various fields with E = 0 [34] and displayed in Fig. 4. For the easy-axis scenario C mag (T ) [ Fig. 4 (a)] shows a single broad maximum at a temperature set by D. In an applied field, the ground state degeneracy is lifted, resulting in the emergence of a second narrow peak at low temperatures. As the field is raised the two peaks merge and move to higher temperatures while increasing in amplitude and breadth. For easy-plane anisotropy [ Fig. 4 (c)] only a single broad peak is apparent, which initially drops in amplitude, gets broader and moves to higher temperatures as the field is applied. The shift in magnetic entropy from low to high temperatures caused by the field-induced splitting of energy levels can been appreciated from the form of the C mag /T curves shown in Fig. 4(b) and 4(d).
Another estimate of the size of D in systems with E = 0 can be obtained from the field dependence of the position, T max , of the broad maximum observed in C mag (T ). The values obtained from the simulated data are plotted in dimensionless units in Fig. 4(e), and can be parametrised as follows, with two field-dependent parameters, δ-the local gradient and γ-the local, extrapolated zero-field intercept. The correspondence between δ and γ is shown in Fig. 4(f). As will be shown later for experimental data, an estimate of γD and δ can be found from a linear fit to the measured values of T max vs gµ 0 µ B H, while the pre-factor γ can be uniquely determined for a particular sign of D by using the fitted value of δ and Fig. 4   are separated by approximately 9.2Å along the [101] crystal direction. As a result, the exchange interactions between the S = 1 Ni(II) ions are expected to be negligible and the low-temperature magnetic properties should be dominated by single-ion anisotropy. The breaking of four-fold rotational symmetry by the equatorial ligands suggests that E will be non-zero. The magnetic heat capacity C mag (T ) at various fields, obtained by subtracting the zero-field lattice contribution from the measured data.

Thermodynamic measurements
The magnetic susceptibility of a polycrystalline sample of [Ni(H 2 O) 2 (3,5lutidine) 4 ](BF 4 ) 2 measured at µ 0 H = 0.1 T is shown in Figure 6(a) and resembles the data for an ensemble of S = 1 moments with single-ion anisotropy, but no significant exchange interactions. The data are fitted with the function χ(T ) = χ av (T ) + χ 0 , under the constraint 0 < 3E < |D|. Here χ av is the simple polycrystalline average 1 3 (χ x + χ y + χ z ) of the susceptibility components defined by Equation 3, and χ 0 is a temperature-independent contribution. Successful fitting of the data requires using the approximation that g x = g y = g z = g in Equation 3. In reality this will not be the case; the presence of single-ion anisotropy suggests that the components of g are unequal, as the same effects give rise to both. Perturbation theory predicts that g z − g xy = 2D/λ and g x − g y = 4E/λ, where λ is the spin-orbit coupling, which is ∼ −500 K for Ni(II) in octahedral environments [31]. However, in the types of system we consider here, typical values are D ∼ 10 K and E ∼ 1 K, so that ∆g is expected to be ∼ 0.01. Thus the uniform g approximation is reasonable within the errors of the thermodynamic measurements. The parameters resulting from the fit to χ(T ) are g = 2.24(1), an easy-plane D = 8.7(2) K, E = 1.2(2) K and χ 0 = −8(1) × 10 −9 m 3 mol −1 . ‡ M (H) data measured at various temperatures are shown in Fig. 6(b). The data increase smoothly towards to the saturated value of 2.21(2)µ B per Ni(II), which is consistent with the polycrystalline averaged value of the g factor resulting from fitting χ(T ). The lowest temperature curves show a kink in M (H) close to 5 T. This is more clear on differentiating the 0.4 K data (inset), where it appears as a small bump resembling the feature discussed earlier that arises from a ground state energy level crossing. The position of the bump is µ 0 H c = 6.0(6) T. Using the easy-plane model (Eq. 5) an estimate of √ D 2 − E 2 = 9.0(9) K is obtained, which is in agreement with the susceptibility results. § Fig. 6(c). The data exhibit two peaks, one around 40 K due to the phonons, and a second at ≈ 3 K which is attributed to single-ion anisotropy. The proximity of lattice and magnetic contributions to the heat capacity mean that dealing with each separately is not possible. Instead we fit the data to a model C/T = C latt /T + C mag /T , where C latt approximates the lattice contribution using a model with one Debye and three Einstein phonon modes [34,35], and C mag is given in Eq. 6. The fit is shown in the figure as a solid red line and is seen to account well for the data across the whole temperature range. Also shown are the separate lattice (dashed line) and magnetic contributions to the fit. The anisotropy parameters resulting from the fit are easy-plane D = 10.4(1) K and E = 2.6(2) K. These values are in agreement with √ D 2 − E 2 = 9.0(9) K estimated from the magnetization data. The size of D is similar to that obtained from the fits of the susceptibility to the simple polycrystalline average model, while E differs by 50% from the susceptibility value.

Zero-field heat capacity measurements of polycrystalline [Ni(H
Heat capacity measurements were also performed at fixed values of applied field 0 ≤ µ 0 H ≤ 9 T. The lattice contribution determined from the fit to the zero-field data is subtracted from each trace to yield C mag (T ) at different fields, which are shown in Fig. 6(d). It is seen that the broad hump due to energy level splitting initially drops in amplitude as the field is turned on, and at higher fields broadens and shifts to higher temperatures. This is very similar to the behaviour of the simulated data shown in ‡ Note that a fit performed without constraining the upper limit of E can also yield D = −6.3(5) K and E = 3.7(1) K, with g and χ 0 similar to the constrained fit. These values for D and E do not obey 3E < |D| and are the result of the permutation of the coordinate axes described earlier. They reduce to the values obtained from the constrained fit under the transformations of Eq. A4.
§ We note that, while these M (H) data are collected using pulsed magnetic fields, the location of H c in this case is within the field and temperature range of more conventional magnetometers equipped with a 3 He refrigerator. The fit of the lattice contribution yields the following characteristic amplitudes, A i (J K −1 mol −1 ), and temperatures, θ i (K), of the Debye (i = D) and Einstein (i = E) phonon modes:

Electron-spin resonance
In order to check the validity of our proposed methodology for determining anisotropy parameters from thermodynamic measurements on polycrystalline S = 1 systems, we also investigated [Ni(H 2 O) 2 (3,5-lutidine) 4 ](BF 4 ) 2 using high-frequency ESR. Arguably, ESR is the technique best suited to the evaluation of single-ion anisotropy in a powdered sample. However, for all but the smallest zero-field splittings, the frequency and field regimes needed to observe the crucial transitions requires highly specialised, nonstandard equipment. Our measurements were performed at the National High Magnetic Field Laboratory, Tallahassee, Florida [34].
In a polycrystalline S = 1 sample with single-ion anisotropy, multiple ESR transitions are expected between the split triplet energy levels. At a given frequency, for the field applied parallel to the local i-axis (i = x, y or z), there are two transitions possible that obey the ESR selection rule, ∆m s = ±1: one at low field and one at high field, which we label β i and γ i , respectively. In addition, it is also possible to see an excitation with ∆m s = ±2, this so-called half-field transition is labelled α i . Formally such transitions are forbidden, but when the Zeeman energy is comparable to the zerofield splitting strong mixing between m s states occurs and the selection rule is relaxed. Examples of the α x , β x and γ x ESR transitions are indicated in Fig. 1(d). Additional lines may also be observed at positions that do not correspond to one of the Cartesian axes. These off-axis resonances may be present at the half-field transitions and have been known to dominate the polycrystalline spectra [36][37][38].
ESR spectra were recorded in first derivative mode at frequencies in the range 100 < ν < 630 GHz at 5 K and representative data are shown in Figure 7(a). A broad feature is observed around 1 T in the 156.0 GHz data that drops in field as the frequency is raised and is attributed to the γ z transition. At higher frequencies a large double resonance is observed (e.g. near 5 T at 412.8 GHz), which sharpens and moves to high fields with increasing frequency. The larger of these two peaks is attributed to the off-axis half-field transition α off , and is expected to lie very close to the α x transition. The broad hump at slightly higher fields is ascribed to β y transitions. At yet higher fields smaller features are seen that can be attributed to the β x and γ x transitions.
Transitions are labelled in the temperature-dependent spectra recorded at 412.8 GHz and shown in Fig. 7(b). Both the α off line and the β x line increase in amplitude as temperature is reduced, identifying them as excitations from the ground state. The γ x transition is seen to be smaller than β x at low temperatures and diminishes in amplitude when the temperature falls, which is expected for a transition between excited states. These observations identify the energy-level splitting as easy plane.
The frequency and field positions of the transitions are modelled with an easyplane energy-level scheme in Figure 7(c). Fitting is performed as described in Ref. [39] and the best fit to the data is found for the parameters g x = 2.21(1), g y = 2.17(1), g z = 2.16(3), D = 10.40(1) K and E = 2.11(4) K. The fit successfully reproduces most of the peak positions and, as shown in panel (b), a simulation of the 412.8 GHz spectrum arising from these parameters compares reasonably well with the measured data. Note that g x > g z , which is consistent with the results of perturbation theory for a Ni(II) ion with easy-plane anisotropy [31]. The observed splitting of the offaxis half-field transition is not explained by the simulations. One possible reason for the extra peak could be the presence of a second Ni(II) site in the ESR sample with slightly different anisotropy parameters, ¶ however no evidence of a significant impurity fraction is observed in the synchrotron x-ray diffraction measurements. Another possible explanation is the presence of small spin-spin couplings between the Ni(II) ions. The crystal structure does not show any evidence for significant exchange pathways and any magnetic interactions must be less than ∼ 1 K or their effect would be observed in the magnetometry measurements, but very small couplings have previously be found to lead to ESR peak splitting in molecular systems, even at elevated temperatures [40].
Whatever the explanation for the extra peak, the best fit D and E parameters account for the majority of the ESR resonances and are in excellent agreement with the values derived from the heat capacity analysis and the position of H c in the low temperature M (H) data.

Discussion
In light of the high-frequency ESR data we can judge the effectiveness of the analysis methodology for magnetometry and heat capacity. The thermodynamic measurements are all strongly indicative of an easy-plane anisotropy in this material, which is confirmed by ESR. The D and E parameters derived from fitting to the zero-field heat capacity agree closely with those obtained from ESR. The agreement is less good for the parameters deduced from fitting susceptibility data. The fitting function in this case makes use of the elementary polycrystalline average χ av = 1 3 (χ x + χ y + χ z ), which simplifies the fitting procedure, but perhaps does not sample enough field angles to fully account for the data. Nevertheless, the estimate of the size of the parameters obtained from the susceptibility matches the ESR and heat capacity values to within less than 20% for D and 50% for E.
Owing to the polycrystalline nature of the sample, it is not possible to identify the easy plane of [Ni(H 2 O) 2 (3,5-lutidine) 4 ](BF 4 ) 2 from the thermodynamic or ESR measurements alone. However the symmetry of the NiN 4 O 2 octahedra would strongly suggest that the hard z-axis is perpendicular to the NiN 4 equatorial plane.

Crystal structure
Having introduced the analysis methods with the previous easy-plane material, we now test them on a Ni(II) material with a different local environment. Ni(H 2 O) 2 (acetate) 2 (4picoline) 2 crystallizes in the orthorhombic space group P cab. Fig. 5(a) shows the coordination environment deduced from single-crystal x-ray diffraction performed at 150 K on a microcrystal, and structural parameters are found in Table 1. The material contains distorted NiO 4 N 2 octahedra, as compared to the NiN 4 O 2 octahedra in the system discussed above. In the present case, the local environment is made up of two axial nitrogen atoms donated by 4-picoline and four equatorial oxygen atoms, two donated by acetate and two from water. The three bond angles between opposite donor atoms in the nickel octahedra are all 180 • , and the cis O-Ni-N angle ranges between 87.3 • and 92.7 • .
The individual Ni(H 2 O) 2 (acetate) 2 (4-picoline) 2 molecular units are well-separated in the c-direction by the 4-picoline molecules. The closest Ni-Ni distance is approximately 7.6Å within the ab-plane, but with no apparent exchange pathway between nearest neighbours. Hence the magnetic properties are again expected to be that of an ensemble of magnetically isolated S = 1 moments with single-ion anisotropy.

Thermodynamic measurements
The polycrystalline magnetic susceptibility of Ni( The line is a fit to the simple polycrystalline average described in the text. (b) Magnetization measured at the temperatures shown. (c) Differential susceptibility dM/dH measured at 0.6 K. (d) Measured heat capacity divided by temperature (circles). The solid red line is a fit to the lattice plus magnetic model described in the text. The dashed green line is the lattice contribution and the dotted blue line is the magnetic part. (e) The magnetic heat capacity C mag (T ) at various fields, obtained by subtracting the zero-field lattice contribution from the measured data.
at µ 0 H = 0.1 T is shown in Fig. 8(a) and resembles that of an S = 1 anisotropic magnet with negligible interactions between the spins. Similar to the previous case, the data are fitted to χ av = 1 3 (χ x + χ y + χ z ) using the expressions in Eq. 3, with an isotropic g and under the constraint 0 < 3E < |D|. The fitted line reproduces the data well and yields the parameters g = 2.20 (1) The magnetization data measured up to 7 T at various temperatures are shown in Fig. 8(b). All traces show a smooth rise towards saturation, with the 0.6 K data approaching a moment of 2.12µ B per Ni(II) ion by 7 T. There is no clear sign of any feature due a ground state energy level crossing either in M (H) or dM/dH (inset) at the lowest temperatures. This is consistent with the expectation of an easy-axis system. Although, as mentioned earlier, a level crossing occurs for fields parallel to x for easyaxis materials, a feature in the polycrystalline magnetization data is only expected to + Note, without the 3E < |D| constraint, the fit also works well with g = 2.20(1), D = 4.9(2) K and E = 2.2(2) K. This is accounted for by the reverse permutation of crystallographic axes described earlier and the parameters can be mapped back on to those above via the inverse relations of Eq. A4.
. The estimate of D and E obtained from susceptibility suggests this condition is not met in our measurements.
Figure 8(c) shows the zero-field heat capacity of polycrystalline Ni(H 2 O) 2 (acetate) 2 (4-picoline) 2 . On cooling, the data exhibit a broad hump between 40 and 50 K due to phonons followed by a steep rise at low temperatures caused by single-ion anisotropy. To extract estimates of the anisotropy parameters it is necessary to fit the data below 18 K to the sum of a Debye phonon mode [34] and the magnetic term given in Equation 6. The resulting fit is displayed as a solid red line in the figure and is seen to compare well with the data at low temperatures. The separate Debye and magnetic terms in the fit are shown as dashed green and dotted blues lines, respectively. The anisotropy parameters resulting from the fit are easy-axis D = −6.7(1) K and E = 1.54(1) K, * and are within 10-15% of the values obtained from fitting the magnetic susceptibility.
Additional low-temperature heat capacity measurements are made in fixed magnetic fields. The fitted zero-field lattice term is subtracted from these data and the results are plotted as C mag (T ) in Figure 8(d). In small applied fields the data exhibit the lowtemperature rise due to the anisotropy. At higher fields this feature moves to higher temperatures and reveals itself to be a peak whose amplitude, width and position increase with increasing fields. This is consistent with the simulated data shown in Fig. 4(a), further confirming the easy-axis nature of this material.

Discussion
Low-field magnetic susceptibility, magnetization and heat capacity measurements all indicate the presence of easy-axis anisotropy in Ni(H 2 O) 2 (acetate) 2 (4-picoline) 2 where the Ni(II) ion is surrounded by four equatorial oxygens and two axial nitrogens. Judging from the local structure, it would be expected that the easy axis lies parallel or close to the axial N-Ni-N bond direction. The parameters taken from the heat capacity analysis, which on the evidence of the previous material offers the most accurate results, are D = −6.7(1) K and E = 1.54(1) K. In contrast, [Ni( where Ni(II) is surrounded by four equatorial nitrogens and two axial oxygens, is an easy-plane system with D = 10.4(1) K and E = 2.6(2) K (also from heat capacity).

Systems with significant exchange
As detailed above, it is possible from polycrystalline thermodynamic measurements alone to obtain good estimates for the parameters governing the magnetic properties of S = 1 systems in the absence of effective exchange pathways. Now we turn to systems containing antiferromagnetic interactions between the spins. The Hamiltonian in this * The fit of the lattice contribution yields the following characteristic amplitude and Debye  case isĤ where the sum in the first term is over unique nearest-neighbour exchange pathways with Heisenberg exchange strength J ij . In the two extreme cases, where the exchange term is much stronger than the anisotropy term or vice versa, then polycrystalline data can be used to parameterize the system. However, in the case where the two are similar in size then interpretation of the data can be problematic, as some of the present authors have discussed previously [18,22]. In this situation, there is a paucity of theoretical models that can be used in fitting either low-field magnetic susceptibility or zero-field heat capacity to obtain reliable estimates of the magnetic parameters. To some extent the application of magnetic field can help. The field both suppresses antiferromagnetism and shifts the features in heat capacity due to the anisotropy to higher temperatures, permitting them to be analysed. Features that provide useful information can also be discerned in polycrystalline measurements of M (H), specifically the spin-flop field (in easy-axis systems) and the saturation fields in the easy and hard directions. Here we illustrate these methods using experimental data collected on a Ni(II) coordination polymer. Fig. 9 shows the structure of this material as determined at 300 K using powder synchrotron x-ray diffraction. The coordination environment consists of NiN 4 O 2 octahedra with a small axial compression, but no distortion in the octahedral bond angles. The equatorial nitrogens are from the pyrazine molecules, which bridge the Ni(II) ions in the ab-plane forming a square planar array. The axial oxygens are provided by the water molecules that tether adjacent nickel-pyrazine sheets along c via a network of H· · · F bonds with the charge-balancing BF − 4 counter ions. The nearest neighbour Ni· · · Ni distance is 6.98Å through the pyrazine molecules, and metal ions in adjacent planes are separated by 7.40Å. Both pyrazine and H· · · F bonds have been shown to be mediators of antiferromagnetic exchange strengths of the order of 1-10 K in Ni(II) complexes [21,22]. It is not possible to tell from the structure alone which pathways will support significant exchange interactions, therefore we define the average nearest-neighbour exchange strength, J , as a sum of the exchange through pyrazine, J pyz , and water, J H 2 O , such that n J = 4J pyz + 2J H 2 O , where n is the total number of effective nearest-neighbour exchange pathways.

Experimental results for [Ni(H
The comparison of structural parameters in Table 1  has four distinct equatorial bond lengths and bond angles that depart somwhat from octahedral symmetry. Thus, in contrast to the Ni-lutidine system, E is expected to be zero in the high-symmetry Ni-pyrazine material. Nevertheless, the comparison suggests that D should have the same sign in the two systems with a similar order of magnitude, i.e. we anticipate that [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 has easy-plane anisotropy and D ∼ 10 K.

Thermodynamic measurements
The susceptibility of polycrystalline [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 taken at µ 0 H = 0.1 T is shown in Fig. 10(a). The data rises smoothly on cooling and exhibits a broad maximum around 4 K, followed by a cusp and a reduction down to 1.8 K. The inverse susceptibility (inset) is fit to a Curie-Weiss model across the range 100 < T < 300 K, yielding g = 2.19(1) and a temperature-independent contribution χ 0 = 1.3(1) × 10 −9 m 3 mol −1 . The same data are plotted as d(χT )/dT in Fig. 10(b). This quantity is known to resemble the behaviour of the heat capacity of simple antiferromagnets in the region of a transition to long-range order [41]. The data show a lambda-like peak close to the cusp observed in χ, indicative of an antiferromagnetic transition at 3.0(1) K. Fig. 10(c) shows pulsed-field magnetization data taken at various fixed values of temperature. As the field is swept, the data display a slightly concave rise followed by a rounded approach to saturation, distinctive of an S = 1 antiferromagnet with single-ion anisotropy. Above 15 T at the lowest temperatures the moment approaches a saturated value of 2.10(1)µ B per Ni(II), which suggests a low-temperature value of g = 2.10(1). There is no indication of a spin flop in the data, which is consistent with the expectation of easy-plane anisotropy in this material. Following Reference [22], we expect to see two characteristic fields in a polycrystalline measurement of M (H) of an easy-plane system: one at the point where moments saturate for fields lying in the easyplane, and the other where moments saturate for fields parallel to the hard axis . These occur at µ 0 H x sat = 2n J /gµ B and µ 0 H z sat = 2(n J + D)/gµ B , respectively. A change in curvature is observed in the low-temperature data between 5 and 6 T, which appears as a peak in dM/dH [as shown in Fig. 10(d)]. We associate this feature with fields within the easy-plane and find µ 0 H x sat = 5.7(3) T. Furthermore, we define the hard-axis saturation as the point at which d 2 M/dH 2 first approaches zero, hence µ 0 H z sat = 15.7(5) T. From these two values we estimate n J = 4.0(2) K and D = 7.1(6) K.
Zero-field heat capacity measurements performed on [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 [plotted as C/T in Fig. 11(a)] reveal a broad hump due to phonons located around 50-70 K and a lambda peak indicating a transition to long-range antiferromagnetic order centred at T N = 3.0(1) K, in agreement with the susceptibility value. To extract the lattice contribution, a fit is made to the data in the range 24 ≤ T ≤ 300 K using a model of one Debye and three Einstein phonon modes [34,35]. The result, shown as a red line, agrees well with the data across the fitted temperature range. The magnetic part of the heat capacity, C mag , is isolated by subtracting the lattice contribution. The magnetic entropy is calculated by integration and found to approach the expected value of R ln 3 at temperatures in excess of 10 K. C mag (T ) measured in fixed applied fields is shown in Fig. 11(c). The lambda peak associated with antiferromagnetic ordering is seen to be suppressed as the field is increased, while a broad shoulder appears to the high-temperature side of the peak and shifts to higher temperatures with increasing field. This broad feature is associated with the single-ion anisotropy and its temperature evolution with field is shown in more detail in the left-hand inset. The right-hand inset shows the position of the hump, which can only be discerned at the highest measured fields, plotted against gµ 0 µ B H/k B . Following the discussion of Eq. 7, we perform a linear fit of these data to find δ = 0.20 (7) and γ|D|/k B = 2.5(8) K. Using these values together with Fig. 4(f), and assuming an easy-plane scenario, we estimate D = 7(2) K, which agrees with the result from the magnetization data.

Muon-spin relaxation
In order to confirm the presence of long-range magnetic order suggested by the heat capacity data, muon-spin relaxation (µ + SR) measurements were made on Example spectra are shown in Fig. 12(a). At temperatures T < 3.2 K the asymmetry shows heavily damped oscillations at two distinct frequencies whose magnitudes decrease with increasing temperature. At temperatures T > 3.2 K, oscillations are seen at lower frequency, but show little variation as the temperature is further increased. The oscillations measured for T < 3.2 K are characteristic of a quasistatic local magnetic field at the muon stopping site usually attributable to long-range magnetic order, which causes a coherent precession of the spins of those muons with a component of their spin polarization perpendicular to this local field. The frequencies of the oscillations are given by f i = γ µ B i /2π, where γ µ is the muon gyromagnetic ratio (= 2π × 135.5 MHz T −1 ) and B i is the average magnitude of the local magnetic field at the ith muon site. The oscillations for T > 3.2 K are caused by dipole-dipole coupling between the muon and fluorine nuclei and are typically resolved in the paramagnetic regime. Detailed modelling of these two regimes is described in [34]. This allows us to conclude that the material undergoes a transition to long-range order throughout its bulk at T N = 3.2(1) K, which is in excellent agreement with the heat capacity data.
A full quantitative structural refinement of the data is made difficult by the dynamics of the water molecules. A LeBail fit of the nuclear Bragg peaks observed at 10 K is fully consistent with the reflection conditions of the space group I4/mcm and yields lattice parameters a = 9.8859(2)Å, b = 9.8859(2)Å, c = 14.6625(4)Å, which are in good agreement with the results of the structural refinement of the room-temperature x-ray diffraction data taken on the non-deuterated material. To account for the fact that the reflections with a sizeable projection on the c * reciprocal lattice vector were found to be broader than others, the fit to the neutron data includes a strain model that represents a small degree of decoherence along the crystalline c-axis.
Taking the difference in scattered neutron intensity obtained at 1.5 and 10 K reveals three magnetic diffraction intensities [34], the positions of which can be indexed by the propagation vector k = (0, 0, 0) with respect to the reciprocal lattice of the paramagnetic unit cell [see Fig. 12(b)]. The three peaks correspond to the following families of reciprocal lattice vectors: {1, 0, 1}, {1, 0, 3}, and {2, 1, 1} and are attributed to long-range magnetic order of Ni(II) ions as observed using µ + SR. None of the observed magnetic reflections violate the I-centring reflection condition, which means that magnetic moments related by I-centring (corresponding to a translation The magnetic moment directions are determined by fitting the relative intensities of the magnetic diffraction peaks in the subtracted data. Because of the difficulties with the nuclear refinement, the magnetic intensity could not be calibrated and the scale of the magnetic phase was left free to refine. As can be seen in the insets to Fig. 12(b), the relative diffraction intensities were consistent with magnetic moments aligned in an undetermined direction perpendicular to c, confirming the presence of easy-plane anisotropy consistent with the results of thermodynamic measurements.
The square root of the integrated intensity of the {1, 0, 1} magnetic peak, which is proportional to the ordered moment, is plotted as a function of temperature in Fig. 12(c). A fit to a power-law dependence, m(T )/m(0) = (1 − T /T c ) β , yields T c = 2.91(2) K, which is consistent with the position of the ordering peak in the heat capacity. While the sparseness of the data in the vicinity of T c limits the sensitivity of the fit to the exponent β, we note that the fitted value of 0.25 (2) is in excellent agreement with values found for experimental realizations of the 2D XY model [43].

Discussion
Thermodynamic measurements on [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 are all consistent with easy-plane anisotropy, which is verified by neutron diffraction. Isothermal magnetization suggests D = 7.1(6) K and n J = 4J pyz +2J H 2 O = 4.0(2) K. The value of D agrees with the analysis of the high-field heat capacity and is within 32% of the value found above for the non-interacting system [Ni(H 2 O) 2 (3,5-lutidine) 4 ](BF 4 ) 2 . Previously measured values of J pyz in related Ni(II) compounds are ∼ 1 K [22], which might suggest that J H 2 O J pyz and that our material is a highly two-dimensional antiferromagnet. This would be consistent with the analysis of the critical exponent extracted from the neutron data.

Using density functional theory to obtain single-ion parameters
To model the Hamiltonian parameters in [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 we performed a sequence of density functional theory (DFT) total energy calculations. Calculations were performed within the DFT plane wave formalism as implemented in the Castep code [44,45]. The exchange-correlation interactions were described with the PBE generalised gradient functional [46], and ultrasoft pseudopotentials [47] were used for the core-valence interactions. Numerical convergence of the plane wave basis set (plane wave cut-off and k-point sampling) was set at a tight tolerance such that total energy differences were converged to better than 0.01 meV/cell to obtain accurate results for coupling constants [48]. Geometry optimisations were performed using a BFGS energy minimisation algorithm until the maximim residual force on atoms were all below 0.05 eV/Å. Spin-orbit coupling, implemented in Castep with the formalism of Dal Corso, et al. [49] was also used, where j = l ± 1/2-resolved pseudopotentials are obtained from a fully relativistic radial atomic Dirac-like equation. This is needed because the strongest part of the spin-orbit interaction is within the core and so, in a plane-wave calculation, it must be dealt with via the construction of a j-dependent pseudopotential. We then use a 4-component spinor as a pseudowavefunction, rather then the usual 2 component spin up/down formalism. The 4-component wavefunction allows for local spin orientations and permits inclusion of spin-orbit coupling, which is closely related to non-collinear magnetism. At each point in space there is a local direction to the spin polarisation and this is used to evaluate the exchange-correlation interaction using standard functionals. The magnetic structure is not the same as the crystallographic structure and hence in the electronic structure calculations we do not impose a predetermined symmetry on the electronic charge densities.
To extract the Heisenberg coupling constants, we employ the method described in [34] involving comparing several collinear spin configurations. After geometry optimisation of the system, there are two different pairs of Ni-Ni interactions. The structure therefore suggests two exchange constants: J 1 within the Ni-pyz planes and J 2 between them. In spin-polarised systems each energy minimum in an electronic structure calculation corresponds to a magnetic structure. To investigate likely magnetic structures (collinear in the first instance) the electronic structure is initialized with various spin configurations and energy minimized to the nearest local minimum electronic magnetic state. The spin structure of the lowest energy state is shown in Fig. 13. It forms an antiferromagnetic state within each of the Ni-pyz planes, offset by (1/2, 1/2) in neighbouring planes. To evaluate the intraplane J 1 coupling and the interplane (next nearest neighbour) J 2 coupling, differences in energy of spin configurations are taken giving J 1 = 0.64(1) meV and J 2 = 0.65(1) meV, where the uncertainty is that of numerical noise in the calculation [50]. Our calculations therefore predict an antiferromagnetically ordered ground state with an isotropic exchange J 1 ≈ J 2 ≈ 8 K. The calculations overestimate the J-couplings compared to the experimental results. This is likely due to the use of the PBE functional which may underestimate the localisation of the Ni d electrons, allowing slightly more neighbour-neighbour overlap and increasing the apparent strength of the magnetic coupling. (Such a systematic effect has been noted previously in GGA+U calculation in Ni-pyz-based systems [18].) It is also worth noting that in previous calculations of exchange effects in a coordination polymer magnet [23], we found that a similar overestimate resulted from the neglect in the calculations of the effects of structural disorder, which acted to strongly reduce the exchange coupling. As noted above (Fig. 9), the water molecules that mediate the interplane exchange in this material exhibit a degree of positional disorder. This then could act to suppress the interplane coupling, leading to a quasi-two-dimensional magnet, which would be consistent with the critical exponent extracted from the neutron data.
In addition to exchange, the energy scales of single-ion anisotropy effects can be investigated by examining the dependence of the energy on the direction of the spin configurations. This requires that the spin-orbit interaction is taken into account in the electronic structure calculations and that the spins are allowed to adopt non-collinear configurations. Spin anisotropy, D, of the system can be found by examining energy differences for the spin configurations that are possible in various orientations. For this the atomic anisotropic energy expression H = DS 2 x + E S 2 x − S 2 y is used. Electronic structure calculations are carried out, initialising the spins to be aligned along the x, y and z directions. We find D = 8.5(2) K and E ≤ 0.2 K. The prediction for D is in good agreement with the value of 7.1(6) K established using magnetometry above. The value of E falls within the limits of resolution of the calculation itself [50], and so can be considered zero within the errors, which agrees with the expectation that E = 0 in this tetragonal system. Lastly, we note that, in computing these values, it is important to include spin-orbit coupling since this contributes significantly to the anisotropy coefficient. Similar calculations without spin-orbit coupling greatly underestimate D, giving |D| 0.2 K.

Conclusions
In this paper we have presented an experimental method for extracting the anisotropy parameters of polycrystalline S = 1 magnets from thermodynamic data and applied it to the situation of magnetically-isolated, exchange-free systems as well as an extended material with antiferromagnetic exchange pathways. We have sought to determine to what extent quantitative information can be achieved using readily accessible measurement techniques.
More generally, in the case of the exchange-free magnets, we have found that fitting zero-field heat capacity data yields reliable values for D and E, and confirmed this using high-frequency ESR. The field-dependence of the heat capacity is a useful check on the sign of the D and the magnitudes of the parameters can be further confirmed by features in the isothermal magnetization. Low-field magnetic susceptibility is a relatively quick and available technique. We find that fitting the results of such measurements using the expressions described can provide rough estimates of D and E. For the antiferromagnetic system we were unable to extract quantitative information about the anisotropy from low-field susceptibility and heat capacity data. Instead, low-temperature magnetization in fields up to the hard-axis saturation is required to find this information. The results can be checked by measuring heat capacity in fields sufficiently high to separate the antiferromagnetic ordering peak and the anomalies that arise due to energy level splittings.
In all cases, the experiments require temperatures low compared to the anisotropy energy. For the exchange-free systems, successful magnetization and fixed-field heat capacity measurements also depend upon applying magnetic fields which are on the scale of the anisotropy energy. The values of anisotropy found above for our materials are representative of those in octahedrally coordinated Ni(II) compounds [31], suggesting that the measurements can be performed in standard, commercially available equipment. On the other hand, the magnetization measurements needed to find useful information in the case of the antiferromagnet materials requires fields ∼ (n J + D)/µ B . In the case of [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 , J is small and so fields < 16 T were sufficient. In other molecule-based Ni(II) systems, fields in the range of pulsed magnets (∼ 70 T) may be required.
For the six-coordinate Ni(II) complexes described here, we have observed a correlation between the Pauling electronegativity (EN) value of the ligand donor atoms and the magnetic ground state. In the case of [Ni(H 2 O) 2 (3,5-lutidine) 4 ](BF 4 ) 2 and Ni(H 2 O) 2 (acetate) 2 (4-picoline) 2 , trans -NiO 2 N 4 and NiO 4 N 2 octahedra are found for which N and O donor atoms have EN values of 3.04 and 3.44, respectively. The difference in EN values determines the resultant Ni(II) magnetic moment direction. Thus, the Ni(II) moment lies in the direction that includes the donor atoms of lower EN; i.e., either along the N-Ni-N axis (Ising-like) or within the NiN 4 plane (XY -like). This observation is fully consistent with all of the thermodynamic data. Finally, we have described density functional theory calculations that incorporate spin-orbit coupling in order to estimate single-ion anisotropy parameters. The calculated values for [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 agree well with the experimentally derived results. Further calculations on other materials are underway to verify the wider applicability of this approach.

Acknowledgements
PAG acknowledges that this project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 681260). This work is also supported by EPSRC (including EP/N023803/1 and EP/N024028/1). We acknowledge Durham for Hamilton supercomputer time and also EPSRC Grant EP/P022782/1 for Archer compute resource. RDJ acknowledges financial support from the Royal Society. A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by National Science Foundation Cooperative Agreement No. DMR-1157490 and the State of Florida, as well as the Strongly Correlated Magnets thrust of the DoE BES "Science in 100 T" program. JAS acknowledges support from the Independent Research/Development (IRD) program while serving at the National Science Foundation. We are grateful to SµS for the provision of beamtime and to A. The work at EWU was supported by the NSF through grant no. DMR-1703003. Data presented in this paper resulting from the UK effort will be made available at XXXXXXX.

Appendix A. Single-ion anisotropy calculations
In zero field, the Hamiltonian of an S = 1 exchange free system may be written where we have used the S = 1 spin matrices. As mentioned in the manuscript, equivalent ways to write this Hamiltonian can be obtained by permutation of the coordinate axes. For example, we can make the transformation (x, y, z) → (z , x , y ).
In a powder measurement, all information about the identity of z with respect to the crystallographic axes is lost. In this case we can shift the Hamiltonian by a constant energy and write it in the standard form, where I is the 3 × 3 identity matrix. The two sets of anisotropy parameters can thus be interconverted via the relations Only one set of parameters will fulfil the constraint 0 < 3E i < |D i |.
The other cyclic permutation (x, y, z) → (y , z , x ) is also possible. In this case the Hamiltonian must be shifted by (D 1 + E 1 ) and the transformed parameters are given by We note that the full derivation of the Hamiltonian (Equation 1) treats the spinorbit and Zeeman interactions as perturbations to the ground state of a single ion [51]. This leads to an expression for the anisotropy Hamiltonian written in terms of the matrix elements of the orbital angular momentum operatorL µ . We define where |0 is the unperturbed ground state of the system with energy E 0 and |n are the excited states with energies E n . The anisotropy Hamiltonian is then given bŷ where λ is the spin-orbit interaction constant. On making the permutations of the coordinate axes it can be verified that the original matrix elements are recovered by making the substitutions given above.

Appendix B. Curie-Weiss temperatures
The expressions for the magnetic susceptibility along the principal axes of the exchangefree S = 1 system are given in Equation 3. Evaluating the expressions in the limit E → 0 and g x = g y = g z = g yields as found elsewhere (e.g. [32]). Inverting these expressions and expanding in the limit of k B T D gives the high-temperature susceptibility in the planar and axial directions: which resemble a Curie-Weiss behaviour with apparent Weiss temperatures Θ xy ≈ D/6 and Θ z ≈ −D/3, respectively. Taking a simple approximation to a powder average χ av = 1 3 (χ x + χ y + χ z ) of these expressions again yields a Curie-Weiss form in the high-temperature limit, but with a zero Weiss temperature. The simulations shown in Figure 2 also indicate that θ av ≈ 0 for systems with non-zero E.

Synthesis of Ni(H
Nickel(II) acetate tetrahydrate (0.2506 g, 1.00 mmol) was dissolved in 2 mL of distilled H 2 O to afford a blue-green solution. To this solution was added 5 mL of neat 4-picoline to give a deep blue reaction mixture. No precipitate formed during initial mixing. Upon standing at room temperature for about 6 weeks, small dark blue crystals formed which were carefully collected by suction filtration (0.3247 g). The crystals were stored in a refrigerator to preserve their quality. The substance formed was identified as Ni(acetate) 2 (H 2 O) 2 (4-picoline) 2 by single crystal x-ray diffraction.

Synthesis of [Ni(H
The most efficient procedure entails a mechanochemical process. Ni(BF 4 ) 2 ·6H 2 O (0.4006 g, 1.18 mmol) was mixed with neat pyrazine (0.1882 g, 2.34 mmol) and the reagents ground together using a mortar and pestle. During the grinding procedure, a moist blue solid was obtained, transferred to a stainless steel autoclave, and placed in an oven set at 110 • C. The reaction mixture was held at 110 • C for 3 days and then removed from the oven to slowly cool back to room temperature after which the vessel was opened. A homogeneous pale purple powder was obtained. Under magnification, no indication for unreacted starting materials was evident. A single crystal of ∼ 50 × 50 × 50 µm 3 in size was mounted on a glass fibre through use of Paratone N oil. Data collection was carried out at the 15ID ChemMatCARS beamline of the Advanced Photon Source, Argonne National Laboratory. The data were collected at 100 K with a Bruker 6000 CCD detector. The structural refinement was carried out through use of Bruker APEX3 software v2015.5-2 [1]. All non-hydrogen atoms were refined with anisotropic displacement parameters. Hydrogen atoms, other than those associated with the water molecule, were placed in calculated positions and refined isotropically using a riding model. Hydrogen atoms associated with the water molecules were located on a difference map and their location restrained by the geometry of an isolated water molecule. The Ni(II) local environment is shown in Fig. 5 of the main text and the stacking of the molecular units is shown in Fig. S1.

1.2.2.
Single crystal x-ray diffraction of Ni(H 2 O) 2 (acetate) 2 (4-picoline) 2 A blue prism-shaped crystal 0.33 × 0.23 × 0.15 mm 3 in size was mounted on a glass fibre with traces of viscous oil and then transferred to a Nonius KappaCCD diffractometer equipped with Mo K α radiation (λ = 0.71073Å). Ten frames of data were collected at 150(1) K with an oscillation range of 1 deg/frame and an exposure time of 20 sec/frame [2]. Indexing and unit cell refinement based on all observed reflections from those ten frames, indicated an orthorhombic P lattice. A total of 4088 reflections (θ max = 27.484 • ) were indexed, integrated and corrected for Lorentz, polarization and absorption effects using DENZO-SMN and SCALEPAC [3] . Post refinement of the unit cell gave a = 8.8996(3)Å, b = 12.3995(4)Å, c = 17.6516(7)Å, and V = 1947.86(12)Å 3 . Axial photographs and systematic absences were consistent with the compound having crystallized in the orthorhombic space group P cab.
The structure was solved using SIR 97 [4]. All of the non-hydrogen atoms were refined with anisotropic displacement coefficients. Hydrogen atoms were assigned isotropic displacement coefficients U (H) = 1.2U (C) or 1.5U (C methyl ), and their coordinates were allowed to ride on their respective carbons using SHELXL [5]. The weighting scheme employed was w = 1/[σ 2 (F 2 0 ) + (0.0635P ) 2 + 3.6101P ] where P = (F 2 0 + 2F 2 c )/3. The refinement converged to R1 = 0.0513, wR2 = 0.1356, and S = 1.145 for 1515 reflections with 1 > 2σ(I), and R1 = 0.0844, wR2 = 0.1515, and S = 1.145 for 2225 unique reflections and 121 parameters. ‡ The maximum ∆/σ in the final cycle of the least-squares was 0, and the residual peaks on the final difference-Fourier map ranged from −0.741 to 0.951 e/Å 3 . Scattering factors were taken from the International Tables for Crystallography, Volume C [6,7]. The resulting Ni(II) local environment is shown in Fig. 5 of the main text and the stacking of the molecular units is shown in Fig. S1.

Synchrotron x-ray powder diffraction of [Ni(H
High-resolution room temperature data were collected using beamline 11-BM at the Advanced Photon Source (APS) and 0.413841Å x-rays. Discrete detectors are scanned at a speed of 0.01 • /s over 34 • in 2θ-range with data points collected every 0.001 • . Data were collected while continually scanning the diffractometer 2θ-arm. Peak indexing was ‡ R1 = (||F 0 | − |F c ||)/ |F 0 |, wR2 = (w(F 2 0 − F 2 c ) 2 )/ (F 2 0 ) 2 , and S = goodness-of-fit on where n is the number of reflections and p is the number of parameters refined. performed using TOPAS Academic. Hydrogen atoms were placed in calculated positions using a fixed C-H bond length of 0.95Å. Each water hydrogen atom is positionally disordered over four equivalent sites.

Magnetometry
Pulsed-field magnetization measurements were performed at the National High Magnetic Field Laboratory in Los Alamos; fields of up to 40 T with typical rise times ≈ 10 ms were used. Powdered samples are mounted in 1.3 mm diameter PCTFE ampules (inner diameter 1.0 mm) that can be moved into and out of a 1500-turn, 1.5 mm bore, 1.5 mm long compensated-coil susceptometer, constructed from 50 gauge high-purity copper wire [8]. When the sample is within the coil and the field pulsed the voltage induced in the coil is proportional to the rate of change of magnetization with time, (dM/dt). Accurate values of the magnetization are obtained by numerical integration of the signal with respect to time, followed by subtraction of the integrated signal recorded using an empty coil under the same conditions [8]. The magnetic field is measured via the signal induced within a coaxial 10-turn coil and calibrated via observation of de Haasvan Alphen oscillations arising from the copper coils of the susceptometer [8]. The susceptometer is placed inside a 3 He cryostat, which can attain temperatures as low as 500 mK.
Pulsed-field data were calibrated using measurements of the sample moment, m, made in a Quantum Design MPMS 7 T SQUID magnetometer fitted with DC transport mode. Powdered samples of a known mass are loaded into a gelatin capsule, placed in a plastic drinking straw and attached to a sample rod. The samples were cooled in zero-field to the measured temperature, then the field was increased in increments between µ 0 H = 0 and 7 T and the field stabilised before each measurement was taken. The pulsed-field measurements were then scaled onto the SQUID data. The same experimental set up was used to obtain DC susceptibility measurements, with the sample zero-field cooled to 1.8 K. The field was then set to 0.1 T and temperature dependent data was taken on warming to 300 K. The molar magnetic susceptibility (χ mol ) is deduced from this measurement using χ mol = m/nH, where n is the number of moles of sample.

Heat capacity
A Quantum Design 9 T PPMS instrument was used to measure temperature-dependent data in fields between 0 and 9 T. The sample was pressed into a thin pellet and mounted on a sapphire sample platform using a small amount of Apiezon-N grease to create a good thermal contact. The platform houses an electric heater and Cernox thermometer and is connected to the thermal bath by thermally conducting wires. The sample was cooled to 1.8 K in a high vacuum and the heat capacity was measured using a thermal relaxation technique. The sample platform and grease addenda was measured over the same field and temperature range beforehand and subtracted from the measurement to obtain the heat capacity of the sample.
To extract the magnetic heat capacity, the lattice contribution is determined by fitting to the high-temperature data. We use the following model [9], which has, as adjustable parameters, the characteristic amplitudes, A i (J K −1 mol −1 ), and temperatures, θ i (K), of the Debye (i = D) and Einstein (i = E) phonon modes.
where x D = θ D /T . The first term is the response from a Debye mode and the second term accounts for n Einstein modes.

Electron-spin resonance
High-field, high-frequency ESR spectra of powdered samples at temperatures ranging from ≈ 3 to 80 K were recorded on a home-built spectrometer at the EMR facility of the National High Magnetic Field Laboratory with the microwave frequencies 52-626 GHz. The instrument is a transmission-type device and uses no resonant cavity. The microwaves were generated by a phase-locked Virginia Diodes source, generating frequency of 13±1 GHz, and equipped with a cascade of frequency multipliers to generate higher harmonic frequencies. A superconducting magnet (Oxford Instruments) capable of reaching a field of 17 T was employed.

Muon-spin relaxation
Zero-field muon spin relaxation (µ + SR) measurements on [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 were carried out on polycrystalline samples using the GPS spectrometer at the Swiss Muon Source (SµS), Paul Scherrer Institut, Switzerland. For the measurements, samples were packed in Ag foil envelopes (foil thickness 25 µm), attached to a silver fork [10] and mounted in a He cryostat. In a µ + SR experiment [10] spin-polarized positive muons are stopped in a target sample, where the muon usually occupies an interstitial position in the crystal. The observed property in the experiment is the time evolution of the muon spin polarization, the behaviour of which depends on the local magnetic field at the muon site. Each muon decays, with an average lifetime of 2.2 µs, into two neutrinos and a positron, the latter particle being emitted preferentially along the instantaneous direction of the muon spin. Recording the time dependence of the positron emission directions therefore allows the determination of the spin-polarization of the ensemble of muons. In our experiments positrons are detected by detectors placed forward (F) and backward (B) of the initial muon polarization direction. Histograms N F (t) and N B (t) record the number of positrons detected in the two detectors as a function of time following the muon implantation. The quantity of interest is the decay positron asymmetry function, defined as where α exp is an experimental calibration constant. A(t) is proportional to the spin polarization of the muon ensemble. The asymmetry A(t) was initially fitted to a model with three components (two oscillatory and one non-oscillatory) across all temperatures: where A rel corresponds to the total relaxing amplitude, p 1 and p 2 are the weights of the oscillatory components and λ 1 and λ 2 are the respective relaxation rates. The constant term A bg accounts for the non-relaxing contribution from those muons that stop at the sample holder/cryostat tail. Throughout the fitting procedure, p 1 , p 2 and the ratio between the oscillatory frequencies α were fixed at 0.66, 0.10 and 0.29 respectively. Fitted values of λ 1 and f are plotted against temperature in the top left panel of Figure S2. A sharp change is observed in both λ 1 and f at around T = 3.2 K, strongly indicative of a magnetic phase transition.
For T < 3.2 K, the oscillatory components are therefore characteristic of quasistatic long-range magnetic ordering, with the two frequencies corresponding to two, magnetically distinct, muon stopping sites. The frequencies fall as the temperature nears the transition at T N [11]. However, unlike the usual case where the frequency vanishes above T N , in [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 the oscillation in A(t) persists up to T = 20 K. Therefore, although Model I describes the data at all measured temperatures, this behaviour of the parameters revealed at T > 3.2 K suggests that the model is not applicable in this regime.
A second model was therefore required to explain the T > 3.2 K behaviour in [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 . Rather than the muons being subject to a quasi-static magnetic field, here it is probable that the persistent oscillations in asymmetry arises from a muon-fluorine entanglement state, often observed in this class of materials in the disordered regime [12]. This occurs in fluorine-containing materials when the disordered electronic spins fluctuate rapidly on the muon timescale, motionally narrowing them from the spectrum. The muon is then relaxed by its dipolar interaction with fluorine nuclear spins [12]. Thus, in the temperature range where λ 1 and f are found to be constant using Model I, the asymmetry is better fitted to the function A(t) = A rel p 1 e −λ F−µ t D z (ω F−µ , t) + p 2 e −Λt + A bg , where p 1 and p 2 were fixed 0.43 and 0.57 respectively. The function describing the muon-fluorine dipolar interaction D z (ω F−µ , t) [12], models the interaction of the muon with a single I = 1/2 fluorine nucleus. Again, the component A bg accounts for the relaxing contribution from the muons stopping at the sample holder/cryostat tail. The extracted relaxation rates λ F−µ and frequency ω F−µ are shown in the bottom left panel of Figure S2. We see that λ F−µ diverges as temperature approaches the magnetic transition, as often observed in similar systems [13] and providing additional evidence for our interpretation. The value of ω corresponds to a muon site with a F-µ + separation of 0.11 nm, which is also typical of muon stopping sites in this class of material [13].

Neutron diffraction
Magnetic diffraction patterns were recorded on the WISH diffractometer (ISIS, Rutherford Appleton Laboratory, UK) [14]. A deuterated polycrystalline sample, [Ni(D 2 O)(d 4 -pyz) 2 ]( 11 BF 4 ) 2 was loaded into a cylindrical vanadium can and placed in an Oxford Instruments cryostat with a base temperature of 1.5 K. Diffraction data were collected over the temperature interval 1.5-10 K, with long counting times (3 hrs) at 1.5 K and 10 K. Intermediate temperature points were measured with an exposure time of 0.5 hrs. Scattered neutron intensity (circles) at 1.5 K as a function of d-spacing fits shown in the right-hand panel of Figure S2.
400 different field orientations with respect to the z-axis, the average C mag (T ) at a given field and temperature is estimated from The field dependence of the heat capacity is examined by repeating the calculation at different fixed values of H. The results for the easy-axis and easy-plane cases with E = 0 are shown in Fig. 4 of the main text.

Electronic structure calculations
Our computational approach involves determining the energy differences between ordered spin states that differ by a number of reversed spins [16,17]. An underlying physical assumption, therefore, is that upon magnetically ordering, the magnetic structure is constrained to be collinear. In other words, the value of the magnetic exchange derived assumes that nearest neighbor spins S at sites i and j obey S i · S j = ±1. If this is not the case, then the error in this assumption is absorbed into the value of the exchange constant that is derived. More specifically, we map the magnetic centres of the system to an Ising Hamiltonian where the sum is over neighbours and I is the Ising exchange energy. In [Ni(H 2 O) 2 (pyrazine) 2 ](BF 4 ) 2 the spins are (nominally) located on the Ni atoms. There are two different pairs of Ni-Ni interactions. The first is within x − y planes with Ni-Ni distance of 7.012Å (coupling I 1 ) and the second aligned along z where Ni atoms are separated by 7.436Å (coupling I 2 ). Each Ni interacts with four neighbours with strength I 1 within the x − y plane, and two neighbours along z with strength I 2 . (Here we count the couplings per Ni-Ni bond, not per Ni atom.) Therefore the I-coupling model simplifies to a nearest x − y neighbour and nearest z neighbour energy: (S10) The spin structure of the lowest energy state is shown in the left panel of Figure S3. It forms an antiferromagnetic state within each of the x − y planes and offset by (1/2, 1/2) in going from the z = 1/4 to z = 3/4 plane.
To evaluate the intraplane I 1 coupling and the interplane (next nearest neighbour) I 2 coupling, differences in energy of the spin configurations are considered. Details are as follows: the second lowest energy state is shown in the middle panel of Fig. S3 which is 15.49 meV higher in energy than the ground state. A single spin is flipped on the central Ni atom in one plane which means that there is an I 1 contribution associated Figure S3. Left: The lowest energy spin state; blue shows the spin density of one spin channel, yellow the other. Middle: Spin flip within one x − y plane with respect to lowest energy spin configuration. Right: The third lowest energy configuration is shown, which is antiferromagnetically stacked ferromagnetic planes.
with each of the four nearest neighbours in the x − y plane. Additionally this atom has two I 2 neighbours along z, hence 2S (4I 1 + 2I 2 ) = 15.49 meV, (S11) where S = 1 for these Ni spins. The next highest energy state is shown in the righthand panel of Fig. S3, which is 20.52 meV above the ground state. Comparing what is flipped, the central atom on both of the x − y planes have changed. Each of these atoms have 4 I 1 (x − y plane) neighbours. However since both atoms have flipped then the I 2 interaction (z-direction) is unchanged, it is still antiferromagnetic. Therefore we have 2S(8I 1 + 0I 2 ) = 20.52 meV. (S12) Solving equations (S11) and (S12) gives I 1 = 1.283 meV and I 2 = 1.308 meV. In order to convert these Ising couplings to Heisenberg couplings, we employ the method from Ref. [17], where the estimated Heisenberg coupling J is related to the Ising coupling I via J = I 2S , where 2S is the spin flip factor again, that we saw above. We conclude that the Heisenberg couplings are J 1 = 0.642 meV and J 2 = 0.652 meV, so the exchange is fairly isotropic, with magnitude around J ≈ 8 K.