Magnetic order in 2D antiferromagnets revealed by spontaneous anisotropic magnetostriction

The temperature dependent order parameter provides important information on the nature of magnetism. Using traditional methods to study this parameter in two-dimensional (2D) magnets remains difficult, however, particularly for insulating antiferromagnetic (AF) compounds. Here, we show that its temperature dependence in AF MPS3 (M(II) = Fe, Co, Ni) can be probed via the anisotropy in the resonance frequency of rectangular membranes, mediated by a combination of anisotropic magnetostriction and spontaneous staggered magnetization. Density functional calculations followed by a derived orbital-resolved magnetic exchange analysis confirm and unravel the microscopic origin of this magnetization-induced anisotropic strain. We further show that the temperature and thickness dependent order parameter allows to deduce the material’s critical exponents characterising magnetic order. Nanomechanical sensing of magnetic order thus provides a future platform to investigate 2D magnetism down to the single-layer limit.

Layered two-dimensional (2D) magnetic materials offer an emerging platform for fundamental studies of magnetism in the 2D limit.Their stackability into van der Waals heterostructures opens pathways to non-trivial magnetic phases and technological applications, including sensors, memories and spintronic logic devices [1].In addition to ferromagnetism, first observed in CrI 3 [2] and Cr 2 Ge 2 Te 6 [3], antiferromagnetism in 2D materials has also been studied in FePS 3 [4] and CrSBr [5].Antiferromagnetic (AF) materials are of particular technological interest due to their high spin-wave propagation speed and lack of macroscopic stray fields, making them strong candidates for spintronic and magnonic applications [6][7][8][9][10].
For insulating, thin AF materials, such as MPS 3 (M(II) = Fe, Co, Ni), few methods are available to study their intrinsic magnetism.Conventional techniques, such as neutron scattering, magnetization measurement by a superconducting quantum interference device (SQUID) or vibrating sample magnetometry are challenging, due to the small volumes of exfoliated 2D materials.Other methods, suited to 2D materials, require electrical conductance, the presence of specific optical modes or ferromagnetic order; they are therefore difficult to apply [1].In contrast, strain applied to 2D magnetic materials was shown to be extremely powerful, inducing magnetization reversal [11], reorientating the easy-axis [12], or reversing the exchange interaction [13].In addition, the direct coupling between strain, resonance frequency and magnetization in membranes of 2D magnets, makes nanomechanical resonance a sensitive method for studying their * These authors contributed equally.
phase transitions [14][15][16].Here, we show, guided by density functional theory (DFT), that the magnetic order parameter of MPS 3 AF membranes can be quantified through the anisotropy in their magneto-elastic response; from its temperature dependence the critical exponents are determined, and their thickness dependence is investigated.

