$J_1$-$J_2$ square lattice antiferromagnetism in the orbitally quenched insulator MoOPO$_4$

TL;DR

This study investigates MoOPO$_4$, a $4d^1$ insulator, revealing a $J_1$-$J_2$ Heisenberg model on a square lattice with a magnetic transition at 16 K, negligible spin-orbit effects, and collinear antiferromagnetic order.

Contribution

It provides detailed experimental and theoretical evidence that MoOPO$_4$ hosts a pure spin-$1/2$ magnetic moment with a $J_1$-$J_2$ model, highlighting $4d$ molybdates as quantum magnet candidates.

Findings

Magnetic transition at 16 K confirmed by multiple techniques.

Mo$^{5+}$ ions have a spin-$1/2$ moment with negligible spin-orbit coupling.

Neutron diffraction shows collinear Néel order and a spin-flop transition.

Abstract

We report magnetic and thermodynamic properties of a $4d^1$ (Mo$^{5+}$) magnetic insulator MoOPO$_4$ single crystal, which realizes a $J_1$-$J_2$ Heisenberg spin-$1/2$ model on a stacked square lattice. The specific-heat measurements show a magnetic transition at 16 K which is also confirmed by magnetic susceptibility, ESR, and neutron diffraction measurements. Magnetic entropy deduced from the specific heat corresponds to a two-level degree of freedom per Mo$^{5+}$ ion, and the effective moment from the susceptibility corresponds to the spin-only value. Using {\it ab initio} quantum chemistry calculations we demonstrate that the Mo$^{5+}$ ion hosts a purely spin-$1/2$ magnetic moment, indicating negligible effects of spin-orbit interaction. The quenched orbital moments originate from the large displacement of Mo ions inside the MoO$_6$ octahedra along the apical direction. The ground state is shown by neutron diffraction to support a collinear N\'eel-type magnetic order, and a spin-flop transition is observed around an applied magnetic field of 3.5 T. The magnetic phase diagram is reproduced by a mean-field calculation assuming a small easy-axis anisotropy in the exchange interactions. Our results suggest $4d$ molybdates as an alternative playground to search for model quantum magnets.

Tables

Table 1. Table 1: Basis functions of irreducible representation Γ ν subscript Γ 𝜈 \Gamma_{\nu} for k = ( 100 ) 𝑘 100 k=(100) separated into real and imaginary components and resolved along the crystallographic axes. The two equivalent Mo 1 and Mo 2 ions are related by an inversion through the origin.

$\nu$		Mo1	Mo2
1	Re	$(0,0,1)$	$(0,0,1)$
2	Re	$(0,0,1)$	$(0,0,\overline{1})$
3	Re	$(1,0,0)$	$(1,0,0)$
3	Im	$(0,\overline{1},0)$	$(0,\overline{1},0)$
4	Re	$(1,0,0)$	$(\overline{1},0,0)$
4	Im	$(0,\overline{1},0)$	$(0,1,0)$
5	Re	$(1,0,0)$	$(1,0,0)$
5	Im	$(0,1,0)$	$(0,1,0)$
6	Re	$(1,0,0)$	$(\overline{1},0,0)$
6	Im	$(0,1,0)$	$(0,\overline{1},0)$

Table 2. Table 2: Relative energies of d 𝑑 d -level states of Mo 5+ ion obtained from CASSCF/NEVPT2 calculations. The corresponding wavefunctions without (coefficients) and with (weights) SOC at the CASSCF level are also provided. At NEVPT2 level, the wavefunction would also contain contributions from the inactive and virtual orbitals. For simplicity only the weights of the SOC wavefunction are provided as they are complex.

$t_{2g}^1$ states	Relative	Wavefunction (CASSCF)
	energies (eV)
