Magnetic structure and spin dynamics of the quasi-2D antiferromagnet Zn-doped copper pyrovanadate

The powder sample of ZnCVO shows a pure phase with an identical structure with $\beta$-CVO, as shown by the Rietveld refinement on the x-ray diffraction patterns in Fig. 1 (a). This result is consistent with the previous work by Pommer et. al., [12] where the Zn_xCu_2-xV₂O₇ compound completely transformed to the $\beta$ phase at $x$ = 0.15. At lower doping concentrations, the samples show the mixed $\alpha-\beta$ phases, and the Zn concentration of $x$ = 0.15 is expectedly at the transition point. From the pure phase powder ZnCVO, the single-crystals of ZnCVO with the largest size of approximately 1 $\times$ 1 $\times$ 1 cm${}^{\text{3}}$ ($m\approx$ 1.5 g) were obtained using the vertical gradient furnace. The natural cleaved facet can be identified as the crystallographic $a$-axis similar to the $\beta$-CVO single-crystals [13]. The $\omega$-scan around $(020)$ Bragg peak using neutron scattering at BT7 (Fig. 5, inset) with a Gaussian fit yields a full-width-at-half-maximum (FWHM) equal to 0.38(4)^∘, indicative of good crystallinity. Rietveld refinements on the powder x-ray diffraction pattern obtained from the ground single-crystals, as shown in Fig. 1 (b), can also be fitted well with the reported $\beta$-CVO crystal structure [10]. In addition, powder neutron diffraction on ZnCVO and $\beta$-CVO powder samples were also performed at 30 K and 2.5 K for crystal structure and magnetic structure determination, respectively. At 30 K, the powder neutron diffraction patterns of both ZnCVO and $\beta$-CVO were refined against the reference $\beta$-CVO crystal structure as shown in Fig. 2. Despite the presence of Zn, the diffraction pattern shows a pure $\beta$-CVO phase without any trace of other phases. The refined occupancy of the Cu site from powder neutron data in Table 2 yields 0.97(1) suggesting that the doping concentration of Zn is approximately 3%, much lower than the stoichiometric ratio of 7.5%. The powder neutron diffraction pattern of $\beta$-CVO, on the other hand, shows some impurity peaks which can be indexed with CuV₂O₆ [30] ($\approx$ 9%) and Cu_0.63V₂O₅ [31] ($\approx$ 6%). The refined parameters obtained from both x-ray and neutron diffractions are summarized in Table 1 and 2, respectively.

Now we discuss the magnetic structure determinations on ZnCVO and $\beta$-CVO using powder neutron diffraction. As mentioned earlier in Section III.1, we prepared both ZnCVO and $\beta$-CVO to confirm that they share not only crystal structure but also magnetic structure. We start with the irreducible representation analysis using the program basirreps in the fullprof [33] suit. According to the crystallographic space group $C2/c$ with commensurate magnetic translation vector $\vec{k}=(0,0,0)$, there are four possible magnetic irreducible representations (IR) as described in Table 3. The corresponding Shubnikov magnetic space groups for $\Gamma_{1}$, $\Gamma_{2}$, $\Gamma_{3}$, and $\Gamma_{4}$ are $C2/c$, $C2/c^{\prime}$, $C2^{\prime}/c^{\prime}$, and $C2^{\prime}/c$, respectively [34]. With the assumption that ZnCVO and $\beta$-CVO have the same magnetic structure, we know that these systems undergo a paramagnetic to antiferromagnetic transition at the Néel temperature of $T_{\text{N}}\simeq$ 26 K which will be discussed in Section III.3. When we start considering the exchange couplings along with the nearest-neighboring pairs $J_{1}$, there are two equivalent bonds between the Cu²⁺ ions i.e., Cu1-Cu3 and Cu2-Cu4 (see Table 3 and Fig. 4). It was originally believed that this system was the antiferromagnetic spin chain with alternating $J_{1}-J_{2}$ bonds (not shown here). However, it has been later proposed using the DFT calculations [15] that this system could be better described by the complex anisotropic honeycomb network. In their proposed model, the leading antiferromagnetic exchange interactions are along two $J_{5}$ and one $J_{6}$ (Fig. 4) i.e., three bonds per site. However, there are still weak but non-negligible antiferromagnetic $J_{1}$ as well as the interplane $J_{14}$ couplings, making the spin network more complex than the simple honeycomb structure. It is therefore presumed that the Cu²⁺ atoms must align antiparallel with their neighbors through the most prominent exchange interactions, here $J_{1}$, $J_{5}$, and $J_{6}$. In addition, the previous magnetization measurements by He et. al., [35] on the single-crystals of $\beta$-CVO strongly suggested that the magnetic easy axis of this system was along the crystallographic $c$-axis. This suggests that the magnetic moment $m_{a}$ and $m_{b}$, despite their possible nonzero values, could be discarded.

