# Two types of alternating spin-$\frac{1}{2}$ chains
and their field-induced transitions in $\varepsilon$-LiVOPO4

Prashanta K. Mukharjee  School of Physics, Indian Institute of Science
Education and Research, Thiruvananthapuram-695551, India  K. M. Ranjith  M.
Baenitz  Max Planck Institute for Chemical Physics of Solids,
N$\ddot{o}$thnitzer Str. 40, 01187 Dresden, Germany  Y. Skourski  Dresden High
Magnetic Field Laboratory (HLD-EMFL), Helmholtz-Zentrum Dresden-Rossendorf,
01314 Dresden, Germany  A. A. Tsirlin<EMAIL_ADDRESS>Experimental
Physics VI, Center for Electronic Correlations and Magnetism, Institute of
Physics, University of Augsburg, 86135 Augsburg, Germany  R. Nath
<EMAIL_ADDRESS>School of Physics, Indian Institute of Science Education
and Research, Thiruvananthapuram-695551, India

###### Abstract

Thermodynamic properties, 31P nuclear magnetic resonance (NMR) measurements,
and density-functional band-structure calculations for $\varepsilon$-LiVOPO4
are reported. This quantum magnet features a singlet ground state and
comprises two types of alternating spin-$\frac{1}{2}$ chains that manifest
themselves by the double maxima in the susceptibility and magnetic specific
heat, and by the two-step magnetization process with an intermediate
$\frac{1}{2}$-plateau. From thermodynamic data and band-structure
calculations, we estimate the leading couplings of $J_{1}\simeq 20$ K and
$J_{2}\simeq 60$ K and the alternation ratios of
$\alpha_{1}=J_{1}^{\prime}/J_{1}\simeq 0.6$ and
$\alpha_{2}=J_{2}^{\prime}/J_{2}\simeq 0.3$ within the two chains,
respectively. The zero-field spin gap $\Delta_{0}/k_{\rm B}\simeq 7.3$ K
probed by thermodynamic and NMR measurements is caused by the
$J_{1}$-$J_{1}^{\prime}$ spin chains and can be closed in the applied field of
$\mu_{0}H_{\rm c1}\simeq 5.6$ T, giving rise to a field-induced long-range
order. The NMR data reveal predominant three-dimensional spin-spin
correlations at low temperatures. Field-induced magnetic ordering transition
observed above $H_{c1}$ is attributed to the Bose-Einstein condensation of
triplons in the sublattice formed by the $J_{1}$-$J_{1}^{\prime}$ chains with
weaker exchange couplings.

##  I Introduction

Field-induced quantum phase transitions in magnets set a link between
fermionic spin systems and lattice boson gas Matsubara and Matsuda (1956);
Giamarchi  _et al._ (2008); Zapf  _et al._ (2014). In this context, spin-dimer
compounds possessing a gap in the excitation spectrum are extensively studied
Giamarchi  _et al._ (2008); Zapf  _et al._ (2014). Their triplet excitations
(triplons) are equivalent to lattice bosons and can be tuned by applying
magnetic field. Gap closing in the applied field will usually lead to magnetic
ordering that can be understood as Bose-Einstein condensation (BEC) of
triplons Mukhopadhyay  _et al._ (2012); Sebastian  _et al._ (2006); Giamarchi
_et al._ (2008); Nikuni  _et al._ (2000). If interactions between the dimers
are one-dimensional (1D) in nature, one further expects a non-Fermi-liquid-
type Tomonaga-Luttinger Liquid (TLL) state realized at intermediate
temperatures before the BEC state is reached Klanjšek  _et al._ (2008);
Matsushita  _et al._ (2017); Willenberg  _et al._ (2015); Thielemann  _et al._
(2009).

The ground-state spin configuration in the applied field may also depend on
the delicate balance between the kinetic energy and repulsive interaction of
the triplons or $S_{\rm z}=+1$ bosons Rice (2002). The dominance of repulsive
interaction would lead to the formation of superlattices, which result in
magnetization plateaus Narumi  _et al._ (1998); Shiramura  _et al._ (1998).
This has been experimentally verified in the celebrated Shastry-Sutherland
compound SrCu2(BO3)2 Kodama  _et al._ (2002); Kageyama  _et al._ (1999). On
the other hand, when the kinetic energy dominates over repulsive interactions,
the triplons become delocalized, and the ground state is a superposition of
singlet-triplet states, which can be approximated as a BEC of triplons. The
phenomenon of BEC has been studied in detail for spin-dimer compounds
TlCuCl3Rüegg  _et al._ (2003), BaCuSi2O6Jaime  _et al._ (2004), (Ba,Sr)3Cr2O8
Aczel  _et al._ (2009a, b), etc. The transition from TLL as a 1D quantum
critical state to the three-dimensional (3D) BEC state has often been observed
in quasi-1D spin systems e.g. spin-$1/2$ ladders (C7H10N)2CuBr4 Jeong  _et
al._ (2013, 2017); Möller  _et al._ (2017) and (C5H12N)2CuBr4 Thielemann  _et
al._ (2009); Rüegg  _et al._ (2008); Klanjšek  _et al._ (2008) and spin-$1/2$
alternating spin chains Cu(NO3)${}_{2}\cdot$2.5D2O and F5PNN Willenberg  _et
al._ (2015); Matsushita  _et al._ (2017). Thus, quasi-1D spin gap materials
provide ample opportunities to tune the spin gap and study the field-induced
quantum phase transitions.

The V4+-based compounds offer an excellent playground to study gapped quantum
magnets and related phenomena. Several of these compounds were already
reported in the past in the context of spin gap physics Ueda (1998); Yamauchi
_et al._ (1999); Johnston  _et al._ (1987); Ghoshray  _et al._ (2005).
Recently, we studied magnetic properties of the $A$VO$X$O4 series ($A$ = Na,
Ag; $X$ = P, As), where all compounds were found to host alternating
spin-$\frac{1}{2}$ chains not matching the structural chains of the VO6
octahedra Mukharjee  _et al._ (2019); Arjun  _et al._ (2019); Ahmed  _et al._
(2017); Tsirlin  _et al._ (2011). In these systems, long-range superexchange
couplings via two oxygen atoms play a central role. They are highly tunable
and allow the variation of the zero-field spin gap $\Delta_{\rm 0}/k_{\rm B}$
from 21 K in NaVOAsO4 to $\sim 2$ K in NaVOPO4. External magnetic field closes
the gap and leads to a field-induced magnetic transition, which is explained
in terms of the triplon BEC Mukharjee  _et al._ (2019); Weickert  _et al._
(2019).

Herein, we report ground-state properties of the chemically similar, but
structurally different LiVOPO4. Unlike the $A$VO$X$O4 compounds, which are all
monoclinic ($P2_{1}/c$), LiVOPO4 crystallizes in several polymorphs with
different symmetries and atomic arrangements Harrison and Manthiram (2013);
Hidalgo  _et al._ (2019). We shall focus on the triclinic
$\varepsilon$-LiVOPO4 ($P\bar{1}$) that can be seen as a distorted version of
monoclinic $A$VO$X$O4 111Triclinic crystal structure is derived from
$\varepsilon$-VOPO4. Alternatively, triclinic LiVOPO4 is sometimes referred to
as the $\alpha$-phase, because it was the earliest discovered LiVOPO4
polymorph. and has been actively studied in the context of battery research
Yang  _et al._ (2008); Quackenbush  _et al._ (2015); Lin  _et al._ (2016); Shi
_et al._ (2018); Chung  _et al._ (2019), but not as a quantum magnet. In its
crystal structure, each of the Li, V, and P reside at two nonequivalent sites,
whereas the O atoms have ten nonequivalent sites. The magnetic V4+ ions form
chains of the VO6 octahedra with the alternation of V1 and V2, leading to two
alternating V–V distances of 3.599 and 3.629 Å along these structural chains
(Fig. 1).

Assuming that strongest magnetic interactions run along the structural chains,
one expects the magnetic behavior of alternating spin-$\frac{1}{2}$ chains
that was indeed proposed by Onoda and Ikeda Onoda and Ikeda (2013) who
reported, among others Yang  _et al._ (2008); Chung  _et al._ (2019), magnetic
susceptibility of $\varepsilon$-LiVOPO4. On the other hand, our recent results
for the monoclinic $A$VO$X$O4 compounds suggest that spin chains may not
coincide with the structural chains, because leading interactions occur
through the double bridges of the $X$O4 tetrahedra. In this case,
$\varepsilon$-LiVOPO4 should feature two types of alternating
spin-$\frac{1}{2}$ chains, one formed by V(1) and the other one formed by
V(2), each with different interactions and different spin gaps (Fig. 1).
Below, we report experimental fingerprints of these two nonequivalent spin
chains and thus directly confirm the proposed microscopic scenario. Moreover,
we detect a field-induced phase transition, which is reminiscent of triplon
BEC.

Table 1: Atomic distances and bond angles along the superexchange paths involving P(1) and P(2) in the two alternating spin chains of $\varepsilon$-LiVOPO4 at $300$ K to highlight the coupling of P atoms with V atoms. Site | Bond Length (Å) | Angle (deg.)
---|---|---
P(1) |  | V(1)-O(1) = 1.96
---
O(1)-P(1) = 1.53
P(1)-O(2) = 1.54
O(2)-V(1) = 1.97
V(2)-O(7) = 1.98
O(7)-P(1) = 1.53
P(1)-O(8) = 1.53
O(8)-V(2) = 2.01
|  $\angle$V(1)-O(1)-O(2) = 152.16
---
$\angle$V(1)-O(2)-O(1) = 109.37
$\angle$V(2)-O(8)-O(7) = 140.38
$\angle$V(2)-O(7)-O(8) = 124.64
Average value $\simeq$ 131.63
P(2) |  | V(1)-O(3) = 1.96
---
O(3)-P(2) = 1.53
P(2)-O(4) = 1.55
O(4)-V(1) = 2.02
V(2)-O(9) = 1.98
O(9)-P(2) = 1.53
P(2)-O(10) = 1.52
O(10)-V(2) = 1.94
|  $\angle$V(1)-O(4)-O(3) = 118.59
---
$\angle$V(1)-O(3)-O(4) = 149.06
$\angle$V(2)-O(9)-O(10) = 134.87
$\angle$V(2)-O(10)-O(9) = 132.66
Average value $\simeq$ 133.79

