# Deformed spin-$\frac{1}{2}$ square lattice in antiferromagnetic
NaZnVOPO4(HPO4)

S. Guchhait  School of Physics, Indian Institute of Science Education and
Research Thiruvananthapuram-695551, India  D. V. Ambika  Qing-Ping Ding  Ames
Laboratory and Department of Physics and Astronomy, Iowa State University,
Ames, Iowa 50011, USA  M. Uhlarz  Dresden High Magnetic Field Laboratory (HLD-
EMFL), Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden, Germany  Y.
Furukawa  Ames Laboratory and Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA  A. A. Tsirlin<EMAIL_ADDRESS>Felix
Bloch Institute for Solid-State Physics, Leipzig University, 04103 Leipzig,
Germany  R. Nath<EMAIL_ADDRESS>School of Physics, Indian Institute of
Science Education and Research Thiruvananthapuram-695551, India

###### Abstract

We report the structural and magnetic properties of a new spin-$\frac{1}{2}$
antiferromagnet NaZnVOPO4(HPO4) studied via x-ray diffraction, magnetic
susceptibility, high-field magnetization, specific heat, and 31P nuclear
magnetic resonance (NMR) measurements, as well as density-functional band-
structure calculations. While thermodynamic properties of this compound are
well described by the $J_{1}-J_{2}$ square-lattice model, ab initio
calculations suggest a significant deformation of the spin lattice. From fits
to the magnetic susceptibility we determine the averaged nearest-neighbor and
second-neighbor exchange couplings of $\bar{J}_{1}\simeq-1.3$ K and
$\bar{J}_{2}\simeq 5.6$ K, respectively. Experimental saturation field of 15.3
T is consistent with these estimates if 20% spatial anisotropy in $J_{1}$ is
taken into account. Specific heat data signal the onset of a magnetic long-
range order at $T_{\rm N}\simeq 2.1$ K, which is further supported by a sharp
peak in the NMR spin-lattice relaxation rate. The NMR spectra mark the
superposition of two P lines due to two noneqivalent P sites where the broad
line with the strong hyperfine coupling and short $T_{1}$ is identified as the
P(1) site located within the magnetic planes, while the narrow line with the
weak hyperfine coupling and long $T_{1}$ is designated as the P(2) site
located between the planes.

##  I Introduction

Low-dimensional spin systems augmented with strong frustration reveal
suppression of a conventional magnetic order and may lead to a quantum
disordered ground state like quantum spin liquid (QSL) [1, *Balents199]. Two-
dimensional (2D) spin-$1/2$ frustrated square-lattice (FSL) is a well-known
example where frustration appears because of the competition between nearest-
neighbor (NN, $J_{1}$) and next-nearest-neighbor (NNN, $J_{2}$) exchange
interactions along the edges and diagonals of a square, respectively
($J_{1}-J_{2}$ model). Theoretical studies have determined a global phase
diagram with different ground states depending on the sign and relative
strength of the exchange couplings ($\alpha=\frac{J_{2}}{J_{1}}$) [3, 4].
Ferromagnetic (FM), Néel antiferromagnetic (NAF), and degenerate columnar
antiferromagnetic (CAF) ordered phases are stabilized and extended in the
regimes of $-0.5<\alpha\leq-\infty$, $-\infty\leq\alpha<0.5$, and
$0.5<\alpha<-0.5$, respectively. Different quantum disordered phases with
novel order parameters are predicted to emerge at the phase boundaries or the
quantum critical regimes. For instance, QSL [5, 6, 7], plaquette valence-bond
solid (PVBS) [8, 9] or columnar valence-bond solid (CVBS) [10] states are
expected around $\alpha\simeq 0.5$ while a spin nematic phase is predicted for
$\alpha\simeq-0.5$ in the phase diagram [11].

In the recent past, a handful number of Cu2+ (3$d^{9}$), V4+ (3$d^{1}$), and
Mo5+ (4$d^{1}$) based spin-$1/2$ FSL magnets have been studied as
$J_{1}-J_{2}$ model compounds. The vanadates are the most celebrated ones and
include Li2VO$X$O4 ($X$= Si, Ge) with $J_{1},J_{2}>0$ and $J_{2}\gg J_{1}$
[12, 13], along with V4+ phosphates that typically show $J_{1}<0$ and
$J_{2}>0$, also resulting in the CAF ground state [14, 15, 16, 17, 18, 19]. A
few other vanadates, such as Zn2VO(PO4)2 [20], VOMoO4 [21, 22, 23], and PbVO3
[24], feature $J_{1}\gg J_{2}$ and lie in the NAF region of the phase diagram.
A further variability becomes possible with Cu2+ as well as $4d$ Mo5+
compounds that span both $J_{2}\gg J_{1}$ and $J_{2}\ll J_{1}$ limits [25, 26,
27, 28, 29, 30, 31, 32, 33]. Unfortunately, none of these compounds fall in
the quantum critical regimes of the phase diagram around $\alpha=\frac{1}{2}$.
Moreover, some of them show intricate deformations of the magnetic square
lattice, because the underlying crystal symmetry is lower than tetragonal
[34]. This deformation was verified experimentally in several V4+ phosphates,
including Pb2VO(PO${}_{4})_{2}$ [35], SrZnVO(PO${}_{4})_{2}$ [36], and
BaCdVO(PO${}_{4})_{2}$ [37]. Each of them showed an interesting pre-saturation
phase [38, 39, 40, 41, 36], which is reminiscent of the nematic phase of the
$J_{1}-J_{2}$ model [11], but deviations from the ideal square lattice put
into question the applicability of this theoretical scenario, and alternative
interpretations were indeed proposed recently [42, 43].

Figure 1: (a) Crystal structure of NaZnVOPO4(HPO4) with the hydrogen position
and hydrogen bonds (dashed lines) determined in this work. (b) Magnetic
[VOPO4] layer with the square-lattice-like arrangement of the V4+ ions. (c)
Deformed square lattice with 5 nonequivalent interactions. VESTA software [44]
was used for crystal structure visualization.

In this context, finding V4+ square-lattice compounds with different
magnitudes of the deformation is important. Here, we report low-temperature
magnetic behavior of the hitherto unexplored NaZnVOPO4(HPO4) as a
spin-$\frac{1}{2}$ square-lattice candidate with a different spacer separating
the magnetic layers and, potentially, a different regime of exchange couplings
compared to the widely studied $AA^{\prime}$VO(PO${}_{4})_{2}$ phosphates
($AA^{\prime}$ = Pb2, SrZn, BaCd). The monoclinic crystal structure of
NaZnVOPO4(HPO4) features V4+ ions in the square-pyramidal coordination. They
are joined into layers via P(1)O4 tetrahedra (Fig. 1), while the interlayer
space is filled with the Na+ and Zn2+ ions, as well as the P(2)O4 tetrahedra
that represent the HPO4 groups.

##  II Methods

Pale-blue coloured powder of the titled compound was synthesized by the
conventional hydrothermal method. 0.150 g Na2CO3 (Aldrich, 99.995%), 0.154 g
V2O5 (Aldrich, 99.99%), and 0.115 g Zn (Aldrich, 99%) powders were mixed with
5 ml of a 1.5 M aqueous solution of H3PO4, sealed in a 23 ml Teflon lined
bomb, and heated at 240oC for 8 days followed by slow cooling
(10$\degree$C/hour) to room temperature. The obtained blue color product was
washed carefully with distilled water and dried in an oven maintained at
$100\degree$C for 24 hours. To check the phase purity of the compound, powder
x-ray diffraction (XRD) experiment was performed using a PANalytical powder
diffractometer with Cu $K_{\alpha}$ radiation ($\lambda_{\rm avg}\simeq
1.5418$ Å) at room temperature. We have also performed temperature-dependent
powder XRD measurements on the pure phase powder sample in the temperature
range 15 K$\leq T\leq 300$ K, using a low-temperature attachment (Oxford
PheniX) to the x-ray diffractometer. Le-Bail fit of the powder XRD patterns
was performed using the FULLPROF software package [45], taking the initial
structural parameters from the previous report [46].

Temperature variation of magnetization ($M$) was measured in the temperature
range 1.8–350 K in different magnetic fields ($H$) using a SQUID magnetometer
(MPMS3, Quantum Design). The isothermal magnetization was measured at $T=1.8$
K from 0 to 7 T. The high-field magnetization measurement was performed at
$T=1.4$ K in pulsed magnetic fields up to 30 T using the facility at the
Dresden High Magnetic Field Laboratory [47, 48]. The temperature dependent
specific heat of this sample was measured on a sintered pellet in a large
temperature range (0.5 K$\leq T\leq 200$ K) using the Physical Property
Measurement System (PPMS, Quantum Design) and adopting the thermal relaxation
technique. For measurements below 2 K, 3He attachment to the PPMS was used.

The Nuclear Magnetic Resonance (NMR) measurements were carried out on the 31P
nuclei (gyromagnetic ratio $\frac{\gamma_{\rm N}}{2\pi}=17.237$ MHz/T and
nuclear spin $I=\frac{1}{2}$) in the temperature range 1.8 K$\leq T\leq 250$
K. The NMR spectra at different temperatures were obtained by varying the
magnetic field at a constant frequency of 121 MHz. The spin-lattice relaxation
rate ($1/T_{1}$) was measured by the single saturation pulse method at two
frequencies (30.2 MHz and 121 MHz). The temperature-dependent NMR shift,
$K(T)=[H_{\rm ref}/H(T)-1]$, was calculated from the resonance field of the
sample $H$ with respect to the resonance field of a non-magnetic reference
sample $(H_{\rm ref})$.