Therefore among the four possible magnetic IRs, where we take into account the first nearest-neighbor couplings Cu1-Cu3 and Cu2-Cu4, we can rule out $\Gamma_{1}$ and $\Gamma_{3}$ where all spins align ferromagnetically along $m_{a}$. The reason that we pay attention to the first nearest neighbor is due to its strongest interaction as we will show later in Section III.4. This leaves us with the two most probable magnetic IRs i.e., $\Gamma_{2}$ and $\Gamma_{4}$. It is obvious that only $\Gamma_{4}$ yields antiferromagnetic interaction on all neighboring bonds whereas $\Gamma_{2}$ gives ferromagnetic coupling on the fifth nearest neighbor. This assumption is based primarily on the DFT results by Tsirlin et. al., [15] and Bhowal et. al., [16] (the citations will be omitted afterward when we mention the DFT results) where the predominant $J_{1}$, $J_{5}$, and $J_{6}$ bonds are all antiferromagnetic. With this initial analysis, we refined the powder neutron diffraction data at 2.5 K with $\Gamma_{1}$, $\Gamma_{2}$, and $\Gamma_{4}$ except for $\Gamma_{3}$ where the symmetry results in the ferromagnetic spin direction along $c$-axis. The magnetic structure of each IR is shown in Fig. 4. The refined patterns from ZnCVO data with $\Gamma_{1}$, $\Gamma_{2}$, and $\Gamma_{4}$ are shown in Fig. 3 (b) for comparison along with their corresponding refined parameters summarized in Table 4.

It should be noted that the powder neutron diffraction patterns of both ZnCVO and $\beta$-CVO samples show very weak magnetic intensities, especially in $\beta$-CVO, and most of them are on top of the structural peaks. It is therefore very difficult to precisely extract the magnetic moment from the refinement. Furthermore, there could be large uncertainties in the refined values of the magnetic moment, and the exact magnetic structure could deviate from our proposed model. In order to present the magnetic intensities from the powder samples, we subtract the 30 K patterns from that of 2.5 K patterns, on both raw data and on the refined results, as shown in the low 2$\theta$ range in Fig 3 where the magnetic scattering is the most intense. It can be seen that the magnetic Bragg peak positions of both samples are consistent with the fitted model. Despite the dilution of the Cu sites by Zn, the magnetic intensities of ZnCVO are more pronounced than those of $\beta$-CVO where the intensities are most likely within the statistical error. Although we attempted to refine the magnetic structure on the $\beta$-CVO data we could not extract the magnetic moment with a reliable value. We could only obtain the magnetic moment from the ZnCVO data. The best fit is obtained from $\Gamma_{4}$ with the refined magnetic moment $m_{c}$ = 0.72(9) $\mu_{B}$, the best among all three IRs. The refined pattern of ZnCVO with $\Gamma_{4}$ along with the plot of $\left|F^{\text{M}}_{\text{Cal}}\right|^{2}$ vs $\left|F^{\text{M}}_{\text{Obs}}\right|^{2}$ are shown in Fig. 3 (a). This magnetic structure will be further used in the spin-wave dispersion analysis in Section III.4.