Figure 1: (a) Crystal structure of $\varepsilon$-LiVOPO4 projected onto the
$ac$-plane. The deep blue and light blue solid circles represent the V(1) and
V(2) sites, respectively. Chain 1 with the couplings $J_{1}$ and
$J_{1}^{\prime}$ is formed by V(1) and chain 2 is formed by V(2) atoms with
the couplings $J_{2}$ and $J_{2}^{\prime}$ via the extended V-O$\ldots$O-V
paths. These chains are nearly orthogonal to each other, whereas the
structural chains are parallel comprising both V1 and V2 atoms at the same
time. These chains run perpendicular to the $ac$-plane. (b) A segment of the
chain 1 formed by the V(1)O6 octahedra with the intrachain couplings $J_{1}$
and $J_{1}^{\prime}$. The distances $d_{1}$ and $d_{1}^{\prime}$ are the
lateral displacements of the VO6 octahedra in chain 1. (c) A section of the
structural chain with the V(1)-O-V(2) paths along the $b$-axis. (d) An
empirical/qualitative sketch of the spin model with all possible exchange
interactions involving chain 1 ($J_{1}$, $J_{1}^{\prime}$) and chain 2
($J_{2},J_{2}^{\prime}$).

##  II Methods

Polycrystalline sample of $\varepsilon$-LiVOPO4 was prepared by the
conventional solid-state reaction method from stoichiometric mixtures of LiPO3
and VO2 (Aldrich, 99.995%). LiPO3 was obtained by heating LiH2PO4.H2O
(Aldrich, 99.995%) for 4 h at 400 ∘C in air. The reactants were ground
thoroughly, pelletized, and fired at 740 ∘C for two days in flowing argon
atmosphere with two intermediate grindings. Phase purity of the sample was
confirmed by powder x-ray diffraction (XRD) recorded at room temperature using
a PANalytical powder diffractometer (CuKα radiation, $\lambda_{\rm avg}\simeq
1.5418$ Å). Rietveld refinement of the acquired data was performed using
FULLPROF software package Carvajal (1993) taking the initial cell parameters
from Ref. [A. Mba  _et al._ , 2012]. The low-temperature XRD data down to $15$
K were recorded using a low-temperature attachment (Oxford Phenix) to the
x-ray diffractometer.

Magnetization ($M$) was measured as a function of temperature (2 K $\leq
T\leq$ 380 K) using the vibrating sample magnetometer (VSM) attachment to the
Physical Property Measurement System [PPMS, Quantum Design]. A 3He attachment
to the SQUID [MPMS-7T, Quantum Design] magnetometer was used for magnetization
measurements in the low-temperature range ($0.5$ K $\leq T\leq$ $2$ K).
Specific heat ($C_{\rm p}$) as a function of temperature was measured down to
$0.35$ K using the thermal relaxation technique in PPMS under magnetic fields
up to $14$ T. For $T\leq 2$ K, measurements were performed using an additional
3He attachment to PPMS. High-field magnetization was measured in pulsed
magnetic field up to $60$ T at the Dresden High Magnetic Field Laboratory
Tsirlin  _et al._ (2009); Skourski  _et al._ (2011).

The NMR experiments on the 31P nucleus (nuclear spin $I=1/2$ and gyromagnetic
ratio $\gamma/2\pi=17.235$ MHz/T) were carried out using pulsed NMR technique
in the temperature range 1.6 K $\leq T\leq$ 230 K. The 31P NMR spectra as a
function of temperature were obtained either by sweeping the field at a fixed
frequency or by taking the Fourier transform (FT) of the echo signal, keeping
the magnetic field fixed. The NMR shift $K(T)$ = [$H_{\rm ref}-H(T)$]/$H(T)$
was determined by measuring the resonant field $H(T)$ of the sample with
respect to the standard H3PO4 sample (resonance frequency $H_{\rm ref}$). The
31P nuclear spin-lattice relaxation rate ($1/T_{1}$) was measured using the
inversion recovery technique at different temperatures. Details of the
positions of both P sites with respect to the magnetic V4+ centers are given
in Table 1.

Density-functional (DFT) band-structure calculations were performed in the
FPLO code Koepernik and Eschrig (1999) using experimental crystal structure
from Ref. [Lavrov  _et al._ , 1982] and the Perdew-Burke-Ernzerhof (PBE)
flavor of the exchange-correlation potential Perdew  _et al._ (1996). Exchange
couplings were obtained within superexchange theory or by a mapping procedure
Tsirlin (2014) using total energies of collinear spin configurations
calculated within DFT+$U$, where the on-site Coulomb repulsion $U_{d}=5$ eV
and Hund’s coupling $J_{d}=1$ eV Nath  _et al._ (2008a); Tsirlin  _et al._
(2008) were used to account for strong correlations in the V $3d$ shell. All
calculations are performed for the $4\times 4\times 4$ $k$-mesh. To resolve
all pertinent exchange couplings, the
$(\mathbf{a}_{t}+\mathbf{c}_{t})\times\mathbf{b}_{t}\times(\mathbf{a}_{t}-\mathbf{c}_{t})$,
$\mathbf{a}_{t}\times 2\mathbf{b}_{t}\times\mathbf{c}_{t}$, and
$(\mathbf{a}_{t}+\mathbf{b}_{t})\times 2\mathbf{b}_{t}\times\mathbf{c}_{t}$
supercells were used, where $\mathbf{a}_{t}$, $\mathbf{b}_{t}$, and
$\mathbf{c}_{t}$ are lattice vectors of the triclinic unit cell given in Ref.
[Lavrov  _et al._ , 1982].

##  III Results

###  III.1 X-ray Diffraction

Figure 2: X-ray powder diffraction patterns of $\varepsilon$-LiVOPO4 at two
different temperatures ($T=300$ K and $15$ K). The solid lines denote the
Rietveld refinement fit of the data. The Bragg-peak positions are indicated by
green vertical bars and bottom blue line indicates the difference between the
experimental and calculated intensities.
Figure 3: Lattice parameters (a) $a$, (b) $b$, (c) $c$, (d) $\alpha$, (e)
$\beta$, and (f) $\gamma$ vs temperature. The unit cell volume ($V_{\rm
cell}$) is plotted as a function of temperature in right $y$-axis of (c) and
the solid line represents the fit using Eq. (1).

In order to confirm the phase purity and to study the temperature variation of
the crystal structure, the powder XRD patterns are analyzed for $15$ K $\leq$
$T$ $\leq$ $300$ K. The XRD patterns at two end temperatures, $T=300$ K and 15
K, along with the refinement are shown in Fig. 2. At room temperature, all the
peaks could be indexed based on the triclinic (space group: $P\bar{1}$)
structure, implying phase purity of the sample. The refined lattice parameters
[$a=6.729(1)$ Å, $b=7.194(1)$ Å, $c=7.920(2)$ Å, $\alpha=89.82(2)^{\circ}$,
$\beta=91.2288(2)^{\circ}$, and $\gamma=116.8799(2)^{\circ}$] at room
temperature are in good agreement with the previous report A. Mba  _et al._
(2012). Identical XRD patterns with no extra peaks are observed in the whole
measured temperature range, which excludes any structural phase transition or
lattice deformation in $\varepsilon$-LiVOPO4, unlike other spin gap compounds
NaTiSi2O6 Isobe  _et al._ (2002); Popović  _et al._ (2004), CuGeO3 Hirota  _et
al._ (1994), NaV2O5Fujii  _et al._ (1997), and K-TCNQ Lépine  _et al._ (1978).

The variation of both lattice parameters and unit cell volume ($V_{\rm cell}$)
as a function of temperature are shown in Fig. 3. The cell parameters $b$ and
$c$ decrease systematically, while the rest of the parameters ($a$, $\alpha$,
$\beta$, and $\gamma$) rise with decreasing temperature. However, the overall
unit cell volume shrinks upon cooling. The temperature variation of $V_{\rm
cell}$ was fitted by the equation Kittel (1986); Bag  _et al._ (2018)

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

where $V_{0}$ is the unit cell volume at $T$ = $0$ K, $K_{0}$ is the bulk
modulus, and $\gamma$ is the Gr$\ddot{\rm u}$neisen parameter. The internal
energy $U(T)$ can be expressed in terms of the Debye approximation as

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

In the above, $p$ stands for the total number of atoms in the unit cell and
$k_{\rm B}$ is the Boltzmann constant. The best fit of the data down to 15 K
[see Fig. 3(c)] was obtained with the Debye temperature $\theta_{D}\simeq 530$
K, $\frac{\gamma}{K_{0}}\simeq 5.78\times 10^{-5}$ Pa-1, and $V_{0}\simeq
340.5$ Å3.

###  III.2 Magnetization

Figure 4: Upper panel: $\chi(T)$ measured in $\mu_{0}H=0.5$ T. The two
shoulders of the broad maximum are indicated by the vertical arrows. Inset:
low-$T$ $\chi(T)$ measured in different applied fields. Lower panel: $1/\chi$
vs $T$ and the solid line is the CW fit using Eq. (3).

The temperature-dependent bulk magnetic susceptibility $\chi~{}(\equiv M/H$)
of $\varepsilon$-LiVOPO4 measured in an applied field of $\mu_{0}H=0.5$ T is
shown in the upper panel of Fig. 4. It is consistent with earlier reports
Chung  _et al._ (2019); Yang  _et al._ (2008); Onoda and Ikeda (2013) and
exhibits a very broad maximum with two shoulders at around $T_{\chi}^{\rm
max1}\simeq 11$ K and $T_{\chi}^{\rm max2}\simeq 27$ K. Such a broad maximum
should not be mistaken with a magnetic transition Chung  _et al._ (2019) and
mimics the short-range ordering of a low-dimensional quantum magnet. However,
already the fact that two shoulders are observed in the susceptibility
indicates the presence of two nonequivalent spin chains in
$\varepsilon$-LiVOPO4. Since peak position is related to the intrachain
exchange coupling Johnston  _et al._ (2000), we can estimate relative
strengths of the interactions in the two chains,
$\bar{J}_{2}/\bar{J}_{1}\simeq 2.45$, where
$\bar{J}_{1}=(J_{1}+J_{1}^{\prime})/2$ and
$\bar{J}_{2}=(J_{2}+J_{2}^{\prime})/2$ are average couplings in the chain 1
and chain 2, respectively. The susceptibility alone does not give information
on which of the chains features stronger couplings, but our DFT calculations
(Sec. III.5) suggest that stronger interactions occur within chain 2.