Without SOC :		coefficients
$|{\varphi }_0⟩$	0	$0.95|xy⟩-0.32|x^2-y^2⟩$
$|{\varphi }_1⟩$	1.79	$0.98|yz⟩+0.21|zx⟩$
$|{\varphi }_2⟩$	1.79	$0.21|yz⟩-0.98|zx⟩$
$|{\varphi }_3⟩$	3.68	$0.32|xy⟩-0.95|x^2-y^2⟩$
$|{\varphi }_4⟩$	4.42	$1.00|z^2⟩$
With SOC :		normalized weights (%)
$|{\psi }_0⟩$	0	$86.0|{\varphi }_0,\uparrow ⟩+14.0|{\varphi }_0,\downarrow ⟩$
$|{\psi }_1⟩$	1.75	$50.0|{\varphi }_1,\uparrow ⟩+50.0|{\varphi }_2,\downarrow ⟩$
$|{\psi }_2⟩$	1.82	$46.0|{\varphi }_1,\uparrow ⟩+46.0|{\varphi }_2,\uparrow ⟩$
		$4.0|{\varphi }_1,\downarrow ⟩+4.0|{\varphi }_2,\downarrow ⟩$
$|{\psi }_3⟩$	3.70	$88.0|{\varphi }_3,\uparrow ⟩+12.0|{\varphi }_3,\downarrow ⟩$
$|{\psi }_4⟩$	4.44	$100.0|{\varphi }_4⟩$

Table 3. Table 3: Computed g 𝑔 g factors of MoOPO 4 at the NEVPT2 level of theory. The ground state multiconfiguration wavefunction as shown in Table  2 produces the correct structure for the g 𝑔 g factors.

CASSCF	$g_a$	$g_c$
active orbital space
$t_{2g}$	1.91	1.99
$t_{2g}+e_g$	1.92	1.84

Equations

$$C_D=9n_{D}R{\left(\frac{T}{{\Theta }_D}\right)}^3{\int }_0^{{\Theta }_D/T}\frac{x^4{e}^x}{{\left(e^x-1\right)}^2}dx,$$

$$C_E=3n_{E}R\frac{y^2{e}^y}{{\left(e^y-1\right)}^2},y\equiv {\Theta }_E/T,$$

$${\Theta }_{CW}=-\frac{S(S+1)}{3k_B}\sum _{i=1,2}{z}_i{J}_i,$$

$$T_N=-\frac{S(S+1)}{3k_B}\sum _{i=1,2}{z}_i(-1{)}^i{J}_i,$$

Full text

I Introduction

The $4d$ transition-metal oxides naturally bridge the two different regimes of the strongly correlated $3d$ compounds and the $5d$ compounds with strong spin-orbit coupling (SOC) Khomskii (2014). To what extent the $4d$ compounds represent the either regime or display original properties are largely open questions of current interest Witczak-Krempa et al. (2014). Most notably, for instance, it is intriguing that seemingly similar $\mathrm{Ca_{2}RuO_{4}}$ and $\mathrm{Sr_{2}RuO_{4}}$ display totally different behavior: the former is a Mott insulator Nakatsuji et al. (1997); Cao et al. (1997); Nakatsuji and Maeno (2000); Fatuzzo et al. (2015) while the latter is a metal and becomes superconducting at low temperature Nakatsuji and Maeno (2000); Maeno et al. (1994); Ishida et al. (1998); Fatuzzo et al. (2015). Despite a large interest, however, purely $4d$ quantum (spin-$1/2$) magnets are rather rare de Vries et al. (2010); Aharen et al. (2010); Clark et al. (2014); Ishikawa et al. (2017) as the electronic structure is often complicated by the presence of other types of, e.g., $3d$ or $4f$ magnetic orbitals Gardner et al. (2010).