Magnetic couplings were determined by density-functional-theory (DFT) band-
structure calculations performed in the FPLO code [49] using experimental
structural parameters from Ref. [46], except the hydrogen position that was
optimized as further explained in Sec. III.5. The Perdew-Burke-Ernzerhof (PBE)
approximation for the exchange-correlation potential was employed [50]. We
used superexchange theory in the vein of Kugel-Khomskii model [51, 52], as
well as the mapping approach [53, 54] based on total energies of collinear
spin configurations obtained from DFT+$U$ calculations with the on-site
Coulomb repulsion $U_{d}=4$ eV, Hund’s coupling $J_{d}=1$ eV, and double-
counting correction in the atomic limit [55].

Experimental thermodynamic properties were modeled by high-temperature series
expansion (HTSE) for the $J_{1}-J_{2}$ spin-$\frac{1}{2}$ square lattice [56]
as well as full diagonalization (FD) for the $4\times 4$ finite lattice with
periodic boundary conditions. Additionally, quantum Monte Carlo (QMC)
simulations were performed for the nonfrustrated spin-$\frac{1}{2}$ square
lattice using directed loop algorithm [57] in the stochastic series expansion
(SSE) [58] representation. QMC simulations were performed for $16\times 16$
finite lattices with periodic boundary conditions, using $4\times 10^{4}$
sweeps and $4\times 10^{3}$ thermalization sweeps. The ALPS package [59] was
used for both FD and QMC.

##  III Results

###  III.1 X-ray Diffraction

Figure 2: Powder XRD patterns (open circles) at (a) $T=300$ K and (b) $T=15$
K. The solid line is the Le-Bail fit, the vertical bars mark the expected
Bragg peak positions, and the lower solid line corresponds to the difference
between the observed and calculated intensities. The goodness-of-fit is
achieved to be $\chi^{2}\sim 5.4$ and $\sim 6.2$ for $T=300$ K and $15$ K,
respectively. Figure 3: (a) Lattice constants ($a$, $b$, and $c$) as a
function of temperature. (b) Monoclinic angle ($\beta$) along with the unit
cell volume ($V_{\rm cell}$) are plotted as a function of temperature from 15
K to 300 K. The solid line represents the fit of $V_{\rm cell}$ using Eq. (1).

The powder XRD patterns of NaZnVOPO4(HPO4) are analyzed by Le Bail fits.
Figure 2 presents the data at two end temperatures ($T=300$ K and 15 K). The
entire XRD pattern down to 15 K could be indexed using the monoclinic crystal
structure with the space group $P2_{1}/c$. Neither any indication of
structural transition nor lattice distortion is found down to 15 K. The
refined lattice parameters and unit cell volume ($V_{\rm cell}$) are
[$a=8.5418(4)$ Å, $b=8.9937(5)$ Å, $c=9.0765(5)$ Å, $\beta=91.186(3)^{\circ}$,
and $V_{\rm cell}\simeq 696.46$ Å3] and [$a=8.5190(2)$ Å, $b=8.9281(4)$ Å,
$c=9.1029(4)$ Å, $\beta=92.013(5)^{\circ}$, and $V_{\rm cell}\simeq 691.73$
Å3] for $T=300$ K and 15 K, respectively. The obtained lattice parameters at
room temperature are in close agreement with the values reported earlier [46].
The temperature dependence of the lattice parameters ($a$, $b$, $c$, and
$\beta$) and $V_{\rm cell}$ are presented in Fig. 3. The lattice constants $a$
and $b$ are found to decrease in a systematic way with decreasing temperature
while $c$ and monoclinic angle $\beta$ increase with decreasing temperature
and then rearch a plateau. These lead to an overall thermal contraction of the
$V_{\rm cell}$ with temperature.

The variation of $V_{\rm cell}$ with temperature can be expressed in terms of
the internal energy [$U(T)$] of the system [60, 61, *Guchhait224415],

$V_{\rm cell}(T)=\frac{\gamma U(T)}{K_{0}}+V_{0}.$ (1)

Here, $V_{0}$ is the unit-cell volume at $T=0$ K, $\gamma$ is the Grüneisen
parameter, and $K_{0}$ is the bulk modulus of the system. According to the
Debye model, $U(T)$ can be written as,

$U(T)=9Nk_{\rm B}T\left(\frac{T}{\theta_{\rm
D}}\right)^{3}\int_{0}^{\frac{\theta_{\rm D}}{T}}\frac{x^{3}}{(e^{x}-1)}dx,$
(2)

where $N$ is the total number of atoms per unit cell, $k_{\rm B}$ is the
Boltzmann constant, and $\theta_{\rm D}$ is the Debye temperature [63]. The
variable $x$ inside the integration stands for the quantity
$\frac{\hbar\omega}{k_{\rm B}T}$ with phonon frequency $\omega$ and Planck
constant $\hbar$. Here, $\theta_{\rm D}=\frac{\hbar\omega_{\rm D}}{k_{\rm B}}$
and $\omega_{\rm D}$ is the upper limit of $\omega$. The best fit of the
$V_{\rm cell}(T)$ data using Eq. (1) [solid line in Fig. 3(b)] yields the
parameters: $\theta_{\rm D}\simeq 312$ K, $V_{0}\simeq$ 691.72 Å3, and
$\frac{\gamma}{K_{0}}\simeq 8.84\times 10^{-12}$ Pa-1.

###  III.2 Magnetization

Figure 4: (a) $\chi(T)$ measured in the magnetic field of $H=0.5$ T. The lines
(solid and dashed) are the HTSE fits using the isotropic spin-$\frac{1}{2}$
FSL model, with two different solutions (solution-$a$ and solution-$b$).
Inset: The simulated $\chi(T)$ using the FD method taking the $J_{1}$ and
$J_{2}$ values from solutions-$a$ and $b$, respectively. (b) $1/\chi$ vs $T$.
The red solid line is the CW fit. Inset: $\chi(T)$ in the low-$T$ region
measured in different fields.

The magnetic susceptibility [$\chi(T)\equiv M/H$] of NaZnVOPO4(HPO4) measured
in an applied field of $H=0.5$ T is shown in Fig. 4(a). In the high-
temperature region, $\chi(T)$ follows a typical Curie-Weiss (CW) behaviour and
shows a rounded maximum at $T_{\chi}^{\rm max}\simeq 4.6$ K. Such a maximum
represents the short-range AFM order in the low-dimensional spin systems.
Below about 2.6 K, a small upturn is likely due to paramagnetic impurities
and/or defects present in the powder sample. No trace of magnetic long-range-
order (LRO) is detected down to 2 K. We have also measured $\chi(T)$ in
different applied fields [inset of Fig. 4(b)] but no obvious features
associated with magnetic LRO is found, except the suppression of a broad
maximum towards low temperatures.

The inverse susceptibility, $1/\chi(T)$, is shown in Fig. 4(b). In the
paramagnetic regime, $1/\chi(T)$ was fitted by the CW law

$\chi(T)=\chi_{0}+\frac{C}{T-\theta_{\rm CW}},$ (3)

where the first term ($\chi_{0}$) represents the combination of temperature-
independent diamagnetic and Van-Vleck paramagnetic susceptibilities. The
second term is the CW law where $C$ is the Curie constant and $\theta_{\rm
CW}$ is the CW temperature. The fit above 25 K returns the parameters:
$\chi_{0}\simeq 4.26\times 10^{-5}$ cm3/mol-V4+, $C\simeq 0.393$ cm3K/mol-V4+,
and $\theta_{\rm CW}\simeq-4.3$ K. Using the value of $C$, the effective
moment is calculated to be $\mu_{\rm eff}$ $[=(3k_{\rm B}C/N_{\rm A}\mu_{\rm
B}^{2})^{\frac{1}{2}}$, where $N_{\rm A}$ is the Avogadro’s number and
$\mu_{\rm B}$ is the Bohr magneton] $\simeq 1.77~{}\mu_{\rm B}$/V4+. This
value of $\mu_{\rm eff}$ is close to the actual value $1.73~{}\mu_{\rm B}$ for
a spin-$\frac{1}{2}$ transition-metal ion with $g=2$. The negative value of
$\theta_{\rm CW}$ indicates the dominant AFM exchange coupling between the V4+
ions. The core diamagnetic susceptibility ($\chi_{\rm core}$) of the compound
caused by the core orbital electrons was calculated to be $-1.37\times
10^{-4}$ cm3/mol by adding the $\chi_{\rm core}$ of Na+, Zn2+, V4+, P5+, and
O2- ions [64]. The Van-Vleck paramagnetic susceptibility ($\chi_{\rm vv}$) was
calculated to be $\sim 1.8\times 10^{-4}$ cm3/mol by subtracting $\chi_{\rm
core}$ from $\chi_{0}$, which is very close to the value reported for other
V4+ based compounds [65].

As evident from the structural data (Fig. 1), the system deviates from the
isotropic square lattice. This spatial anisotropy is mainly due to five non-
equivalent exchange couplings, three between nearest neighbors ($J_{1}$,
$J_{1}^{\prime}$, and $J_{1}^{\prime\prime}$) and two between next-nearest
neighbors ($J_{2}$ and $J_{2}^{\prime}$). Previous studies suggest that this
anisotropy has only a minor effect on thermodynamic properties at higher
temperatures [34]. Therefore, fits with the isotropic $J_{1}-J_{2}$ model
return averaged values of the NN and NNN couplings,
$\bar{J}_{1}=(2J_{1}+J_{1}^{\prime}+J_{1}^{\prime\prime})/4$ and
$\bar{J}_{2}=(J_{2}+J_{2}^{\prime})/2$, respectively. We fitted the data with

