Lattice Distortions and Magnetic Interactions in Single-Layer VOCl

Atomically thin layers exfoliated from magnetic van der Waals layered materials are currently of high interest in solid state physics. VOCl is a quasi-two-dimensional layered antiferromagnet which was recently synthesized in monolayer form. Previous theoretical studies have assumed the high-temperature orthorhombic lattice symmetry also in the low temperature range, where the bulk system is known to be monoclinic due to a strong magnetoelastic coupling. We demonstrate from \textit{ab-initio} calulations that this monoclinic distortion is prevalent also in monolayers, which is in line with recent experimental indications of monoclinic symmetry. Our calculations also show that competing ferromagnetic and antiferromagnetic interactions give rise a frustrated two-fold magnetic superstructure where higher-order magnetic interactions play a key role to stabilize the observed magnetic ground state.


I. INTRODUCTION
The recent discovery of spontaneous long-range ferromagnetic order in the two-dimensional (2D) material CrI 3 [1] has lead to a surge in the search for such materials by experiments and theoretical calculations alike.Stable long-range ordering in low dimension that prospectively could be combined with various tunable properties make them appealing for next generation spintronics devices and functional materials [2][3][4][5][6][7][8].Yet, the fundamental understanding of magnetic interactions in such 2D magnets is still a developing field in solid state theory [9,10], as ferro-or antiferromagnetic order in a 2D spin array with isotropic interactions is forbidden at non-zero temperature by the Mermin-Wagner theorem [11].The observed long-range ordering is commonly attributed to magnetic anisotropy introducing a spin-wave excitation gap [12].In this pursuit, the magnetic van der Waals (vdW) layered materials receive considerable attention, as the weakly bonded layers may be easily exfoliated, and they can be expected to retain the magnetic properties of the bulk material [12][13][14].
VOCl is a layered vdW material consisting of V-O bilayers connected by Cl ions on each side, as illustrated in Fig. 1.In its bulk form, these bilayers are separated by a large vdW gap, taking orthorhombic P mmn symmetry (space group No. 58) at ambient conditions [15].The crystal structure is common to all the so-called transition metal oxychlorides, M OCl, where M ∈ {Ti, V, Cr, Fe}.
At room temperature, bulk VOCl is a paramagnetic insulator, but a twofold antiferromagnetic (AFM) superstructure develops below the Néel temperature, T N ≈ 80 K [16,17].The large vdW gap makes VOCl suitable for intercalation applications, and it is currently being considered for novel transistors [18] and battery archi-tectures, with a demonstrated stability to air exposure and cyclic ion shuttling [19][20][21].
Single crystals of VOCl with a thickness of only a few atomic layers were first synthesized by Wang et al. [22], and were shown to retain the crystal symmetry of the bulk form at room temperature.However, detailed measurements of the magnetic order of single-layers are challenging and scant.An ab-initio study by Marouche et al. [23] found ferromagnetic ordering to be the most favorable configuration on a single bilayer with orthorhombic symmetry.Subsequent theoretical studies [24,25] suggested that the corresponding AFM configuration observed in bulk VOCl (see Fig. 1) would constitute the magnetic ground state of the single-layers.
Theoretical studies have so far assumed the orthorhombic lattice structure experimentally observed at room temperature [23,25] The system is then highly frustrated, as each V ion is connected to two V ↑ and two V ↓ ions.Nevertheless, in bulk VOCl, the development of magnetic order below T N is accompanied by a monoclinic distortion of the crystal structure, lowering the symmetry to P 2 /n (space group No. 13) [16,17].Ab-initio calculations have shown that this distortion is related to magnetoelastic coupling, reducing the V ↑ -V ↓ distance [26].This monoclinic distortion appears to lift the apparent magnetic frustration that would prevail for the orthorhombic lattice, as the bonds, and the exchange interactions, would be equivalent by symmetry [16].Indeed, recent low temperature measurements on VOCl single-layers have inferred a monoclinic lattice symmetry, although the detailed lattice geometry and magnetic properties require further investigation [27].
In this study, we perform structural relaxation of VOCl monolayers by density functional theory (DFT) [28,29] calculations to show that the AFM configuration leads to a distortion of the lattice that is completely analogous to the monoclinic distortion of bulk VOCl.Assuming this lower lattice symmetry, we derive a magnetic Hamiltonian to study the role of exchange interactions, single-ion anisotropy and the Dzyaloshinskii-Moriya (DM) [30,31] interaction in the monoclinic phase.Monte Carlo simulations recover a T N comparable to the bulk form.Our study shows that, counterintuitively, the non-equivalent nearest and next-nearest neighbor interactions are both ferromagnetic; the AFM configuration is due to more long-ranged exchange interactions.This shows that the system remains frustrated even in the monoclinic phase.
The paper is structured as follows.In Section II we provide details of the electronic structure calculations and Monte Carlo simulations.In Section III we first report on the structural optimization and magnetic order.We then describe the electronic structure before detailing the magnetic interactions.Finally, in Section IV we discuss the implications of our results for VOCl single-layers and in the broader context of magnetic vdW layered materials.