Below $T_{\chi}^{\rm max1}$, $\chi$ decreases rapidly suggesting the opening
of a spin gap. However, below $2$ K, a large upturn is observed, which can be
attributed to the effect of extrinsic paramagnetic contributions. As shown in
the inset of the upper panel of Fig. 4, with the application of magnetic field
this low-temperature Curie tail gets suppressed. Moreover, our powder XRD
suggests high purity of the sample. Therefore, this low-temperature upturn in
$\chi(T)$ may be due to the uncorrelated V4+ free spins or chain-end effects
that largely affect $\chi(T)$ at low temperatures Eggert and Affleck (1995).

To extract the magnetic parameters, we first fitted the $1/\chi(T)$ data (see
the lower panel of Fig. 4) above 150 K by the modified Curie-Weiss (CW) law,

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

where $\chi_{\rm 0}$ represents the temperature-independent susceptibility,
which includes the Van Vleck paramagnetic and core diamagnetic contributions,
$C$ is the Curie constant, and $\theta_{\rm CW}$ is the CW temperature. The
resulting fitting parameters are $\chi_{0}\simeq 7.76\times 10^{-5}$
cm3/mol-V4+, $C\simeq 0.383$ cm3K/mol-V4+, and $\theta_{\rm CW}\simeq-13.4$ K.
Negative value of $\theta_{\rm CW}$ indicates that the dominant interactions
between the V4+ spins are antiferromagnetic (AFM) in nature. The effective
moment was calculated by using the experimental value of $C$ in the relation
$\mu_{\rm eff}=\sqrt{3k_{\rm B}C/N_{\rm A}}$, where $N_{\rm A}$ is the
Avogadro’s number. The calculated value of $\mu_{\rm eff}\simeq 1.72\mu_{\rm
B}$/V4+ is very close to the theoretical spin-only value [$\mu_{\rm
eff}=g\sqrt{S(S+1)}\simeq 1.73\mu_{\rm B}$] for $S=1/2$, assuming the
Land$\acute{e}$ $g$ factor $g=2$. The experimental value of $\mu_{\rm eff}$
corresponds to a slightly reduced $g$ value of $g\simeq 1.98$ which is
identical to the case of other V4+-based compounds Mukharjee  _et al._ (2019);
Tsirlin  _et al._ (2011); Yogi  _et al._ (2015). The core diamagnetic
susceptibility ($\chi_{\rm core}$) of $\varepsilon$-LiVOPO4 was estimated to
be $-6.8\times 10^{-5}$ cm3/mol by adding the $\chi_{\rm core}$ of the
individual ions Li1+, V4+, P5+, and O2- Selwood (2013). The Van Vleck
paramagnetic susceptibility ($\chi_{\rm VV}$) of $\sim 14.6\times 10^{-5}$
cm3/mol is obtained by subtracting $\chi_{\rm core}$ from $\chi_{\rm 0}$.

Figure 5: Magnetization vs field measured at $T=1.5$ K using pulsed magnetic
field. Inset: $dM$/$dH$ vs $H$ to highlight the critical fields $H_{\rm c1}$,
$H_{\rm c2}$, and $H_{\rm c3}$.

Two alternating spin chains with different exchange parameters should also
manifest themselves in the magnetization process. Indeed, experimental
magnetization curve measured at the base temperature of $T=1.5$ K (Fig. 5)
shows a kink due to the gap closing at $\mu_{0}H_{\rm c1}\simeq 5.6$ T that
corresponds to the spin gap of $\Delta_{\rm 0}^{\rm M}/k_{\rm
B}=g\mu_{0}\mu_{\rm B}H_{\rm c1}/k_{\rm B}\simeq 7.3$ K. At $\mu_{0}H_{\rm
c2}\simeq 25$ T, half of the spins are saturated, with the
$\frac{1}{2}$-plateau observed up to $\mu_{0}H_{\rm c3}\simeq 35$ T. At even
higher fields, the remaining spins are gradually polarized with nearly 80 % of
the saturation magnetization reached at 58 T, the highest field of our
experiment.

The two-step increase of the magnetization with the spin gap and intermediate
$\frac{1}{2}$-plateau is common for dimers built by two spin-1 ions Samulon
_et al._ (2009). However, $\varepsilon$-LiVOPO4 contains spins-$\frac{1}{2}$,
so this behavior should have a different origin. We infer that at
$\mu_{0}H_{\rm c1}$ the spin gap corresponding to chain 1 is closed. This
chain is fully polarized at $\mu_{0}H_{\rm c2}$, whereas at $\mu_{0}H_{\rm
c3}$ the gap corresponding to chain 2 is closed, and the magnetization
increases further. Such a scenario is confirmed by our quantitative analysis
in Sec. III.5.

###  III.3 Specific heat

Figure 6: Upper panel: specific heat $C_{\rm p}$ vs $T$ of
$\varepsilon$-LiVOPO4 in zero applied field. The dashed line stands for the
phonon contribution to the specific heat ($C_{\rm ph}$) using Debye fit [Eq.
(4)]. The solid line is the magnetic contribution to the specific heat
($C_{\rm mag}$). Inset: $C_{\rm mag}$ vs $T$ in the low-temperature region.
The solid line is an exponential fit as described in the text. Lower panel:
$C_{\rm mag}/T$ and $S_{\rm mag}$ vs $T$ in the left and right $y$-axes,
respectively. Inset: $C_{\rm p}/T$ vs $T$ in different magnetic fields in the
low-temperature regime.

The specific-heat ($C_{\rm p}$) data measured under zero field are shown in
the upper panel of Fig. 6. No sharp anomaly/peak was noticed down to $T=0.35$
K, thus ruling out the possibility of any magnetic or structural transition. A
broad hump associated with low-dimensional short-range order is observed
around $T_{C}^{\rm max}$ $\simeq 6$ K which moves weakly toward low
temperatures with increasing magnetic field (see the inset of the lower panel
of Fig. 6), reflecting the closing of the spin gap. The gap was estimated by
fitting the data below 4 K with the activated behavior, $C_{\rm
mag}\propto\exp(-\Delta_{0}^{\rm C}/k_{\rm B}T)$. The fit, as illustrated in
the inset of the upper panel of Fig. 6, returns the zero-field spin gap of
$\Delta_{0}^{\rm C}/k_{\rm B}\simeq 7.3$ K that matches nicely the value
obtained from the high-field magnetization data.

Typically, in magnetic insulators, the high-temperature part of $C_{\rm p}$ is
dominated by the phonon contribution, whereas the magnetic contribution
becomes prominent at low temperatures. To estimate the phonon contribution
($C_{\rm ph}$), the experimental data at high temperatures (60 K $\leq T\leq$
110 K) were fitted by a linear combination of three Debye functions Nath  _et
al._ (2008b)

$C_{\rm ph}(T)=9R\displaystyle\sum\limits_{\rm n=1}^{3}c_{\rm
n}\left(\frac{T}{\theta_{\rm Dn}}\right)^{3}\int_{0}^{\frac{\theta_{\rm
Dn}}{T}}\frac{x^{4}e^{x}}{(e^{x}-1)^{2}}dx.$ (4)

In the above, $R$ is the universal gas constant, the coefficients $c_{\rm n}$
stand for the groups of different atoms present in the crystal, and
$\theta_{\rm Dn}$ are the corresponding Debye temperatures. The $C_{\rm ph}$
was extrapolated down to low temperatures and subtracted from the total
specific heat to obtain the magnetic contribution to the specific heat
($C_{\rm mag}$). The obtained $C_{\rm mag}(T)$ is presented in the upper panel
of Fig. 6. The accuracy of the above fitting procedure was further verified by
calculating the magnetic entropy $S_{\rm mag}$ obtained by integrating $C_{\rm
mag}/T$ (see the lower panel of Fig. 6). The value of $S_{\rm mag}$ is
calculated to be $\sim 5.9$ J/mol K at $T\simeq 100$ K, which is close to
$S_{\rm mag}=R\ln 2=5.76$ J/mol K expected for spin $\frac{1}{2}$.

As shown in the upper panel of Fig. 6, $C_{\rm mag}$ develops two broad maxima
at $T_{C}^{\rm max1}\simeq 7$ K and $T_{C}^{\rm max2}\simeq 27$ K, similar to
the two shoulders in the $\chi(T)$ data. The clear separation of these maxima
indicates the different interaction strength in chains 1 and 2.

Figure 7: $C_{\rm p}/T$ vs $T$ measured in different applied magnetic fields
up to 14 T.

With the spin gap closed around $\mu_{0}H_{\rm c1}\simeq 5.6$ T, the system
may enter a long-range order (LRO) state. This state is indeed observed in the
specific-heat data measured above $\mu_{0}H_{\rm c1}$. As shown in Fig. 7, no
anomaly in $C_{\rm p}$ is found down to $T=0.35$ K and up to $\mu_{0}H=5$ T.
However, for $\mu_{0}H>7$ T a clear anomaly appears indicating the onset of a
field-induced magnetic LRO ($T_{\rm N}$). This peak shifts toward higher
temperature with increasing magnetic field.

###  III.4 31P NMR

As mentioned previously, $\chi(T)$ does not show an exponential decrease at
low temperatures anticipated for a gapped spin system, and may be influenced
by extrinsic contributions. To access the intrinsic susceptibility of
$\varepsilon$-LiVOPO4, we performed NMR measurements on the 31P nuclei.

####  III.4.1 NMR Shift

Figure 8:  Typical 31P field-sweep NMR spectra of $\varepsilon$-LiVOPO4
measured at 24.96 MHz at different temperatures. The dashed lines track the
shifting of NMR line positions. The vertical line corresponds to the 31P
reference line ($H_{\rm ref}$) measured for H3PO4. The solid line is the fit
of the spectra at $T=50$ K using a double Gaussian function.

Figure 8 presents the field-sweep 31P NMR spectra measured over a wide
temperature range. At high temperatures, a narrow and symmetric spectral line
typical for a $I=1/2$ nucleus is observed. As the temperature is lowered, the
line shape becomes asymmetric, followed by a complete splitting of the two
spectral lines below about 120 K. This suggests the existence of two
nonequivalent P sites, P(1) and P(2), with a different crystallographic
environment, which is consistent with the structural data (Table 1). Both
lines shift with temperature and merge below about 4 K. The absence of drastic
line broadening and/or line splitting down to 1.6 K, except for the one caused
by the two nonequivalent P sites, rules out the occurrence of any structural
and magnetic transitions. The line with a lower intensity shifts more strongly
than the one with the higher intensity. The former can be assigned to P(2) and
the latter to P(1), because the P(2)O4 tetrahedra mediate stronger exchange
interactions, $J_{1}$ in chain 1 and $J_{2}$ in chain 2, whereas the P(1)O4
tetrahedra mediate weaker interactions $J_{1}^{\prime}$ and $J_{2}^{\prime}$,
respectively (see Sec. III.5). At $T=1.6$ K, the position of the peak is very
close to the zero shift value suggesting that the ground state of
$\varepsilon$-LiVOPO4 is nonmagnetic.