$\chi(T)=\chi_{0}+\chi_{\rm spin}(T)$ (4)

using the temperature-independent term ($\chi_{0}$) and the 9th order HTSE
[$\chi_{\rm spin}(T)$],

$\displaystyle\chi_{\rm spin}(T)=\frac{N_{\rm A}g^{2}\mu_{\rm B}^{2}}{k_{\rm
B}T}\sum_{n}\left(\frac{\bar{J}_{1}}{k_{\rm
B}T}\right)^{n}\sum_{m}c_{mn}\left(\frac{\bar{J}_{\rm 2}}{\bar{J}_{\rm
1}}\right)^{m}.$ (5)

The values of the coefficients $c_{mn}$ are taken from Ref. [56]. Our fit for
$T\geq 7$ K yields two solutions: solution-$a$ with $\chi_{0}\simeq 6\times
10^{-5}$ cm3/mol-V4+, $g\simeq 2.04$, $\bar{J}_{1}\simeq 5.18$ K,
$\bar{J}_{2}\simeq-0.66$ K and solution-$b$ with $\chi_{0}\simeq 6.6\times
10^{-5}$ cm3/mol-V4+, $g\simeq 2.04$, $\bar{J}_{1}\simeq-1.28$ K,
$\bar{J}_{2}\simeq 5.59$ K.

The main difficulty of the HTSE fit is choosing the appropriate $T$-range. The
convergence of the HTSE depends on the $\bar{J}_{2}/\bar{J}_{1}$ ratio, hence
the lower limit of the fitting range ($T_{\rm min}$) should be chosen with
caution depending on the results of the fit [16]. For a precise estimation of
the $J$ values, we varied the lower limit of the fitting $T$-range ($T_{\rm
min}$) between 5 and 9 K and estimated the exchange couplings for both the
solutions. Figures 5(a) and (b) present the variation of $\bar{J}_{1}$ and
$\bar{J}_{2}$ with $T_{\rm min}$ for the solutions-$a$ and $b$, respectively.
To check the convergence of the HTSE, we calculated the ratio of the ninth-
order term to the total susceptibility [$\chi_{9}(T)/\chi_{\rm total}(T)$]
using the parameters obtained from HTSE fits with a different $T_{\rm min}$.
Both solutions are stable above 7 K. We also performed the FD simulations of
$\chi(T)$ using the values of $\chi_{0}$, $g$, $\bar{J}_{1}$, and
$\bar{J}_{2}$ and found a good agreement with the experimental data [inset of
Fig. 4(a)].

Figure 5: Results of the $\chi(T)$ fit using HTSE for isotropic
spin-$\frac{1}{2}$ FSL model by varying the minimum temperature of the fitting
range ($T_{\rm min}$). (a) and (b) Averaged exchange interactions
($\bar{J}_{1}$ and $\bar{J}_{2}$) vs $T_{\rm min}$ for the solutions $a$ and
$b$, respectively. (c) Convergence test of the HTSE fit for both solutions,
$\chi_{9}(T)/\chi_{\rm total}(T)$ vs $T_{\rm min}$. (d) Saturation fields for
the solutions $a$ and $b$ vs $T_{\rm min}$. Figure 6: Magnetization ($M$) vs
field ($H$) at $T=1.4$ K measured using pulse magnetic field and scaled with
respect to the SQUID data. The solid line represents the QMC simulation,
assuming the spatially isotropic nonfrustrated spin-$1/2$ square-lattice model
with $J_{\rm eff}=5.1$ K. The arrow indicates the position of the saturation
field ($H_{\rm S}$).

Discriminating between the solutions-$a$ and $b$ may be possible using
measurements of the saturation field [48]. The $M$ vs $H$ data measured up to
25 T in pulsed magnetic fields at $T=1.4$ K are shown in Fig. 6. We have
scaled the high-field data with respect to the magnetic isotherm at $T=1.8$ K
measured up to 7 T using a SQUID magnetometer. $M$ increases linearly with $H$
in the low-field region, shows a positive curvature in the intermediate
fields, and then saturates at around $H_{\rm S}\simeq 15.3$ T. The positive
curvature is typical for low-dimensional and frustrated spin systems [66].

Saturation field of an FSL magnet depends on the type of magnetic order. In
the NAF case (solution-$a$), $H_{\rm S}=4\bar{J}_{1}k_{\rm B}/(g\mu_{\rm
B})\simeq 15.1$ T in good agreement with the experiment. On the other hand, in
the CAF case $H_{\rm S}=(\mathcal{J}_{1}+2\bar{J}_{2})2k_{\rm B}/(g\mu_{\rm
B})$ where $\mathcal{J}_{1}$ is the weaker of the couplings $J_{1}$ and
$(J_{1}^{\prime}+J_{1}^{\prime\prime})/2$. In the CAF state, spins align
antiferromagnetically along the direction of this coupling and
ferromagnetically along the orthogonal direction. Therefore, only the weaker
coupling enters the saturation field. Using $\mathcal{J}_{1}=\bar{J}_{1}$
leads to $H_{\rm S}^{\rm CAF}\simeq 14.45$ T for the solution-$b$, lower than
in the experiment. On the other hand, the actual value of $\mathcal{J}_{1}$ is
reduced owing to the deformation of the square lattice. Assuming 20% spatial
anisotropy in the NN couplings ($\mathcal{J}_{1}=0.8\bar{J}_{1}$) according to
the DFT results (Sec. III.5), we arrive at $H_{\rm S}^{\rm CAF}\simeq 15.15$
T, which is on par with the result for the solution-$a$ and also matches the
experimental value.

Interestingly, the magnetization curve can be also well described by the
simple NN square-lattice model with an effective coupling $J_{\rm eff}=5.1$ K
(Fig. 6), which is in good agreement with the leading exchange couplings
extracted from the HTSE fits.

###  III.3 Specific Heat

Figure 7: (a) Variation of $C_{\rm p}$ with temperature in the absence of
magnetic field. The red solid line represents the simulated phonon
contribution ($C_{\rm ph}$) taking into account the Debye (dash-dotted line)
and Einstein (dashed line) terms. Inset: $C_{\rm p}/T$ vs $T$ at low
temperatures and in different applied fields. The blue solid line is the
magnetic contribution ($C_{\rm mag}$). (b) $C_{\rm mag}/R$ vs $T$ (left
$y$-axis). The solid line represents the change in the magnetic entropy,
$S_{\rm mag}$ vs $T$ (right $y$-axis). The red (dash-dotted) and blue (solid)
lines are the HTSE fits corresponding to the two solutions, respectively. The
green (dashed) line is the QMC simulation for an isotropic non-frustrated
spin-$1/2$ square lattice model with $J_{\rm eff}=6.3$ K.

In Fig. 7(a), we have plotted the temperature-dependent specific heat [$C_{\rm
p}(T)$] measured from 0.5 to 140 K in zero applied field. It decreases
systematically with temperature and passes through a broad maximum at
$T_{C}^{\rm max}\simeq 3.7$ K, typical for low-dimensional oxides. An
anomalous behavior was observed in a narrow temperature range $T_{\rm
N}=2-2.2$ K where we couldn’t stabilize the temperature. This is a possible
indication of the onset of a magnetic LRO.

In a magnetic insulator, the total specific heat $C_{\rm p}(T)$ is the sum of
two main contributions: one is the phonon/lattice contribution [$C_{\rm
ph}(T)$], which dominates in the high-temperature region, and another one is
the magnetic contribution [$C_{\rm mag}(T)$], which dominates in the low-
temperature region depending upon the strength of the exchange coupling. In
order to bring out the magnetic part of the specific heat, we first quantified
the lattice contribution and then subtracted it from the total specific heat.
We simulated the high-temperature $C_{\rm ph}(T)$ data taking into account the
sum of one Debye [$C_{\rm D}(T)$] and three Einstein [$C_{\rm E}(T)$] terms,
i.e., $C_{\rm ph}(T)=C_{\rm D}(T)+\sum\limits_{i=1}^{3}C_{\rm Ei}(T)$. The
Debye and Einstein terms are expressed as

$C_{\rm D}(T)=9n_{\rm D}R\left(\frac{T}{\theta_{\rm
D}}\right)^{3}\int_{0}^{\frac{\theta_{\rm
D}}{T}}\frac{x^{4}e^{x}}{(e^{x}-1)^{2}}dx$ (6)

and

$C_{\rm Ei}(T)=3n_{\rm
Ei}R\left(\frac{\theta_{Ei}}{T}\right)^{2}\frac{e^{\left(\frac{\theta_{Ei}}{T}\right)}}{\left[e^{\left(\frac{\theta_{Ei}}{T}\right)}-1\right]^{2}}.$
(7)

