Helmholtz-Zentrum Dresden-Rossendorf (HZDR) Extremely well isolated two-dimensional spin-1/2 antiferromagnetic Heisenberg layers with a small exchange coupling in the molecular- based magnet CuPOF

Extremely well isolated two-dimensional spin-1/2 antiferromagnetic Heisenberg layers with a small exchange coupling in the molecularbased magnet CuPOF Opherden, D.; Nizar, N.; Richardson, K.; Monroe, J. C.; Turnbull, M. M.; Polson, M.; Vela, S.; Blackmore, W. J. A.; Goddard, P. A.; Singleton, J.; Choi, E. S.; Xiao, F.; Williams, R. C.; Lancaster, T.; Pratt, F. L.; Blundell, S. J.; Scurschii, I.; Uhlarz, M.; Ponomaryov, O.; Zvyagin, S.; Wosnitza, J.; Baenitz, M.; Heinmaa, I.; Stern, R.; Kühne, H.; Landee, C. P.;


I. INTRODUCTION
The study of critical phenomena, related to phase transitions between exotic ground states that emerge from complex underlying electronic correlations, is a subject of high importance in the research of low-dimensional magnetism. As has been established by extensive theoretical work, in contrast to the cases of one-and three-dimensional magnetic lattices, the critical behavior and ground-state properties of two-dimensional (2D) lattice systems are strongly dependent on the symmetry of the interactions between the magnetic moments, e.g., Ising, XY , and Heisenberg types [1]. For instance, in accordance with Onsager's exact solution [2], the 2D spin- 1 2 Ising antiferromagnet undergoes a Néel-type transition to long-range order at T N = 1.06J/k B [3], where J is the exchange strength between neighboring magnetic moments. In contrast, thermal fluctuations in the 2D quantum Heisenberg antiferromagnet (2D QHAF) prevent long-range order at any finite temperature, as was rigorously proven by Mermin and Wagner [4]. The 2D planar, or XY , magnetic lattice does not exhibit Néel-type order, despite the divergence of the susceptibility at finite temperatures. Instead, an unusual topological order with characteristic algebraic decay of the spin correlations was proposed by Berezinskii [5] and Kosterlitz and Thouless [6]. Here, the formation of bound pairs of topological excitations, called vortex and antivortex states, characterizes the ordered state, where the unbinding of the vortex-antivortex pairs constitutes the Berezinskii-Kosterlitz-Thouless (BKT) transition, which occurs at T BKT = 0.353J/k B for the S = 1 2 case [7]. Advancing the experimental research of phase transitions and critical phenomena relies on the availability of welldefined model systems. In particular, a targeted synthesis and characterization of materials with interaction parameters that closely approximate those of theoretical model systems and yield energy scales of temperature and magnetic field that are accessible by existing experimental infrastructures is required. The effective Hamiltonian to describe a quasi-2D QHAF in an applied magnetic field is where J and J are the intra-and interlayer exchange couplings, and scales the deviation from an ideal Heisenberg interaction to an easy-axis or easy-plane characteristic. The first sum in Eq. (1) is taken over nearest-neighbor (NN) spins in the quasi-2D planes, the second summation is over NN spins in adjacent layers, and the Zeeman term applies to all moments. For external magnetic fields smaller than the anisotropy field H A and T J/k B , the formation of anisotropic magnetic correlations is driven by the intrinsic anisotropy . Conversely, for fields exceeding H A = ×H sat , where H sat denotes the saturation field, the anisotropy of magnetic correlations is mainly determined by the strength and direction of the applied magnetic field. In the Hamiltonian (1), a positive J corresponds to an antiferromagnetic exchange, and 0 < 1 describes the degree of easy-plane, or XY , anisotropy.
Whereas the critical phenomena of several systems were mapped onto the BKT theory, such as the superfluid transition of thin 4 He films [28], the solid-on-solid model [29], the two-dimensional melting [30], the superconducting transition of Josephson junctions [31], the collision physics of 2D atomic hydrogen [32], the loss of quasicoherence of a trapped degenerate quantum gas of rubidium atoms [33], and the magnetic van der Waals antiferromagnet in the atomically thin limit [34], a very clean bulk realization of the square 2D XY Heisenberg lattice is lacking up to now. There are two main obstacles for the experimental realization of a 2D XY antiferromagnet. One challenge relates to the unavoidable existence of finite interlayer interactions. In case of the 2D Heisenberg magnets, even an infinitesimal interlayer interaction is sufficient to perturb the critical behavior of the system and stabilize a conventional type of Néel order below a transition temperature T N [35]. Second, a weak intrinsic easyplane anisotropy constraints the temperature range of XY -type correlations. For the case of crystalline magnetic lattices of Cu 2+ ions, the exchange coupling between spin moments as well as the single-ion properties are almost isotropic.
Another approach to realize a well-defined investigation of the magnetic correlations in low-dimensional spin systems is by tuning the Zeeman terms of the effective Hamiltonian with the application of an external magnetic field. This gives rise to the unique possibility of probing a quasi-2D spin system with well-defined XY anisotropy in the experiment. It was shown by quantum Monte Carlo calculations that the application of a magnetic field to an isotropic 2D Heisenberg antiferromagnet can be mapped onto the anisotropic 2D XY model in zero field, where the strength of the spin-exchange anisotropy can be continuously tuned by the external field [36,37]. This overall context has triggered extensive and ongoing efforts to synthesize novel quasi-2D Heisenberg model materials with highly isolated layers and relatively small antiferromagnetic interaction energies, allowing for the investigation of field-induced effects in moderate applied magnetic fields. Molecular-based materials consisting of 3d transition-metal ions, such as copper, embedded into an organic matrix are of particular interest. The combination of different ligands gives the opportunity to engineer a wide range of materials with well-defined magnetic properties [38][39][40][41][42][43][44][45][46][47][48].
In this paper, we report a comprehensive investigation of the magnetic properties of the newly synthesized compound [Cu(pz) 2 (2-HOpy) 2 ](PF 6 ) 2 (CuPOF), where pz = C 4 H 4 N 2 and 2-HOpy = C 5 H 4 NHO. The material is characterized by magnetometry, ESR, specific heat, μ + SR, and NMR. CuPOF is shown to be a very clean realization of the square 2D spin- 1 2 Heisenberg lattice with moderate intralayer antiferromagnetic NN exchange, J/k B = 6.80(5) K, and highly isolated magnetic layers, with J /J 10 −4 . A weak intrinsic easy-plane anisotropy, revealed by bulk magnetometry, yields a temperature-driven crossover of the spin-correlation anisotropy from isotropic Heisenberg to anisotropic XY -type behavior, which, under the influence of a finite interlayer coupling J , constitutes a driving mechanism for a transition to long-range commensurate order. A strong increase of the transition temperature upon application of a magnetic field from 1.38(2) K at zero field to 2.8 K at 7 T is caused by the field-driven increase of the anisotropy of spin correlations.