The temperature-dependent NMR shift $K(T)$ for both 31P sites was extracted by
fitting each spectrum to a sum of two Gaussian functions. The results are
shown in the upper panel of Fig. 9. One advantage of the NMR shift over the
bulk $\chi(T)$ is that the Curie-Weiss term due to foreign phases and/or
defects does not appear. The randomly distributed defects/impurities only
broaden the NMR line but do not contribute to the NMR shift Walstedt and
Walker (1974). Therefore, $K(T)$ is more favorable than bulk $\chi(T)$ data
for a reliable determination of magnetic parameters.

The $K(T)$’s corresponding to the P(1) and P(2) sites pass through a very
broad maximum, similar to the $\chi(T)$ data. The overall temperature
dependence is similar for P(1) and P(2), but the absolute values differ due to
the different hyperfine couplings. The presence of several distinct P sites in
the structure is reminiscent of the ambient-pressure polymorph of (VO)2P2O7
with its two non-equivalent alternating spin-$\frac{1}{2}$ chains. In that
case, different phosphorous sites probe the behavior of different spin chains
in the structure Yamauchi  _et al._ (1999); Kikuchi  _et al._ (1999). In
contrast, each of the P sites in $\varepsilon$-LiVOPO4 is coupled to both spin
chains, so it is not possible to probe $K(T)$ separately. Two distinct
shoulders are observed in $K(T)$ at $T_{K}^{\rm max1}\simeq 10$ K and
$T_{K}^{\rm max2}\simeq 26$ K and closely resemble bulk magnetic
susceptibility (Fig. 4).

At low temperatures, both shifts decrease rapidly toward zero, suggesting the
opening of a spin gap between the singlet ($S=0$) ground state and triplet
($S=1$) excited states. As NMR shift is insensitive to the impurities and
defects, one can use it to accurately measure the intrinsic spin
susceptibility. In the upper panel of Fig. 9, we show that for both P sites,
$K(T)$ decreases toward zero, which is in contrast to the upturn observed in
the low-temperature $\chi(T)$ data. This confirms the extrinsic nature of the
low-temperature upturn observed in $\chi(T)$. In powder samples, the defects
often break spin chains, with the unpaired spins at the ends of finite chains
giving rise to the staggered magnetization, which also appears as the low-
temperature Curie tail in $\chi(T)$.

Figure 9:  Upper panel: temperature-dependent 31P NMR shift $K(T)$ measured at
24.96 MHz as a function of temperature for both P(1) and P(2) sites. Lower
panel: $K$ vs $\chi$ (measured at 1.5 T) with temperature as an implicit
parameter. The solid lines are the linear fits.

The direct relation between $K(T)$ and spin susceptibility $\chi_{\rm spin}$
can be written as

$K(T)=K_{0}+\frac{A_{\rm hf}}{N_{A}\mu_{\rm B}}\chi_{\rm spin},$ (5)

where $K_{0}$ is the temperature-independent NMR shift and $A_{\rm hf}$ is the
total hyperfine coupling constant between the 31P nuclei and V4+ spins. The
$A_{\rm hf}$ consists of the contributions due to transferred hyperfine
coupling and nuclear dipolar coupling constants. Since both of the
aforementioned couplings are temperature-independent, $K(T)$ is a direct
measure of $\chi_{\rm spin}$. Using Eq. (5), $A_{\rm hf}$ can be calculated
from the slope of the linear $K$ vs $\chi$ plot with temperature as an
implicit parameter. The lower panel of Fig. 9 presents the $K$ vs $\chi$ plots
for both P sites showing linear behavior at high temperatures. From the linear
fit, the hyperfine coupling constants $A_{\rm hf}^{\rm P(1)}\simeq 3290$
Oe/$\mu_{\rm B}$ and $A_{\rm hf}^{\rm P(2)}\simeq 7068$ Oe/$\mu_{\rm B}$ are
estimated for the P(1) and P(2) sites, respectively. Thus, the P(2) site is
coupled with the V4+ ions twice stronger than the P(1) site. These values are
comparable to the $A_{\rm hf}$ values reported for other spin chains with
similar interaction geometry Mukharjee  _et al._ (2019); Nath  _et al._
(2008c, 2005).

Figure 10: Upper panel: $K$ vs $T$ for both P(1) and P(2) sites. The solid
lines are the fits below 5 K by $K=K_{0}+mT^{(d/2-1)}e^{-\Delta/k_{\rm B}T}$,
fixing $\Delta/k_{\rm B}\simeq 5.5$ K. Lower panel: $(K-K_{0})e^{5.5/T}$ vs
$T$ is shown in the low-temperature regime. The solid lines are the fits with
$d\simeq 2.83$ and 2.73 for the P(1) and P(2) sites, respectively.

It is also possible to estimate the value of spin gap and the effective
lattice dimensionality ($d$) by analyzing the low-temperature $K(T)$ data.
Typically, the non-negligible interchain couplings are inevitably present in
real materials and significantly influence the $K(T)$ data at low
temperatures. For a $d$-dimensional system, the susceptibility at $k_{\rm
B}T\ll\Delta$ can be approximated as Taniguchi  _et al._ (1995)

$\chi_{d}\propto T^{(d/2)-1}\times e^{-\Delta/k_{\rm B}T}.$ (6)

Our NMR experiments were carried out in the magnetic field of $\mu_{0}H=1.4$
T. Assuming a linear variation of $\Delta/k_{\rm B}$ with $H$, the zero-field
spin gap $\Delta_{\rm 0}/k_{\rm B}\simeq 7.3$ K determined from the specific
heat and magnetization is expected to be reduced to $\Delta_{\rm
1.4~{}T}/k_{\rm B}\simeq 5.5$ K at $\mu_{0}H=1.4$ T. In the upper panel of
Fig. 10, the $K(T)$ data below 5 K are fitted by
$K=K_{0}+mT^{(d/2-1)}e^{-\Delta/k_{\rm B}T}$, fixing $\Delta/k_{\rm B}\simeq
5.5$ K where $m$ is a proportionality constant. In order to highlight the low-
temperature linear regime and the fit, in the lower panel of Fig. 10, we have
plotted $(K-K_{0})e^{5.5/T}$ vs $T$ in the log-log plot. The fit in this
region returns $K_{0}\simeq 0.04694$%, $m\simeq 0.4654$ %/K1/2, and $d\simeq
2.83$ for the P(1) site and $K_{0}\simeq 0.0462$%, $m\simeq 1.1663$ %/K1/2,
and $d\simeq 2.73$ for the P(2) site. The average value of $d\simeq 2.78$,
which is very close to 3, suggests the dominant role of 3D spin-spin
correlations at low temperatures. Indeed, we find rather strong interchain
couplings from DFT (Sec. III.5).

####  III.4.2 Spin-lattice relaxation rate $1/T_{1}$

Figure 11: Upper panel: recovery of the longitudinal magnetization as a
function of waiting time $t$ at three different temperatures for the P(1)
site. Solid lines are the fits using Eq. (7). Lower panel: temperature
variation of $1/T_{1}$ measured in different magnetic fields. For
$\mu_{0}H=1.4$ T, measurements are done on both P(1) and P(2) sites while for
other fields, only P(1) site is probed. The solid line represents the fit
using $1/T_{1}\propto T^{\alpha}$ to the 10 T data at low temperatures with
$\alpha\simeq-0.6$.

The spin-lattice relaxation rate $1/T_{1}$ measures the dynamic
susceptibility, which provides direct access to the low-energy spin
excitations or spin-spin correlation function Moriya (1956). $1/T_{1}$ was
measured at the central peak position of the spectra at each temperature using
an inversion pulse sequence down to $T=1.6$ K. Since 31P has the nuclear spin
$I=1/2$, the value of $T_{1}$ at each temperature was estimated by fitting the
recovery curve of the longitudinal magnetization to a single exponential
function

$\frac{1}{2}\left[\frac{M(0)-M(t)}{M(0)}\right]=Ae^{-(t/T_{1})}.$ (7)

Here, $M(t)$ is the nuclear magnetization at a time $t$ after the inversion
pulse and $M(0)$ is the equilibrium magnetization. The upper panel of Fig. 11
shows the recovery curves at three different temperatures probed for the P(1)
site at $\mu_{0}H=1.4$ T. The recovery curves show linearity over one and half
decades when the $y$ axis is plotted in log scale. $1/T_{1}$ was estimated by
fitting Eq. (7) in this linear regime.

$1/T_{1}$ estimated from the above fit is shown in the lower panel of Fig. 11.
Our measurements are done at different field values ranging from 1.4 T to 10
T, above and below $H_{\rm c1}$. For $\mu_{0}H=1.4$ T, the measurements are
done at both P(1) and P(2) sites and over the whole temperature range, while
for other fields only the P(1) site is probed and the measurements are
restricted to low temperatures ($T<30$ K). Since there is a large difference
in the magnitude of $A_{\rm hf}$ for the P(1) and P(2) sites, they experience
different local magnetic fields induced by the V4+ spins. Therefore, it is
expected that the resulting temperature-dependent $1/T_{1}$ will have
different values accordingly. Indeed, for $\mu_{0}H=1.4$ T, $1/T_{1}$ of the
P(2) site has larger magnitude than that of the P(1) site, as $A_{\rm hf}$ of
P(2) is larger than that of P(1). For both the P-sites, $1/T_{1}$ follows a
temperature-independent behavior due to the random fluctuation of the
paramagnetic moments at high temperatures Moriya (1956). At lower
temperatures, $1/T_{1}$ starts to decrease and below about 10 K it drops
rapidly towards zero. The 1.6 K value is almost two orders of magnitude lower
than the room-temperature one, indicating the opening of a spin gap. In higher
fields, the low-temperature values of $1/T_{1}$ increase and show an upward
curvature.