Among a few known $4d^{1}$ magnets de Vries et al. (2010); Clark et al. (2014); Ishikawa et al. (2017) a molybdenium phosphate $\mathrm{MoOPO_{4}}$ is reported Kierkegaard and Longo (1970). The $\mathrm{MoO_{6}}$ octahedra with Mo5+ ions are corner shared to form a chain along the crystallographic $c$ axis of the tetragonal structure (Fig. 1(a)), and these chains are further coupled next to each other via corner sharing $\mathrm{PO_{4}}$ tetrahedra (Fig. 1(b)) Kierkegaard and Longo (1970); Canadell et al. (1997). Previous susceptibility data on a powder sample of $\mathrm{MoOPO_{4}}$ shows a Curie-Weiss behavior with antiferromagnetic $\Theta_{CW}=-8$ K and a magnetic transition at 18 K Lezama et al. (1990). 31P NMR on a powder evidences a substantial exchange through the $\mathrm{PO_{4}}$ tetrahedra and a sharp powder ESR line infers a rather isotropic $g$ factor Lezama et al. (1990). However, so far there have not been any studies on the magnetic structure in the ordered state nor magnetic properties of a single crystal. Moreover, any discussion on possible interplay between the crystal electric field and SOC is absent.

Here we report the magnetic and thermodynamic properties of a $\mathrm{MoOPO_{4}}$ single crystal using specific heat, susceptibility, magnetization, ESR, and neutron diffraction experiments. We also elucidate the electronic states and magnetic aspects in light of SOC and crystal field effects, with the help of ab initio quantum chemistry calculations.

II Experimental Details

High-quality single crystals of $\mathrm{MoOPO_{4}}$ were grown following the procedure described in Ref. Kierkegaard and Longo (1970). $\mathrm{H_{2}MoO_{4}}$ was mixed with concentrated phosphoric acid and heated up to 1000 *∘*C for reaction in an open platinum crucible. After being cooled to room temperature, the resulting dark-blue solid was dissolved in a large amount of hot water. The yellow transparent crystals were obtained in a plate-like shape (Fig. 1(c)). Large-sized crystals have a typical dimension of $3\times 2\times 0.4$ mm3 with the $c$ axis normal to the plate. The crystal belongs to the space group of $P4/n$ with the lattice parameters of $a=b=6.2044$ Å and $c=4.3003$ Å, obtained by a single-crystal x-ray diffraction, in agreement with Ref. Kierkegaard and Longo (1970).

Specific heat was measured using a physical properties measurement system (PPMS, Quantum Design, Inc.), and magnetization was measured using a magnetic properties measurement system SQUID (MPMS, Quantum Design, Inc.). ESR measurements were performed using a Bruker X-band spectrometer with a TE102 resonant cavity around 9.4 GHz. Neutron diffraction experiments were performed on TRICS and D23 beamlines at PSI and ILL, respectively. An incident neutron wavelength of 2.3109 Å was employed.

III Results

III.1 Specific heat

Figure 2(a) shows the specific heat $C_{p}$ measured from 2 to 150 K in zero field and in a magnetic field of 14 T. The $C_{p}$ above 25 K for the both fields is essentally the same, increasing monotonically with increasing temperature. In zero field a pronounced peak is found at 16.1 K, while the peak is shifted to a slightly lower temperature of 15.4 K in 14 T. These peaks correspond to a transition into a magnetically long-range-ordered phase, as evidenced by other experimental measurements discussed in later sections.

In order to extract the magnetic part of the specific heat, $C_{\mathrm{mag}}$, and to deduce the corresponding entropy $S_{\mathrm{mag}}$, we simulate the lattice contribution from the high temperature data by taking into account the Debye and Einstein contributions. We fit the $C_{p}$ data above 30 K by a lattice-only model $C_{p}=C_{D}+\sum_{i}C_{E,i}$, where $C_{D}$ and $C_{E,i}$ represent the Debye and Einstein terms, respectively. The Debye term is expressed as

[TABLE]

and the Einstein term as

[TABLE]

where $R$ denotes the gas constant, $\Theta_{D}$ and $\Theta_{E}$ are the Debye and Einstein temperatures, and $n_{D}$ and $n_{E}$ are the numbers of the corresponding modes, respectively; the sum $n_{D}+n_{E}$ is the total number of atoms per formula unit. The best fit for the zero-field, using one Debye and two Einstein terms, yields the characteristic temperatures $\Theta_{D}$ = 1177 K, $\Theta_{E,1}=372$ K, and $\Theta_{E,2}=154$ K, and the numbers $n_{D}$ = 4, $n_{E,1}$ = 2, $n_{E,2}$ = 1. Solid line in Fig. 2(a) is the best fit result for the total lattice contribution while dash-dotted and dashed lines are the corresponding Debye and Einstein contributions, respectively. While the parameters in the phonon fit may not be directly physical, they provide a parametrisation of the lattice contribution to the specific heat, which can be substracted to estimate the magnetic specific heat.