II. EXPERIMENTAL
The compound [Cu(pz) 2 (2-HOpy) 2 ](PF 6 ) 2 was synthesized by using conventional solution chemistry techniques, described in detail in the Supplemental Material (SM) [49]. Slow evaporation of methanol solutions of the product, [Cu(pz) 2 (2-HOpy) 2 ](PF 6 ) 2 , (pz = pyrazine, 2-HOpy = 2pyridone) produces thin, rectangular, green plates, typically 3 to 5 mm on a side and about 1 mm thick. The crystals extinguish well under polarized light and are dichroic. 064431-2 X-ray data were obtained by using an Agilent Technologies Gemini Eos CCX-ray diffractometer using Cu-Kα radiation (λ = 1.5418 Å) with ω scans using CrysAlisPro software [50] refined cell parameters and SCALE3 ABSPACK [51] scaling algorithm defined absorption corrections. Data were collected at 120 K. SHELXS97 [52] was used to solve the structures, which were refined via least-squares analysis using SHELXL-2016 [53]. All non-hydrogen atoms were refined anisotropically. Hydrogen atoms bonded to nitrogen atoms were located in the difference Fourier maps and their positions refined with fixed isotropic thermal parameters. The remaining hydrogen atoms were geometrically located and refined by using a riding model with fixed isotropic thermal parameters. The structure has been deposited with the Cambridge Crystallographic Data Centre (CCDC) (1553982). A Bruker D8 powder x-ray diffractometer was used to confirm the purity and phase of powdered samples prior to magnetic measurements.
Measurements of magnetic bulk properties between 1.8 and 310 K were carried out by using a Quantum Design MPMS-XL SQUID magnetometer with a 5 T magnet, as well as a vibrating sample magnetometer (VSM) Quantum Design PPMS with a 14 T magnet. Corrections were made to the data for the background signal of the sample holder, as well as diamagnetic contributions. Studies below 2 K were performed by employing a 3 He cooling stage.
The high-field magnetization of CuPOF single crystals in pulsed fields up to 35 T and at temperatures of 0.37 and 1.4 K were performed at the Dresden High Magnetic Field Laboratory (HLD) at the Helmholtz-Zentrum Dresden-Rossendorf. Additional measurements of a polycrystalline sample in DC fields up to 35 T were done by using a vibrating sample magnetometer at the National High Magnetic Field Laboratory (NHMFL) in Tallahassee, as well as in pulsed fields up to 25 T at the NHMFL facility in Los Alamos. The results are fully consistent with the data from the HLD and can be found in the SM [49].
Room-temperature electron-spin resonance (ESR) studies were performed on polycrystalline material and single crystals of CuPOF, using a commercially available X-Band Bruker ESR spectrometer operating at 9.8 GHz at Clark University. EASYSPIN [54] was used to determine the g factors and linewidths. Additional ESR measurements of angulardependent spectra, as well as temperature-dependent spectra for field parallel to the c axis, were performed at the HLD between 3 and 300 K at 9.4 GHz, employing an X-band Bruker ELEYSYS E500 ESR spectrometer. The obtained values of the anisotropic g factor are consistent with the results measured at Clark University and are presented in the SM [49]. High-frequency ESR measurements along the crystallographic c axis at 1.5 K and fields up to 16 T were performed at the HLD by using a home-built transmission-type tunable-frequency ESR spectrometer, similar to that described in Ref. [55], with a probe in Faraday configuration. These results can also be found in the SM [49].
Heat-capacity measurements between 1.8 and 300 K were performed by using a Quantum Design PPMS system. Furthermore, a 3 He insert was used to record the heat capacity at temperatures down to 0.4 K. Powdered samples with masses of 1.065(5) and 1.832(5) mg, for measurements at 4 He and 3 He temperatures, respectively, were pressed into pellets and attached to the sample platforms by using Apiezon N grease. The addenda was determined from measurements with an empty sample holder and subtracted from the data to obtain the heat capacity of the sample.
Zero-field muon-spin relaxation (μ + SR) measurements on a polycrystalline sample were carried out by using the EMU spectrometer at the ISIS facility at the Rutherford Appleton Laboratory. The sample was mounted on a Ag backing plate and covered with a 12.5-μm-thick Ag foil mask before being inserted into a dilution refrigerator. Further technical details are provided in the SM [49]. 1 H nuclear magnetic resonance (NMR) spectra were recorded using a commercial solid-state spectrometer. A standard Hahn spin-echo pulse sequence with stepped sweep of the carrier frequency was employed to record the broadbandwidth spectra. The NMR probe was equipped with a single-axis goniometer for precise orientation of the magnetic field parallel to the crystallographic c axis. The measurements at 1.6 K and above were performed in a 8 T high-resolution magnet equipped with a 4 He flow cryostat.