The spin gap should be closed at $H_{\rm c1}$ and thereafter an AFM LRO sets
in. Therefore, we measured $1/T_{1}$ at different fields above $H_{\rm c1}$.
The increase in the low-temperature values of $1/T_{1}$ confirms the closing
of the spin gap and the growth of 3D AFM correlations due to the field-induced
LRO Klanjšek  _et al._ (2008). Since our measurements are limited down to 1.6
K, we are unable to detect the field-induced LRO from the $1/T_{1}$ data.
Nevertheless, the systematic increase of $1/T_{1}$ with field implies that
$T_{\rm N}$ shifts toward higher temperatures with increasing $H$, in good
agreement with our $C_{\rm p}(T)$ data, where field-induced LRO is detected
above $H_{\rm c1}$.

The data sets for the P(1) site at various field values exhibit a fanlike
pattern in the low-temperature regime. Similar behavior has been reported for
spin-$1/2$ ladder compound (C5H12N)2CuBr4 (BPCB) and spin-1 chain compound
NiCl2-4SC(NH2)2 (DTN) Mukhopadhyay  _et al._ (2012). It is expected that
$1/T_{1}(T)$ data for different fields, around the quantum critical point
(QCP) (i.e., around $H_{\rm c1}$) should follow a power-law behavior,
$1/T_{1}\propto T^{\alpha}$. The exponent $\alpha$ should vary across $H_{\rm
c1}$, resulting in a fanlike pattern of the $1/T_{1}$ data. For instance, in
the gapless TLL region ($H>H_{\rm c1}$), the power law has the form
$1/T_{1}\propto T^{1/(2K)-1}$, with $K$ being the TLL exponent Klanjšek  _et
al._ (2008); Giamarchi and Tsvelik (1999). The value of $K$ increases
gradually as one approaches QCP from the TLL regime and reaches the maximum
value 1 that corresponds to $\alpha[=1/(2K)-1]=-0.5$. In order to test this
formalism, we fitted the $1/T_{1}$ data below 5 K measured at 10 T by the
power law. As shown in the lower panel of Fig. 11, our fit yields
$\alpha\simeq-0.6$, which is very close to the predicted value ($\alpha=-0.5$)
in the TLL regime. This indicates the pronounced 1D character of the compound.

Figure 12: $1/T_{1}e^{5.5/T}$ vs $T$ in the log-log plot. The solid lines are
the fits, as described in the text with $d\simeq 2.74$ and 3.2, for the P(1)
and P(2) sites, respectively.

On the other hand, from analyzing the $K(T)$ data in magnetic fields below
$H_{c1}$ we already inferred that 3D spin-spin correlations caused by the
interchain coupling may be relevant. According to Mukhopadhyay et al.
Mukhopadhyay  _et al._ (2012), in the gapped region ($H\leq H_{c1}$) with 3D
correlations $1/T_{1}$ follows an activated behavior of the form,

$1/T_{1}\propto T^{\alpha_{\rm 0}}\exp\left[\frac{g\mu_{B}(H-H_{\rm
c1})}{k_{\rm B}T}\right].$ (8)

The exponent $\alpha_{0}$ in the above equation depends on the effective
dimension of the magnon dispersion relation as set by the thermal fluctuations
$k_{\rm B}T$. With increasing temperature, $\alpha_{0}$ slowly decreases from
2 for $k_{\rm B}T<J_{\rm 3D}$ (3D regime) to 0 for $J_{3D}<k_{\rm B}T<J_{1D}$
(1D regime). In order to detect the effective dimensionality of the spin-spin
correlations, Eq. (8) was further simplified as $1/T_{1}\propto
T^{d-1}e^{-\Delta_{\rm 1.4\,T}/k_{\rm B}T}$ by putting $\alpha_{\rm 0}=d-1$
and $\Delta_{\rm 1.4~{}T}/k_{\rm B}=g\mu_{\rm B}(H_{\rm c1}-H)/k_{\rm B}$ and
fitted to the low temperature $1/T_{1}$ data. The fit for $T\leq 5$ K, taking
$\Delta_{\rm 1.4\,T}/k_{\rm B}\simeq 5.5$ K results in $d\simeq 2.74$ and 3.2
for the P(1) and P(2) sites, respectively. The average value of $d\simeq 2.97$
is consistent with the value obtained from the $K(T)$ analysis. This further
confirms the importance of interchain couplings at low temperatures, where
activated behavior is observed. Figure 12 presents the $1/T_{1}e^{5.5/T}$ vs
$T$ plot for the data measured at $\mu_{0}H=1.4$ T along with the fit using
Eq. (8). Both the $x$ and $y$ axes are shown in log scale to highlight the
power-law prefactor to the activated behavior, at low temperatures.

###  III.5 Microscopic magnetic model

Similar to Refs. Tsirlin  _et al._ (2011); Arjun  _et al._ (2019); Mukharjee
_et al._ (2019), we use two complementary computational methods to derive
exchange couplings in $\varepsilon$-LiVOPO4. For a single magnetic orbital of
V4+, superexchange theory yields antiferromagnetic exchange couplings
$J_{i}^{\rm AFM}=4t_{i}^{2}/U_{\rm eff}$, where $t_{i}$ are V–V hoppings
extracted from the uncorrelated (PBE) band structure, and $U_{\rm eff}$ is an
effective Coulomb repulsion in the V $3d$ bands. On the other hand, exchange
couplings $J_{i}$ can be obtained from DFT+$U$ by a mapping procedure, where
both ferromagnetic and antiferromagnetic contributions are taken into account.

Table 2:  Exchange couplings in $\varepsilon$-LiVOPO4. The $t_{i}$ values are the V–V hoppings extracted from the uncorrelated band structure, and show relative strengths of the AFM contributions to the exchange couplings $J_{i}^{\rm AFM}\sim t_{i}^{2}$. The $J_{i}$ are exchange interactions obtained by the mapping procedure within DFT+$U$.  |  $d_{\rm V-V}$ (Å) |  |  $t_{i}$ (meV) |  $J_{i}$ (K)
---|---|---|---|---
$J_{1}$ | 5.250 | V1–V1 | $-72$ | 33
$J_{1}^{\prime}$ | 5.101 | V1–V1 | $-55$ | 23
$J_{2}$ | 5.275 | V2–V2 | $-117$ | 63
$J_{2}^{\prime}$ | 5.303 | V2–V2 | $-78$ | 22
$J_{c1}$ | 3.599 | V1–V2 | 0 | $-15$
$J_{c2}$ | 3.629 | V1–V2 | 0 | $-15$
$J_{a1}$ | 6.018 | V1–V2 | $-21$ | 12
$J_{a2}$ | 6.070 | V1–V2 | $-32$ | 7

In Table 2, we list the $t_{i}$ values for the uncorrelated band structure and
the exchange couplings $J_{i}$ obtained from DFT+$U$. The two methods are in
excellent agreement and consistently show the stronger couplings within chain
2. Moreover, we find that within each spin chain the stronger couplings
involve the P(2) bridges and the weaker couplings involve the P(1) bridges. On
the structural level, this difference should be traced back to the lateral
displacements $d_{i}$ of the VO6 octahedra within the spin chain [Fig. 1(b)],
where smaller displacement leads to a stronger coupling Mukharjee  _et al._
(2019); Roca  _et al._ (1998). Indeed, we find $d_{1}=0.71$ Å for $J_{1}$ vs
$d_{1}^{\prime}=0.97$ Å for $J_{1}^{\prime}$ and $d_{2}=0.08$ Å for $J_{2}$ vs
$d_{2}^{\prime}=0.23$ Å for $J_{2}^{\prime}$. The smaller lateral
displacements $d_{2}$ and $d_{2}^{\prime}$ could also explain the overall
stronger couplings in chain 2, although in this case other geometrical
parameters Roca  _et al._ (1998) are relevant as well, because $J_{1}$ is
about as strong as $J_{2}^{\prime}$, despite the fact that
$d_{1}>d_{2}^{\prime}$.

Regarding the interchain couplings, the microscopic scenario is very similar
to that of monoclinic $A$VO$X$O4 with $A$ = Ag, Na and $X$ = P, As Tsirlin
_et al._ (2011); Arjun  _et al._ (2019); Mukharjee  _et al._ (2019). Shorter
V–O–V bridges render $J_{c1}$ and $J_{c2}$ ferromagnetic, whereas the long-
range couplings $J_{a1}$ and $J_{a2}$ are weakly antiferromagnetic. These
ferromagnetic and antiferromagnetic interactions compete and make the spin
lattice of $\varepsilon$-LiVOPO4 frustrated.

Table 3:  Exchange couplings (in K) extracted from the $\chi(T)$ and $M(H)$ data using the fits shown in Fig. 13. The susceptibility fit using Eq. (9) returns $\chi_{0}=2.7\times 10^{-5}$ cm3/mol, $C_{\rm imp}=0.013$ cm3K/mol (3.5% of paramagnetic impurities), $\theta_{\rm imp}=0.9$ K, and $g=2.03$. This $g$ value is slightly higher than 1.98 obtained from the Curie-Weiss fit, probably because the interchain couplings were not taken into account. For magnetization data, $g=1.98$ has been used as a fixed parameter.  | $J_{1}$ | $J_{1}^{\prime}$ | $J_{2}$ | $J_{2}^{\prime}$
---|---|---|---|---
$\chi(T)$ | 20 | 12 | 70 | 20
$M(H)$ | 19 | 12 | 63 | 22

Our DFT results suggest that chain 1 shows only a moderate degree of
alternation ($\alpha_{1}=J_{1}^{\prime}/J_{1}\simeq 0.6$) that, together with
the lower energy scale of the couplings, leads to a relatively small spin gap
closed at $H_{\rm c1}$. In contrast, the alternation ratio of
$\alpha_{2}=J_{2}^{\prime}/J_{2}\simeq 0.3$ renders chain 2 strongly dimerized
with a larger spin gap that is closed at the much higher field $H_{\rm c3}$.
The model of two alternating spin-$\frac{1}{2}$ chains was further used to
calculate temperature-dependent magnetic susceptibility and field-dependent
magnetization of $\varepsilon$-LiVOPO4 (Fig. 13). For the susceptibility, we
used the analytical expression from Ref. Johnston  _et al._ (2000) augmented
by additional terms that account for the temperature-independent ($\chi_{0}$)
and Curie-type impurity contributions,

$\chi=\chi_{0}+\frac{C_{\rm imp}}{T+\theta_{\rm imp}}+\chi_{\rm ch1}+\chi_{\rm
ch2},$ (9)