Figure 2(b) shows the resulting $C_{\mathrm{mag}}$ in zero field (circle, left axis) obtained by subtracting the lattice contribution from the measured $C_{p}$. Solid line in Fig. 2(b) plots the $S_{\mathrm{mag}}(T)$ obtained by integrating $C_{\mathrm{mag}}/T$ over temperature (right axis). $S_{\mathrm{mag}}(T)$ is found to reach and stay at $R\ln 2$ at high temperatures, indicating two-level degrees of freedom. The thin colored band in Fig. 2(b) represents the entropy range obtained when fitting the $C_{p}$ data by varying the lower bound of temperature between 25 K to 35 K, to confirm the little dependence of the result on the chosen fit range. The similar analysis for the 14 T data (not shown) indicates negligible field effects.

III.2 Susceptibility and magnetization

Figure 3(a) shows the DC magnetic susceptibility $\chi=M/H$, where $M$ is magnetization, in a field of $H=0.1$ T applied parallel and perpendicular to the $c$ axis. For both cases, $\chi(T)$ shows almost identical behavior from 300 K down to 20 K. However, for $H\parallel c$ the $\chi(T)$ exhibits a sharp drop toward zero as temperature is decreased across 17 K, while the one for $H\perp c$ remains only weakly temperature dependent. This is indicative of antiferromagnetic transition where the ordered moments at low temperatures are collinear to each other, and parallel to the $c$ axis.

The nearly isotropic, high-temperature part of $\chi(T)$ could be well fit by the Curie-Weiss formula, $\chi(T)=C/(T-\Theta_{CW})+\chi_{0}$, where $\Theta_{CW}$ is the Curie-Weiss temperature and $\chi_{0}$ a temperature-independent diamagnetic and background term. The best fit yields the effective moment $\mu_{\mathrm{eff}}=1.67(1)\,\mu_{\mathrm{B}}$ per Mo5+ ion, $\Theta_{CW}=-6(1)$ K, and $\chi_{0}=2.2(1)\times$ $10^{-4}$ emu/mol for $H\parallel c$, and $\mu_{\mathrm{eff}}=1.69(1)\,\mu_{\mathrm{B}}$, $\Theta_{CW}=-4(1)$ K, and $\chi_{0}=4.6(1)\times 10^{-4}$ emu/mol for $H\perp c$. The best fit for $H\parallel c$ is shown as dashed line in Fig. 3(a). The negative $\Theta_{CW}$ indicates that antiferromagnetic interactions are dominant. The effective moments indicate a spin-only value consistent with the specific-heat results.

The isothermal magnetization $M(H)$ for $H\parallel c$ and $H\perp c$ at several temperatures are shown in Fig. 3(b). At 5 K, $M(H)$ increases slowly with the field $H\parallel c$ up to 3 T, but then sharply increases in a narrow field range of $3-4$ T until it eventually converges to the high temperature $M(H)$ data obtained at 16 K or 20 K. This stepwise increase of $M(H)$ becomes smeared out as temperature is increased. On the other hand, no such stepwise behavior was observed at any temperatures for $H\perp c$. These are typical signatures of a spin-flop transition which occurs when the field is applied along an easy axis, along which the ordered moments align.

The magnetic phase diagram is thus mapped out by combining the above bulk magnetic and specific-heat results, as shown in Fig. 3(c). The antiferromagnetic transition temperatures in different fields are obtained from the peaks in $\chi(T)$ and $C_{p}(T)$, and the spin-flop transition fields at different temperatures are obtained from the peak positions in the $dM/dH$ versus $H$ plot (inset of Fig. 3(b)).