First principles analysis of spontaneous magnetostriction in MPS3
Transition-metal phosphorus trisulphides, with general formula MPS 3 , are layered materials stacked in a monoclinic lattice with symmetry group C2/m [17], as shown in the top view of a single-layer in the paramagnetic phase, Fig. 1a, top panel.The spins of FePS 3 point outof-plane, whereas both CoPS 3 and NiPS 3 are in-plane systems with their spins preferentially aligned along the a axis.The intralayer AF order forms a zigzag configuration, as shown in bottom panel of Fig. 1a, leading to two opposite aligned magnetic sub-latices.The difference of the magnetisation between these sub-latices is the Néel vector.In bulk CoPS 3 and NiPS 3 , these layers with this staggered magnetism are stacked in a ferromagnetic (FM) fashion with Néel transition temperatures, T N , around 119 and 155 K, respectively [18,19].The interlayer magnetic interactions in FePS 3 are AF with a transition around 118 K [20].
To analyse the effect of magnetic ordering on the lattice, we performed first principles structural optimiza- tions of FePS 3 , CoPS 3 and NiPS 3 based on density functional theory (DFT).For the ground state zigzag magnetic configuration, the calculations predict a compression of the a lattice parameter with respect to the crystallographic, non-magnetic structure of 2.545% and 1.328% for the Co and Fe derivatives respectively (see Table 1).
In addition, the b axis expands by 0.402% (Co) and 0.359% (Fe).In contrast, in NiPS 3 the lattice parameters remain almost unchanged.The crystal and magnetic structures are strongly connected in these compounds, which is further corroborated by simulations of different spin configurations (see Supplementary Note 1).
The microscopic mechanism governing the spontaneous magnetostriction in these materials is studied using orbital-resolved magnetic exchange analyses based on maximally localized Wannier functions, (see Supplementary Note 1).The analysis shows that the spontaneous magnetostriction calculated in FePS 3 and CoPS 3 arises from isotropic magnetic exchange interactions between t 2g -t 2g orbitals.Specifically, for FePS 3 the main magnetic exchange channels, substantially affected by the compression of the a and expansion of the b lattice parameters, are the ones involving t 2g -t 2g interactions of FM nature.The changes in the lattice parameters result in an increase in J 1 and J 1 due to a decrease in distance between the d yz -d yz and d xz -d xz orbitals, respectively (Fig. 1c).Simultaneously, these changes cause a decrease of J 2 due to a larger separation of the d xy -d xy orbitals (Fig. 1d).This is compatible with the electron configuration of Fe 2+ (d 6 ), which has these orbitals partially filled and allows FM hopping between them (Fig. 1e,f).
This hopping effect also occurs for Co 2+ (d 7 ) although the additional electron present for Co blocks the d xyd xy pathway (Supplementary Note 1, Fig. S2).This results in a stronger effect along J 1 and J 1 for the optimized structure, maximizing FM interactions in the zigzag chain, which involve the d yz -d yz and d xz -d xz orbitals, respectively.For the Ni 2+ derivative (d 8 ), the t 2g energy levels are fully occupied (Supplementary Note 1, Fig. S3), which results in a blocking of the t 2g -t 2g magnetic super-exchange channels.This leads to an almost negligible modification in the lattice parameters of the optimized structure with respect to the crystallographic non-magnetic one.

Resonance frequency changes due to spontaneous magnetostrictive strain
The predicted anisotropic change of lattice parameters when going from the paramagnetic to the AF phase, causes compressive stress, σ a , and tensile stress, σ b , along the a axis and b axis respectively, as illustrated in Fig. 1a, bottom pannel.To quantify this anisotropy appearing at the phase transition, we use rectangular membranes, shown in Fig. 2b, to nanomechanically probe stress varia-SUPPLEMENTARY TABLE 1. CoPS3, FePS3 and NiPS3 lattice parameters of the crystallographic non-magnetic (NM) and fully optimized zigzag antiferromagnetic (AF-zigzag) configurations, as calculated by DFT (see Supplementary Note 1).

CoPS3
FePS3 NiPS3 tions, along a specific crystallographic axis [21] (see Supplementary Note 2).In the following analysis, we neglect the stress contribution from the thermal expansion of the substrate, as this is small compared to that of the MPS 3 compounds [14].
The resonance frequency of the fundamental mode of a rectangular membrane, f res , is approximately given by [22]: where ρ is the mass density, w and l are respectively the width and length of the membrane, as indicated in Fig. 2b, and σ w,l are the stresses parallel to these directions.For high-aspect-ratio membranes (w l), the mechanical resonance frequency is mostly determined by the stress along the shortest direction, σ w .
We study the resonance frequency of thin MPS 3 flakes suspended over star-shaped cavities with 30 • angular resolution, as shown in an example device in Fig. 2b.When the longest side of the cavity is aligned along a crystallographic axis (a or b) and w l, its fundamental resonance frequency (f a or f b ) is determined by the stress along the perpendicular axis (σ b or σ a ): On cavities oriented at an intermediate angle, θ, (defined with respect to the b axis), the resonance frequency is: where we have used the constitutive equations for a magnetostrictive membrane with plane stress [23], while only keeping the anisotropy in the magnetostriction coefficient, Supplementary Note 2. Here, E is the Young's modulus and ν is Poisson's ratio of the material.Moreover, we have ¯ = fab − th , with fab the residual fabrication strain and th the phononic thermal expansion induced strain variation.The magnetostrictive strain along the a and b-axes is given by ms,a,b = λ a,b L 2 , respectively (see Supplementary Note 3 for a detailed derivation of Eq. ( 3)), where λ a,b are magnetostriction coefficients and L 2 is the AF order parameter squared.
The temperature dependence of the resonance frequency comprises two contributions: one due to the phononic thermal expansion coefficient α, given by th (T ) = T T0 α( T )d T, where T 0 is a reference temperature and T the integration variable, and the magnetostrictive contribution ms,a,b (T ) = λ a,b L 2 (T ).The former contribution is a slowly varying function of T , while the latter term contains the staggered magnetization, which increases abruptly near the phase transition; it thus can be used to determine L(T ), as we will show below.We assume λ a,b to be T independent, as its temperature dependence will be negligible when compared to that of L(T ).