where $\chi_{\rm ch1}$ and $\chi_{\rm ch2}$ are susceptibilities of two
nonequivalent alternating spin-$\frac{1}{2}$ chains with the interaction
parameters $J_{1}$, $J_{1}^{\prime}$ and $J_{2}$, $J_{2}^{\prime}$,
respectively (same $g$ factor is used for both chains). Magnetization curve of
$\varepsilon$-LiVOPO4 was modeled by the sum of simulated magnetization curves
for the V(1) and V(2) sublattices obtained by the method described in Ref.
Tsirlin  _et al._ (2011).

Figure 13:  Temperature-dependent susceptibility (top) and field-dependent
magnetization (bottom) of $\varepsilon$-LiVOPO4 fitted using the model of two
non-equivalent alternating spin-$\frac{1}{2}$ chains, as explained in the
text. See Table 3 for the fitting parameters.

The fitting results listed in Table 3 show good agreement between the fits to
the susceptibility and magnetization. They are also consistent with the
exchange parameters calculated by DFT (Table 2). We chose not to include
interchain couplings into the susceptibility fit, which became ambiguous when
more than four exchange parameters were involved. The effect of the interchain
couplings can be seen from the fact that for an isolated chain 1 one finds,
using $J_{1}$ and $J_{1}^{\prime}$ from the susceptibility fit, the zero-field
spin gap of 11.3 K, which is somewhat larger than 7.3 K obtained
experimentally 222Similar to Ref. Tsirlin  _et al._ (2011), magnetization
curve was simulated by including a weak interchain coupling of
$J_{\perp}/J_{1}=-0.05$ in order to reproduce the experimental data around
$H_{\rm c1}$. This interchain coupling is much lower than estimated by DFT
(Table 2), possibly because $J_{\perp}$ reflects a cumulative effect of the
frustrated couplings $J_{a1}$, $J_{a2}$ and $J_{c1}$, $J_{c2}$. .

##  IV Discussion and summary

$\varepsilon$-LiVOPO4 features two types of alternating spin-$\frac{1}{2}$
chains that manifest themselves in the double maxima of the susceptibility and
magnetic specific heat and in the two-step magnetization process. This unusual
microscopic scenario is reminiscent of the ambient-pressure polymorph of
(VO)2P2O7 Johnston  _et al._ (2001), where two spin gaps corresponding to two
types of spin chains were directly observed by NMR Yamauchi  _et al._ (1999)
and inelastic neutron scattering Garrett  _et al._ (1997). On the other hand,
large size of these gaps (35 K and 68 K, respectively) and high critical
fields associated with them preclude experimental access to field-induced
transitions, where triplon excitations of the spin chains consecutively
condense leading to a long-range magnetic order.

Figure 14: $H-T$ phase diagram of $\varepsilon$-LiVOPO4 obtained using the
data points from $C_{\rm p}(T)$ measurements. The solid line corresponds to
Eq. (10). Inset: $T_{\rm N}$ vs ($H-H_{\rm c1})^{1/1.5}$ to highlight the low-
temperature linear regime and the agreement of the simulated curve with the
experimental data points.

With its lower critical fields, $\varepsilon$-LiVOPO4 offers a much better
platform for studying these transitions experimentally. Indeed, we observed
field-induced magnetic order already above $\mu_{0}H_{\rm c1}\simeq 5.6$ T.
Transition temperature systematically increases with field and tracks the
$H$-$T$ phase boundary shown in Fig. 14.

The field-induced transition in gapped quantum magnets is often understood as
triplon BEC. In this case, the phase boundary close to $H_{\rm c1}$ should
follow the universal power law Zapf  _et al._ (2014); Giamarchi and Tsvelik
(1999); Nohadani  _et al._ (2004),

$T_{\rm N}\propto(H-H_{\rm c1})^{\frac{1}{\phi}},$ (10)

where $\phi=d/2$ is the critical exponent reflecting the universality class of
the quantum phase transition at $H_{\rm c1}$, and $d$ is the lattice
dimensionality. In the absence of the low-temperature data immediately above
$H_{\rm c1}$, we simulate the anticipated critical behavior by choosing $d=3$,
$\phi=\frac{3}{2}$, and $\mu_{0}H_{\rm c1}\simeq 5.6$ T in Eq. (10). The
resulting curves match our data up to 1.2 K (Fig. 14), although this agreement
may be partly accidental, because the temperature of 1.2 K is rather high
compared to the temperature scale of the BEC transition given by
$T_{N}^{\max}$ at the tip of the BEC dome. Since 14 T is roughly the midpoint
between $H_{\rm c1}$ and $H_{\rm c2}$, we expect $T_{N}^{\max}\simeq 1.7$ K,
and 1.2 K is more than half of this temperature. Further measurements of the
phase boundary around $H_{\rm c1}$ and below 0.7 K would be clearly
interesting to confirm the tentative BEC critical behavior in
$\varepsilon$-LiVOPO4. Impurity effects witnessed by the low-temperature
upturn in the magnetic susceptibility (Fig. 4) may also be important.

The fate of the ordered phase in fields above 14 T is of significant interest
too. The two-step increase in the magnetization suggests that the transition
at $H_{\rm c1}$ corresponds to chain 1 and will lead to a BEC dome between
$H_{\rm c1}$ and $H_{\rm c2}$, while chain 2 remains in the singlet state up
to $H_{\rm c3}$, where another BEC dome should appear. Phenomenologically,
this behavior would be similar to the two-dome $H$-$T$ phase diagrams of
spin-1 dimer magnets Samulon  _et al._ (2009), although in
$\varepsilon$-LiVOPO4 with local spin $\frac{1}{2}$ it should have a very
different microscopic origin. It does in fact arise from the coexistence of
the nonequivalent magnetic sublattices. We also point out that the
$\frac{1}{2}$-magnetization plateau in this compound is not caused by Wigner
crystallization and is thus dissimilar to the magnetization plateaus in
SrCu2(BO${}_{3})_{2}$. Because $H_{\rm c3}$ lies above $H_{\rm c2}$, magnetic
orders related to chain 1 and chain 2 should be mostly decoupled. This is
different from other systems with multiple spin gaps, where intertwined BEC
transitions lead to an unusual critical behavior, as in the dimer magnet
BaCuSi2O6 Mazurenko  _et al._ (2014); Allenspach  _et al._ (2020).

In summary, we have shown that $\varepsilon$-LiVOPO4 is a gapped quantum
magnet that features singlet ground state in zero field. With two non-
equivalent alternating spin-$\frac{1}{2}$ chains, it shows double maxima in
the susceptibility and magnetic specific heat and a two-step increase in the
magnetization. Chain 1 features weaker couplings and a weaker alternation
($J_{1}\simeq 20$ K, $\alpha_{1}\simeq 0.6$), whereas chain 2 reveals stronger
couplings and lies closer to the dimer limit ($J_{2}\simeq 60$ K,
$\alpha_{2}\simeq 0.3$). The zero-field spin gap of $\Delta_{0}/k_{B}\simeq
7.3$ K is closed at $\mu_{0}H_{\rm c1}\simeq 5.6$ T. The magnetization
increases up to $\mu_{0}H_{\rm c2}\simeq 25$ T, flattens out within the
$\frac{1}{2}$ plateau, and increases again above $\mu_{0}H_{\rm c3}\simeq 35$
T. The gap closing above $H_{\rm c1}$ leads to a field-induced LRO that can be
understood as Bose-Einstein condensation of triplons.

###### Acknowledgements.

P.K.M. and R.N. acknowledge BRNS, India for financial support bearing sanction
No.37(3)/14/26/2017-BRNS. We also thank C. Klausnitzer (MPI-CPfS) for the
technical support. A.A.T was funded by the Federal Ministry for Education and
Research through the Sofja Kovalevskaya Award of Alexander von Humboldt
Foundation. We also acknowledge the support of the HLD at HZDR, member of
European Magnetic Field Laboratory (EMFL).