III.3 Electron spin resonance

In order to gain a microscopic insight into the magnetic properties, we have performed ESR measurements as a function of field orientation and temperature. Figure 4(a) plots the obtained room-temperature $g$ factor as the field direction is rotated by $\phi$ in the $ab$ and the $ac$ planes. The $g$ factor in the $ac$ plane shows a $\phi$-variation as large as 2 $\%$ with characteristic $\cos^{2}\phi$ angular dependence. On the other hand, the $g$ factor in the $ab$ plane remains essentially constant as expected from the tetragonal symmetry, within the error of 0.08 % which might have arisen from a slight misorientation of the crystal. We obtain the $g$ factor along the principal axes as $g_{a}=1.926(2)$ and $g_{c}=1.889(2)$. The average value $g=(2g_{a}+g_{c})/3=1.913(2)$ agrees with the one previously obtained by powder ESR Lezama et al. (1989). These $g$ values correspond to the effective moments of 1.64 and 1.66 $\mu_{B}$ for spin-$1/2$, for $H\parallel c$ and $H\perp c$, respectively, which are very close to the effective moment values obtained from the Curie-Weiss fit in the previous section.

For a system with tetragonal symmetry with short distances between the transition metal and ligand ions, one would expect $g_{a}<g_{c}$ Abragam and Bleaney (1970). However, we find an opposite structure for the $g$ factor in MoOPO4, even though the orbital energy diagram for the Mo5+ ion is expected to be similar to that of tetragonally compressed octrahedron with stabilized $d_{xy}$ orbital (see Fig. 7). As explained in Sect. E, the multi-orbital character of the ground state in MoOPO4 results in the observed $g$ values.

Figure 4(b) shows the temperature dependence of the resonance field $B_{0}$ and the line width $\Delta B_{0}$ of the ESR spectrum. $B_{0}$ slowly decreases as temperature is lowered from 300 K down to 24 K, which may be attributed to a lattice contraction. As temperature is further lowered below 24 K, $B_{0}$ starts increasing sharply, which indicates that a magnetic transition is approached. Similarly, $\Delta B_{0}$ slowly decreases as temperature is lowered down to 25 K, but then starts broadening significantly as temperature is further lowered down to 15 K due to critical spin fluctuations. Inset of Fig. 4(b) plots temperature dependence of the local spin susceptibility which is obtained from the spectral area at each temperature normalized by the one at 300 K, $\chi_{s}(T)/\chi_{s}(\mathrm{300\,K)}$. The data could be fit to the Curie-Weiss formula with $\Theta_{CW}=-8.9$ K, which is in reasonable agreement with the bulk susceptibility result shown in Fig. 3(a).

Across the transition, the ESR line changes in shape and intensity as shown in Fig. 4(c). The line sustains a perfect Lorenztian shape down to 16 K. On the other hand, the line below 16 K close to the transition fits better to a sum of two Lorentzian: one corresponds to the intrinisic sample susceptibility while the other may correspond to some defects. Indeed, the ESR signal at the paramagnetic resonance field position below 15 K corresponds to about 0.1 % concentration of paramagnetic impurities. The response below 15 K represents the summation of the possible defect contribution and the intrinsic susceptibility. The tiny intrinsic response below the transition temperature may represent clusters of which spins remain fluctuating within the ESR time window, which essentially disappears at lower temperatures below 14 K. At 5 K, a broad hump of weak signal is observed around 0.27 T for $\mathrm{B}\parallel c$, which is absent for $B\perp c$. This signal may correspond to an antiferromagnetic resonance.

III.4 Neutron diffraction