Magnetic susceptibility of single-crystal ZnCVO was measured along two crystallographic axes i.e., $\chi_{\parallel a}$ with $H\parallel a$ (along the cleaved surface), and $\chi_{\perp a}$ with $H\perp a$ (parallel to the cleaved surface). The results, as shown in Fig. 6(a), reveal a broad peak at $T\approx$ 50 K indicating short-range correlations among the Cu²⁺ spins. The paramagnetic upturn below $T\approx$ 20 K can be observed. This upturn, which corresponds to approximately 0.006$\mu_{\text{B}}$ at the field of 1 T and at the base temperature, is 13 times smaller than the ferromagnetism observed in $\alpha$-Cu₂V₂O₇ [7] and most likely a result of the presence of defective magnetic sites where Cu²⁺ ions were substituted by Zn²⁺ and thus the free Cu²⁺ spins were produced [12, 36]. There is a large anisotropy between $\chi_{\parallel a}$ and $\chi_{{\perp a}}$ up to $T$ = 300 K similar to that observed in $\beta$-CVO by He et. al., [35]. This unusual anisotropy was suggested as a result of the Jahn-Teller distortion [37]. The similarity of the magnetic susceptibility behavior between ZnCVO in this work and $\beta$-CVO by the previous works, as well as our powder neutron diffraction data analysis, strongly suggest that both systems share the same magnetic properties.

The plot of inverse magnetic susceptibility versus temperature, shown in Fig. 6(b), can be fitted well with the Curie-Weiss law ($\chi=C/(T-\Theta)$) at $T>$ 100 K. The fit yields the Curie-Weiss temperature of $\Theta$ = -79(1) K (-89(1) K) with $H\perp a$ ($H\parallel a$) indicating the dominant antiferromagnetic exchange interactions, and the Curie-Weiss constant $C$ = 0.429(1) cm${}^{\text{3}}$K/molCu and 0.593(3) cm${}^{\text{3}}$K/molCu for $H\perp a$ and $H\parallel a$, respectively. The effective magnetic moment can be estimated to $\mu_{\text{eff}}=\sqrt{3k_{B}C/N_{A}}=1.852(4)\mu_{B}$ for $H\perp a$ and $2.17(1)\mu_{B}$ for $H\parallel a$. These values are slightly larger than the spin-only value of $\mu_{\text{eff}}=g\mu_{B}\sqrt{s(S+1)}=1.73\mu_{B}$ for $g=2$ and $S=1/2$. The Néel temperature, $T_{\text{N}}\simeq$ 26 K is obtained from the exponent fit to the order parameter scans as a function of temperature on the magnetic Bragg peaks using elastic neutron scattering as shown in Fig. 5. The fits were done in the range $15~{}\text{K}<T<30~{}\text{K}$, close to the phase transition temperature, using equation $I=I_{0}\left(1-T/T_{\text{N}}\right)^{2\beta}$. This value is consistent with the observed $\lambda$-like transition at around 26 K from the heat capacity measurement on $\beta$-CVO single-crystal [35]. The obtained critical exponent of $\beta\sim$ 0.2 is comparable to its cousin phase $\alpha$-Cu₂V₂O₇ [38].

It should be noted that the doping of Zn on Cu sites results in a dilution of magnetic spin and typically decreases the Néel temperature [39, 22, 36]. However, in this case the value of $T_{\text{N}}$ is nearly the same as that of $\beta$-CVO [13, 35]. In addition, the finite magnetic susceptibility below $T_{\text{N}}$ does not fit the paramagnetic impurity upturn as that observed in the powder sample by the previous works [12]. The lattice parameters obtained from the Rietveld refinements (Table 1) on the powder samples reveal that the lattice parameter $a$ slightly decreases while $b$, $c$, and the angle $\beta$ slightly increase upon the presence of Zn compared to pure $\beta$-CVO. Those lattice parameters on the ground single-crystals are also consistent with the powder ZnCVO sample. This suggests that the refined value of approximately 3% Zn substitution on Cu sites only slightly alters the overall lattice parameters and does not affect the macroscopic magnetic properties.