## References

  * Matsubara and Matsuda (1956) T. Matsubara and H. Matsuda, “A Lattice Model of Liquid Helium, I,” Prog. Theor. Phys. 16, 569 (1956).
  * Giamarchi  _et al._ (2008) T. Giamarchi, Ch. Rüegg, and O. Tchernyshyov, “Bose-Einstein condensation in magnetic insulators,” Nat. Phys. 4, 198 (2008).
  * Zapf  _et al._ (2014) V. Zapf, M. Jaime, and C. D. Batista, “Bose-Einstein condensation in quantum magnets,” Rev. Mod. Phys. 86, 563 (2014).
  * Mukhopadhyay  _et al._ (2012) S. Mukhopadhyay, M. Klanjšek, M. S. Grbić, R. Blinder, H. Mayaffre, C. Berthier, M. Horvatić, M. A. Continentino, A. Paduan-Filho, B. Chiari, and O. Piovesana, “Quantum-critical spin dynamics in quasi-one-dimensional antiferromagnets,” Phys. Rev. Lett. 109, 177206 (2012).
  * Sebastian  _et al._ (2006) S. E. Sebastian, N. Harrison, C. D. Batista, L. Balicas, M. Jaime, P. A. Sharma, N. Kawashima, and I. R. Fisher, “Dimensional reduction at a quantum critical point,” Nature 441, 617 (2006).
  * Nikuni  _et al._ (2000) T. Nikuni, M. Oshikawa, A. Oosawa, and H. Tanaka, “Bose-Einstein condensation of dilute magnons in TlCuCl3,” Phys. Rev. Lett. 84, 5868 (2000).
  * Klanjšek  _et al._ (2008) M. Klanjšek, H. Mayaffre, C. Berthier, M. Horvatić, B. Chiari, O. Piovesana, P. Bouillot, C. Kollath, E. Orignac, R. Citro, and T. Giamarchi, “Controlling luttinger liquid physics in spin ladders under a magnetic field,” Phys. Rev. Lett. 101, 137207 (2008).
  * Matsushita  _et al._ (2017) T. Matsushita, N. Hori, S. Takata, N. Wada, N. Amaya, and Y. Hosokoshi, “Direct three-dimensional ordering of quasi-one-dimensional quantum dimer system near critical fields,” Phys. Rev. B 95, 020408(R) (2017).
  * Willenberg  _et al._ (2015) B. Willenberg, H. Ryll, K. Kiefer, D. A. Tennant, F. Groitl, K. Rolfs, P. Manuel, D. Khalyavin, K. C. Rule, A. U. B. Wolter, and S. Süllow, “Luttinger liquid behavior in the alternating spin-chain system copper nitrate,” Phys. Rev. B 91, 060407(R) (2015).
  * Thielemann  _et al._ (2009) B. Thielemann, Ch. Rüegg, K. Kiefer, H. M. Rønnow, B. Normand, P. Bouillot, C. Kollath, E. Orignac, R. Citro, T. Giamarchi, A. M. Läuchli, D. Biner, K. W. Krämer, F. Wolff-Fabris, V. S. Zapf, M. Jaime, J. Stahn, N. B. Christensen, B. Grenier, D. F. McMorrow, and J. Mesot, “Field-controlled magnetic order in the quantum spin-ladder system (Hpip)2CuBr4,” Phys. Rev. B 79, 020408(R) (2009).
  * Rice (2002) T. M. Rice, “To condense or not to condense,” Science 298, 760 (2002).
  * Narumi  _et al._ (1998) Y. Narumi, M. Hagiwara, R. Sato, K. Kindo, H. Nakano, and M. Takahashi, “High field magnetization in a $S=1$ antiferromagnetic chain with bond alternation,” Physica B: Cond. Mat. 246, 509 (1998).
  * Shiramura  _et al._ (1998) B. Shiramura, K. Takatsu, B. Kurniawan, H. Tanaka, H. Uekusa, Y. Ohashi, K. Takizawa, H. Mitamura, and T. Goto, “Magnetization plateaus in NH4CuCl3,” J. Phys. Soc. Jpn. 67, 1548 (1998).
  * Kodama  _et al._ (2002) K. Kodama, M. Takigawa, M. Horvatić, C. Berthier, H. Kageyama, Y. Ueda, S. Miyahara, F. Becca, and F. Mila, “Magnetic superstructure in the two-dimensional quantum antiferromagnet SrCu2(BO3)2,” Science 298, 395 (2002).
  * Kageyama  _et al._ (1999) H. Kageyama, K. Yoshimura, R. Stern, N. V. Mushnikov, K. Onizuka, M. Kato, K. Kosuge, C. P. Slichter, T. Goto, and Y. Ueda, “Exact dimer ground state and quantized magnetization plateaus in the two-dimensional spin system SrCu2(BO3)2,” Phys. Rev. Lett. 82, 3168 (1999).
  * Rüegg  _et al._ (2003) C. Rüegg, N. Cavadini, A. Furrer, H. U. Güdel, K. Krämer, H. Mutka, A. Wildes, K. Habicht, and P. Vorderwisch, “Bose-Einstein condensation of the triplet states in the magnetic insulator TlCuCl3,” Nature 423, 62 (2003).
  * Jaime  _et al._ (2004) M. Jaime, V. F. Correa, N. Harrison, C. D. Batista, N. Kawashima, Y. Kazuma, G. A. Jorge, R. Stern, I. Heinmaa, S. A. Zvyagin, Y. Sasago, and K. Uchinokura, “Magnetic-field-induced condensation of triplons in han purple pigment BaCuSi2O6,” Phys. Rev. Lett. 93, 087203 (2004).
  * Aczel  _et al._ (2009a) A. A. Aczel, Y. Kohama, C. Marcenat, F. Weickert, M. Jaime, O. E. Ayala-Valenzuela, R. D. McDonald, S. D. Selesnic, H. A. Dabkowska, and G. M. Luke, “Field-induced Bose-Einstein condensation of triplons up to 8 K in Sr3Cr2O8,” Phys. Rev. Lett. 103, 207203 (2009a).
  * Aczel  _et al._ (2009b) A. A. Aczel, Y. Kohama, M. Jaime, K. Ninios, H. B. Chan, L. Balicas, H. A. Dabkowska, and G. M. Luke, “Bose-Einstein condensation of triplons in Ba3Cr2O8,” Phys. Rev. B 79, 100409(R) (2009b).
  * Jeong  _et al._ (2013) M. Jeong, H. Mayaffre, C. Berthier, D. Schmidiger, A. Zheludev, and M. Horvatić, “Attractive Tomonaga-Luttinger liquid in a quantum spin ladder,” Phys. Rev. Lett. 111, 106404 (2013).
  * Jeong  _et al._ (2017) M. Jeong, H. Mayaffre, C. Berthier, D. Schmidiger, A. Zheludev, and M. Horvatić, “Magnetic-order crossover in coupled spin ladders,” Phys. Rev. Lett. 118, 167206 (2017).
  * Möller  _et al._ (2017) J. S. Möller, T. Lancaster, S. J. Blundell, F. L. Pratt, P. J. Baker, F. Xiao, R. C. Williams, W. Hayes, M. M. Turnbull, and C. P. Landee, “Quantum-critical spin dynamics in a Tomonaga-Luttinger liquid studied with muon-spin relaxation,” Phys. Rev. B 95, 020402(R) (2017).
  * Rüegg  _et al._ (2008) Ch. Rüegg, K. Kiefer, B. Thielemann, D. F. McMorrow, V. Zapf, B. Normand, M. B. Zvonarev, P. Bouillot, C. Kollath, T. Giamarchi, S. Capponi, D. Poilblanc, D. Biner, and K. W. Krämer, “Thermodynamics of the spin Luttinger liquid in a model ladder material,” Phys. Rev. Lett. 101, 247202 (2008).
  * Ueda (1998) Y. Ueda, “Vanadate family as spin-gap systems,” Chem. Mater. 10, 2653 (1998).
  * Yamauchi  _et al._ (1999) T. Yamauchi, Y. Narumi, J. Kikuchi, Y. Ueda, K. Tatani, T. C. Kobayashi, K. Kindo, and K. Motoya, “Two gaps in (VO)2P2O7: Observation using high-field magnetization and NMR,” Phys. Rev. Lett. 83, 3729 (1999).
  * Johnston  _et al._ (1987) D. C. Johnston, J. W. Johnson, D. P. Goshorn, and A. J. Jacobson, “Magnetic susceptibility of (VO)2P2O7: A one-dimensional spin-1/2 Heisenberg antiferromagnet with a ladder spin configuration and a singlet ground state,” Phys. Rev. B 35, 219 (1987).
  * Ghoshray  _et al._ (2005) K. Ghoshray, B. Pahari, B. Bandyopadhyay, R. Sarkar, and A. Ghoshray, “${}^{51}\rm{V}$ NMR study of the quasi-one-dimensional alternating chain compound ${\mathrm{BaCu}}_{2}{\mathrm{V}}_{2}{\mathrm{O}}_{8}$,” Phys. Rev. B 71, 214401 (2005).
  * Mukharjee  _et al._ (2019) P. K. Mukharjee, K. M. Ranjith, B. Koo, J. Sichelschmidt, M. Baenitz, Y. Skourski, Y. Inagaki, Y. Furukawa, A. A. Tsirlin, and R. Nath, “Bose-Einstein condensation of triplons close to the quantum critical point in the quasi-one-dimensional spin-$\frac{1}{2}$ antiferromagnet NaVOPO4,” Phys. Rev. B 100, 144433 (2019).
  * Arjun  _et al._ (2019) U. Arjun, K. M. Ranjith, B. Koo, J. Sichelschmidt, Y. Skourski, M. Baenitz, A. A. Tsirlin, and R. Nath, “Singlet ground state in the alternating spin-$\frac{1}{2}$ chain compound NaVOAsO4,” Phys. Rev. B 99, 014421 (2019).
  * Ahmed  _et al._ (2017) N. Ahmed, P. Khuntia, K. M. Ranjith, H. Rosner, M. Baenitz, A. A. Tsirlin, and R. Nath, “Alternating spin chain compound AgVOAsO4 probed by ${}^{75}\mathrm{As}$ NMR,” Phys. Rev. B 96, 224423 (2017).
  * Tsirlin  _et al._ (2011) A. A. Tsirlin, R. Nath, J. Sichelschmidt, Y. Skourski, C. Geibel, and H. Rosner, “Frustrated couplings between alternating spin-$\frac{1}{2}$ chains in AgVOAsO4,” Phys. Rev. B 83, 144412 (2011).
  * Weickert  _et al._ (2019) F. Weickert, A. A. Aczel, M. B. Stone, V. O. Garlea, C. Dong, Y. Kohama, R. Movshovich, A. Demuer, N. Harrison, M. B. Gamża, A. Steppke, M. Brando, H. Rosner, and A. A. Tsirlin, “Field-induced double dome and Bose-Einstein condensation in the crossing quantum spin chain system AgVOAsO4,” Phys. Rev. B 100, 104422 (2019).
  * Harrison and Manthiram (2013) K.L. Harrison and A. Manthiram, “Microwave-assisted solvothermal synthesis and characterization of various polymorphs of LiVOPO4,” Chem. Mater. 25, 1751 (2013).
  * Hidalgo  _et al._ (2019) M.F.V. Hidalgo, Y.-C. Lin, A. Grenier, D. Xiao, J. Rana, R. Tran, H. Xin, M. Zuba, J. Donohue, F.O. Omenya, I.-H. Chu, Z. Wang, X.G. Li, N.A. Chernova, K.W. Chapman, G. Zhou, L. Piper, S.P. Ong, and M.S. Whittingham, “Rational synthesis and electrochemical performance of LiVOPO4 polymorphs,” J. Mater. Chem. A 7, 8423 (2019).
  * Note (1) Triclinic crystal structure is derived from $\varepsilon$-VOPO4. Alternatively, triclinic LiVOPO4 is sometimes referred to as the $\alpha$-phase, because it was the earliest discovered LiVOPO4 polymorph.
  * Yang  _et al._ (2008) Y. Yang, H. Fang, J. Zheng, L. Li, G. Li, and G. Yan, “Towards the understanding of poor electrochemical activity of triclinic LiVOPO4: Experimental characterization and theoretical investigations,” Solid State Sci. 10, 1292 (2008).
  * Quackenbush  _et al._ (2015) N. F. Quackenbush, L. Wangoh, D. O. Scanlon, R. Zhang, Y. Chung, Z. Chen, B. Wen, Y. Lin, J. C. Woicik, N. A. Chernova, S. P. Ong, M.S. Whittingham, and L. F. J. Piper, “Interfacial effects in $\varepsilon$-LixVOPO4 and evolution of the electronic structure,” Chem. Mater. 27, 8211 (2015).
  * Lin  _et al._ (2016) Y.-C. Lin, B. Wen, K. M. Wiaderek, S. Sallis, H. Liu, S. H. Lapidus, O. J. Borkiewicz, N. F. Quackenbush, N. A. Chernova, K. Karki, F. Omenya, P. J. Chupas, L. F. J. Piper, M. S. Whittingham, K. W. Chapman, and S. P. Ong, “Thermodynamics, kinetics and structural evolution of $\epsilon$-LiVOPO4 over multiple lithium intercalation,” Chem. Mater. 28, 1794 (2016).
  * Shi  _et al._ (2018) Y. Shi, H. Zhou, I.D. Seymour, S. Britto, J. Rana, L.W. Wangoh, Y. Huang, Q. Yin, P.J. Reeves, M. Zuba, Y. Chung, F. Omenya, N.A. Chernova, G. Zhou, L.F.J. Piper, C.P. Grey, and M.S. Whittingham, “Electrochemical performance of nanosized disordered LiVOPO4,” ACS Omega 3, 7310 (2018).
  * Chung  _et al._ (2019) Y. Chung, E. Cassidy, K. Lee, C. Siu, Y. Huang, F. Omenya, J. Rana, K.M. Wiaderek, N.A. Chernova, K.W. Chapman, L.F.J. Piper, and M.S. Whittingham, “Nonstoichiometry and defects in hydrothermally synthesized $\varepsilon$-LiVOPO4,” ACS Appl. Energy Mater. 2, 4792 (2019).
  * Onoda and Ikeda (2013) M. Onoda and S. Ikeda, “Crystal structure and spin-singlet state of the LixVOPO4 insertion electrode system with alternating-bond chain,” J. Phys. Soc. Jpn. 82, 053801 (2013).
  * Carvajal (1993) J. R. Carvajal, “Recent advances in magnetic structure determination by neutron powder diffraction,” Physica B: Cond. Mat. 192, 55 (1993).
  * A. Mba  _et al._ (2012) J. A. Mba, C. Masquelier, E. Suard, and L. Croguennec, “Synthesis and crystallographic study of homeotypic LiVPO4F and LiVOPO4,” Chem. Mater. 24, 1223 (2012).
  * 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).
  * 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).
  * 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).
  * Lavrov  _et al._ (1982) A.V. Lavrov, V.P. Nikolaev, G.G. Sadikov, and M.A. Poray-Koshits, “Synthesis and crystal structure of mixed vanadyl and lithium orthophosphate LiVOPO4,” Dokl. Akad. Nauk SSSR 266, 343 (1982).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Tsirlin (2014) A.A. Tsirlin, “Spin-chain magnetism and uniform Dzyaloshinsky-Moriya anisotropy in BaV3O8,” Phys. Rev. B 89, 014405 (2014).
  * Nath  _et al._ (2008a) R. Nath, A.A. Tsirlin, E.E. Kaul, M. Baenitz, N. Büttgen, C. Geibel, and H. Rosner, “Strong frustration due to competing ferromagnetic and antiferromagnetic interactions: Magnetic properties of M(VO)2(PO${}_{4})_{2}$ (M = Ca and Sr),” Phys. Rev. B 78, 024418 (2008a).
  * Tsirlin  _et al._ (2008) A.A. Tsirlin, R. Nath, C. Geibel, and H. Rosner, “Magnetic properties of Ag2VOP2O7: An unexpected spin dimer system,” Phys. Rev. B 77, 104436 (2008).
  * Isobe  _et al._ (2002) M. Isobe, E. Ninomiya, A. N. Vasil’ev, and U. Yutaka, “Novel phase transition in spin-$\frac{1}{2}$ linear chain systems: NaTiSi2O6 and LiTiSi2O6,” J. Phys. Soc. Jpn. 71, 1423 (2002).
  * Popović  _et al._ (2004) Z. S. Popović, Z. V. Šljivančanin, and F. R. Vukajlović, “Sodium pyroxene NaTiSi2O6: Possible Haldane spin-1 chain system,” Phys. Rev. Lett. 93, 036401 (2004).
  * Hirota  _et al._ (1994) K. Hirota, D. E. Cox, J. E. Lorenzo, G. Shirane, J. M. Tranquada, M. Hase, K. Uchinokura, H. Kojima, Y. Shibuya, and I. Tanaka, “Dimerization of CuGeO3 in the spin-Peierls state,” Phys. Rev. Lett. 73, 736 (1994).
  * Fujii  _et al._ (1997) Y. Fujii, H. Nakao, T. Yosihama, M. Nishi, K. Nakajima, K. Kakurai, M. Isobe, Y. Ueda, and H. Sawa, “New inorganic spin-Peierls compound NaV2O5 evidenced by x-ray and neutron scattering,” J. Phys. Soc. Jpn. 66, 326 (1997).
  * Lépine  _et al._ (1978) Y. Lépine, A. Caillé, and V. Larochelle, “Potassium-tetracyanoquinodimethane (K-TCNQ): A spin-Peierls system,” Phys. Rev. B 18, 3585 (1978).
  * Kittel (1986) Charles Kittel,  _Introduction to Solid State Physics_ , 6th ed. (John Wiley & Sons, Inc., New York, 1986).
  * Bag  _et al._ (2018) P. Bag, P. R. Baral, and R. Nath, “Cluster spin-glass behavior and memory effect in Cr0.5Fe0.5Ga,” Phys. Rev. B 98, 144436 (2018).
  * Johnston  _et al._ (2000) D.C. Johnston, R. K. Kremer, M. Troyer, X. Wang, A. Klümper, S. L. Budko, A. F. Panchula, and P. C. Canfield, “Thermodynamics of spin $S$= $\frac{1}{2}$ antiferromagnetic uniform and alternating-exchange Heisenberg chains,” Phys. Rev. B 61, 9558 (2000).
  * Eggert and Affleck (1995) S. Eggert and I. Affleck, “Impurities in $S=\frac{1}{2}$ Heisenberg antiferromagnetic chains: Consequences for neutron scattering and knight shift,” Phys. Rev. Lett. 75, 934 (1995).
  * 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).
  * Selwood (2013) P. W. Selwood,  _Magnetochemistry_ (Read Books Ltd, 2013).
  * Samulon  _et al._ (2009) E.C. Samulon, Y. Kohama, R.D. McDonald, M.C. Shapiro, K.A. Al-Hassanieh, C.D. Batista, M. Jaime, and I.R. Fisher, “Asymmetric quintuplet condensation in the frustrated $S=1$ spin-dimer compound Ba3Mn2O8,” Phys. Rev. Lett. 103, 047202 (2009).
  * Nath  _et al._ (2008b) 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 (2008b).
  * Walstedt and Walker (1974) R. E. Walstedt and L. R. Walker, “Nuclear-resonance line shapes due to magnetic impurities in metals,” Phys. Rev. B 9, 4857 (1974).
  * Kikuchi  _et al._ (1999) J. Kikuchi, K. Motoya, T. Yamauchi, and Y. Ueda, “Coexistence of double alternating antiferromagnetic chains in (VO)2P2O7:NMR study,” Phys. Rev. B 60, 6731 (1999).
  * Nath  _et al._ (2008c) R. Nath, D. Kasinathan, H. Rosner, M. Baenitz, and C. Geibel, “Electronic and magnetic properties of K2CuP2O7: A model $S=\frac{1}{2}$ Heisenberg chain system,” Phys. Rev. B 77, 134451 (2008c).
  * Nath  _et al._ (2005) R. Nath, A. V. Mahajan, N. Büttgen, C. Kegler, A. Loidl, and J. Bobroff, “Study of one-dimensional nature of $S=\frac{1}{2}$ of (Sr,Ba)2Cu(PO4)2 and BaCuP2O7 via 31P NMR,” Phys. Rev. B 71, 174436 (2005).
  * Taniguchi  _et al._ (1995) S. Taniguchi, T. Nishikawa, Y. Yasui, Y. Kobayashi, M. Sato, T. Nishioka, M. Kontani, and K. Sano, “Spin gap behavior of $S$= $\frac{1}{2}$ quasi-two-dimensional system CaV4O9,” J. Phys. Soc. Jpn. 64, 2758 (1995).
  * Moriya (1956) T. Moriya, “Nuclear magnetic relaxation in antiferromagnetics,” Prog. Theor. Phys. 16, 23 (1956).
  * Giamarchi and Tsvelik (1999) T. Giamarchi and A. M. Tsvelik, “Coupled ladders in a magnetic field,” Phys. Rev. B 59, 11398 (1999).
  * Roca  _et al._ (1998) M. Roca, P. Amorós, J. Cano, M. Dolores Marcos, J. Alamo, A. Beltrán-Porter, and D. Beltrán-Porter, “Prediction of magnetic properties in oxovanadium(IV) phosphates: The role of the bridging PO4 anions,” Inorg. Chem. 37, 3167 (1998).
  * Note (2) Similar to Ref. Tsirlin  _et al._ (2011), magnetization curve was simulated by including a weak interchain coupling of $J_{\perp}/J_{1}=-0.05$ in order to reproduce the experimental data around $H_{\rm c1}$. This interchain coupling is much lower than estimated by DFT (Table 2), possibly because $J_{\perp}$ reflects a cumulative effect of the frustrated couplings $J_{a1}$, $J_{a2}$ and $J_{c1}$, $J_{c2}$.
  * Johnston  _et al._ (2001) D.C. Johnston, T. Saito, M. Azuma, M. Takano, T. Yamauchi, and Y. Ueda, “Modeling of the magnetic susceptibilities of the ambient- and high-pressure phases of (VO)2P2O7,” Phys. Rev. B 64, 134403 (2001).
  * Garrett  _et al._ (1997) A. W. Garrett, S. E. Nagler, D. A. Tennant, B. C. Sales, and T. Barnes, “Magnetic excitations in the $S$ = $1/2$ alternating chain compound (VO)2P2O7,” Phys. Rev. Lett. 79, 745 (1997).
  * Nohadani  _et al._ (2004) O. Nohadani, S. Wessel, B. Normand, and S. Haas, “Universal scaling at field-induced magnetic phase transitions,” Phys. Rev. B 69, 220402(R) (2004).
  * Mazurenko  _et al._ (2014) V.V. Mazurenko, M.V. Valentyuk, R. Stern, and A.A. Tsirlin, “Nonfrustrated interlayer order and its relevance to the Bose-Einstein condensation of magnons in BaCuSi2O6,” Phys. Rev. Lett. 112, 107202 (2014).
  * Allenspach  _et al._ (2020) S. Allenspach, A. Biffin, U. Stuhr, G. S. Tucker, S. Ohira-Kawamura, M. Kofu, D. J. Voneshen, M. Boehm, B. Normand, N. Laflorencie, F. Mila, and Ch. Rüegg, “Multiple magnetic bilayers and unconventional criticality without frustration in BaCuSi2O6,” Phys. Rev. Lett. 124, 177205 (2020).