To determine the microscopic magnetic structure, we have performed neutron diffraction measurements. Magnetic intensity appears at the position of the $k=(100)$ wave vector at 5 K, as shown in the rotation scan in Fig. 5(a). No appreciable change in scattering is found close to $(001)$ between 5 and 25 K as shown in Fig. 5(b). A small shoulder of the $(001)$ reflection is likely to originate from a closely oriented secondary grain. A non-zero $(100)$ reflection would be consistent with Mo spins related by a spatial inversion being antiparallel. Due to the dipolar nature of the magnetic interaction, only magnetization perpendicular to the scattering wavevector gives a non-zero structure factor. As no change is observed for the $(001)$ reflection upon cooling below $T_{N}$, we can conclude that the moments are parallel to the $c$ axis. To verify that this is consistent with the symmetry of the lattice and rule out any other magnetic structures, we utilize BasiReps Rodríguez-Carvajal (1993) and outline the results here.

The magnetic representation is decomposed into six one-dimensional irreducible representations $\Gamma_{\nu}$ whose resulting basis functions are shown in Table 1. Examining the results of the irreducible magnetic representations, we find that only $\Gamma_{2}$ is consistent with our observations. These results are in contrast to the closely related $A$MoO(PO4)Cl ($A=$ K and Rb) materials. Unlike the tilted arrangement of MoO6 octahedra and PO4 tetrahedra in MoOPO4, $A$MoO(PO4)Cl possesses a higher symmetry where the octahedra and tetrahedra are arranged untilted in the $ab$ plane Borel et al. (1998). Powder neutron diffraction measurements on $A$MoO(PO4)Cl reveal an antiferromagnetic structure where Mo moments are instead confined to the $ab$ plane Ishikawa et al. (2017).

Figure 5(c) shows the temperature dependence of the (100) Bragg peak integrated intensity. By fitting a powerlaw dependence to the intensity, we find $T_{N}=16.17\pm 0.06$ K – consistent with the magnetization and specific-heat measurements. The order parameter exponent is found to be $\beta=0.23$, corresponding to the 2D XY universality class, However, dedicated measurements with better resolution and separating critical scattering would be needed before any conclusions are drawn from this. In Fig. 5(d) we show the magnetic Bragg peak intensity as a function of applied field along the $c$ axis recorded at 2 K. Above 3 T, we find a sharp decrease in intensity which then appears to saturate above 5 T. The change in the Bragg peak intensity is consistent with a spin-flop transition that is observed in the magnetization measurements shown in Fig. 3(b). This corresponds to a tilt of the moments by approximately 35∘ away from the $c$ axis for the fields above 5 T.

III.5 Model calculations

In order to gain insight into the magnetic interactions, we have fit the experimental susceptibility shown in Fig. 3(a) using a high-temperature series expansion Schmidt et al. (2011) assuming a $J_{1}$-$J_{2}$ spin-$1/2$ Heisenberg model on a square lattice. The best fit (solid line in Fig. 3(a)) returns $J_{1}=11.4(0.4)$ K and $J_{2}=-5.2(1.0)$ K corresponding to $J_{2}/J_{1}=-0.46$. This ratio supports a collinear Néel order for the ground state (see Fig. 6) in agreement with the neutron diffraction result. Using the mean-field expression for the Curie-Weiss temperature,

[TABLE]

where $z_{i}$ is the number of neighbors for the corresponding couplings (4 both for $J_{1}$ and $J_{2}$ in the present case), the high-temperature expansion fit yields $\Theta_{CW}=-6.2$ K, which agrees with the value obtained from the simple Curie-Weiss fit. Next, we simulate the phase diagram using a mean-field calculation. The results are presented by the colored background in Fig. 3(c). A slight exchange anisotropy, $\Delta=0.02$, has been introduced in the Hamiltonian, i.e., for a pairwise interaction $\mathcal{H}=J[S^{x}S^{x}+S^{y}S^{y}+(1+\Delta)S^{z}S^{z}]$, to account for the spin-flop transition in a spin-$1/2$ system where single-ion anisotropy is not expected to be present. We note that the mean-field calculation reproduces the temperature dependence of the spin-flop field. From the mean-field expression for the Néel temperature,

[TABLE]

we obtain $T_{N}=16.6(1.4)$ K, which is in excellent agreement with the actual value from the experiments.