III. RESULTS
a. Crystal structure. [Cu(pz) 2 (2-HOpy) 2 ](PF 6 ) 2 crystallizes in the orthorhombic space group Cmca. The asymmetric crystallographic unit comprises one Cu 2+ ion, two half-pyrazine molecules, one 2-pyridone molecule and two PF − 6 ions. The local coordination sphere of the Cu 2+ ion is presented in Fig. S1 in the SM [49]. The Cu 2+ ion sits on a twofold axis (parallel to b), one coordinated, dissymmetric (N11/N14) pyrazine sits on the same twofold axis while the other pyrazine molecule (N21) lies across a mirror plane normal to the a axis. The two PF − 6 ions also lie on mirror planes such that there are five independent fluoride ions in each. The Cu 2+ coordination sphere exhibits a classic distorted Jahn-Teller octahedral geometry with four bridging pyrazine molecules in the equatorial plane [Cu-N = 2.05(1) Å] and elongated Cu-O bonds [2.285(1) Å] in the axial sites. The copper ion and all four bound nitrogens are coplanar as required by symmetry. The Cu-O1 bond is nearly perpendicular to the plane, making an angle of 1.0 • with the normal to the plane. Crystal data and structure refinement parameters, as well as selected bond lengths and angles of CuPOF, are presented in Tables S1 and S2 in the SM [49].
The bridging pyrazine ligands link the Cu 2+ ions into a nearly square layer, see Fig. 1(a). The pyrazine rings exhibit a propeller twist relative to the Cu plane, with the N11 rings canted 66.9 • and the N21 rings canted 53.1 • out of the plane. This results in slightly different Cu-Cu distances of 6.676 Å through the N11 ring and 6.680 Å through the N21 ring. The layers are further separated by the PF − 6 anions, which are located in the pockets between the 2-pyridone ligands, resulting in a minimum Cu-Cu distance in adjacent layers of 13.097 Å.
The structure of CuPOF possesses a hidden canting, see Fig. 1(b). The Cu-pyrazine layers lie in the ab plane (into the page and horizontal) while the Cu-oxygen bonds are nearly parallel to the c axis but canted by ±1.0 • towards b. Within each layer, the canting is in the same direction but adjacent layers are canted in the opposite direction. A similar canting is seen in the orientation of the pyridone rings; they are alternately tilted by 8.7 • away from the normal to the Cu-pyrazine planes. The chemical equivalence of the 2-pyridone molecules is confirmed by our 13 C magic-angle spinning (MAS) NMR spectroscopy (see Fig. S11 in the SM [49]).
b. Magnetometry. The magnetic susceptibility of a polycrystalline sample of CuPOF is shown in Fig. 2 and yields a rounded maximum near 6.8 K. At higher temperatures (T > 40 K), the data are well described by a Curie-Weiss law with a Curie constant of 0.440(5) emuG −1 mol −1 K and a Curie- Weiss temperature CW = −5.2(6) K, indicating a small antiferromagnetic interaction. Accordingly, the data were compared with the susceptibility of a 2D S = 1 2 Heisenberg antiferromagnet [56] in which the Curie constant, exchange strength J, and small paramagnetic impurity fraction were adjusted. The best agreement between data and the model calculations, denoted by the red curve in Fig. 2, is obtained with C = 0.445(5) emuG −1 mol −1 K, J/k B = 6.80(5) K, and 1.1(1)% paramagnetic contribution. The agreement between the Curie constants from the Curie-Weiss model and the 2D QHAF model is excellent, with their values corresponding to an average g value of 2.17. The average g factor has been determined at room temperature by using ESR and was found to be g = 2.15.
The existence of two distinct pyrazine molecules in the unit cell allows for the possibility of a rectangular magnetic lattice in which the exchange strengths (J and αJ, 0 α 1) along the a and b axes are different. This possibility has been tested by comparing the susceptibility data to the susceptibilities of a rectangular 2D QHAF [57,58] for various values of J and α. The square-lattice (α = 1) case gives by far the best fit and it is possible to rule out any rectangular contribution with α < 0.96.
The low-temperature static susceptibility of a single crystal of CuPOF at 0.1 T is shown in the inset of Fig. 3. The out-of-plane DC susceptibility (H c) has a minimum at 1.86(5) K, whereas, for a field applied in the ab plane (H ⊥ c), the DC susceptibility steadily decreases to the lowest measured temperature. The minimum in the out-of-plane susceptibility at T co indicates the presence of XY anisotropy, where, with decreasing temperature, the correlation of spin moments crosses from isotropic to easy-plane behavior. The temperature-dependent out-of-plane susceptibility at different magnetic fields is presented in the main panel of Fig. 3. With increasing magnetic field, the broad minimum of the static susceptibility shifts to higher temperatures, as indicated by the triangles. At fields above around 4 T, the minimum broadens and cannot be observed anymore.
The magnetization of a single crystal of CuPOF has been measured up to 1 T at 0.5 K, see the insets of Figs. 4(a) and 4(b). For a field parallel to the c axis, the magnetization monotonically increases over that range. In contrast, when the field is applied in the ab plane (H ⊥ c), a small spin-flop anomaly is observed at around 0.36(1) T.
The relative magnetization of single crystals of CuPOF was determined at several temperatures for fields up to 35 T both parallel to the c axis (H c) and within the layers (H ⊥ c).
The experimental values at 0.37 K are presented in Figs. 4(a) and 4(b). The absolute values of the magnetization were obtained from a direct comparison with the magnetization recorded at 4 K for H c (at 2 K for H ⊥ c) up to 14 T, by using a VSM magnetometer; see Fig. S3 in the SM [49]. For both magnetic-field orientations, the saturated moment is about 1 μ B per formula unit, as expected for the Cu 2+ ion. The experimental data (symbols) are compared with quantum Monte Carlo (QMC) simulations (lines) for the 2D QHAF at the relative temperature k B T /J = 0.05. Excellent agreement between the experimental data and the QMC simulations is found for both field orientations. Over the full field range, the deviation between experimental and theoretical data is below ±1% for H c and below ±2% for H ⊥ c. The inand out-of-plane saturation fields were determined as 17.6 and 19.6 T for H c and H ⊥ c, respectively. In the meanfield approximation, the saturation field H sat is defined by the exchange strength J, where z = 4 is the number of nearest-neighbor magnetic moments. The exchange strength can be computed from Eq. by using the experimentally determined saturation fields and the respective g values g c = 2.29(1) at 1.5 K and g ab = 2.070(7) at room temperature; see the ESR section below. The obtained values, J/k B = 6.75(5) and 6.78(5) K, are in very good agreement with each other and with those determined from the susceptibility and specific-heat measurements.
c. Electron-spin resonance. The anisotropy of the roomtemperature ESR spectrum of a single-crystalline CuPOF sample was investigated at 9.8 GHz. A single ESR line was observed for each field orientation; see Fig. S6 in the SM [49]. The resonance field, used to extract the electronic g factor and the spectral linewidth, was determined by modeling a nearly Lorentzian function to the experimentally obtained spectra [54,59,60]. The angular dependence of the g factor in the ac, bc, and ab planes indicates a strong planar-like anisotropy ( Fig. S7 in the SM [49]), with g a = 2.073(2), g b = 2.066(4), and g c = 2.298(2). 1 1 The average values of the g factors from two independent ESR measurements on the single-crystalline CuPOF samples are presented here; for details see the SM [49]. The temperature-dependent integrated ESR-line intensity, scaling as the bulk susceptibility of the sample [61] is shown in Fig. 2 for temperatures between 3 and 50 K with the field applied along the c axis. Analogously to the modeling of the DC susceptibility, the integrated ESR intensity was modeled, varying only the Curie constant, with fixed exchange coupling of J/k B = 6.80(5) K and a paramagnetic impurity percentage of 1.1(1)% (black dash-dotted line in Fig. 2).
High-frequency ESR spectroscopy at 1.5 K in the frequency range between 52 and 500 GHz revealed a single resonance mode with a linear frequency-field dependence ( Fig. S9 in the SM [49]), yielding g c = 2.29(1).
d. Specific heat. The specific heat was measured between 0.4 and 300 K. The data increase smoothly from about 0.2 J mol −1 K −1 at 1 K, approaching 640 J mol −1 K −1 at 300 K. No sharp anomalies, corresponding to structural changes or ordering transitions, were observed in this range (Fig. S10 in the SM [49]). The data between 1 and 9 K are shown in Fig. 5(a) (black open circles) revealing a broad hump exceeding the normal phonon contribution. The data in this temperature range were analyzed as a sum of the magnetic specific heat of a 2D QHAF and a phononic contribution. The low-temperature lattice contribution to the specific heat, stem-ming from a complex phononic spectrum in the molecularbased material CuPOF, is best approximated by choosing a three-term polynomial with C pho = AT 3 + BT 5 + CT 7 . The magnetic specific heat was represented as a ratio of polynomials, similar to the approach used in a previous study [62], but is based on recent QMC simulations of the magnetic specific heat [63] that extended to lower relative temperatures. The range of validity for the square lattice is 0.15 k B T /J 5.0 (see SM and Table S4 therein [49]).
The resulting best fits to C pho and C mag are shown in Fig. 5(a). The sum of the two individual contributions appears as the solid black line and shows very good agreement with the experimental data. The modeling parameters are J/k B = 6.75 (5) As seen in Fig. 5(a), the magnetic contribution dominates at lower temperatures. Including data at higher temperatures in the fitting process did not change the value of the exchange strength, since C mag is a rapidly decreasing fraction of the total specific heat. Figure 5(b) shows a direct comparison of the calculated magnetic specific heat (red line) to the difference between the total and estimated lattice specific heat (black triangles). The obtained exchange strength of J/k B = 6.75 (5) K is in excellent agreement with that obtained from the susceptibility and magnetization results. Similar to the analysis of the bulk susceptibility, the specific-heat data were investigated in terms of the possibility of a rectangular magnetic lattice [63]. The results are consistent with a magnetic square lattice.
e. Muon-spin relaxation. Zero-field muon-spin relaxation (μ + SR) spectra [64] for CuPOF are shown in Fig. 6(a). Spontaneous oscillations in the asymmetry function A(t ), which is proportional to the spin polarization of the muon ensemble, see SM [49], were observed at low temperatures. An oscillating behavior of the muon-spin polarization is characteristic of the presence of quasistatic long-range magnetic order (LRO). The local magnetic field that results from LRO causes those muons with spin perpendicular to the local field to precess coherently at the frequency ν i , where ν i is proportional to the magnitude of the local field B i at the ith muon site. Changes in the spectra are observed in measurements across the temperature range 0.1 < T < 2 K, which is parametrized by modeling the asymmetry A(t ) with the relaxation function where A 0 is the initial asymmetry at t = 0, the parameters ν and λ 1 are the precession frequency and relaxation rate of the oscillatory component, respectively. The parameters A bg and λ 2 account for the background muons and those with spin initially parallel to the local magnetic field. The evolution of the parameters A 0 , λ 1 , and ν are shown in Figs. 6(b)-6(d). The precession frequency, see Fig. 6(d), which is proportional to the magnitude of internal magnetic field probed by the muon spins, shows a monotonic decrease from base temperature up to 1.4 K, but starts to rise again before becoming roughly constant above 1.5 K. A sharp maximum of A 0 and discontinuity of λ 1 are observed just below T = 1.4 K, coincident with the minimum of ν.
The behavior of ν, together with the peak of A 0 and the discontinuity of λ 1 , suggest that CuPOF undergoes a magnetic phase transition around 1.4 K. However, unlike many magnetic systems where the precession vanishes above T N , the asymmetry still appears to show oscillatory behavior. Such an oscillatory signal is common in similar materials containing fluorine nuclei in the paramagnetic phase [65]. It arises because the electronic moments fluctuate outside of the muon time window and are consequently removed from the spectrum, leaving the muon sensitive to the nuclear magnetic moments. In fluorinated materials, there are frequently muon sites where the positive muon sits close to the electronegative fluorine and enters a dipole-dipole coupled entangled state, leading to heavily damped, low-frequency oscillations [65]. A similar scenario for CuPOF is suggested, allowing identification of the antiferromagnetic transition temperature with the discontinuities in the modeled parameters.
To study the critical behavior close to the phase transition, we exclude the data above 1.4 K. To extract the critical exponent of the order parameter, the temperature dependence of the precession frequency was modeled by The best fit with this model gives an ordering temperature T N = 1.38(2) K, a critical exponent β = 0.37 (2), and a phenomenological parameter α = 1.4(2). The resulting curve is shown in Fig. 6(d). Fixing α to unity and fitting the data between 0.85 and 1.40 K, i.e., closely below the transition temperature, results in very similar values for the critical exponent β = 0.344 (30) and the transition temperature T N = 1.382(10) K. Note that the critical temperature is consistent with the discontinuities found in Figs. 6(b) and 6(c). f. Nuclear magnetic resonance. Selected 1 H-NMR spectra, recorded for an out-of-plane field of 7 T and temperatures in the regime of long-range order, are presented in Fig. 7(a). Due to several nonequivalent hydrogen sites in the crystallographic unit cell, the resulting 1 H-NMR spectrum is composed of many resonance peaks, and, therefore, rather complicated. Since the same qualitative temperature dependence was observed for all 1 H lines of the spectrum, we consider in the following only selected lines with comparably little overlap. 064431-7  1 2 Heisenberg square-lattice antiferromagnets with relevant exchange and anisotropy parameters. The exchange interaction J, ordering temperature T N , g factor, anisotropy field H A , and saturation field H sat are experimentally determined. The inter-to intralayer coupling ratio J /J is estimated by the use of Eq. (5). The easy-plane anisotropy parameter is calculated from the out-of-plane static susceptibility minimum by the use of Eq. (7). A direct estimate of from ESR measurements in the ordered state [42] is denoted by ( ) . At temperatures above T N and an out-of-plane field of 2 T, a single, slightly nonsymmetric Gaussian-like line is observed.
The larger field of 7 T allows one to assign two nonequivalent hydrogen sites with Gaussian lineshape in the paramagnetic regime, as exemplified by the blue and red fits to the spectrum at 3.5 K in the inset of Fig. 7(c). As shown in Figs. 7(b) and 7(c), with decreasing temperature, the 1 H-NMR line splits into two sets of doublets, revealing a phase transition to long-range order at 2.3 and 2.8 K for μ 0 H c = 2 and 7 T, respectively. The splitting of the NMR spectrum is a clear signature of commensurate antiferromagnetic order, where each of the two lines represents a sublattice magnetization of opposite local spin polarization. The observation of multiple doublets is caused by several nonequivalent hydrogen sites in the lattice, with coincidental overlap of the nuclear resonance frequency in the paramagnetic temperature regime. Due to the different hyperfine coupling constants of the corresponding 1 H sites, this overlap is lifted in the ordered state. Considering the quasi-2D, almost-square structure of the Cu 2+ ions in CuPOF, the commensurate antiferromagnetic order is presumably of checkerboard type.
We note that the temperature dependence of the sublattice magnetization curves deviates from the mean-field-type behavior probed by the μ + SR precession frequency at zero field. The details of this field-induced behavior will be a subject of future, more detailed investigations by local-probe techniques.
As part of a thorough characterization of CuPOF, additional room-temperature 13 C magic-angle spinning (MAS) NMR and 31 P cryo-MAS NMR studies [66] were performed, and are presented in the SM (Figs. S11-S15) [49].