Here, the Einstein temperatures $\theta_{\rm
Ei}=\frac{\hbar\omega_{Ei}}{k_{\rm B}}$, $\omega_{Ei}$ are the respective
Einstein frequencies, and $R$ denotes the universal gas constant. The values
of $n_{\rm D}$ and $n_{\rm Ei}$ are chosen in such a way that the sum $n_{\rm
D}+\sum\limits_{i=1}^{3}n_{\rm Ei}$ matches with the total number of atoms per
formula unit. The best fit of the $C_{\rm p}$ data in the high-temperature
regime using one Debye and three Einstein branches yields the characteristic
temperatures: $\theta_{\rm D}\simeq 122$ K, $\theta_{\rm E1}\simeq 307$ K,
$\theta_{\rm E2}\simeq 850$ K, and $\theta_{\rm E3}\simeq 161$ K with $n_{\rm
D}=1$, $n_{\rm E1}=5$, $n_{\rm E2}=7$, and $n_{\rm E3}=2$, respectively [67].
The red solid line in Fig. 7(a) is the total phononic contribution to the
specific heat ($C_{\rm ph}$) extrapolated down to low temperature, while the
dashed and dash-dotted lines are the Einstein and Debye contributions,
respectively. The value of $\theta_{\rm D}$ and the average value of
$\theta_{\rm E}$ match well with $\theta_{\rm D}$ estimated from the $V_{\rm
cell}$ vs $T$ analysis.

$C_{\rm mag}$ estimated after subtracting $C_{\rm ph}$ from $C_{\rm p}$ is
plotted in Fig. 7(b) as a function of temperature. The pronounced broad
maximum at $T_{C}^{\rm max}\simeq 3.7$ K mimics short-range antiferromagnetic
correlations. The change in magnetic entropy [$S_{\rm mag}(T)$] is obtained by
integrating $C_{\rm mag}/T$ over temperature. It saturates to a value of $\sim
5.6$ J.mol-1.K-1 at around $T\simeq 25$ K, which is close to the expected
value $R\ln 2=5.76$ J.mol-1.K-1 for a two-level (spin-$\frac{1}{2}$) system,
thus justifying our subtraction procedure and the evaluation of $C_{\rm mag}$.

The magnetic LRO is highlighted by plotting $C_{\rm p}/T$ vs $T$ in the inset
of Fig. 7(a) in different applied fields. While the anomaly associated with
the magnetic LRO is not pronounced in the zero-field data, the peak at $T_{\rm
N}\simeq 2.1$ K becomes more pronounced with increasing the field. This is due
to the transfer of entropy from the broad maximum to the transition anomaly.
Surprisingly, no visible shift in $T_{\rm N}$ is perceived even with the field
change of 9 T, indication of a robust AFM transition.

One can estimate the exchange couplings by analyzing $C_{\rm mag}(T)$ using
the HTSE of a spin-$\frac{1}{2}$ FSL model [56],

$\displaystyle\frac{C_{\rm mag}(T)}{R}=\frac{\bar{J}_{1}}{k_{\rm
B}T}\sum_{n}(-n)\left(\frac{\bar{J}_{1}}{k_{\rm
B}T}\right)^{n}\sum_{m}e_{mn}\left(\frac{\bar{J}_{\rm 2}}{\bar{J}_{\rm
1}}\right)^{m}$ (8)

We have fitted the $C_{\rm mag}(T)/R$ data at $T>6$ K [see Fig. 7(b)] and
arrived at two solutions that strongly resemble the two solutions from the
susceptibility fit: solution-$a$ ($\bar{J}_{1}/k_{\rm B}\simeq 5.9$ K,
$\bar{J}_{2}/k_{\rm B}\simeq-0.2$ K) and solution-$b$ ($\bar{J}_{1}/k_{\rm
B}\simeq-1.6$ K, $\bar{J}_{2}/k_{\rm B}\simeq 6.1$ K). Alternatively, we can
compare our experimental data with the QMC simulation for the nonfrustrated
square-lattice model [see Fig. 7(b)]. The position of the maximum in $C_{\rm
mag}$ is reproduced with $J_{\rm eff}=6.3$ K, which is notably higher than in
the $M(H)$ fit (Fig. 6). Moreover, the maximum value of $C_{\rm mag}(T)/R$ is
higher than in the experiment. The simple NN square-lattice model is thus
insufficient to describe the magnetic behavior of NaZnVOPO4(HPO4). Both
frustration and deformation of the square lattice affect thermodynamic
properties of this compound.

###  III.4 31 P NMR

NMR is a convenient local probe of both static and dynamic properties. In
NaZnVOPO4(HPO4), the P(1) site is strongly coupled to the V4+ ions within the
VOP(1)O4 layer, while the P(2) site which is located between the adjacent
layers is weakly coupled to the V4+ ions [see Fig. 1(a)]. This difference
allows a useful comparison and highlights the magnetic behavior of the square-
lattice planes.

####  III.4.1 31P NMR Spectra

Figure 8: Temperature evolution of the 31P NMR spectra of NaZnVOPO4(HPO4)
measured at 121 MHz. The dashed line indicates the reference field position.
Inset: 31P NMR spectrum at $T=12$ K with the dashed and dotted lines are the
fits of the P(1) and P(2) sites, respectively, and the solid line (final fit)
is the superposition of the P(1) and P(2) fits. The NMR shift values along the
$x$, $y$, and $z$ directions, obtained from the fitting are $K_{x}^{(1)}\simeq
0.48\%,K_{y}^{(1)}\simeq-0.68\%,$ and $K_{z}^{(1)}\simeq 0.24\%$ for P(1) site
and $K_{x}^{(2)}\simeq 0.13\%,K_{y}^{(2)}\simeq-0.1\%$, and
$K_{z}^{(2)}\simeq-0.09\%$ for P(2) site, respectively.

The field-sweep 31P NMR spectra above $T_{\rm N}$ measured in a radio
frequency of 121 MHz are shown in Fig. 8. Each spectrum is normalized by its
maximum amplitude and offset vertically by adding a constant. At high
temperatures, the line is found to be narrow but asymmetric and the central
peak appears at the zero-shift position. As the temperature is lowered, the
line width increases drastically and becomes more anisotropic with two
shoulders on either side of the central peak. This abnormal spectral shape can
be attributed to two nonequivalent P sites in the crystal structure [46].
Remarkably, the complete shape of the spectra could be reproduced considering
the superposition of two spectral lines. The inset of Fig. 8 portrays the
spectral fit at $T=12$ K. We have a narrow central line with weak anisotropy
and a broad asymmetric background with two distinct shoulders.

The narrow central peak shifts weakly, whereas the shoulders move
significantly with decreasing temperature. Thus, the narrow central line with
a weak shift can be assigned to the P(2) site, which is weakly coupled while
the broad line with the strong temperature-dependent behavior corresponds to
the in-plane P(1) site, which is strongly coupled to the V4+ spins. The
asymmetric shape of both P(1) and P(2) sites is likely due to the anisotropy
in $\chi(T)$ or asymmetry in the hyperfine coupling constant between the P
nuclear spins and the V4+ electronic spins. The overall spectral shape matches
exactly with the spectral shape reported for other V4+-based FSL compounds,
$AA^{\prime}$VO(PO4)2 ($AA^{\prime}$ = Pb2, SrZn), on the polycrystalline
sample [14, 17].

Figure 9: Temperature-dependent 31P NMR spectra of NaZnVOPO4(HPO4) measured at
121 MHz around $T_{\rm N}$. The dashed line indicates the reference field
position.

The nature of the magnetic ordering can also be gleaned from the analysis of
NMR spectra below $T_{\rm N}$. Figure 9 shows the normalized 31P NMR spectra
measured at 121 MHz around $T_{\rm N}$. Neither significant change in the line
shape nor any visible line broadening are observed below $T_{\rm N}$. This
indicates that the 31P site experiences only a weak static field in the
ordered state. This observation is quite opposite to that reported for
Pb2VO(PO4)2 [14] but similar to the 29Si NMR results on Li2VOSiO4 where local
field on the Si site cancels out [68]. All these compounds show a similar type
of magnetic layers and a similar mutual arrangement of the V4+ ions and
PO4/SiO4 tetrahedra. However, only in Pb2VO(PO${}_{4})_{2}$ the magnetic
[VOPO4] layers are strongly buckled [34]. This may explain why the hyperfine
couplings to V4+ spins with the opposite alignment do not lead to the
cancellation of the local field. On the other hand, in NaZnVOPO4(HPO4) the
layers are almost flat (Fig. 1), and the absence of the 31P NMR line
broadening below $T_{\rm N}$ can be still ascribed to the filtering of the
hyperfine fields due to the AFM spin alignments in the ground state.

####  III.4.2 31P NMR Shift

Figure 10: Anisotropic components of $K$ along the $x$, $y$, and
$z$-directions as a function of temperature for (a) P(1) and (b) P(2) sites,
respectively. $K$ vs $\chi$ measured at 7 T are plotted for all the three
orientations for (c) P(1) and (d) P(2) sites, respectively. The solid lines
are the straight line fits.

The spectrum at each temperature was fitted following the same procedure as
described above for $T=12$ K. The estimated temperature-dependent NMR shift
[$K(T)$] along different orientations ($K_{\rm x}$, $K_{\rm y}$, and $K_{\rm
z}$) for both P(1) and P(2) sites are plotted in Fig. 10(a) and (b),
respectively. All the components of $K(T)$ exhibit a broad maximum/minimum at
around 4.3 K, an indication of the 2D AFM short-range-ordering. Further, the
magnitude of $K(T)$ for the P(2) site is weaker than the P(1) site, as
expected. The isotropic NMR shift was calculated as $K_{\rm
iso}=(K_{x}+K_{y}+K_{z})/3$, which is found to be almost temperature-
independent for both the 31P sites. This also suggests that the isotropic part
of hyperfine coupling at the P site from the four neighboring V4+ spins is
nearly averaged out.