With strong ferromagnetic second nearest neighbor interactions, MoOPO4 populates a region of the $J_{1}$-$J_{2}$ phase diagram, which has so far seen rather few investigations (see Fig. 6). In the context of $(\pi,0)$ zone-boundary anomalies on the square lattice, linear spin-wave theory would for MoOPO4 predict a dispersion with significantly higher energy at $(\pi,0)$ than at $(\pi/2,\pi/2)$, opposite to the case of weak antiferromagnetic $J_{2}$ in $\mathrm{Cu(pz)_{2}(ClO_{4})_{2}}$ Tsyrulin et al. (2010). Compared to the 39 % reduction in ordered moment due to quantum fluctuations for the nearest-neighbour Heisenberg model, the estimate for $J_{2}/J_{1}=-0.46$ is only a 24 % reduction of the ordered moment. Adding the weak anisotropy for MoOPO4 yields 21 % reduction in ordered moment. Hence quantum fluctuations are likely much weaker in MoOPO4 than in e.g. CFTD Rønnow et al. (1999, 2001); Christensen et al. (2007); Dalla Piazza et al. (2015) or Sr2CuTeO6 Babkevich et al. (2016), and it would be interesting in future investigations to examine whether this leads to a similar suppression of the quantum dispersion and continuum around $(\pi,0)$.

III.6 Ab initio calculations

An interesting feature in $\mathrm{MoOPO_{4}}$ is that the axial position of the Mo4+ ion inside MoO6 octahedron is heavily shifted such that the short and long distances to the apical oxygens are 1.652 Å and 2.641 Å, respectively. As a consequence the octahedral symmetry around the Mo ion is reduced, resulting in the removal of octahedral orbital degeneracies and an orbitally mixed electronic ground state. To elucidate the electronic levels of the Mo5+ ion in low-symmetry crystal fields in $\mathrm{MoOPO_{4}}$, we performed ab initio quantum chemistry calculations using the cluster-in-embedding formalism Hozoi and Fulde (2011). A cluster of a single active MoO6 octahedron along with surrounding nearest neighbor (NN) PO4 tetrahedra within the plane and the out-of-plane MoO6 octahedra embedded in an array of point charges that reproduces the effect of the solid enviroment Klintenberg et al. (2000) was considered for many-body calculations. The NN polyhedra were included within the cluster region to better describe the charge density within the active MoO6 region. Such calculations have provided excellent insights into the interplay of crystal field and SOC effects for several $4d$ and $5d$ transition metal compounds Katukuri et al. (2014, 2012); Bogdanov et al. (2012); Liu et al. (2012).

A perfect octahedral arrangement of the oxygen ligands around the transition-metal ion splits the $d$ levels into high-energy $e_{g}$ and low-energy $t_{2g}$ manifolds (see Fig. 7). In $\mathrm{MoOPO_{4}}$, the low-symmetry crystal fields further split the $t_{2g}$ and $e_{g}$ levels of the Mo5+ ion resulting in an orbital singlet ground state. In Table 2 the ground-state wavefunction and the $d$-$d$ excitations of the Mo5+ ion are summarized. These have been obtained from many-body multi-configurational self consistent field (MCSCF) Helgaker et al. (2000) and $N$-electron valence-state perturbation theory (NEVPT2) Angeli et al. (2001) calculations for the atoms in the active cluster region. All-electron DKH (Douglas-Kroll-Hess) basis sets of triple-zeta quality Pantazis et al. (2008) were used to represent the Mo and oxygen ions in the central MoO6 octahedron, and for the Mo and P ions in the NN polyhedra we employed effective core potentials Peterson et al. (2007); Igel-Mann et al. (1988) with valence triple-zeta Peterson et al. (2007) and a single basis function, respectively. The oxygen ions corresponding to the NN MoO6 and PO4 polyhedra were expanded in atomic natural orbital (ANO) type Pierloot et al. (1995) two $s$ and one $p$ functions. All the caluations were performed using the orca quantum chemistry package Neese (2012)