II. COMPUTATIONAL DETAILS
Calculations were performed with the Quantum Espresso [32,33] code using the GBRV ultra-soft pseudopotentials [34], and the all-electron FLEUR [35] code.In Quantum Espresso calculations, we used the cutoffs 50 Ry and 550 Ry when expanding wave functions and charge density in plane waves, respectively.
In FLEUR-based calculations, the wave function expansion cut-off in the interstitial region was set to k max = 4.2 a.u.−1 .The muffin-tin radius of V, Cl, and O atoms were set to 2.28, 2.13, and 1.29 a.u., respectively.We have included the 3s and 3p V-orbitals as semicore states.
For the exchange-correlation energy functional, we have employed the Perdew-Burke-Ernzerhof (PBE) parametrization of the generalized gradient approximation (GGA) [36] with the on-site Coulomb repulsion (DFT+U) [37,38] applied to the V-3d orbitals.In FLEUR calculations, the on-site Hund's exchange Jparameter was set to J = 1 eV [39], and the on-site Coulomb repulsion, U , was varied.In Quantum Espresso calculations we used the Dudarev parametrization, which requires only the on-site effective Coulomb repulsion, U eff = U − J [40].
Magnetic interactions were obtained by fitting a model Hamiltonian to total energy calculations for various magnetic configurations, as described in the Supplemental Material [15].To simulate an isolated bilayer, we increased the c lattice parameter to over 30 Å.For primitive cell calculations (6 atoms), we used a 20×20×1 optimized Monkhorst-Pack [41] k-mesh.The AFM structures require a 2×2×1 cell, and we used a 10×10×1 Monkhorst-Pack k-point mesh.
Monte Carlo simulations were performed for a simulation cell containing 12800 spins, using the replica exchange method [42].We performed 2 × 10 6 steps for each spin at each temperature.To reduce correlation between successive data, statistics were collected every 10 Monte Carlo steps.Figures of the crystal structures were created with the VESTA software [43].