As $K(T)$ is an intrinsic measure of the spin susceptibility $\chi_{\rm
spin}(T)$ and is free from extrinsic contributions, one can write

$K(T)=K_{0}+\frac{A_{\rm hf}}{N_{\rm A}}\chi_{\rm spin}(T).$ (9)

Here, $K_{0}$ is the temperature-independent chemical (orbital) shift and
$A_{\rm hf}$ is the hyperfine coupling constant between the 31P nucleus and
V4+ spins. In order to estimate $A_{\rm hf}$, we have plotted $K$ vs $\chi$,
assuming temperature as an implied variable in Fig. 10(c) and (d) for P(1) and
P(2) sites, respectively. In every case, $K$ vs $\chi$ plot is linear in the
whole temperature range. A straight line fit results $A^{(1)}_{\rm
hf}=(1134\pm 19)$ Oe/$\mu_{\rm B}$, $(-1716\pm 20)$ Oe/$\mu_{\rm B}$, and
$(582.59\pm 11)$ Oe/$\mu_{\rm B}$ along the $x$, $y$ and $z$-directions for
the P(1) site and $A^{(2)}_{\rm hf}=(343\pm 13)$ Oe/$\mu_{\rm B}$, $(-248\pm
4)$ Oe/$\mu_{\rm B}$, and $(-192\pm 3)$ Oe/$\mu_{\rm B}$ along the $x$, $y$
and $z$-directions for the P(2) site, respectively. Clearly, the magnitude of
$A^{(1)}_{\rm hf}$ is almost one order of magnitude larger than $A^{(2)}_{\rm
hf}$ in all the three directions proving the stronger coupling for P(1) than
P(2).

For the P1 site, the transferred hyperfine coupling mainly arises from the
interactions with the four nearest-neighbor V4+ spins in the plane. The
isotropic and anisotropic transferred hyperfine couplings originate from
P(3s)-O(2p)-V(3d) and P(3p)-O(2p)-V(3d) covalent bonds, respectively. Since P1
is surrounded by four V4+ ions forming a nearly square lattice in the plane,
the experimentally observed asymmetry in hyperfine field indicates
inequivalent P(3p)-O(2p)-V(3d) bonds for the four NN V4+ ions and hence a
distortion in the square lattice, consistent with the low symmetry of the
crystal structure as pointed out earlier.

####  III.4.3 Spin-lattice relaxation rate ${}^{31}1/T_{1}$

Figure 11: (a) Longitudinal recovery curves at three selective temperatures
and the solid lines are fits using Eq. (10). (b) 31P spin-lattice relaxation
rates for P(1) ($1/T_{11}$) and P(2) ($1/T_{12}$) sites as a function of
temperature measured in 30.2 MHz and 121 MHz. The data are shown in log-log
scale in order to highlight the peak at $T_{\rm N}$. The vertical dash-dotted
lines indicate the magnetic LRO at $T_{\rm N}\simeq 2.04$ K and $T_{\rm
N}\simeq 2.2$ K for 30.2 MHz and 121 MHz data, respectively.

We have performed 31P spin-lattice relaxation rate ($1/T_{1}$) measurements as
a function of temperature down to 1.7 K in two different frequencies, 30.2 MHz
and 121 MHz. The recovery of the longitudinal nuclear magnetization after a
saturation pulse could be fitted by a double exponential function

$1-\frac{M(t)}{M(\infty)}=Ae^{-t/T_{11}}+Be^{-t/T_{12}},$ (10)

where $M(t)$ is the nuclear magnetization at a time $t$ after the saturation
pulse and $M(\infty)$ is the equilibrium nuclear magnetization. As each NMR
spectrum is a superposition of two P-sites, which are inseparable, a double
exponential function is used to fit the recovery curves. In Eq. (10),
$1/T_{11}$ and $1/T_{12}$ are the spin-lattice relaxation rates for P(1) and
P(2) sites, respectively, and $A$ and $B$ account for their respective weight
factors. At 121 MHz, the total spectral width was large and we were not able
to saturate the whole spectrum using a single saturation pulse. Therefore,
$1/T_{1}$ is also measured at a lower frequency of 30.2 MHz where the spectral
width is reduced significantly and we were able to saturate the whole spectrum
above $T_{\rm N}$.

The extracted $1/T_{11}$ and $1/T_{12}$ as a function of temperature at 30.2
MHz and 121 MHz and for both the P-sites are presented in Fig. 11(b). In the
high-temperature ($T\geq 4.3$ K) region, $1/T_{1}$ for both the P-sites are
almost constant, typically expected in the paramagnetic regime [69]. At low
temperatures, $1/T_{11}$ and $1/T_{12}$ show sharp peaks at around $T_{\rm
N}\simeq 2.04$ K and 2.36 K for 30.2 MHz and 121 MHz, respectively, indicating
the slowing down of fluctuating moments as we approach the magnetic LRO. Below
$T_{\rm N}$, both $1/T_{11}$ and $1/T_{12}$ decrease toward zero due to the
scattering of magnons by the nuclear spins [70].

###  III.5 Microscopic magnetic model

Before discussing the magnetic model, we determine the position of hydrogen,
which may be crucial for the correct evaluation of exchange couplings [71]. In
hydrophosphates, one expects a deformation of PO4 tetrahedra because one of
the oxygens is linked to hydrogen and should thus weaken its bond to
phosphorous in order to keep the overall bond valence unchanged. Such a
deformation is observed in the P(2)O4 tetrahedra with the P–O bond distances
of 1.490 [O(1)], 1.530 [O(7)], 1.537 [O(2)], and 1.585 Å [O(9)] [46]. We thus
considered different hydrogen positions in the vicinity of O(9) and found the
lowest energy for hydrogen located at (0.7658,0.2823,0.7588) and separated
from oxygen by 0.988 Å. The O–H bond is directed toward the [VOPO4] layer,
with the H atom forming three hydrogen bonds of $2.2-2.4$ Å to oxygen atoms of
the VO5 pyramids and P(1)O4 tetrahedra (Fig. 1).

Table 1:  Interatomic V–V distances (in Å) and exchange couplings (in K). The FM and AFM contributions are obtained from the superexchange model, Eq. (11), whereas total exchange couplings $J_{i}$ are calculated via the DFT+$U$ mapping procedure and thus deviate from $J_{i}^{\rm FM}+J_{i}^{\rm AFM}$.  |  $d_{\rm V-V}$ | $J_{i}^{\rm AFM}$ |  $J_{i}^{\rm FM}$ | $J_{i}$
---|---|---|---|---
$J_{1}$ | 4.616 | 1.4 |  $-5.6$ | $-5.1$
$J_{1}^{\prime}$ | 4.542 | 1.4 |  $-7.4$ | $-7.9$
$J_{1}^{\prime\prime}$ | 4.785 | 2.6 |  $-3.1$ | $-1.0$
$J_{2}$ | 6.380 | 18.7 |  $-0.4$ | 8.1
$J_{2}^{\prime}$ | 6.395 | 9.1 |  $-0.4$ | 4.6

We now adopt this hydrogen position in DFT calculations and evaluate exchange
couplings using two complementary methods. The first one is superexchange
model based on electron hoppings extracted by Wannier fits to the PBE band
structure, resulting in [51]

$J_{i}=\frac{4t_{xy\rightarrow xy}^{2}}{U_{\rm
eff}}-\sum_{\alpha}\frac{4t_{xy\rightarrow\alpha}^{2}J_{\rm eff}}{(U_{\rm
eff}+\Delta_{\alpha})(U_{\rm eff}+\Delta_{\alpha}-J_{\rm eff})},$ (11)

where $d_{xy}$ is the half-filled orbital of V4+, $\alpha$ labels unoccupied
$d$-orbitals, and $\Delta_{\alpha}$ is the crystal-field splitting. The first
term stands for AFM superexchange arising from electron hoppings between the
half-filled orbitals, whereas the second term is FM superexchange due to
electron hoppings between the half-filled and empty orbitals. Using the
effective Coulomb repulsion of $U_{\rm eff}=4$ eV and Hund’s coupling of
$J_{\rm eff}=1$ eV 111Note that $U_{\rm eff}$ and $J_{\rm eff}$ apply to
V-based Wannier functions that include both V $3d$ and O $2p$ orbitals, while
the DFT+$U$ parameters $U_{d}$ and $J_{d}$ mentioned in Sec. II describe,
respectively, the Coulomb repulsion and Hund’s coupling in V $3d$ orbitals
only. from Refs. [52, 54], we arrive at the exchange couplings $J_{i}^{\rm
AFM}$ and $J_{i}^{\rm FM}$ listed in Table 1. Additionally, we compute total
exchange couplings $J_{i}$ via the DFT+$U$ mapping approach [53, 54].

Both methods arrive at qualitatively similar results. The NN couplings are FM,
whereas the NNN couplings are AFM in nature, thus favoring solution-$b$. The
largest anisotropy is observed between $J_{1}^{\prime}$ and
$J_{1}^{\prime\prime}$. However, on average these couplings – both running
along $c$ – differ from $J_{1}$ (along $b$) by 20% only. The spatial
anisotropy of the NNN couplings is more significant with
$J_{2}/J_{2}^{\prime}\simeq 1.8$. The overall energy scale of the NNN
couplings, $\bar{J}_{2}\simeq 5.6$ K, is similar to Li2VOSiO4 ($J_{2}\simeq
5.9$ K [56]) and Na1.5VOPO4F0.5 ($\bar{J}_{2}\simeq 6.6$ K) and much lower
than in Pb2VO(PO${}_{4})_{2}$ with $\bar{J}_{2}\simeq 9.3$ K [14] or
SrZnVO(PO${}_{4})_{2}$ with $\bar{J}_{2}=8.6$ K [17]. This difference can be
traced back to the buckling of the magnetic layers in
$AA^{\prime}$VO(PO${}_{4})_{2}$ ($AA^{\prime}$ = Pb2, SrZn), while the layers
are almost or even perfectly flat in NaZnVOPO4(HPO4) (Fig. 1), Li2VOSiO4, and
Na1.5VOPO4F0.5.