In the complete active space formalism of the MCSCF (CASSCF) calculation, a self-consistent wavefunction was constructed with an active space of one electron in five Mo $d$ orbitals. On top of the CASSCF wavefunction, NEVPT2 was applied to capture the dynamic electronic correlation. It can be seen in Table 2 that the ground state is predominantly of $d_{xy}$ character but has significant contributions from the $d_{x^{2}-y^{2}}$ orbital. The first orbital exciations are nearly degenerate at 1.79 eV and are composed of $d_{yz}$ and $d_{zx}$ atomic orbitals. This scenario contrasts the situation in other $t_{2g}$ class of compounds with regular transition-metal oxygen octahedra where the $t_{2g}$ manifold remains degenerate with an effective orbital angular momentum ${\tilde{l}}=1$. In the latter scenario the spin-orbit interaction admixes all the $t_{2g}$ states to give rise to a total angular momentum $J_{\mathrm{eff}}$ ground state Kim et al. (2009, 2008). Due to the large non-cubic crystal field splittings in the $t_{2g}$ manifold in $\mathrm{MoOPO_{4}}$, the spin-orbit interaction has negligble effect on the Mo5+ ground state $\psi_{0}$, see ‘with SOC’ results in Table 2. However, the orbital angular momentum is unquenched in the $d_{zx}$ and $d_{yz}$ and hence the SOC results in the splitting of the high-energy $\psi_{1}$ and $\psi_{2}$. Our calculations result in excitation energies of 3.68 eV and 4.42 eV into the $e_{g}$ states.

To understand the unusual structure of $g$ factors deduced from the ESR experiments, we computed the same from the ab initio wavefunction as implemented in orca Neese (2007). In Table 3, the $g$ factors obtained from CASSCF calculations with two different active orbital spaces, $t_{2g}$ only and $t_{2g}+e_{g}$, are presented. With only $t_{2g}$ orbitals in the active space, we find $g_{a}<g_{c}$ as expected for tetragonal symmetry with $d_{xy}$ ground state. By enlarging the active space, the multiconfiguration wavefunction now contains configurations involving the $e_{g}$ orbitals and this is crucial to produce the experimentally observed $g$ factors with $g_{a}>g_{c}$.

IV Conclusion

We have shown by a variety of experimental and computational techniques that $\mathrm{MoOPO_{4}}$ realizes a spin-$1/2$ magnetic system of $4d^{1}$ electrons, with the quenched orbital moment due to the large displacement of the Mo ions inside the MoO6 octahedra. The magnetic ground state supports a Néel-type collinear staggered order on the square-lattice with the moments pointing normal to the plane, while the moments align ferromagnetically along the stacking axis. The compound likely realizes a spin-$1/2$ Heisenberg model on a $J_{1}$-$J_{2}$ square lattice, with an unfrustrated configuration of antiferromagnetic $J_{1}$ and ferromagnetic $J_{2}$. The spin-flop transition suggests a small easy-axis anisotropy in the dominant antiferromagnetic exchange, and the mean-field calculation reproduces the experimental magnetic phase diagram. The small anisotropy in the $g$ factor observed in ESR, which is reproduced by the quantum chemistry calculations, points to that the ground state involves the higher-energy $e_{g}$ orbitals in addition to the $t_{2g}$ orbitals. Our results suggest that $4d$ molybdates provide an alternative playground to search for model quantum magnets other than $3d$ compounds.

Acknowledgment

We thank R. Scopelliti and O. Zaharko for their help with the x-ray and neutron diffraction, respectively. We also thank V. Favre and P. Huang for their help with the specific heat analysis. V.M.K is grateful to H. Stoll for discussions on effective core potentials. This work was supported by the Swiss National Science Foundation, the MPBH network, and European Research Council grants CONQUEST and TopoMat (No. 306504). M.J. is grateful to support by European Commission through Marie Skłodowska-Curie Action COFUND (EPFL Fellows). M.K. is supported by a Grants-in-Aid for Scientific Research (C) (JSPS, KAKENHI No 15K05140). The ab initio calculations have been performed at the Swiss National Supercomputing Centre (CSCS) under project s675.