III. RESULTS AND DISCUSSION
A. Crystal structure and magnetic order Using the DFT+U method we have optimized the crystal structure while adopting the AFM order previously established for the bulk (see Fig. 1), for various values of the parameter U eff = U −J.As a first step, we constrained the lattice symmetry to orthorhombic, which yields the lattice constants in Table I.These values are in agreement with previous calculations [24,25].Ref. 23 reported similar lattice constants for FM ordering.Lifting the orthorhombic symmetry constraint of the unit cell we find a monoclinic distortion of the crystal lattice for all considered values of the U eff -parameter.This is in agreement with the recent experimental results by Villalpando et al. [27], who reported a monoclinic lattice symmetry.The distortion is due to magnetoelastic coupling and is induced by the two-fold AFM superstructure, which breaks the translational symmetry of FM order.It is completely analogous to what is seen in the bulk; the V ↑ -V ↓ distance is decreased at the expense of the V ↑ -V ↑ distance.
Table II accounts for the monoclinic angle γ and the lattice parameters obtained with various U eff -values.A larger U eff -value will reduce the monoclinic angle, while expanding the lattice.In the bulk, the value U eff = 2 eV has been shown to simultaneously reproduce structural, electronic and magnetic properties reasonably well [20,26].We have calculated the U eff -parameter with density functional perturbation theory (DFPT) [44], which resulted in the value U eff =5.67 eV for both monolayers as well as the bulk.Consequently, we have taken U eff = 5 eV as an upper limit while considering U eff = 2.0 eV a reasonable value.
For U eff = 2 eV, the distortion lowers total energy by 2.2 meV / atom, and the difference in V ↑ -V ↓ and V ↑ -V ↑ distances is 0.019 Å.The values of a and b more or less the same, and they are within ∼ 0.02 Å of what was previously found for bulk VOCl for the same U effvalue [26].These values may in turn be compared to the experimental a = 3.30 Å and b = 3.78 Å reported in Ref. 22 for single crystal results at room temperature.The local V spin magnetic moment is 1.5µ B and the orbital moment is −0.079µB , which is quite insensitive to the particular choice of U eff .We find the magnetic easy axis to be along b, which is agreement with the orthorhombic structure, [16,22,23,25], and is also analogous to the bulk [26].
These results clearly demonstrate that the magnetoelastic properties seen in bulk VOCl carry over to isolated single-layers as well.All calculations indicate a monoclinic ground state, which is induced by the AFM magnetic order.An orthorhombic lattice symmetry would indicate magnetic disorder.
Nevertheless, by expanding the lattice we may recover an AFM orthorhombic structure.Fig. 2 shows the resulting γ-angle as a function of the a lattice constant.At each point, the basis coordinates and the b/a-ratio has been optimized with U eff = 2 eV.Above a = 3.37 Å (b = 3.85 Å) the lattice symmetry abruptly changes from monoclinic to orthorhombic, with γ = 90.0• .As the interionic distances are increased, the magnetic interactions are diminished until there is no elastic energy gain in the distortion, whereupon the lattice changes its symmetry accordingly.Compressing the lattice has the effect of slightly decreasing the γ-angle, but no transition is seen in the examined range.
Before discussing further details of the magnetic interactions in Sec.III C, we will outline how the electronic structure of the single-layers compare to the bulk in Sec.III B.
B. Electronic structure Fig. 3(a) shows the total density of states (DOS) for a VOCl single-layer compared with that of the bulk, calculated with U eff = 2 eV and the AFM order shown in Fig. 1.The electronic structure is very similar in the two cases, with an insulating gap of 1.2 eV.The similarity underlines the two-dimensional aspects of bulk VOCl, as the single-layers are seen to be quite independent.Consequently, VOCl single-layers can be expected to retain the electronic and magnetic properties of the bulk.In Fig. 3(b) we show the site-projected DOS, indicat-ing that the V-3d electrons dominate the valence states.These are in turn separated by a gap of 1.05 eV from a manifold of Cl and O states hybridizing with a single V electron.
The character of the valence V states are seen in Fig. 4(a) to be of d zx and d x 2 −y 2 character.The lowest unoccupied orbital is of d zy character.Thus, the degeneracy of the 3d-levels is completely lifted by the crystal field of the strongly distorted VO 4 Cl 2 octahedron.Referring to Fig. 1(d), the d x 2 −y 2 lobes are directed along the a and c axes, between the V-O and V-Cl bonds of the ac-plane.The d zx lobes would be most pronounced in the ab-plane pointing towards V next-nearest neighbors.
The band structure, shown in Fig. 4(b), reveals several indirect band gaps of approximately the same size.However, this is highly dependent on the U eff -value.At the large value of U eff = 5 eV, shown in Fig. 4(c) and 4(d), the conduction bands hybridize with the high-binding energy manifold, and the delicate balance between the top and bottom of the conduction and valence bands will shift.In Ref. [24] it was reported an indirect Γ-X gap, which we cannot reproduce.
Having established the crystal symmetry, magnetic order and electronic structure, we will detail the magnetic interactions in the following section.