##  IV Discussion

Our data suggest that NaZnVOPO4(HPO4) is well described by the
spin-$\frac{1}{2}$ FSL model if deformation of the square lattice is taken
into account. Although individual data sets such as field-dependent
magnetization may be consistent even with a nonfrustrated square lattice, the
reduced size of the specific heat maximum, $C_{\rm max}^{\rm mag}/R\simeq
0.41$, indicates the presence of frustration. Indeed, the value of the
frustration parameter $f=\frac{|\theta_{\rm CW}|}{T_{\rm N}}\simeq 2.1$
reflects a moderate frustration in the compound.

In a spin system, the spin-lattice relaxation rate carries information on the
low-lying excitations or spin dynamics in the momentum space. Typically,
$\frac{1}{T_{1}T}$ can be expressed in terms of the dynamic susceptibility
$\chi_{M}(\vec{q},\omega_{0})$ as [73]

$\frac{1}{T_{1}T}=\frac{2\gamma_{N}^{2}k_{B}}{N_{\rm
A}^{2}}\sum\limits_{\vec{q}}\mid
A(\vec{q})\mid^{2}\frac{\chi^{{}^{\prime\prime}}_{M}(\vec{q},\omega_{0})}{\omega_{0}},$
(12)

where the sum is over the wave vector $\vec{q}$ within the first Brillouin
zone, $A(\vec{q})$ is the form-factor of the hyperfine interaction, and
$\chi^{{}^{\prime\prime}}_{M}(\vec{q},\omega_{0})$ is the imaginary part of
the dynamic susceptibility at the nuclear Larmor frequency $\omega_{0}$. The
dominance of different $q$-components [$\vec{q}=0$ and
$\vec{q}=(\pm\pi/a,\pm\pi/b)$] is often visible in the $1/T_{1}$ data when
plotted against temperature, especially for the low-dimensional spin systems
with strong exchange coupling. At very high temperatures ($T>J/k_{\rm B}$),
$1/T_{1}$ is almost temperature-independent due to uncorrelated moments and
can be expressed as [14]

$\displaystyle\begin{split}\lx@intercol\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\left(\frac{1}{T_{1}}\right)_{T\to\infty}=\frac{(\gamma_{N}g\mu_{\rm
B})^{2}\sqrt{2\pi}z^{\prime}S(S+1)}{3\omega_{\rm
ex}}$\mbox{}\hfil\lx@intercol\\\
&\qquad&\qquad&\qquad\times\frac{\left(\frac{A_{x}}{z^{\prime}}\right)^{2}+\left(\frac{A_{y}}{z^{\prime}}\right)^{2}+\left(\frac{A_{z}}{z^{\prime}}\right)^{2}}{3},\end{split}$
(13)

where $\omega_{\rm ex}=[\max(|J_{1}|,|J_{2}|)k_{\rm
B}/\hbar]\sqrt{\frac{2zS(S+1)}{3}}$ is the Heisenberg exchange frequency, $z$
is the number of NN spins of each V4+ ion, and $z^{\prime}$ is the number of
NN V4+ spins of the P(1) site. In the above expression, the hyperfine
couplings along different directions are divided by $z^{\prime}$ in order to
account for the coupling of P site with the individual V4+ ion. As the
measurements are carried out on the powder sample, we have taken the rms
average of the couplings along three directions.

For a tentative estimation of the in-plane exchange coupling between the V4+
ions we took the high temperature value of $1/T_{11}$ for the strongly coupled
P(1) site. Using the experimental parameters obtained for this site
($A^{1}_{x}\simeq 1134$ Oe/$\mu_{\rm B}$, $A^{1}_{y}\simeq-1716$ Oe/$\mu_{\rm
B}$, $A^{1}_{z}\simeq 582$ Oe/$\mu_{\rm B}$, $\gamma_{N}=108.303\times 10^{2}$
rad.sec-1/Oe, $z^{\prime}=4$, $z=4$, $g=2.04$, $S=1/2$, and $1/T_{11}\simeq
150$ sec-1), the magnitude of the maximum exchange coupling strength between
V4+ ions is estimated to be $\max(|J_{1}|,|J_{2}|)\simeq 4.2$ K. This value is
indeed very close to the dominant in-plane AFM exchange coupling for both the
solutions ($a$ and $b$), obtained from the $\chi(T)$ analysis.

It is also instructive to compare the transition temperatures of V4+-based FSL
magnets. NaZnVOPO4(HPO4) and Li2VOSiO4 feature almost the same
$\bar{J}_{2}\simeq 6$ K but different signs of $J_{1}$ (FM and AFM,
respectively), while the former compound has a much lower $T_{N}\simeq 2.1$ K
than the latter ($T_{N}\simeq 2.8$ K [12]). This indicates that FM couplings
$J_{1}$ together with the deformation of the square lattice lead to a visible
reduction in $T_{N}$ in NaZnVOPO4(HPO4). On the other hand,
Pb2VO(PO${}_{4})_{2}$ shows an even higher $T_{N}\simeq 3.65$ K [14] because
of the 30% increase in the magnitude of $\bar{J}_{2}$. This comparison
illustrates that buckling of the [VOPO4] layers caused by variable spacers
between the magnetic layers effectively tunes magnetic interactions in FSL-
like compounds.

##  V Summary

We have studied the magnetism of NaZnVOPO4(HPO4) using wide variety of
experimental techniques and complementary ab initio calculations. With almost
flat magnetic layers, this compound shows weaker next-nearest-neighbor
couplings ($\bar{J}_{2}$) than other V4+ phosphates. This reduction in
$\bar{J}_{2}$ along with the deformation of the spin-$\frac{1}{2}$ square
lattice lead to a lower magnetic transition temperature than in most of the
other FSL candidates.

##  VI Acknowledgments

SG and RN would like to acknowledge SERB, India for financial support bearing
sanction Grant No. CRG/2019/000960. SG is supported by the Prime Minister’s
Research Fellowship (PMRF) scheme, Government of India. We also acknowledge
the support of the HLD at HZDR, member of European Magnetic Field Laboratory
(EMFL). Work at the Ames Laboratory was supported by the U.S. Department of
Energy, Office of Science, Basic Energy Sciences, Materials Sciences and
Engineering Division. The Ames Laboratory is operated for the U.S. Department
of Energy by Iowa State University under Contract No. DEAC02-07CH11358.