Nanomechanical determination of the order parameter
To quantify the anisotropy in the magnetic membranes, a laser interferometry technique is used to measure their resonance frequency as a function of temperature [24].A MPS 3 flake, suspended over holes in a patterned Si/SiO 2 chip, Fig. 2b, is placed inside a cryostat with optical access as shown in Fig. 2a.Both actuation and detection are done optically, by means of a powermodulated blue laser which opto-thermally excites the fundamental resonance, and a constant red laser which measures the change in the reflected signal resulting from the membrane's motion [14].A typical resonance is shown in Fig. 2c, along with the damped harmonic oscillator model fit defining the resonance frequency.Figure 2d shows that in CoPS 3 f a and f b exhibit a similar temperature dependence for T > T N , while opposite behaviour below the phase transition is visible, namely an increase of f a and a relative decrease of f b .This sudden change in f (T ) for the perpendicular cavities, occurring near T N , constitutes, in accordance with the DFT calculations, the central result of this work as it shows that the magnetic ordering in MPS 3 leads to anisotropic strain and thus spontaneous magnetostriction.We further note that strictly speaking, T N should be replaced by T * N which includes the effects of strain (see Supplementary Note 2).For simplicity, we here use the notation T N for the measured transition temperatures.
The anisotropic behavior of CoPS 3 in the AF state is even more evident in Fig. 2e, where f res (T ) − f res (140 K) for the different cavities of the star-shaped sample are plotted as a function of θ and temperature.The po-lar plot in Fig. 2f shows the data along the red dashed line at T = 70 K in Fig. 2e and results in a characteristic dumbbell-shape.Similar anisotropic behaviour is observed in FePS 3 as shown in Supplementary Note 4. On the contrary, for NiPS 3 negligible anisotropy is observed in the angle-resolved magnetostriction data in Fig. 2g-i   To obtain L(T ) from the data, we first subtract the pretension contribution from the resonance frequency , for each angle, where T 0 = 150 K is the highest temperature in our measurements.The resulting values of f 2 θ (T ) along the crystalline axes a and b are shown in Fig. 3a,d,g for the three MPS 3 compounds.With Eq. (3), we then calculate the We can now use Eq. ( 4) to access the critical behaviour of L below T N by plotting f 2 b − f 2 a as a function of temperature.As shown in Fig. 3b,e,h, the trend presents the typical critical behaviour with a non-zero order pa-rameter appearing in the ordered state for T < T N .Figures 3c,f,i show the same critical curve as Fig. 3b,e,h respectively, plotted on a logarithmic scale against the reduced temperature (1 − T /T N ).Note that the difference f 2 b − f 2 a for NiPS 3 , is substantially smaller than that of the Fe/CoPS 3 membranes indicative of a weaker anisotropic magnetostrictive behaviour.
The angle dependence of the resonance frequencies allows us to estimate the ratio r ab = λ a /λ b between the magnetostriction parameters, λ a,b , (see Supplementary Note 3).This ratio we directly compare to DFT calculations: Experimentally, we find for FePS 3 , r ab = −2.3± 0.3 while from the DFT calculations we estimate r ab = −3.70.For CoPS 3 (taking [25] ν CoPS3 = 0.293), the experimental value is −1.42 ± 0.07 and the DFT one −6.33.We conclude that although both the sign and order of magnitude of the magnetostrictive anisotropy in these compounds are well reproduced in the current work, more detailed studies will be needed to obtain full quantitative correspondence with theory.