IV. DISCUSSION
The ideal 2D QHAF is described by the Hamiltonian (1) for the case J = = 0. Without applied field, the thermal fluctuations at arbitrarily low temperatures prevent a semiclassical order, regardless of the strength of the intralayer interaction J [4]. Any small perturbations, such as a finite interlayer coupling or spin-exchange anisotropy, give rise to quasi-long-range order below a nonzero transition temperature. Hence, the ratio k B T N /J can be treated as a measure of perturbations to the pure 2D QHAF [35]. It varies between zero for the 2D and 0.946(1) for the three-dimensional (3D) isotropic spin-1 2 Heisenberg model [67]. For typical molecular-based Cu-pyrazine materials, the reported values of k B T N /J are between about 0.2 and 0.3, see Table I. The lowest ratio k B T N /J = 0.176 yet found applies to the inorganic material Sr 2 CuO 2 Cl 2 , with a strong intralayer coupling of J/k B 1450 K.
The nature of the transitions induced by finite J or are fundamentally different. Whereas the interlayer coupling J induces a Néel state [35], a finite easy-plane or XY -like anisotropy yields a temperature-driven crossover of the spin correlations from isotropic Heisenberg to anisotropic XY -type behavior at a finite temperature T co , predicted to lead to a topological BKT transition at T BKT < T co [25,26,68].
In all bulk materials, finite interlayer couplings and magnetic anisotropies are present; see Table I for some selected cases. However, until now, there has been no theoretical framework that allows for an accurate experimental determination of J in the presence of a nonzero anisotropy [40]. Still, the challenges associated with a clean experimental realization of the 2D QHAF can be discussed as follows: Predicted critical temperatures for the two limiting cases of a finite J with = 0, as well as that of a finite with J = 0 are presented in Fig. 8 as functions of the perturbation parameter, J /J or , respectively. The black curve represents the normalized Néel temperature, k B T N /J, of a 3D array of isotropic square 2D spin- 1 2 Heisenberg planes with a coupling J between adjacent layers as a function of the exchange ratio  [71,72]. The values indicated by are for CuPOF, are from Ref. [40], is from Ref. [42], and are from Ref. [47]. The numbers denote 2D QHAF compounds according to the labeling in Table I. The two distinct values for CuPOF and Sr 2 CuO 2 Cl 2 denote slightly different values of , determined from the DC susceptibility and low-field magnetization measurements, respectively. J /J from Ref. [35]: where ρ s = 0.183J is the renormalized spin-stiffness constant. Note that the vertical axis is linear, whereas the horizontal axis is logarithmic and spans three orders of magnitude. The very weak decrease of the ordering temperature with reduction of the ratio J /J results from the exponential divergence of the correlation length of the 2D QHAF at low temperatures [17,69].
In the other limiting case with J = 0 and finite > 0, quantum Monte Carlo calculations showed that, even for anisotropies as small as 10 −3 , the critical behavior of the magnetic lattice resembles that of the Berezinskii-Kosterlitz-Thouless universality class. A slow logarithmic decrease of the T BKT temperature with reduction of the spin-exchange anisotropy was determined as [68] as depicted by the red dashed line in Fig. 8. From the comparison of these results, it is apparent that the interlayer interaction determines the ordering process for equal magnitudes of J /J and . Upon cooling from T > J/k B , the onset of the 3D long-range order occurs before any signatures of the exchange anisotropy can be observed. The influence of the spin anisotropy is only relevant if the interlayer interaction J /J is significantly smaller than . Therefore, minimizing the interlayer spin interaction is of key importance for developing any material approximation to the ideal 2D QHAF.
In contrast to a strongly anisotropic spin system with 1, where the topological excitations of unpaired vortices and antivortices are formed well above T BKT , a qualitatively different behavior is expected for a weakly anisotropic system with 1 [26]. At high temperatures, the spin correlations can be well approximated as isotropic. With decreasing temperatures, the XY anisotropy becomes relevant and stabilizes a planar spin configuration. The formation of vortices and antivortices starts in the regime of the temperature T co , which indicates the crossover between isotropic and XY behavior. Quantum Monte Carlo calculations revealed that the uniform susceptibility is very sensitive to this crossover phenomenon [26,68,70]. Whereas the in-plane susceptibility monotonically decreases with temperature below about J/k B , a characteristic minimum of the out-of-plane susceptibility marks the crossover from isotropic to anisotropic behavior at T co . The dependence of the crossover temperature on the spin anisotropy can be described by the empirical expression [26] and is depicted by the dash-dotted blue line in Fig. 8. Furthermore, in the presence of XY spin anisotropy, the field-dependent magnetization is expected to yield a qualitatively different behavior for in-and out-of-plane orientations of the field at T < T co . Below T co , the antiferromagnetically coupled magnetic moments preferably fluctuate in the easy plane. The application of a magnetic field suppresses longitudinal spin fluctuations, effectively inducing an easy-plane anisotropy of spin correlations perpendicular to the field direction. For a field applied perpendicular to the intrinsic easy plane, this yields a monotonic increase of the magnetization. On the other hand, for a magnetic field applied in plane, when the Zeeman term becomes larger than the intrinsic anisotropy energy, the total spin polarization in the field direction is enhanced, yielding a slope increase of the magnetization curve at the anisotropy field H A . Therefore, H A represents a measure of the spin-anisotropy parameter and can be estimated as = H A /H sat . It is worth noting that, for most bulk realizations of the 2D QHAF model, the weak intrinsic XY anisotropies do not stem from anisotropic exchange interactions. For the prevalent case of 2D QHAF materials based on Cu 2+ ions, the exchange between spins is almost isotropic. Typically, a weak anisotropy of a few percent of the isotropic superexchange interactions stems from crystal electric-field effects and related nonquenched orbital contributions. Commonly, the magnetic properties are then described by an effective-spin formalism, yielding an anisotropic g factor. Equivalently, this anisotropy can be treated in terms of an anisotropic exchange of isotropic effective spin moments, as described by the Hamiltonian (1) [73].
Since the impact of a finite interlayer interaction J and spin anisotropy on the critical behavior differ significantly [35,68], both parameters need to be experimentally determined for a newly synthesized compound to be accurately described as a quasi-2D QHAF. The strength of the antiferromagnetic exchange, J, as well as confirmation of lattice and exchange geometry being square rather than rectangular,  are of key importance. Secondary considerations include the possible existence of any next-nearest-neighbor interactions and spin-canting terms.
For the present case of CuPOF, from the comparison of theoretical modeling to our experimental results of DC susceptibility, specific heat, and high-field magnetization, a leading intralayer exchange constant J/k B = 6.80(5) K is consistently determined. All magnetic properties are in excellent agreement with theoretical predictions for the 2D squarelattice spin- 1 2 Heisenberg antiferromagnet; see Figs. 2, 4, and 5. The resulting values of J, independently obtained by means of the aforementioned techniques, are in excellent agreement within experimental error. Additionally, the results of DC susceptibility and magnetic specific heat were analyzed in terms of a possible rectangular rather than a square magnetic in-plane structure [57,58] and are both fully consistent with the square-lattice case.
To further characterize CuPOF, the anisotropy of the electronic g factor was investigated by means of ESR spectroscopy. A very weak in-plane anisotropy was found, with g a = 2.073(1) and g b = 2.066(3), in contrast with the outof-plane g factor g c = 2.298 (2). Overall, the electronic g factor anisotropy evinces a planar-like magnetic structure, with half filled d x 2 −y 2 orbitals oriented in the crystallographic ab plane. The integrated ESR intensity, recorded for the out-of-plane field orientation at 9.4 GHz, is in excellent agreement with the bulk magnetic susceptibility, see Fig. 2, and confirms the intralayer antiferromagnetic coupling as J/k B = 6.80(5) K. Moreover, the linear frequency-field dependence of the ESR spectrum at 1.5 K, observed in the frequency range between 52 and 500 GHz, indicates the absence of a notable energy gap upon approaching the zero-field limit.
The absence of any specific-heat anomaly, see Fig. 5, which would be associated with the transition to long-range order at T N = 1.38(2) K, as determined by μ + SR, indicates a high isolation of the magnetic layers and related small change of the residual entropy at T N [74]. This sets the upper limit of the interlayer coupling as J /J < 10 −2 [74]. Due to its pronounced low dimensionality, the thermodynamic properties of CuPOF show a very good agreement with those for an ideal isotropic 2D QHAF, where both the magnetic susceptibility and the specific heat exhibit a broad maximum at temperatures of the order of the nearest-neighbor exchange interaction J/k B . On the other hand, the local-probe techniques μ + SR and NMR are highly sensitive to the onset of static internal-field components that are associated with long-range order. μ + SR was successfully used to determine the long-range-ordering temperature T N for a wide range of quasi-2D square-lattice quantum Heisenberg antiferromagnets [44,46,47,62].
In the case of CuPOF, the temperature evolution of the asymmetry parameter, precession frequency, and the relaxation rate of the oscillatory component of the μ + SR asymmetry relaxation function (3) revealed a zero-field transition to long-range order at T N = 1.38(2) K, see Fig. 6. Thus, the ratio k B T N /J is found to be 0.203 for CuPOF, which is the smallest among all yet-characterized molecular-based materials with a magnetic lattice of Cu 2+ ions, compare Table I. By use of the empirical formula (5), the inter-to intralayer exchange ratio is obtained as J /J = 1.4×10 −4 , with J 1 mK.
Since strictly isotropic exchange interactions, = 0, are assumed in the derivation of Eq. (5), J /J = 1.4×10 −4 sets an upper limit to the effective interlayer coupling. The dipolar interaction between nearest-neighbor Cu 2+ ions in adjacent layers, with a Cu-Cu distance of 13.097 Å, is estimated as about 1 mK. Therefore, the interlayer interaction likely stems from dipole-dipole coupling, rather than superexchange via the interlayer molecules. A planar magnetic structure with highly isolated layers is further supported by our densityfunctional theory (DFT) calculations [75,76]; see the SM [49]. Both the interlayer interactions and the next-nearest-neighbor intralayer interactions are negligible in comparison with the antiferromagnetic in-plane coupling mediated by the pyrazine molecules.
As discussed above, the anisotropy parameter can be estimated from the low-temperature magnetization or lowfield susceptibility measurements [26,40,47,68,72,77]. From the anisotropy field μ 0 H A = 0.36(1) T, see inset of Fig. 4(b), the exchange anisotropy in CuPOF is evaluated as = H A /H sat = 1.85(5)×10 −2 . This results in the largest ratio H A /H sat yet found for Cu-pyrazine-based layered materials; compare Table I. A slightly smaller anisotropy parameter is estimated from the out-of-plane DC susceptibility at 0.1 T; see inset of Fig. 3. Employing the empirical Eq. (7), derived for the case of J = 0, the anisotropy parameter = 0.9(2)×10 −2 is determined. The qualitative behavior of the field-dependent magnetization (Fig. 4) and temperaturedependent DC susceptibility (Fig. 3) are in very good agreement with calculations for a 2D QHAF with an XY anisotropy of 1%-2% of the intralayer coupling [26,68], and resembles the characteristic behavior of previously studied Cu-pyrazinebased compounds [40,47,72].
To evaluate the relevant perturbations with respect to the ideal 2D QHAF and to shed light on the driving mechanism of the long-range order in CuPOF, the normalized Néel temperature k B T N /J is shown as a function of the evaluated spin anisotropy in Fig. 8 and compared with other molecularbased quasi-2D QHAFs. Where possible, the parameter was determined from experimental values of H A /H sat for the compounds labeled by black solid circles and stars. The two values for the inorganic compound Sr 2 CuO 2 Cl 2 are based on H A /H sat = 1.8×10 −4 [19] and the minimum of the out-ofplane DC susceptibility k B T co /J = 0.221 [72], which corresponds to = 8.3×10 −4 , see Table I. The experimental values are compared with the theoretical expectation of the BKT transition temperature for a weakly anisotropic 2D QHAF, described by the empirical function (6).
For highly isolated quasi-2D materials, such as Sr 2 CuO 2 Cl 2 with J /J = 3×10 −5 , the experimentally observed value of k B T N /J is very close to the prediction of k B T BKT /J (red dashed line in Fig. 8). Although the transitions to long-range and topological order at T N and T BKT , respectively, are of a different nature, it was argued that, due to the exceptional low dimensionality of Sr 2 CuO 2 Cl 2 , the long-range order is triggered by the inceptive BKT-type topological transition at T BKT T N [25,26]. Due to the finite interlayer coupling, a Néel-type antiferromagnetic order is stabilized at T N , before the topological transition is completed. A very similar scenario is proposed for CuPOF, motivated by the close agreement between the experimentally obtained value of T N and the theoretically predicted value of T BKT . Moreover, the larger ratio k B T N /J = 0.203 as compared with k B T N /J = 0.176 for Sr 2 CuO 2 Cl 2 , is attributed only to the stronger intrinsic spin anisotropy in CuPOF.
Furthermore, the field-driven increase of the crossover temperature T co , observed by DC susceptibility, see Fig. 3, reveals the corresponding increase of the effective spin anisotropy. This field-induced spin anisotropy is evinced as well by the strong increase of the ordering temperature T N in applied magnetic field, as found by 1 H-NMR spectroscopy, see Fig. 7. This strong increase of T N is attributed to the exceptional two-dimensionality of CuPOF, as compared with other less-isolated quasi-2D Cu-pyrazine-layered compounds, for which weaker field-induced changes of T N were reported [42,45,78,79]. The splitting of the 1 H-NMR spectrum below the transition to long-range order is a clear signature of commensurate antiferromagnetic order.
In conclusion, we present a comprehensive characterization of the newly synthesized molecular-based compound [Cu(pz) 2 (2-HOpy) 2 ](PF 6 ) 2 (CuPOF). Employing bulk magnetometry, specific heat, density-functional theory calculations, ESR, μ + SR, and NMR spectroscopy, CuPOF is characterized as an excellent realization of the 2D square-lattice spin-1 2 Heisenberg model with a moderate nearest-neighbor exchange interaction of J/k B = 6.80(5) K, and well-separated magnetic layers. The intralayer interaction is about four orders of magnitude larger than the estimated upper limit on the interlayer interaction, J 1 mK. A weak intrinsic easyplane anisotropy, revealed by bulk magnetometry, yields a temperature-driven crossover of the spin correlations from isotropic Heisenberg to anisotropic XY -type behavior and constitutes the driving mechanism of a transition to magnetic long-range order at T N = 1.38(2) K, as revealed by μ + SR spectroscopy. The application of a magnetic field normal to the easy plane yields a field-driven increase of the magnetic anisotropy, as shown by the evolution of the crossover temperature T co in the DC susceptibility data. The application of magnetic fields of several tesla leads to a strong increase of T N , as revealed by 1 H-NMR spectroscopy, in agreement with the pronounced two dimensionality of the magnetic lattice in CuPOF. As an outlook, our comprehensive characterization of [Cu(pz) 2 (2-HOpy) 2 ](PF 6 ) 2 as a clean realization of a 2D square-lattice spin- 1 2 Heisenberg antiferromagnet with moderate intralayer coupling and highly isolated magnetic layers calls for further studies of the field-induced effects on the anisotropy of the magnetic correlations [37], in particular by scattering and local-probe techniques.
Data presented in this paper resulting from the UK effort will be made available [80].
Extremely well-isolated 2D spin-1/2 antiferromagnetic Heisenberg layers with small exchange coupling in the molecular-based magnet CuPOF: Synthesis and Crystal Structure. Copper(II) bromide (0.58 g, 2.5 mmol) and silver hexafluorophosphate (1.261 g, 5.0 mmol) were separately dissolved in 5 ml of methanol each. The copper-bromide solution was added to the other solution while stirring; an immediate precipitate resulted. The combined solution was stirred for thirty minutes, then filtered. The precipitate yielded a dry weight of 0.820 g, equivalent to 4.36 mmol of AgBr (87% yield). Pyrazine (0.40 g, 5.0 mmol) and 2-hydroxypyridine (0.475 g, 5.0 mmol) were dissolved separately in 5 ml of MeOH. These solutions were combined and added dropwise to the stirred filtrate; the green cupric solution gradually changed to an olive-green color but no precipitate appeared. The solution was stirred for one hour, filtered, then partially covered and set to evaporate. Within several days, flat green crystals appeared. The product was collected by vacuum filtration, rinsed once with MeOH, then dried under vacuum. The crystals obtained were pale green, pleochroic plates, several mm on a side. Under one orientation of polarized light, the crystals are emerald green, and sky blue along a perpendicular orientation. The crystals extinguish between crossed polarizers when the growth edges are parallel to the axes of polarization. Total product (0.504 g, 0.716 mmol) for a yield of 28.6%. The infrared spectrum by the attenuated total reflection (IR (ATR)) and the results of the elemental analysis (CHN) appear in the following. IR (ATR): 3411w, 3134w, 1646s, 1595s,  1543m, 1455m, 1425s, 1401m, 1367w, 1258w, 1231w,  1158m, 1170w, 1155w, 1120m, 1094w, 1086w, 1066m (11.76, 11.90).
The local coordination sphere of the Cu 2+ ion in [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 (CuPOF) is presented in Fig. S1. The copper ion and the four pyrazine nitrogen atoms lie in the ab plane, whereas the coordinated 2pyridone molecule is nearly normal to the plane. Symmetry equivalent atoms A are generated via a mirror plane perpendicular to the a axis, whereas symmetry equivalent atoms B are generated by a two-fold rotation axis parallel to b.
Magnetometry. The single-crystal susceptibilities of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 between 1.8 and 12 K are shown in Fig. S2, where they are plotted as the ratios χ i /C i to demonstrate the equality of the temperature dependence for the out-of-plane and in-plane directions c and ab, respectively, where C i is the corresponding Curie constant. C c and C ab are 0.496 and 0.402 emuG −1 mol −1 K, respectively. This equivalence demonstrates that the ideal Heisenberg model is appropriate to describe the data in this temperature range. In order to obtain the absolute values of the magnetization recorded in the pulsed-field experiments, the data were directly compared with the magnetization recorded for a CuPOF single crystal at 4 K for H c (at 2 K for H ⊥ c) and up to 14 T, using a vibrating sample magnetometer (VSM) in a superconducting magnet and 4 He cryostat. These results are presented in Fig. S3.
Additional measurements of the magnetization of a polycrystalline sample in DC fields up to 35 T were done using a VSM magnetometer at the National High Magnetic Field Laboratory in Tallahassee, as well as in pulsed fields up to 25 T at the NHMFL facility in Los Alamos. Representative examples of these data are presented in Fig. S4. The results are fully consistent with the data from the HLD (see Fig. 3

in the main text).
Electron Spin Resonance. The room-temperature ESR powder spectrum of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 , measured at 9.8 GHz, is shown in Fig. S5. A least-squares fit to the spectrum performed with EasySpin is shown as the red line. The fit to the spectrum indicates a slightly rhombic g-tensor with g a = 2.071, g b = 2.056, and g c = 2.322. The lineshape is primarily Lorentzian with a smaller Gaussian contribution, as indicated by the fit parameters 1.207 mT (Lorentzian contribution) and 0.524 mT (Gaussian contribution). The parameters of the chosen pseudo-Voigt function represent coefficients in the linear combination of Lorentzian and Gaussian contributions to the lineshape [1]. The mostly Lorentzian lineshape indicates a dominant superexchange interaction [2]. A small value of G strain = 0.027 mT was permitted in fitting the lineshape contribution of g c , where G strain accounts for anisotropic broadening due to a local distribution of the g-factor, by adding a Gaussian envelope to the line [1]. The average g-factor determined from the powder ESR spectrum, g pow = 2.150, agrees well with the values determined from the single crystal measurements, g SC = 2.146(3), as well as from the magnetic DC susceptibility, g χ = 2.167. The anisotropy of the room-temperature ESR spectrum of a single-crystalline CuPOF sample was investigated at 9.8 GHz. A single ESR line of nearly Lorentzian lineshape was observed for each field orientation. As a representative example, the ESR spectrum obtained for a field along the crystallographic a axis is shown in Fig. S6.
178.02 (7)  formed at the Dresden High Magnetic Field Laboratory. The anisotropy of the g-factor obtained from these measurements is in very good agreement with those determined from the measurements at Clark University. A comparison of all experimentally determined g-factors and linewidths is presented in Table S3. The angular dependences ∆H pp (θ) of the room-    FIG. S7. Angular variation of the g-factor of [Cu(pyz)2(2-HOpy)2](PF6)2 at room temperature measured at 9.8 GHz for three orthogonal planes. The anisotropies of the g-factor were modeled using a cos 2 (θ)-type angular dependence, as represented by the red solid lines. Note the fine scale for the ab-plane variation.

Hamiltonian
where H Z is the Zeeman interaction, H J is the isotropic exchange interaction, and D and d are the anisotropic and antisymmetric exchange interactions, represented by the term H D/d . Due to the relatively strong isotropic exchange, J/k B = 6.80(5) K, the hyperfine interaction was neglected. The anisotropy of the g-factor likely accounts for some of the broadening observed in the c direction. Based on the powder spectrum of CuPOF, a minor Gaussian broadening of about 0.267 Oe is expected from the anisotropy of the g-factor. Due to the large Cu-Cu distance, the dipolar interaction in the c direction is expected to be small. A possible spin diffusive behavior, expected for 2D systems [3], is not observed. In the ab plane, no angular variation of the linewidth is observed within experimental uncertainty. The angular variation of the linewidth in the ac and bc planes was interpreted as the sum of the isotropic scaling factor A and the variation of the linewidth due to the antisymmetric and anisotropic The ESR frequency vs. magnetic field diagram at 1.5 K and H c is presented in Fig. S9. A single resonance line was observed at all frequencies up to 500 GHz. A linear fit of the frequency-field dependence of the ESR mode (red solid line in Fig. S9) yields a slope, and, thus, gfactor of g c = 2.29(1). At intermediate fields, when crossing to the quasi long-range ordered regime (compare Fig. 7), the resonance is expected to change its character from paramagnetic to antiferromagnetic. In general terms, the development of 3D correlations at the transition to a quasi long-range ordered state results in the opening of an antiferromagnetic excitation gap that can be probed in ESR experiments [4]. The gap in CuPOF appears to be small, as no signatures of it are detected within the accessible frequency range of our ESR spec-trometer, with a lower limit of about 50 GHz. On the other hand, ESR spectroscopy is well established as a sensitive probe for enhanced spin correlations in the vicinity of the transition to magnetically ordered states [5]. For the case of CuPOF, more detailed investigations of the temperature-and field-driven evolution of spin correlations as probed by ESR spectroscopy will be the subject of future work.  The specific heat between 1 to 10 K reveals a broad hump exceeding the phonon contribution, see main panel of Fig. S10. The 3 He data were analyzed as the sum of the magnetic specific heat of a 2D QHAF and a phononic contribution. The magnetic specific heat capacity was represented as a ratio of polynomials where R is the gas constant and the values of the coefficients N i and D i are given in Table S4. This new polynomial is similar in form to one used in a previous study [6], but is based on recent quantum Monte Carlo simulations of the magnetic specific heat [7] that extended to lower relative temperatures and can also represent the specific heat of rectangular lattices in which the exchange strengths (J and αJ) along the a and b axes are different. The range of validity for the square lattice is between 0.15 ≤ T /J ≤ 5.0. Muon Spin Relaxation. In a µ + SR experiment [8], 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 behavior 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.
Magic-Angle Spinning NMR spectroscopy. 13 C magicangle spinning (MAS) NMR spectra were recorded employing a Bruker AVANCE-II-600 spectrometer in a 14.1 T magnetic field, using a home-built MAS probe for 4×25 mm Si 3 N 4 rotors. 137 mg of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 powder was packed into the rotor and spun at 13.5 kHz. The temperature dependence of the spectrum was recorded between 276 and 332 K. Cooling of the sample was maintained using a gas cooling unit BCUX from Bruker. The actual temperature of the rotating sample was determined by 207 Pb chemical shift of lead nitrate rotating under identical conditions [9]. The temperature variation within the sample was typically ±1.5 K, and ±3 K for temperatures above 300 K.
Room-temperature 31 P spectra have been recorded on a Bruker AVANCE-II-360 spectrometer in an 8.45 T external magnetic field with a 31 P resonance frequency of 145.56 MHz, using a home-built MAS probe for 1.8×15 mm rotors. Variable-temperature 31 P MAS NMR spectra were studied at temperatures between 27 and 310 K. Low-temperature spectra have been recorded with the Bruker AVANCE-II spectrometer at an external field of 4.7 T with a 31 P resonance frequency of 80.985 MHz, using a home-built cryoMAS probe [10]  where the spectrum was recorded at 25 kHz spinning rate. This slowing down of the spinning rate is necessary to stabilize such a low temperature. The temperature of the spinning sample was measured with a temperature sensor at the spinner assembly and corrected for a given spinning frequency with the temperature dependent 207 Pb chemical shift of Pb(NO 3 ) 2 [9], which was rotating at the same conditions. The temperature dependence (276 -332 K) of the 13 C MAS NMR spectrum is shown in Fig. S11. All five resonances of the 2-pyridone molecules are resolved and appear at frequency shifts between 110 and 150 ppm, see left panel in Fig. S11. Spinning sidebands are also observed and marked by asterisks in the right panel of Fig. S11. The frequency shifts of the 2-pyridone carbons are close to the regular shift values and are independent of temperature, except for the carbon site C2. There are four independent carbons on the two independent pyrazine molecules. However, only three MAS NMR spectral lines have been resolved at the lowest temperature measured. These lines appear between 350 to 420 ppm with monotonous increase of the frequency shift with decreasing temperature, see Table S5. The hyperfine shift of the carbon sites can be evaluated as K = δ OBS − δ CS , where δ OBS and δ CS are the observed frequency shift and the chemical shift value in solution, respectively.
The hyperfine shift K is caused by the magnetic moments of the copper d-shell electrons. It is related to the molar magnetic susceptibility χ mol by H hf χ mol /N A µ B , where H hf is the hyperfine field, µ B is the Bohr magneton and N A is Avogadro's number. The hyperfine field can be determined from the plot K vs χ mol . Fig. S12 shows that the hyperfine shift of the 2-pyridone C2 carbon (b) and that of the pyrazine carbons (a) linearly depend on the magnetic susceptibility. The slopes in the plots give  The room-temperature 31 P-MAS NMR spectrum is shown in Fig. S13. Due to the cancellation of local dipole fields by spinning the sample with a frequency of 40 kHz at the magic-angle orientation, the 31 P-MAS spectrum is revealed as J resolved with 7 lines, separated by a scalar spin-spin coupling of J = 712 Hz between the 19 F and 31 P nuclei. The chemical shift of -143.2 ppm and J are in good agreement with previously reported values for compounds containing PF 6 molecules [11,12]. As the hyperfine fields at the phosphorus sites and chemical shift variation of the PF6 anions are small, the two phosphorus sites are not resolved in the MAS experiments. The absence of considerable spinning sidebands of the resonance indicates an isotropic surrounding of the phosphorus. The temperature-dependent 31 P-MAS NMR spectra are shown in Fig. S14. At high temperatures, the 31 P-MAS NMR spectrum appears J resolved with 7 components, as described above for the room-temperature spectrum. Towards lower temperatures, the components become broader and the resonance becomes gradually a single slightly asymmetric line. Usually, the main broadening mechanism in paramagnetic rotating solids is the influence of the magnetic susceptibility of the powder particles [13]. The frequency shift of the 31 P resonance line shows almost no variation with temperature, see Fig. S15. Usually, the paramagnetic shift is proportional to the magnetic susceptibility. With decreasing temperature from 300 to 27 K the latter increases significantly (see Fig. 2 in main text), whereas there is nearly no variation of the isotropic value of the 31 P NMR frequency shift in that temperature range.
Density Functional Theory (DFT) Calculations. The magnetic exchange couplings J between nearest-neighbor (NN) and next-nearest-neighbor (NNN) copper pairs were analyzed using the isotropic Heisenberg Hamiltonian The broken-symmetry (BS) approach [14,15] has been used to properly describe the open-shell low-spin states (LS, BS), while the equation proposed by Yamaguchi and co-workers [16] has been employed to account for the different weight of the open-shell and closed-shell states in the BS solution. All calculations have been performed with GAUSSIAN 09 [17] at the UB3LYP [18][19][20] level, and using a triple-zeta polarization basis set (TZVP). The magnetic exchange couplings (J 1 , J 2 , J 3 , and J 4 ) have been computed for four possible types of dimers in the unit cell of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 . Pairs separated by d 1 (6.676Å) and d 2 (6.680Å) correspond to the NN interaction through the pyrazine molecules, while d 3 (9.727Å) is the diagonal next-nearest-neighbor interaction within the magnetic layer. In turn, d 4 (13.10Å) is the nearest-neighbor distance between Cu units in adjacent layers. Pairs at larger distances have not been considered. The cluster models selected to compute J 1 to J 3 have been designed to incorporate the potentially- relevant counter ions in an explicit manner. Ideally, the treatment of the whole unit cell of the crystal, including periodic boundary conditions, would be necessary to properly account for effects of the environment in the J values. However, the small magnitude of the magnetic interactions make them incompatible with the approximations adopted in this type of calculations to keep the computational cost under control. In the present case, with the aim at balancing cost and accuracy, we have used a tetramer cluster model, i.e., four Cu centers, for the J 1 to J 3 intralayer interactions, rather than the common dimer model [21]. For comparison purposes, the dimer model has also been used to compute J 1 to J 3 and the interlayer value J 4 . The cluster models are shown in Fig. S17.
The computed magnetic interactions confirm the 2D magnetic topology of the system, see Table S6. The dominant magnetic interactions are between adjacent copper centers in the layers, with J 1 /k B of 18.7 K and J 2 /k B of 16.4 K. The diagonal interaction, J 3 /k B , is much weaker with only 0.3 K, whereas the interlayer interaction is negligible, with J 4 /k B < |0.15| K. Notice that the strengths of J 1 and J 2 are slightly different despite the near equality of the respective copper-copper distances. Such a difference was also observed in copper pyrazine perchlorate Cu(pz) 2 (ClO 4 ) 2 [22] and was ascribed to (i) the different relative disposition of the counterions and (ii) to a change in the Cu-N and N-N distances involved in the magnetic pathways.
The calculations using the dimer cluster model yield J 1 /k B = J 2 /k B = 29.3 K, and J 3 /k B = J 4 /k B < |0.15| K. Note that J 1 and J 2 are equivalent within this model, and J 3 vanishes in comparison with the calculations using a tetramer model. Finally, the dimer cluster model  in the absence of counter ions leads to the weakening of J 1 /k B = J 2 /k B = 11.8 K. This suggests, similar to the case of Cu(pz) 2 (ClO 4 ) 2 [22], that the counterions have a strong influence on the magnetic pathway, enhancing the strength of the antiferromagnetic interactions. It is conceivable that further enlargement of our cluster model, by increasing either the number of Cu centers or the neighboring PF 6 molecules, would change our computational estimation of J 1 . Generally, the DFT calculations confirm the magnetic topology of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 , as consisting of magnetically isolated 2D layers.
Using the computational estimates of J 1 -J 3 by the tetramer cluster model, the macroscopic properties of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 can be computed by using well-known statistical-thermodynamics expressions within the first-principles bottom-up procedure (FPBU) [23,24]. A subset of the magnetic topology, called model space, is selected in a way that, when the Heisenberg Hamiltonian is applied on the basis of a regionally reduced density-matrix approach, the resulting set of eigenvalues reproduces those that result from the application of the Heisenberg Hamiltonian to the full, infinite crystal. In the present study, we have selected a model space of sixteen Cu atoms arranged in a 4 × 4 plane, which aims at reproducing the infinite 2D layers of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 . Sixteen sites are usually sufficient to achieve convergence of the magnetic properties, especially considering the simplicity of the magnetic topology of [Cu(pyz) 2 (2-HOpy) 2 ](PF 6 ) 2 . Finally, the matrix representation of the Heisenberg Hamiltonian is built and fully diagonalized by using as a basis set the space of m s spin functions of a magnetic model of the selected 16 sites. The resulting energy and spin multiplicity of all possible magnetic states (up to 12870) are then used to calculate the macroscopic properties of the system, such as the magnetic susceptibility and specific heat. As shown in Fig. S18, the experimental and calculated values are qualitatively similar, with temperature dependences that scale with the values of the interaction strengths.