To further estimate the average exchange interactions, we performed QMC simulations and fit the resulting simulated data to the broad peak of the magnetic susceptibility, provided that the spin network model and the values of leading exchange interactions were predicted. We proceed with the very first report on the DFT results (here we label the couplings according to the order of nearest-neighbor distances. The notation used by Tsirlin et. al., in Ref [15] will be recalled in the parentheses). Among their various models, they suggested that the best realization of the spins network in $\beta$-CVO can be described by the fifth $J_{5}$ ($J_{1}$) and sixth $J_{6}$ ($J^{\prime}_{1}$) neighboring bonds, represented by the blue and green bonds, respectively, in Fig. 4 and the inset in Fig. 6 (a). These two bonds connect the Cu²⁺ ions into the irregular honeycomb network, i.e., three bonds per site, spanning the $bc$-plane when viewed along the $a$-axis. These honeycomb planes, according to the DFT results, are however not the perfect 2D since there are non-zero $J_{14}$ ($J_{\perp}$), represented by the magenta bonds in Fig. 4, that connect between the adjacent honeycomb planes. There is also the suspicious $J_{1}$, formerly believed to be the leading exchange interaction, that appeared to be non-negligible from the DFT making the spin network in this system to be topologically the anisotropic magnetic 2D lattice (four bonds per site) with weak interplane couplings.

In our QMC simulation, we, therefore, construct the 2D spin network with anisotropic exchange interactions $J_{1}$, $J_{5}$, and $J_{6}$ as shown in the inset of Fig. 6 (a). The values of these couplings were obtained from the spin-wave dispersion fit on our inelastic neutron scattering data which will be discussed in Section III.4. We simplify our spin network model by truncating the interplane fourteenth neighboring bond $J_{14}$ in the QMC simulation due to its very weak value. Although we fit the spin-wave dispersion based on the DFT model, the fitted parameters were obtained differently. Here we use the ratio $J_{1}:J_{5}:J_{6}$ = 1 : 0.61 : 0.25 for the QMC model. With $J_{1}-J_{5}-J_{6}$ interactions, the spin network resembles the irregular 2D edge-sharing trapezoid shape. We then conducted the QMC with the loop algorithm [40] using the simulation package alps [41]. The obtained QMC simulation result and the experimental magnetic susceptibility data were fitted using the equations,

with

where $N_{\text{A}}$, $\mu_{\text{B}}$, and $k_{\text{B}}$ are the Avogadro constant, Bohr magneton, and Boltzmann constant, respectively. The function $\chi^{*}(t)$ is the susceptibility as a function of reduced temperature $t=k_{\text{B}}T/J_{\text{max}}$ which was obtained by fitting the simulated QMC to the Padé approximant [42]. Here $J_{\text{max}}$ is $J_{1}$, the leading exchange interaction. The fitting parameters are the background $\chi_{0}$, the Landé $g$-factor, and the leading exchange interaction $J_{\text{max}}$ ($J_{1}$). The results are shown by the solid red lines in Fig. 6 (a) along with the two orthogonal magnetic field directions while the fitted parameters are summarized in Tabel 5. The QMC simulation fits well with the magnetic susceptibility data over the broad maximum from $T\simeq$ 35 K up to 300 K yielding the leading exchange interaction $J_{1}\simeq$ 73 K ($\simeq$ 6.4 meV). Although the fitted values of the Landé $g$-factors are slightly deviated between $H\parallel a$ and $H\perp a$ data due most likely to the anisotropy, their average $g_{\text{av}}$ = 2.09(1) is still very close to the theoretical value of 2.