Thickness dependence of critical behaviour
As follows from Landau's theory of phase transitions (see Supplementary Note 2), L(T ) near T N is given by where A and B are constants and β is a critical exponent representative of the magnetic order.We fit Eq. ( 5) to the data in Fig. 3b,e,h in the region close to T N (indicated by the black dashed line in Fig. 3b,e,h) to extract the critical exponent β and T N for the three materials (see Supplementary Note 6 for more details on the fitting procedure).In the logarithmic plot of the critical curve the fitting of a straight line shows good agreement to the data points, consistent with the result of Eq. ( 5).
The values for β and T N are plotted in Fig. 4 as a function of thickness, t, and listed in Supplementary Note 6, Table 49.
For thicker CoPS 3 samples (t = 40 − 60 nm) we find β = 0.289 ± 0.034 close to what is reported in literature for the bulk (β bulk = 0.3 ± 0.01 [18]) and consistent with the 3D Ising model.For samples with t < 10 nm the measured β, on the other hand, is 0.195 ± 0.045, closer to β 2DXY as shown in the top panel of Fig. 4.This constitutes a noticeable change in β while going from bulk to thinner samples.Similarly, we observe for CoPS 3 a decrease in T N from the bulk value of 118 K down to ∼ 100 K, similar to what was previously reported in Ref. [26].We fit a power law to the dependence of T N on thickness, where C is a non-universal constant related to the interlayer coupling, and ν eff is an effective critical exponent related to the correlation length [31].Fitting the CoPS 3 data points with T 3D N = 118 K [26] yields C = 1.43 ± 0.457 nm and ν eff = 0.84 ± 0.096.This value of ν eff is intermediate between the expected values of ν eff = 0.630 for the 3D Ising and ν eff = 1 for the 2D Ising models, and indicative of a transition regime [32].

CONCLUSIONS
In conclusion, we provide a comprehensive analysis of the anisotropic magnetostriction effect in MPS 3 compounds and its implications to the dynamics of membrane made from them.DFT calculations provide a microscopic explanation for the anisotropic lattice deformation in CoPS 3 , FePS 3 and NiPS 3 which are consistent with our measurements.We further demonstrate the relation between magnetic ordering and anisotropy in the mechanical resonance frequency of suspended MPS 3 resonators, providing a direct measure of the AF order parameter in absence of an external magnetic field.We observe a thickness dependence in the critical behaviour of CoPS 3 resonators [18,26], which is absent in the case of FePS 3 .The presented technique is of particular interest for the study of 2D magnetism given the scarcity of methods available to investigate critical phenomena of van der Waals materials in the atomically thin limit.

Sample fabrication. Substrates consist of thermal
SiO 2 of 285 nm thickness, grown on highly doped (Si ++ ) silicon.The rectangular cavities are defined via e-beam lithography using AR-P 6200 resist.After development, the exposed SiO 2 areas are fully etched via reactive ion etching.The AR-P 6200 resist is stripped in PRS-3000 and the sample is cleaned in an O 2 plasma before stamping.The exfoliation and transfer of multi-layer MPS 3 flakes is done using a polydimethylsiloxaan (PDMS) transfer method.First, MPS 3 crystals are exfoliated onto the PDMS through scotch tape.Selected flakes are then transferred on the star-shaped cavities in the SiO 2 /Si substrate.