C. Magnetic interactions
We have derived the coefficients of the following magnetic Hamiltonian [15]: where the unit vector Ŝi denotes the magnetic spin at site i.J ij and B are the bilinear exchange and the nearest neighbor (nn) biquadratic [47] exchange couplings, respectively.∆ denotes the single-ion anisotropy, which is responsible for aligning the spins along the easy axis, d.
As the monoclinic distortion removes inversion symmetry, the nearest neighbor DM interaction, D, may be nonzero.According to the Moriya rules [30], D should lie in the ab-plane.Our calculations indeed show that D is directed along the easy axis, b.For orthorhombic symmetry, we recover D = 0.
The bilinear Heisenberg exchange interactions, J ij , deserve particular attention and will be discussed in detail in Sec.III C 1 below.In Sec.III C 2 we will present results from Monte Carlo simulations based on the Hamiltonian (1) and discuss the impact of higher order magnetic interactions.

Bilinear exchange interactions
In Fig. 5(a) we highlight the most relevant Heisenberg exchange-couplings on the VOCl lattice.We denote the nearest neighbor interaction by J 1 and the second nearest neighbor J 1 .In the orthorhombic case, J 1 and J 1 are equivalent by symmetry as discussed in Sec.III A. This means that the spin configuration in Fig. 5(b) is degenerate with that of Fig. 5(a).For monoclinic symmetry, this degeneracy has been lifted, and the configuration in Fig. 5(a) is the ground state.
The J 2 and J 3 interactions act between ions separated by a and b (note that a > b), while the J 4 and J 5 act along ±a ± b.
The obtained magnetic interactions are summarized in Table III, which shows that J 4 and J 5 are negligible, and will therefore not be discussed further.The nearest neighbor interactions, J 1 and J 1 are both seen to favor FM order, while the J 2 and J 3 interactions are AFM, regardless of the U eff -parameter.This means that it is not primarily the nearest neighbor interactions which are responsible for the magnetic ordering, as anticipated in the literature [16], but the more distant J 2 and J 3 interactions.It also means that the the monoclinic phase is frustrated: each V ↑/↓ ion has two V ↑ and two V ↓ nearest neighbors, although all four nearest neighbor exchange interactions actually favor FM alignment.
The nearest-neigbor exchange interactions in bulk VOCl have previously been discussed in terms of a combination of direct exchange mediated by the V d zx orbitals and superexchange involving d x 2 −y 2 electrons [51].Indeed, the former indeed point approximately along a line connecting the V-V nearest neighbors, while the latter involves two V-O-V paths.The sign of the TN is the Néel temperature obtained from Monte Carlo simulations [48][49][50].J 1 and J 1 interactions are consistent with the Goodenough-Kanamori-Anderson (GKA) rule [52] for superexchange, applied to the two V-O-V paths: the bond angles are close to 90 • (99.5 • and 100.4 • ) with the same total bond length, which would favor FM ordering.The V-O-V bond angle along b is 147.5 • , closer to 180 • which would favor AFM order for the J 2 interaction, as observed.As noted in Sec.III B the lobes of the d x 2 −y 2 orbitals point along a and b, which indeed would support the V-O-V hopping path.However, the J 3 interaction along a is mediated by a V-Cl-V bond, which forms a 97.3 • angle, together with the V-O-V 103 • bond, yet the J 3 interaction is AFM, contradicting the GKA rule.
It remains an open question how the superexchange mechanism would work in polyvalent materials, such as VOCl, although attempts have been made to construct a theory for CrOCl and FeOCl [53].Although direct overlap may seem unlikely, it cannot be ruled out that the exchange interactions are mediated by a combination of direct and indirect exchange.Most likely, there is a competition between Pauli exchange, Hund's coupling, and dynamical electron correlation [54,55].
Varying the U eff -parameter, we also find that interactions are reduced.The nearest neighbor interaction J 1 is always smaller than J 1 , although it has a shorter V-V distance and also corresponds to the smaller 99.5 • V-O-V bond.It seems as if the forced AFM order between the V-V nearest neighbors leads to a reduction of the FM exchange interactions.
However, the ratio between the AFM J 2,3 and the FM J 1 /J 1 will vary with U eff .In particular, the J 2 interaction which connects V -spins along a is affected most strongly.For U eff = 1 eV, the J 2,3 interactions dominate J 1 and J 1 , and J 2 is by far the strongest.For U eff = 5 eV, J 2 is instead the weakest interaction, and J 3 is comparable to the J 1 /J 1 interactions.
In the orbital-resolved DOS of Fig. 4, the t 2g and e g orbitals are seen to be well separated from the high energy manifold for U eff = 2 eV.But for U eff = 5 eV, the hybridization of these orbitals is significant.The band gap has also been effectively doubled, which reduces the hopping tendency to the unoccupied d zy states along V-V bonds in the bc-plane.The dominating effect seems to be the latter, which reduces the hopping of the largeangle V-O-V bond along b responsible for J 2 .
It is interesting to compare the exchange interactions on the monoclinic lattice with those of the orthorhombic lattice.Our calculations [15] (for U eff = 2 eV) yield a smaller FM J 1 parameter of 1.27 meV, with comparable AFM J 2,3 .These results are in qualitative agreement with Ref. [25], although those interactions were derived from a smaller set of magnetic configurations, and the implications of the results were never discussed.Allowing the lattice to we thus observe an increase in the FM nearest neighbor interactions.
Glawion et al. [51] considered exchange interactions for the orthorhombic bulk system and also reported AFM interactions along a and b, but found the sign of J 1 to depend on the assumed U eff -value.We do not see this effect in the single-layers, which may be due to an assumed FM state in Ref. 51.In any case, all theoretical work agrees on competing FM and AFM interactions in the VOCl system, giving rise to frustration.