All of the experimental data and analysis in the previous sections have led us to believe that the magnetic properties of ZnCVO could be a good realization of the $\beta$-CVO system. In this final section, we investigate the spin dynamics of ZnCVO single-crystals and analyze the obtained dispersion relation using the linear spin-wave theory (LSWT) [43, 44]. We measured inelastic neutron scattering along two directions around the magnetic zone center at (0,2,0) i.e., along $(0,k,0)$ and along $(h,2,0)$. At low energy transfer ($\hbar\omega~{}<~{}5$ meV), we conducted the experiments at SPINS and CTAX whereas at high energy transfer ($5$ meV $<\hbar\omega<15$ meV) the experiments were done at BT7. The intensity map along the two directions at the base temperatures (depending on the spectrometer) is shown in Fig 7 (a) - (b) and Fig 7 (e) - (f). Figure 7 (a) shows the whole extent of the dispersion along $(0,k,0)$ from the magnetic zone center at (0,2,0) to the zone boundary at (0,3,0). The dispersion reaches its maximum at the energy transfer of $\approx$ 11 meV. At the magnetic zone center, we can see an energy gap clearly when using CTAX and SPINS spectrometers in Fig. 7 (b) and (f) respectively. The dispersion is however different from its cousin phase $\alpha$-CVO where we found the splitting of the dispersion into two branches away from the magnetic zone center [7]. This splitting, as mentioned earlier, was due to the presence of the DM interaction. On the other hand, in the $\beta$-CVO system which in this case is the ZnCVO, the crystal is centrosymmetric and thus DM interaction is absent due to the symmetry of the underlying crystal structure. Therefore, as we expected, the magnon dispersion in ZnCVO shows only one symmetric branch without the bidirectional shift. This evidence is a great test that the nonreciprocal magnon vanishes in the centrosymmetric crystal in the Cu₂V₂O₇ system. Along the $(h,2,0)$, on the other hand, the dispersion gradually increases from (0,2,0) up to the magnetic zone boundary at (1,2,0). This suggests that the spin interactions along the reciprocal lattice $b^{*}$ i.e., within the anisotropic lattice plane, are stronger than those along $a^{*}$ between the planes, and that the interaction between the planes should be relatively weak compared to the in-plane interactions.

In order to quantitatively analyze the exchange coupling values, we extract the dispersion relation from the convolute fit to the energy scan at each $Q$. The obtained dispersions from all data sets were plotted altogether as shown by the circle symbol in Fig 7 (c) and (g) for $(0,k,0)$, and $(h,2,0)$ directions, respectively. Since the spin structure has been analyzed in Section III.3, we need to construct spin interactions network for the LSWT fit. Again we start with the predicted models by the DFT calculations. In their work [15], they performed various computational approaches and showed that the leading exchange interactions were $J_{5}$ and $J_{6}$ forming the anisotropic honeycomb network with weak but non-negligible $J_{1}$ and the interplane $J_{14}$. We started with this model by introducing $J_{1}$, $J_{5}$, $J_{6}$, and $J_{14}$ into our spin model. The Hamiltonian that we used in our spin-wave fit is shown in Eq. 3 below.

where $J_{ij}$ is the exchange interaction between spins $S_{i}$ and $S_{j}$, $\beta$ is the angle between the $a$-axis and $c$-axis due to the monoclinic system, and $G_{ij}=GJ_{ij}$, defined to be proportional to the exchange couplings, is the anisotropic parameter which gives rise to the spin gap at the magnetic zone center. We applied $J_{ij}$ and $G_{ij}$ to the first, fifth, sixth, and fourteenth neighboring bonds then fit the spin wave along both $(0,k,0)$ and $(h,2,0)$ directions simultaneously using least-square fitting routine to the modeled Hamiltonian.

Our result is, although qualitatively consistent with the DFT in terms of the representative leading exchange couplings, still quantitively deviated from the proposed honeycomb model. We note that for the data collected at BT7, despite its broad range covering from the magnetic zone center to the zone boundary, the resolution is rather low. It is possible that the exact values of the exchange interactions could slightly deviate from our results. As a result, despite the proposed honeycomb model with $J_{5}$ and $J_{6}$ as the leading exchange interactions, we instead get the largest value of 8.5(6) meV on $J_{1}$ which is also much higher than that in $\alpha$-CVO [7]. The fitted results are shown by the red lines in Fig. 7 whereas the fitted parameters are summarized in Table 6. These couplings yield the average in-plane exchange interactions $(J_{1}+2J_{5}+J_{6})/4$ = 5.3(2) meV. It should be noted that we also failed to fit our data when the second neighbor $J_{2}$ ($J_{a}$ in Ref [15]) was introduced, in agreement with the DFT that this bond is rather weak and hence the previously proposed spin-chain model for this system is unfeasible. Although $J_{14}$ is rather weak compared to $J_{1}$ and $J_{5}$, this bond is non-negligible. This evidence leads us to conclude that the spin network of ZnCVO should be better described by the anisotropic 2D lattice with weak interplane couplings. Lastly, the calculated intensities of the dispersion along both $(0,k,0)$ and $(h,2,0)$ directions using the parameters in Table 6 as shown in Fig. 7 (d) and (h) can well describe the measured intensity maps.