Laser interferometry.
Samples are mounted on a heater stage which is cooled down to 4 K using a Montana Instruments Cryostation s50 cryostat with optical access.A blue diode laser (λ = 405 nm) is used to excite the membrane optothermally via AC powermodulation from a vector network analyzer (VNA) [33].Displacements are detected by focusing a red He-Ne laser beam (λ = 632 nm) on the cavity formed by the membrane and Si substrate.The reflected light, which is modulated by the position-dependent membrane motion, is recorded by a photodiode and processed by a phase-sensitive VNA.Laser spot size is ∼ 1 µm.
DFT calculations.First principles spin-polarized DFT calculations in the plane wave formalism are preformed as implemented in the Quantum ESPRESSO package [34].The exchange-correlation energy is calculated using the generalized gradient approximation using the Perdew-Burke-Ernzerhof functional [35] and standard Ultra-soft (USPP) solid-state pseudopotentials.The electronic wave functions are expanded with well-converged kinetic energy cut-offs for the wave functions (charge density) of 75 (800), 85 (800) and 85 (800) Ry for Fe, Co and Ni, respectively.The crystal structures are fully optimized using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [36] until the forces on each atom are smaller than 1 × 10 −3 Ry/au and the energy difference between two consecutive relaxation steps is less than 1 × 10 −4 Ry.In order to avoid unphysical interactions between images along the non-periodic direction, we add a vacuum of 18 Å in the z direction for the monolayer calculations.The Brillouin zone is sampled by a fine Γ-centered 5 × 5 × 1 k-point Monkhorst-Pack [37].
A tight-binding Hamiltonian derived from first-principles is constructed in the base of Maximally-localized Wannier functions, as implemented in the Wannier90 code [38].For that, we select the d orbitals of the metal centre (Fe, Co, Ni) and the s and p orbitals of P and S to construct the connected subspace.Magnetic interactions are determined using the Green's function method in the TB2J software [39].The orbital resolved analysis is performed after rotating the coordinate system of the crystal to align the metal-sulfur bonds direction of the octahedra with the cartesian axes.
Crystal growth Crystal growth of MPS 3 (M(II) = Ni, Fe, Co) is performed following a solid-state reaction inside a sealed evacuated quartz tube (pressure ∼ 5 × 10 −5 mbar).I 2 was used as a transport agent to obtain large crystals.A three-zone furnace is used, where a tube with the material was placed in the leftmost zone.This side is then heated up to 700 • C in 3 hours so that a temperature gradient of 700/650/675 • C is established.The other two zones are heated up in 24 hours from room temperature to 650 • C and kept at that temperature for one day.The temperature is kept constant for 28 days and cooled down naturally.With this process crystals with a length up to several centimeters are obtained.Detailed description of the crystal growth and characterization can be found in earlier work [14].Magnetostriction is a coupling between magnetic and mechanical parts of our system.This coupling can be described by an energy term in the total free energy of our system [1].We can then write the total free energy as Here F is the total free energy of in the AF phase at zero magnetic field, F 0 is the free energy of paramagnetic phase, U el (z) is the elastic energy of a membrane with deflection z at its centre, T is the temperature of our system and T N is the Néel temperature, L i are the components of the Néel vector, β is a critical exponent, a and B are positive constants, σ ij (z) is the stress tensor and λ ijkl is the magnetostriction tensor.The last term couples the stress to the Néel vector thereby describing the magnetostriction.If we assume the Néel vector to be fully aligned with the easy axis, Eq. ( 7) simplifies to: where L is the magnetic order parameter (i.e. the magnitude of the Néel vector).For notational convenience we write λ ij in dropping the third and fourth index of λ ijkl as only the component where kl corresponds to the easy axis contributes.The elastic energy in a homogeneous membrane is given by [2] U el = S ijkl 2 σ ij (x, y, z)σ kl (x, y, z)dxdy, where the integration runs over the in plane dimensions of the membrane, z is the membrane deflection at its centre and should not be confused with the out of plane coordinate.For ease of notation we will not explicitly write the integration and coordinate dependence from here on.Assuming the membrane thickness does not vary significantly we can take the out of plane component stress component to vanish, σ zx = σ zy = σ zz = 0. Eq. ( 9) then simplifies to Taking our coordinates such that x, y correspond with the principle stress directions all terms containing σ xy vanish.This simplifies the elastic energy further to Assuming the material has isotropic elastic properties the relevant compliance tensor components are Substituting this in to Eq. ( 11) we find By taking the derivative of the free energy with respect to either z or L we find the forces acting on these degrees of freedom, φ L and φ z respectively to be given by Order parameter and critical exponent In order to find an equation that describes the order parameter as a function of temperature, we find a solution for equation (14) for the case where φ L = 0. Aside from the trivial solution L = 0, we find for below the transition the additional solution: which could be rewritten for 1 2β [3] as: This equation now describes the temperature dependence of the order parameter in a critical region near T N with a corresponding critical exponent β.