Spin texture
As a test of the calculated magnetic interactions, we have performed Heisenberg Monte Carlo simulations on the monoclinic lattice.We define an AFM order parameter as m = 1 N N i=1 Ŝi • di , where di = ± b is the ideal direction of the spin at site i, and N is the total number of spins.m = 1 thus corresponds to the AFM ground state and m = 0 indicates complete disorder.Fig. 6(a) shows the order parameter as a function of temperature for various values of U eff , which reaches m = 0.98 at T = 0.5 K.The magnetic heat capacity is plotted in Fig. 6(b) and is seen to reach a finite value in the T → 0 limit.The Néel temperature, T N , is taken from the divergence of the heat capacity, and is listed in Table III.For U eff = 2 eV, we obtain T N = 70 K, which is comparable to the experimental value of 80 K for bulk VOCl [16,56].
T N is mostly determined by the size of the isotropic J ij -parameters.However, the biquadratic exchange interaction, B, is comparable with the J ij -values and gives a non-negligible contribution to T N .Neglecting the biquadratic exchange would lower T N by 13 K for U eff = 2 eV.
The DM-interaction parameter, D, is seen in Table III to be quite small and only weakly dependent on U eff .It does not give any appreciable contribution to T N .The negative values of D indicate that the DM interaction tends to make the spins collinear to each other [57,58].Fig. 7(a) shows a snapshot of the spin texture at 0.5 K on the monoclinic lattice.Although the correct AFM ground state is reached, we observe slight deviations from collinearity.Averaging the deviation angle, θ, of the local spin moments from the global quantization axis over all atoms in 100 individual simulation cells, we find θ = 2.8 • .The corresponding distribution is shown in Specific Heat (a.u.)  .Apart from J ij , the most important term for the alignment of the spins seems to be the biquadratic interaction, B. Removing the biquadratic term by setting B = 0 leads to θ = 3.9 • with a larger variation of θ.
Ref. [48] reported the much smaller value of T N = 23 K from Monte Carlo simulations on the orthorhombic lattice.In the orthorhombic case we find that the system jumps between the two degenerate magnetic solutions below T N (see Fig. 5) leading to a non-monotonic dependence of the order parameter with temperature.This is interesting, as the monoclinic angle is temperature dependent, pointing at the importance of the spin-lattice coupling, that ideally should be taken into account.However, this is beyond the scope of this study, which targets the magnetic ground state.