## References

  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80, 016502 (2016).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Shannon  _et al._ [2004] N. Shannon, B. Schmidt, K. Penc, and P. Thalmeier, Finite temperature properties and frustrated ferromagnetism in a square lattice Heisenberg model, Eur. Phys. J. B 38, 599 (2004).
  * Schmidt and Thalmeier [2017] B. Schmidt and P. Thalmeier, Néel temperature and reentrant $H\text{$-$}T$ phase diagram of quasi-two-dimensional frustrated magnets, Phys. Rev. B 96, 214443 (2017).
  * Zhang  _et al._ [2003] G.-M. Zhang, H. Hu, and L. Yu, Valence-Bond Spin-Liquid State in Two-Dimensional Frustrated Spin-$1/2$ Heisenberg Antiferromagnets, Phys. Rev. Lett. 91, 067201 (2003).
  * Hu  _et al._ [2013] W.-J. Hu, F. Becca, A. Parola, and S. Sorella, Direct evidence for a gapless ${Z}_{2}$ spin liquid by frustrating Néel antiferromagnetism, Phys. Rev. B 88, 060402(R) (2013).
  * Wang and Sandvik [2018] L. Wang and A. W. Sandvik, Critical Level Crossings and Gapless Spin Liquid in the Square-Lattice Spin-$1/2$ ${J}_{1}-{J}_{2}$ Heisenberg Antiferromagnet, Phys. Rev. Lett. 121, 107202 (2018).
  * Gong  _et al._ [2014] S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Plaquette Ordered Phase and Quantum Phase Diagram in the Spin-$\frac{1}{2}$ ${J}_{1}\text{$-$}{J}_{2}$ Square Heisenberg Model, Phys. Rev. Lett. 113, 027201 (2014).
  * Doretto [2014] R. L. Doretto, Plaquette valence-bond solid in the square-lattice ${J}_{1}$-${J}_{2}$ antiferromagnet Heisenberg model: A bond operator approach, Phys. Rev. B 89, 104415 (2014).
  * Haghshenas and Sheng [2018] R. Haghshenas and D. N. Sheng, U(1)-symmetric infinite projected entangled-pair states study of the spin-1/2 square ${J}_{1}\text{$-$}{J}_{2}$ Heisenberg model, Phys. Rev. B 97, 174408 (2018).
  * Shannon  _et al._ [2006] N. Shannon, T. Momoi, and P. Sindzingre, Nematic Order in Square Lattice Frustrated Ferromagnets, Phys. Rev. Lett. 96, 027213 (2006).
  * Melzi  _et al._ [2000] R. Melzi, P. Carretta, A. Lascialfari, M. Mambrini, M. Troyer, P. Millet, and F. Mila, Li2VO(Si,Ge)O4, a Prototype of a Two-Dimensional Frustrated Quantum Heisenberg Antiferromagnet, Phys. Rev. Lett. 85, 1318 (2000).
  * Rosner  _et al._ [2002] H. Rosner, R. R. P. Singh, W. H. Zheng, J. Oitmaa, S.-L. Drechsler, and W. E. Pickett, Realization of a large $J_{2}$ quasi-2D spin-half Heisenberg system: Li2VOSiO4, Phys. Rev. Lett. 88, 186405 (2002).
  * Nath  _et al._ [2009] R. Nath, Y. Furukawa, F. Borsa, E. E. Kaul, M. Baenitz, C. Geibel, and D. C. Johnston, Single-crystal ${}^{31}\text{P}$ NMR studies of the frustrated square-lattice compound Pb2VO(PO4)2, Phys. Rev. B 80, 214430 (2009).
  * Roy  _et al._ [2011] B. Roy, Y. Furukawa, R. Nath, and D. C. Johnston, Low-temperature 31P NMR study of the two-dimensional frustrated square lattice compound BaCdVO(PO4)2, J. Phys. Conf. Ser 320, 012048 (2011).
  * Tsirlin  _et al._ [2010] A. A. Tsirlin, R. Nath, A. M. Abakumov, R. V. Shpanchenko, C. Geibel, and H. Rosner, Frustrated square lattice with spatial anisotropy: Crystal structure and magnetic properties of $\text{PbZnVO}{({\text{PO}}_{4})}_{2}$, Phys. Rev. B 81, 174424 (2010).
  * Bossoni  _et al._ [2011] L. Bossoni, P. Carretta, R. Nath, M. Moscardini, M. Baenitz, and C. Geibel, NMR and $\mu$SR study of spin correlations in SrZnVO(PO4)2: An $S=\frac{1}{2}$ frustrated magnet on a square lattice, Phys. Rev. B 83, 014412 (2011).
  * Tsirlin  _et al._ [2011a] A. A. Tsirlin, R. Nath, A. M. Abakumov, Y. Furukawa, D. C. Johnston, M. Hemmida, H.-A. Krug von Nidda, A. Loidl, C. Geibel, and H. Rosner, Phase separation and frustrated square lattice magnetism of Na1.5VOPO4F0.5, Phys. Rev. B 84, 014429 (2011a).
  * Nath  _et al._ [2008] R. Nath, A. A. Tsirlin, H. Rosner, and C. Geibel, Magnetic properties of $\text{BaCdVO}{({\text{PO}}_{4})}_{2}$: A strongly frustrated spin-$\frac{1}{2}$ square lattice close to the quantum critical regime, Phys. Rev. B 78, 064422 (2008).
  * Yogi  _et al._ [2015] A. Yogi, N. Ahmed, R. Nath, A. A. Tsirlin, S. Kundu, A. V. Mahajan, J. Sichelschmidt, B. Roy, and Y. Furukawa, Antiferromagnetism of Zn2VO(PO4)2 and the dilution with Ti4+, Phys. Rev. B 91, 024413 (2015).
  * Kiani and Pavarini [2016] A. Kiani and E. Pavarini, Electronic correlation and magnetic frustration in Li2VOSiO4 and VOMoO4, Phys. Rev. B 94, 075112 (2016).
  * Bombardi  _et al._ [2005] A. Bombardi, L. C. Chapon, I. Margiolaki, C. Mazzoli, S. Gonthier, F. Duc, and P. G. Radaelli, Magnetic order and lattice anomalies in the $J_{1}-{J}_{2}$ model system VOMoO4, Phys. Rev. B 71, 220406 (2005).
  * Carretta  _et al._ [2002] P. Carretta, N. Papinutto, C. B. Azzoni, M. C. Mozzati, E. Pavarini, S. Gonthier, and P. Millet, Frustration-driven structural distortion in VOMoO4, Phys. Rev. B 66, 094420 (2002).
  * Tsirlin  _et al._ [2008] A. A. Tsirlin, A. A. Belik, R. V. Shpanchenko, E. V. Antipov, E. Takayama-Muromachi, and H. Rosner, Frustrated spin-1/2 square lattice in the layered perovskite PbVO3, Phys. Rev. B 77, 092402 (2008).
  * Ishikawa  _et al._ [2017] H. Ishikawa, N. Nakamura, M. Yoshida, M. Takigawa, P. Babkevich, N. Qureshi, H. M. Rønnow, T. Yajima, and Z. Hiroi, $J_{1}-J_{2}$ square-lattice Heisenberg antiferromagnets with $4{d}^{1}$ spins: $A$MoOPO4Cl ($A$ = K,Rb), Phys. Rev. B 95, 064408 (2017).
  * Takeda  _et al._ [2021] H. Takeda, T. Yamauchi, M. Takigawa, H. Ishikawa, and Z. Hiroi, Pressure-induced phase transition in the ${J}_{1}\text{$-$}{J}_{2}$ square lattice antiferromagnet RbMoOPO4Cl, Phys. Rev. B 103, 104406 (2021).
  * Xu  _et al._ [2017] Y. Xu, S. Liu, N. Qu, Y. Cui, Q. Gao, R. Chen, J. Wang, F. Gao, and X. Hao, Comparative description of magnetic interactions in Sr2CuTeO6 and Sr2CuWO6, J. Phys. Condens. Matter 29, 105801 (2017).
  * Vasala  _et al._ [2014] S. Vasala, M. Avdeev, S. Danilkin, O. Chmaissem, and M. Karppinen, Magnetic structure of Sr2CuWO6, J. Phys. Condens. Matter 26, 496001 (2014).
  * Koga  _et al._ [2016] T. Koga, N. Kurita, M. Avdeev, S. Danilkin, T. J. Sato, and H. Tanaka, Magnetic structure of the $S=\frac{1}{2}$ quasi-two-dimensional square-lattice Heisenberg antiferromagnet Sr2CuTeO6, Phys. Rev. B 93, 054426 (2016).
  * Babkevich  _et al._ [2016] P. Babkevich, V. M. Katukuri, B. Fåk, S. Rols, T. Fennell, D. Pajić, H. Tanaka, T. Pardini, R. R. P. Singh, A. Mitrushchenkov, O. V. Yazyev, and H. M. Rønnow, Magnetic Excitations and Electronic Interactions in Sr2CuTeO6: A Spin-$1/2$ Square Lattice Heisenberg Antiferromagnet, Phys. Rev. Lett. 117, 237203 (2016).
  * Walker  _et al._ [2016] H. C. Walker, O. Mustonen, S. Vasala, D. J. Voneshen, M. D. Le, D. T. Adroja, and M. Karppinen, Spin wave excitations in the tetragonal double perovskite Sr2CuWO6, Phys. Rev. B 94, 064411 (2016).
  * Watanabe  _et al._ [2022] M. Watanabe, N. Kurita, H. Tanaka, W. Ueno, K. Matsui, T. Goto, and M. Hagihala, Contrasting magnetic structures in SrLaCuSbO6 and SrLaCuNbO6: Spin-$\frac{1}{2}$ quasi-square-lattice $J_{1}-J_{2}$ Heisenberg antiferromagnets, Phys. Rev. B 105, 054414 (2022).
  * Guchhait  _et al._ [2021a] S. Guchhait, S. Baby, M. Padmanabhan, A. Medhi, and R. Nath, Quasi-two-dimensional magnetism in spin-1/2 square lattice compound Cu[C6H2(COO)4][H3N-(CH2)2-NH3]$\cdot$3H2O, EPL 133, 57006 (2021a).
  * Tsirlin and Rosner [2009] A. A. Tsirlin and H. Rosner, Extension of the spin-$\frac{1}{2}$ frustrated square lattice model: The case of layered vanadium phosphates, Phys. Rev. B 79, 214417 (2009).
  * Bettler  _et al._ [2019] S. Bettler, F. Landolt, O. M. Aksoy, Z. Yan, S. Gvasaliya, Y. Qiu, E. Ressouche, K. Beauvois, S. Raymond, A. N. Ponomaryov, S. A. Zvyagin, and A. Zheludev, Magnetic structure and spin waves in the frustrated ferro-antiferromagnet $\mathrm{Pb}_{2}\mathrm{VO}{({\mathrm{PO}}_{4})}_{2}$, Phys. Rev. B 99, 184437 (2019).
  * Landolt  _et al._ [2021] F. Landolt, Z. Yan, S. Gvasaliya, K. Beauvois, E. Ressouche, J. Xu, and A. Zheludev, Phase diagram and spin waves in the frustrated ferro-antiferromagnet SrZnVO(PO4)2, Phys. Rev. B 104, 224435 (2021).
  * Bhartiya  _et al._ [2021] V. K. Bhartiya, S. Hayashida, K. Y. Povarov, Z. Yan, Y. Qiu, S. Raymond, and A. Zheludev, Inelastic neutron scattering determination of the spin hamiltonian for $\mathrm{BaCdVO}{({\mathrm{PO}}_{4})}_{2}$, Phys. Rev. B 103, 144402 (2021).
  * Povarov  _et al._ [2019] K. Y. Povarov, V. K. Bhartiya, Z. Yan, and A. Zheludev, Thermodynamics of a frustrated quantum magnet on a square lattice, Phys. Rev. B 99, 024413 (2019).
  * Bhartiya  _et al._ [2019] V. K. Bhartiya, K. Y. Povarov, D. Blosser, S. Bettler, Z. Yan, S. Gvasaliya, S. Raymond, E. Ressouche, K. Beauvois, J. Xu, F. Yokaichiya, and A. Zheludev, Presaturation phase with no dipolar order in a quantum ferro-antiferromagnet, Phys. Rev. Research 1, 033078 (2019).
  * Skoulatos  _et al._ [2019] M. Skoulatos, F. Rucker, G. J. Nilsen, A. Bertin, E. Pomjakushina, J. Ollivier, A. Schneidewind, R. Georgii, O. Zaharko, L. Keller, C. Rüegg, C. Pfleiderer, B. Schmidt, N. Shannon, A. Kriele, A. Senyshyn, and A. Smerald, Putative spin-nematic phase in $\mathrm{BaCdVO}(\mathrm{PO}_{4}{)}_{2}$, Phys. Rev. B 100, 014405 (2019).
  * Landolt  _et al._ [2020] F. Landolt, S. Bettler, Z. Yan, S. Gvasaliya, A. Zheludev, S. Mishra, I. Sheikin, S. Krämer, M. Horvatić, A. Gazizulina, and O. Prokhnenko, Presaturation phase in the frustrated ferro-antiferromagnet ${\mathrm{Pb}}_{2}\mathrm{VO}{({\mathrm{PO}}_{4})}_{2}$, Phys. Rev. B 102, 094414 (2020).
  * Ranjith  _et al._ [2022] K. M. Ranjith, F. Landolt, S. Raymond, A. Zheludev, and M. Horvatić, NMR evidence against a spin-nematic nature of the presaturation phase in the frustrated magnet SrZnVO(PO${}_{4})_{2}$, Phys. Rev. B 105, 134429 (2022).
  * Landolt  _et al._ [2022] F. Landolt, K. Povarov, Z. Yan, S. Gvasaliya, E. Ressouche, S. Raymond, V. O. Garlea, and A. Zheludev, Spin correlations in the frustrated ferro-antiferromagnet SrZnVO(PO${}_{4})_{2}$ near saturation, arXiv preprint arXiv:2204.08778 (2022).
  * Momma and Izumi [2011] K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr. 44, 1272 (2011).
  * Rodríguez-Carvajal [1993] J. Rodríguez-Carvajal, Recent advances in magnetic structure determination by neutron powder diffraction, Physica B Condens. Matter 192, 55 (1993).
  * Le Fur  _et al._ [2001] E. Le Fur, Y. Moreno, and J. Pivan, On the new oxovanadium phosphate NaZnVOPO4(HPO4) obtained at 510 K from hydrothermal treatment, J. Mater. Chem. 11, 1735 (2001).
  * Skourski  _et al._ [2011] Y. Skourski, M. D. Kuz’min, K. P. Skokov, A. V. Andreev, and J. Wosnitza, High-field magnetization of Ho2Fe17, Phys. Rev. B 83, 214420 (2011).
  * Tsirlin  _et al._ [2009] A. A. Tsirlin, B. Schmidt, Y. Skourski, R. Nath, C. Geibel, and H. Rosner, Exploring the spin-$\frac{1}{2}$ frustrated square lattice model with high-field magnetization studies, Phys. Rev. B 80, 132407 (2009).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59, 1743 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Mazurenko  _et al._ [2006] V. V. Mazurenko, F. Mila, and V. I. Anisimov, Electronic structure and exchange interactions of Na2V3O7, Phys. Rev. B 73, 014418 (2006).
  * Tsirlin  _et al._ [2011b] A. A. Tsirlin, O. Janson, and H. Rosner, Unusual ferromagnetic superexchange in CdVO3: The role of Cd, Phys. Rev. B 84, 144429 (2011b).
  * Xiang  _et al._ [2011] H. J. Xiang, E. J. Kan, S.-H. Wei, M.-H. Whangbo, and X. G. Gong, Predicting the spin-lattice order of frustrated systems from first principles, Phys. Rev. B 84, 224429 (2011).
  * Tsirlin [2014] A. A. Tsirlin, Spin-chain magnetism and uniform Dzyaloshinsky-Moriya anisotropy in BaV3O8, Phys. Rev. B 89, 014405 (2014).
  * Weickert  _et al._ [2016] F. Weickert, N. Harrison, B. L. Scott, M. Jaime, A. Leitmäe, I. Heinmaa, R. Stern, O. Janson, H. Berger, H. Rosner, and A. A. Tsirlin, Magnetic anisotropy in the frustrated spin-chain compound $\beta$-TeVO4, Phys. Rev. B 94, 064403 (2016).
  * Rosner  _et al._ [2003] H. Rosner, R. R. P. Singh, W. H. Zheng, J. Oitmaa, and W. E. Pickett, High-temperature expansions for the ${J}_{1}-{J}_{2}$ Heisenberg models: Applications to ab initio calculated models for Li2VOSiO4 and Li2VOGeO4, Phys. Rev. B 67, 014416 (2003).
  * Alet  _et al._ [2005] F. Alet, S. Wessel, and M. Troyer, Generalized directed loop method for quantum monte carlo simulations, Phys. Rev. E 71, 036706 (2005).
  * Sandvik [1999] A. W. Sandvik, Stochastic series expansion method with operator-loop update, Phys. Rev. B 59, R14157 (1999).
  * Bauer  _et al._ [2011] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, _et al._ , The ALPS project release 2.0: open source software for strongly correlated systems, J. Stat. Mech. Theory Exp. 2011, P05001 (2011).
  * Wallace [1998] D. C. Wallace,  _Thermodynamics of Crystals_ (Dover Publications, 1998).
  * Budd  _et al._ [2015] L. E. Budd, R. M. Ibberson, W. G. Marshall, and S. Parsons, The effect of temperature and pressure on the crystal structure of piperidine, Chem. Cent. J. 9, 18 (2015).
  * Guchhait  _et al._ [2021b] S. Guchhait, Q.-P. Ding, M. Sahoo, A. Giri, S. Maji, Y. Furukawa, and R. Nath, Quasi-one-dimensional uniform spin-$\frac{1}{2}$ Heisenberg antiferromagnet ${\mathrm{KNaCuP}}_{2}{\mathrm{O}}_{7}$ probed by ${}^{31}\mathrm{P}$ and ${}^{23}\mathrm{Na}$ NMR, Phys. Rev. B 103, 224415 (2021b).
  * Kittel [2004] C. Kittel,  _Introduction to Solid State Physics_ (Wiley, 2004).
  * Selwood [2013] P. W. Selwood,  _Magnetochemistry_ (Read Books Ltd, 2013).
  * Mukharjee  _et al._ [2020] P. K. Mukharjee, K. M. Ranjith, M. Baenitz, Y. Skourski, A. A. Tsirlin, and R. Nath, Two types of alternating spin-$\frac{1}{2}$ chains and their field-induced transitions in $\varepsilon\text{$-$}{\mathrm{LiVOPO}}_{4}$, Phys. Rev. B 101, 224403 (2020).
  * Thalmeier  _et al._ [2008] P. Thalmeier, M. E. Zhitomirsky, B. Schmidt, and N. Shannon, Quantum effects in magnetization of $J_{1}-J_{2}$ square lattice antiferromagnet, Phys. Rev. B 77, 104441 (2008).
  * Sebastian  _et al._ [2021] S. J. Sebastian, K. Somesh, M. Nandi, N. Ahmed, P. Bag, M. Baenitz, B. Koo, J. Sichelschmidt, A. A. Tsirlin, Y. Furukawa, and R. Nath, Quasi-one-dimensional magnetism in the spin-$\frac{1}{2}$ antiferromagnet ${\mathrm{BaNa}}_{2}\mathrm{Cu}{({\mathrm{VO}}_{4})}_{2}$, Phys. Rev. B 103, 064413 (2021).
  * Melzi  _et al._ [2001] R. Melzi, S. Aldrovandi, F. Tedoldi, P. Carretta, P. Millet, and F. Mila, Magnetic and thermodynamic properties of Li2VOSiO4: A two-dimensional $S=1/2$ frustrated antiferromagnet on a square lattice, Phys. Rev. B 64, 024409 (2001), cond-mat/0101066.
  * Moriya [1956] T. Moriya, Nuclear magnetic relaxation in antiferromagnetics, Prog. Theor. Phys. 16, 23 (1956).
  * Islam  _et al._ [2018] S. S. Islam, K. M. Ranjith, M. Baenitz, Y. Skourski, A. A. Tsirlin, and R. Nath, Frustration of square cupola in Sr(TiO)${\mathrm{Cu}}_{4}({\mathrm{PO}}_{4}{)}_{4}$, Phys. Rev. B 97, 174432 (2018).
  * Lebernegg  _et al._ [2013] S. Lebernegg, A. A. Tsirlin, O. Janson, and H. Rosner, Two energy scales of spin dimers in clinoclase Cu3(AsO4)(OH)3, Phys. Rev. B 87, 235117 (2013).
  * Note [1] Note that $U_{\rm eff}$ and $J_{\rm eff}$ apply to V-based Wannier functions that include both V $3d$ and O $2p$ orbitals, while the DFT+$U$ parameters $U_{d}$ and $J_{d}$ mentioned in Sec. II describe, respectively, the Coulomb repulsion and Hund’s coupling in V $3d$ orbitals only.
  * Moriya [1963] T. Moriya, The effect of electron-electron interaction on the nuclear spin relaxation in metals, J. Phys. Soc. Jpn. 18, 516 (1963).