Magnetostrictive strain and resonance frequency
To assess the magnetostriction contribution to strain and thus the frequency of a rectangular membrane resonator, we need to find stiffness of the membrane from its force-deflection equation.In doing that we analyse equation (15).First, we describe strain equation for the rectangular membrane at its centre as: where c 1 is a geometrical pre-factor that describes the deflection shape of the fundamental mode of vibration [4,5].For w l we can neglect the z dependence of xx (z).Now, we substitute ( 18) to ( 15), using the relation σ ij = C ijkl kl , we find where we used that where k 1 is the elastic linear stiffness and k 3 is the cubic elastic stiffness, given by Assuming small deflections we can neglect the z 3 contribution in eq. ( 19) and find that the linear stiffness is changed with respect to the purely elastic case.If we consider a rectangular cavity with it's long axis is parallel to the crystalline axis b or a respectively we find: where λ a,b are the phenomenological magnetostriction coefficients, chosen such that to couple a and b crystalline directions and L. That leads to a change in the effective linear stiffness k b,a : which can be used to write down the frequency equations using f where m is the mass of the membrane.And define magnetostrictive strain by: Taking a difference of the squares of equations ( 25) and assuming 0,a = 0,b , we arrive at the final equation: which relates antiferromagnetic order parameter L and measured resonance frequencies of orthogonal resonators aligned to crystalline axes f a,b in ordered phase.Finally, one can show that by plugging equation ( 17) to (28): which could be used to fit experimental data to extract critical exponent β near the phase transition temperature T * N .

Supplementary Note 3. DERIVATION OF ANISOTROPIC RESONANCE FREQUENCY
Here, we derive the general equation of the resonance frequency of a rectangular cavity oriented at an angle θ with respect to the crystalline axes, as schematically shown in Fig. 9.The global coordinate system is defined by the crystallographic axes a and b, along which the material deforms resulting in stresses σ aa and σ bb .The longest side of the cavity, with length l, can be oriented at an arbitrary angle θ with respect to b.Let us first consider a cavity oriented parallel to a crystallographic axis.Since the membranes are very thin, we can assume that the stress in the direction perpendicular to the plane is zero, σ cc = 0, the membrane's stress tensor can be expressed as where the subscript ab indicates we expressed the stress tensor in the basis of the crystallographic coordinate system.If we assume that there are no shear forces acting on the crystal lattice, σ ab = σ ba = 0, there will be no shear on cavities oriented along the main crystallographic axes.Now, if we consider a rectangular cavity rotated by θ with respect to the crystallographic axes, we can define a rotated xy-coordinate system oriented along the main axis of the rectangle.To express σ in this coordinate system we use the tensor transformation rule, σ ij = q ki q lj σ kl where q ij are components of the rotation tensor transforming the ab-coordinate system, e, into the xy-coordinate system, e as e i = q ij e j .We then get: Then, the resonance frequency of a rectangular membrane oriented at an angle θ with respect to the crystallographic axis can be expressed as In the case of high-aspect ratio membranes (w l), Eq. 32 can be approximated to which is Eq. 3 of the main text.Now, let us consider the constitutive equations of the material: where fab is residual fabrication strain at T = T 0 , th and ms are respectively the thermal expansion and magnetostriction contributions to strain, α is the thermal expansion coefficient, λ the magnetostriction coefficient and E is the Young's modulus, which is assumed to be isotropic.We can thus write σ aa = c 1 + νσ bb (36) which can be combined in the following expressions for σ aa and σ bb : We can now rewrite Eq. 33 in terms of the different contributions to strain, i.e. residual strain from fabrication ( fab ), thermal expansion (∝ α) and magnetostriction (∝ λ): Which is consistent with Eq. (25).We can eliminate the pretension fab terms by considering fθ (T ) = f θ (T ) − f θ (T 0 ).In the following, we assume that the only anisotropic temperature-dependent contribution to the total strain comes This behaviour is observed for the thicker samples (t > 10 nm) and it is exploited to have a better estimate of the critical parameters β and T N as discussed in Supplementary Note 6.For thinner resonators every irregularity, like wrinkles or tears, can strongly affect their mode shapes.In some cases, these imperfections can drastically change the resonance frequency of the fundamental mode, as well as its temperature dependence.Therefore, when analysing the critical behaviour of thin flakes, we choose only the most pristine and unaffected membranes out of all fabricated out of a single flake to be sure we are not affected by these irregularities.