IV. SUMMARY AND CONCLUSIONS
In summary, using DFT+U calculations we have determined structural and magnetic properties of the sin-gle VOCl bilayer.Our PBE+U calculations show that the system undergoes the same monoclinic distortion as previously observed in bulk VOCl [16,17,26].The monoclinic AFM phase is magnetically frustrated, as the nearest neighbor interactions are all FM, and the observed AFM order is in fact enforced by longer ranged AFM interactions.Thus, the monoclinic distortion does not remove the magnetic frustration.These conclusions are independent of the particular value of the U eff -parameter and are in line with recent experimental reports of a monolinic lattice symmetry [27].
Together with our calculations of the electronic structure, we conclude that the physical properties of the individual layers of bulk VOCl carry over to single-layers.Nevertheless, it should be remembered that the layers of bulk VOCl are not completely magnetically independent, as they do form a well ordered two-fold magnetic superstructure along c as well.
By means of Monte Carlo simulations we have calculated the Néel temperature, which will depend on the U eff -value.With U eff = 2 eV we obtain results in good agreement with the experimental transition temperature for the bulk.In addition, our results underline the importance of higher-order exchange-interactions, such as biquadratic exchange, in line with previous theoretical predictions for layered vdW materials [47].Nevertheless, the spin-phonon coupling is most likely more pronounced in the single layer systems [59] and our calculations also do not include the contribution of low-energy excitations, such as magnons.
We hope that our results can aid in the interpretation of future experiments on atomically thin VOCl layers, as well as the other members of the M OCl family, which also become distorted at low temperature, and which currently receive increasing attention [53-55, 60, 61].

FIG. 1 .
FIG. 1. (a)-(c) Geometry of a single VOCl bilayer viewed along the c-, a-, and b-axes.The green, blue and red spheres denote Cl, V, and O, respectively.The AFM magnetic order corresponding to bulk VOCl is indicated.(d) The distorted VO4Cl2 octahedron and the local (x, y, z) coordinate system, where x = −â, ŷ = ĉ, and ẑ = b.

FIG. 4 .
FIG. 4. Orbital-resolved partial DOS (PDOS) and band structure obtained with U eff =2 eV and U eff =5 eV.In (b) and (d), the indirect band gaps are shown by dotted lines.The path are selected based on crystallography [45, 46].

FIG. 5 .
FIG. 5. (a) Exchange couplings in a VOCl bilayer.The J1 and J 1 interactions act between V ↑ -V ↓ and V ↑ -V ↑ ions, respectively.The J2 and J3 interactions are between nearest neighbors along a and b.The magnetic configuration in (a) and (b) are degenerate for orthorhombic symmetry, where J1 = J 1 .A monoclinic distortion of the lattice reduces the J 1 distance.
FIG. 6.(a) AFM order parameter, m, and (b) magnetic (specific) heat capacity as a function of temperature, T , for different U eff parameters.

Fig. 7 (
Fig. 7(b).Apart from J ij , the most important term for the alignment of the spins seems to be the biquadratic interaction, B. Removing the biquadratic term by setting B = 0 leads to θ = 3.9 • with a larger variation of θ.

FIG. 7 .
FIG. 7. (a) Spin snapshot from Monte Carlo simulations at T = 0.5 K.The black rectangle indicates a 2 × 2 magnetic unit cell.Yellow lines highlight nearest-neighbors and blue lines next-nearest neighbors.(b) Histogram of the angle between the spins and the global quantization axis, θ, obtained with the full Hamiltonian of Eq. (1) (blue), and setting the biquadratic term B = 0 (red).

TABLE I .
Optimized lattice constants for a single VOCl layer with enforced orthorhombic symmetry, obtained with various U eff -values, along with theoretical literature values.

TABLE II .
Optimized lattice constants and monoclinic angle, γ, for a VOCl single-layer, obtained with various U eff -values.