1 'SUPPLEMENTARY FIG. 1 .
Magnetostriction in MPS3 membranes.a, top panel, Crystalline structure of MPS3 in the paramagnetic phase (T > TN).Black hexagons indicate the organisation of magnetic atoms in the lattice.a, bottom panel, Crystalline structure of MPS3 at the AF phase (T < TN) as it elongates in the b and contracts in the a direction.Light blue and red arrows indicate the axial lattice distortion.b, Illustration of the exchange interaction parameters included into the Heisenberg spin Hamiltonian.c-d, Calculated maximally localized Wannier orbitals.Green arrows illustrate the most relevant FM superexchange channels for J1 (J 1 ) (c) and J2 (d), corresponding with the dyz-dyz (dxz-dxz) and dxy-dxy orbitals, respectively.e-f, Electron configuration of the Fe 2+ magnetic ions connected by J1 (e) and J2 (f ), showing parallel and antiparallel spin orientations, respectively.
res (70 K) -f res (140 K) (MHz) f res (5 K) -f res (190 K) (MHz) SUPPLEMENTARY FIG. 2. Angle-resolved mechanical characterization via laser interferometry.a, Schematic illustration of the laser interferometry setup and sample with rectangular cavity array.b, Optical image of the rectangular membranes array for a CoPS3 sample.The a and b axis are determined from the resonance frequency behaviour.Scale bar: 12 µm.Schematic of 0 • and 90 • membranes from the array where w is the width of the membrane and l its length.c, Measured amplitude of the fundamental resonance peak in a CoPS3 drum at T = 10 K and Lorentzian fit used to extract the fundamental resonance frequency, fres, and quality factor, Q. d, Temperature dependence of fres of a CoPS3 rectangular membrane orientated along the a and b axes.The dashed line indicates the transition temperature TN extracted from the data.e, Resonance frequency difference, fres(T ) − fres(140 K), as a function of angle and temperature.The dashed line indicates the transition as in d f, Polar plot of fres(T ) − fres(140 K) taken along the red dashed line in (e).Panels g-i, follow the same structure as (c-e) for NiPS3 resonators with negligible anisotropy, measured between 5 K and 190 K. .

3 .
Anisotropy and critical behaviour in resonance frequency of MPS3 (M(II) = Co, Fe, Ni) membranes.a, Pretension corrected resonance frequency ( f 2 a rectangular membranes of CoPS3 oriented along the b axis (blue) and a axis (red) b, Difference of the corrected frequency squared f 2 b − f 2 a proportional to the order parameter L 2 from Eq. 4. The dashed-dotted line indicates the measured transition temperature TN.The dashed black line is a powerlaw fit through the data close to TN (see Supplementary Note 6).c, Difference of the corrected frequency squared f 2 b − f 2 a as a function of the reduced temperature 1 − T /TN.The dashed black line is the fit from b where the slope defines the critical exponent 2β.d-f, and g-i, follow the same structure as (a-d) for FePS3 and NiPS3 resonators, respectively.

4 .
Thickness dependence of critical behaviour.Average critical exponent, β, and critical temperature, TN, of MPS3 resonators plotted as a function of thickness.The blue stars indicate CoPS3 bulk values from [26].Critical parameters have been determined from power law fits to f 2 b − f 2 θ , as shown in Fig. 3b,e,h, and then taking the average value over the fit parameter for all angles θ = 0. Error bars are calculated from standard deviation of fit results for all θ.The horizontal gray dashed lines in the upper plot indicate the expected values of β for the 3D or 2D versions of the Heisenberg (H), XY or Ising (Is) models.The blue dashed line in the lowe panel indicates a fit to Eq. (6) through the CoPS3 data with ν eff = 0.84 ± 0.13.

9 .
Schematic illustration of rectangular membrane Rectangular membrane of width w and length l oriented with its long side at an angle θ with respect to the crystalline direction b.The x-y direction refer to the main directions of the rectangular membrane.

2 ( 1 −be fitted to b 1 C Debye + b 2 dL 2 dT where b 1 and b 2 2 ( 1 13 .
10. Angle-resolved df 2 dT : Plot of measured (light blue dots) df 2 dT and fit to b1C Debye + b2 dL 2 dT (blue full line) for all angles of a star-shaped array of CoPS3.from magnetostriction, thus we take th,aa = th,bb = th .By definition of θ we have b → θ = 0 • and a → θ = 90 • .From Eq. (41), we then find that f 2 a ν) − ms,aa + ms,bb(42)= − E 4ρw 2 (1 + ν) (λ a − λ b )L 2(43)from which we can directly extract the order parameter.The thermal expansion contribution to strain α is proportional to the integral over the temperature of the thermal expansion coefficient α, which is proportional to the Debye specific heat, C Debye , via the Grünesen parameter.Thus, the derivative with respect to temperature off 2 ν 2 ) (sin 2 θ + ν cos 2 θ) −α − λ a dL 2 dT + (cos 2 θ + ν sin 2 θ) −α − λ b dL 2 dT (44) = −E 4ρw 2 (1 − ν 2 ) α(1 + ν) + sin 2 θ(λ a + νλ b ) dL 2 dT + cos 2 θ(λ b + νλ a )are fit parameters, and dL 2 dT is estimated from Eq. 43.The results of these fits along to measured data of df 2 θ dT are shown in Fig. 10.Polar plot of the resulting b 1 (θ) and b 2 (θ) are shown in Fig. 11, which confirm that the thermal contribution to strain does not exhibit significant anisotropic behavior.From Eq. 42 and 40, the expected angle dependance of the parameter b 2 is b2 (θ) = E 4ρw 2 (1 − ν 2 ) − ν)(λ a − λ b )[(sin 2 θ(λ a + νλ b ) + cos 2 θ(λ b + νλ a )] ,(46)which we use to fit b 2 (θ) in Fig.11to A sin 2 θ + B cos 2 θ whereA B = λ a + νλ b λ b + νλ a .(47)derivation in Supplementary Note 2), f 2 b − f 2 θ is also proportional to L 2 , where f 2 θ is the pretension corrected resonance frequency of a rectangular cavity oriented at an angle θ with respect to the b-axis.We show this quantity for the CoPS 3 and FePS 3 star-cavity resonators in Fig.13a,dby plotting f 2 b − f 2 θ as a function of angle and temperature in Fig. 13b,e.Figure 13c,f shows the polar plot of f 2 b − f 2 θ taken along the red dashed line in 13b,e.Angle-resolved data of f 2 b − f 2 θ ∝ L 2 of CoPS3 and FePS3 membranes.a,c Optical image of the CoPS3 (a) and FePS3 (c) resonators.Scale bar 12 µm.b Resonance frequency difference, f 2 b − f 2 θ ∝ L 2 , as a function of angle θ and temperature of the CoPS3 sample in (a).c Polar plot of f 2 b − f 2 θ taken along the red dashed line in (b).e,f follows the same structure as (b,c) for the FePS3 sample in (c).