The Hard Ferromagnetism in FePS3 Induced by Non‐Magnetic Molecular Intercalation

Abstract

Manipulating the magnetic ground states of 2D magnets is a focal point of recent research efforts. Various methods have demonstrated efficacy in modulating the magnetic properties inherent to van der Waals (vdW) magnetic systems. Herein, the emergence of robust anisotropic ferromagnetism within antiferromagnetic FePS₃ is unveiled via intercalation with non‐magnetic pyridinium ions. A one‐step ion exchange reaction facilitates the formation of energetically favorable B‐phase and metastable P‐phase. Notably, both B‐ and P‐phases manifest hard ferromagnetic behavior, featuring substantial unsaturated coercive fields (>7 T) and high Curie temperatures (72–87 K). First‐principles calculations elucidate the pivotal role of electron transfer from pyridinium ions to FePS₃ in engineering magnetic exchange interactions. Calculated effective spin Hamiltonian corroborates the observed hard ferromagnetism in intercalated FePS₃. This study offers crucial insights into hard magnetism in intercalated vdW materials, thereby presenting promising avenues for 2D vdW magnet‐based magnetic devices.

Introduction

Van der Waals (vdW) crystals have been vastly investigated in the past decades owing to their exotic electrical, optical, magnetic, catalytic, energy storage, and topological properties.^[1–9] Recently, much attention has been given to the vdW/2D magnets after the first observation of ferromagnetism (FM) in atomically thin CrI₃,^[3] leading to the discovery of a wealth of 2D magnets. However, the earlier attempts toward the isolation of 2D magnets were carried out in antiferromagnetic (AFM) transition metal phosphorous trichalcogenides MPX₃ (M = Mn, Fe, Ni, X = S, Se).^[10–13] These material systems exhibited different magnetic properties, though mostly AFM, such as Heisenberg-type AFM in MnPS₃, Ising-type AFM in FePS₃, and XY- or XXZ-type AFM in NiPS_3.^[11,14–18]

Although the magnetic ground states of MPX₃ or any other vdW magnets are determined by the spin Hamiltonian, their magnetic properties can be effectively modified by the electric field,^[19] strain,^[20] electrostatic doping,^[21] and intercalation.^[22] Significantly, the weak vdW forces between covalently bonded 2D magnetic sheets facilitate the insertion of foreign species, including atoms, ions, and molecules, into the vdW gap without disturbing the planar structure yet strongly altering their properties.^[23] Previously, phosphorous trichalcogenides have been subjected to intercalation by ions and molecules, and only powder FePS₃ has indicated the AFM to FM transition after intercalation.^[24–27] Since the vdW FePS₃ powder grains are randomly oriented in all possible directions, the magnetic properties and the origins of ferromagnetism remained mostly unexplored and unanswered.

In this work, we demonstrate the emergence of strong ferromagnetism in pyridinium ion (PyH⁺) intercalated single crystal FePS₃. The one-step ion exchange reaction was utilized to introduce a selected concentration of PyH⁺ molecules into the vdW FePS₃ under different PyH⁺ concentrations, temperatures, and intercalation durations. Results show that two distinct intercalated phases with planar PyH⁺ molecule aligning to the ab-plane of FePS₃ in parallel (Basal or B-phase) and perpendicular (P-phase) geometry, each associated with unique ferromagnetic behaviors, can be realized. The interaction between PyH⁺ and the FePS₃ slab emerges as a crucial factor influencing the emergence of hard ferromagnetism. Our results offer fundamental insights into the engineering of magnetism in vdW FePS₃ systems.

Results

vdW FePS₃ crystalizes into a monoclinic structure with the space group of C2/m. Here, the Fe atoms form a honeycomb lattice in the ab-plane and are surrounded by bipyramidal (P₂S₆)4⁻, where two P atoms (P-P dimer) are covalently bonded to six S atoms (Figure 1a). Each Fe atom is bonded to six S atoms while each P is bonded to three S atoms and one P atom. Within each FePS₃ layer, two sets of FM coupled zigzag Fe chains, passing through the long Fe-Fe bond, antiferromagnetically interacted with each other, making FePS₃ an Ising-type antiferromagnet with Néel temperature (T_N) of ≈120 K.^[28] The broad range of the electromagnetic spectrum (from infrared to ultraviolet) and 1.23 eV bandgap of FePS₃ enable the application of ultraviolet photodetectors and nonlinear optics.^[29,30] Previously, the spontaneous magnetization was reported in powder FePS₃ by intercalation with organics in the range of 77–90 K. The magnetic field dependent magnetization (M–H) curves saturated at temperature slightly lower than the critical temperature and a λ-like zero-field-cooled (ZFC) and field-cooled (FC) magnetization was observed. However, the dramatically weakened saturation of M–H curve when the temperature was far lower than the critical temperature puzzles the real FM ordering in the intercalated FePS₃.

[IMAGE OMITTED. SEE PDF]

Structural Understanding of PyH⁺ Intercalation

Single-crystal FePS₃ is immersed in pyridine hydrochloride (PyH⁺Cl⁻) dissolved in ethanol to facilitate intercalation of PyH⁺ molecules into FePS₃, and fully (all the vdW gap of FePS₃ is filled by molecule) PyH⁺-intercalated single crystal FePS₃ was obtained after 50 h soaking (Figure 1a).^[24] The intercalation process is a direct ion exchange reaction.^[24] Fe²⁺ ions are released into the solution during the process. The elemental analysis (Table S1, Supporting Information) identified the loss of Fe²⁺ ions and the chemical formula as Fe_0.85P_1.06S₃(PyH⁺)_0.30.

After the intercalation, X-ray diffraction (XRD) patterns show a large shift in the (0 0 l) peaks as depicted in Figure 1b. More specifically, the out-of-plane (OP) lattice spacing (d-spacing) increased from 6.45 to 9.60 Å, which marks ≈42% increase in their vdW gap when PyH⁺ molecules are parallelly aligned to the basal plane of FePS₃,^[31] which is referred as B-phase (basal phase). In this geometry, the in-plane (IP) peak (0 6 0) shifts to the lower angle after PyH⁺ intercalation (Figure 1c). The ϕ-scan at peak (0 6 0) along the normal of ab-plane in Figure 1d revealed the same d-spacing of {0 6 0}, $$\{ {3\ 3\ \bar{1}} \}$$, and $$\{ {\bar{3}\ 3\ 1} \}$$ and 60° angle between each of them. After the intercalation with PyH⁺, the ab-plane symmetry doesn't change as shown in Figure 1e, which means the guest molecule only intercalated into the vdW gap of FePS₃ slabs without causing structural phase transition in FePS₃ sheets. The crystal parameters for pristine FePS₃ and B-phase can be obtained as a = 5.9458 Å, b = 10.3023 Å, c = 6.7526 Å, b = 107.028° and a = 6.0308 Å, b = 10.2943 Å, c = 10.3835 Å, b = 109.1216° according the XRD patterns. Raman spectra in Figure 1f at lower frequency exhibit the softening of the P-P and [P₂S₆]⁴⁻ unit vibrations.^[32] The stand-alone PyH⁺Cl⁻ displays characterized Raman peaks from 1150 to 1700 cm⁻¹ in Figure 1g as previously reported.^[33] In comparison, no Raman active modes can be detected in the same spectral range in FePS₃. After intercalation, a broad hump ranging from 1250 to 1700 cm⁻¹ emerges, implying the great suppression of the active modes of PyH⁺ due to the constraint from the FePS₃ sheets, which evidences the successful PyH⁺ intercalation.

Magnetic Properties of B-Phase and Emergence of Strong FM

While the ZFC and FC curves of the pristine (unintercalated) FePS₃ display a typical AFM ordering behavior with the Néel temperature of T_N = 120.8 K (Figure 2a), similar measurements on B-phase show vastly different magnetic response: the B-phase sample shows a typical bifurcation of FM ordering with an external magnetic field (µ₀H_⊥ = 100 Oe) perpendicular to the ab-plane (OP) of FePS₃ in Figure 2b. The FC curve quickly saturated when the temperature decreased below 72 K where the T_C value can be determined by the first derivative of ZFC curve as marked by the blue dash line. To further confirm the long-range FM ordering in B-phase, the hysteresis measurements were carried out between ±7 T at different temperatures. The square hysteresis loop was observed in B-phase comparing with the flat magnetic response of pristine FePS₃ (Figure 2c). These results show the presence of FM hysteresis below T_C and interestingly the saturation magnetization values are on the same order of magnitude as hard magnets (Figure 2e). The squareness (SQR), the ratio of remanent, and saturation magnetization (M_r/M_s) of the hysteresis is >0.82 below 40 K and 1 below 15 K (Figure S4b, Supporting Information). As plotted in Figure 2e, the coercive field (H_c) of B-phase sample developed very quickly upon decreasing temperature and exceeded 7 T (the limit of our system) below 15 K. At 15 K, the external magnetic field sweeping from 7 T to −7 T could not flip all the magnetic domains in the sample even at -7 T. When the external magnetic field swept back to 7 T, a minor loop was observed in the second panel of Figure 2e. When temperature reaches 2 K, no observable hysteresis loop can be detected because the H_c is far >7 T. The M_r and H_c in the OP magnetic field extracted from the hysteresis loop scan at different temperatures (Figure S4a, Supporting Information) are plotted in Figure 2d as a function of temperature. The H_c⊥ increases exponentially as temperature decreases, although the M_r tends to be saturated below 40 K. We also measured the magnetic response of the B-phase sample with the external magnetic field parallel to the ab-plane. The data under IP field show the FM ordering below T_C but it was much harder to align the magnetic moments in B-phase sample to IP (Figure 2f; Figure S5, Supporting Information), suggesting the B-phase sample is magnetically anisotropic with the easy axis aligning to OP direction.

[IMAGE OMITTED. SEE PDF]

It is known that a new Raman magnon mode with a frequency of ≈3.7 THz (122 cm⁻¹) emerges when pristine FePS₃ enters AFM state.^[34,35] Low temperature Raman measurements on B-phase samples show no observable peaks ≈122 cm⁻¹ as shown in Figure S6 (Supporting Information), which elaborates the vanishing of AFM ordering and the formation of FM ordering.

Emergence of Metastable Phases

A metastable phase emerges during the initial stages of B-phase synthesis.^[31] After optimizing the synthesis method, this metastable phase can be obtained by lowering the concentration of PyH⁺Cl⁻ in ethanol (see Experimental Section). This phase, referred to as P-phase, demonstrates a distinct arrangement of guest PyH⁺ molecules perpendicular to the ab-plane of FePS₃. Elemental analysis determines the formula as Fe_0.79P_0.99S₃(PyH⁺)_0.42 (Table S1, Supporting Information). In comparison to phase-B, the XRD peaks in P-phase appear at considerably lower 2θ values, resulting in a d-spacing of 12.0040 Å (2θ = 7.3585(6)°) immediately after extraction from the solution. The lattice parameters can be obtained from the XRD patterns (Figure S3, Supporting Information) as a = 6.0661 Å, b = 10.3567 Å, c = 12.6106 Å, b = 107.9105 °.

Notably, P-phase gradually transitions into a mixture of P-phase and B-phase, eventually transforming into B’-phase (2θ = 9.05(1)°) over a 10-day period, regardless of the conditions (vacuum, inert, or low temperature) to which it is exposed (Figure 3a). Here, the finally relaxed B’-phase is differentiated from B-phase since its magnetic behavior is slightly different from B-phase as discussed later in the article.

[IMAGE OMITTED. SEE PDF]

FC M-T curves in Figure 3b and Figure S7a,d,g (Supporting Information) show the differences in magnetic properties among P-phase, B’-phase, and the intermediate state. The P-phase exhibits typical FM behavior with a T_c of 87 K which is higher than that of B-phase (72 K). Similar to B-phase, however, H_c increases as the temperature is lowered, exceeding our system limit (7 T) below 10 K, as depicted in Figure 3d (lower panel). At the same range, M_r shows an upward trend below T_c, reaching a peak at 30 K before declining.

The influence of PyH⁺ structural evolution on magnetic properties is evident in Figure 3b. Following the initial VSM measurement, the sample was stored in a glovebox, and a subsequent measurement was conducted 4 days after its removal from the solution. The XRD pattern indicates the coexistence of P- and B-phase, suggesting partial relaxation from perpendicular to parallel orientation (Figure 3a). In this regime, both M–T and M–H curves (Figure 3b,c) exhibit a significant enhancement in total magnetization in the mixed phase, while the T_c, M_r-T, and H_c-T (Figure 3d) curves remain similar. After 9 days, P-phase fully transforms into B’-phase, and the magnetization remains unchanged. The IP magnetic data of P-phase (Figure S8, Supporting Information) still reveals an easy axis along the normal direction of the ab-plane.

The B’-phase transforming from P-phase can be treated as degraded B-phase directly obtained from the direct ion exchange reaction. The XRD patterns in Figures 1b and 3a show larger full width at half maximum (FWHM) in B’-phase compared with B-phase, which indicates a less alignment of (0 0 l) planes after relaxing in inert gas. During the relaxation of P- to B- phase, the solution provides a good support in the vdW gap and lets the PyH⁺ completely react with FePS₃. When the relaxation happens in inert gas, the FePS₃ sheets collapse randomly without any support and result in a much-disordered structure as observed by the broadening of (0 0 l) peaks. The less perfect crystallinity finally leads to a deterioration compared with the B-phase.

The Raman measurement in Figure 3e reveals an identical red shift of the peaks in P-phase as B-phase. This indicates that the softening of these Raman modes is not directly linked or sensitive to the orientation of the PyH⁺. However, the differences appear at much higher frequencies. A broad peak at 1435 cm⁻¹ was observed on the base of the broad hump in the P-phase, implying a partial relief of the constraint due to the less contact/interaction between the PyH⁺ molecule and FePS₃ in the perpendicular orientation.^[33,36] We argue that this peak may be related to the changes in the Raman scattering cross-section and thus the Raman sensitivity of the PyH⁺ molecule with respect to their molecular orientation under the polarized incoming laser beam and warrants further investigations in future studies.

Theoretical Insights

Given all phases of intercalated FePS₃ exhibiting similar magnetic behaviors, we expect that the mechanism of engineering magnetic properties by PyH⁺ intercalation should be the same and the different orientations of molecules only lead to quantitative difference. Therefore, without loss of generality, here we focus on the stable B-phase with a chemical formula of FePS₃(PyH⁺)_0.25 as an explicit example. By inserting one molecule into the conventional cell of pristine FePS₃ and fully relaxing the structure, we obtained the optimized structure as shown in Figure 4a. The calculated crystal parameters are a = 6.0591 Å, b = 10.3595 Å, c = 9.4185 Å and β = 101.0531°.

[IMAGE OMITTED. SEE PDF]

To determine the magnetic ground state of the intercalated FePS₃, we calculated the relative energies of FM and three possible AFM configurations (Figure 4b). The energies of zigzag-, Néel-, and stripy-AFM state are 7.28, 27.75, and 4.63 meV/FePS₃(PyH⁺)_0.25 higher than FM state, respectively, indicating the FM ground state and AFM to FM transition of FePS₃ after PyH⁺ intercalation. Since no significant structural distortion or additional chemical bonding formation was observed, we deduced that the AFM to FM transition could be induced by charge transfer between the inserted molecules and the host matrix. Bader charge analysis^[37] illustrates that one molecule would denote ≈0.74 electrons to the FePS₃ framework. This corresponds to doping 0.185 electrons per FePS₃ formula and the doping concentration of ≈1.18 × 10¹⁴ cm⁻², which is sufficient to induce AFM to FM transition in monolayer FePS₃ as predicted.^[38] Moreover, the electron transfer was also evidenced by our X-ray photoelectron spectroscopy (XPS) measurement (Figure S9, Supporting Information). All results point to the AFM to FM transition rooted in the electron transfer from the intercalants to the FePS₃ framework.

Another striking property in the intercalated FePS₃ is the giant coercive field (>7 T) at low temperature. Magnetic anisotropy energy (MAE) was calculated to understand this phenomenon. By varying the magnetization axis of the FM state, we calculated MAE (Figure 4c) and found that FePS₃(PyH⁺)_0.25 is an OP easy-axis magnet with a large MAE of −0.83 meV per Fe atom. To directly calculate the coercive field, we performed atomistic spin dynamics simulations on a reduced Fe-site spin lattice whose Hamiltonian was derived from the first-principles calculations. Considering the negligible interlayer coupling in intercalated FePS₃ (supplementary information F), we focus on the intralayer spin interactions and treat the system as a single layer. We adopted a simplified effective spin Hamiltonian to capture the essential physics in the intercalated FePS₃, written as: 1 $$$\begin{eqnarray}H &=& - A\sum_i S_{iz}^2 - {{J}_{{\mathrm{NN}}}}\sum_{\left\langle {i,j} \right\rangle } {{S}_i} \cdot {{S}_j} - {{J}_{2{\mathrm{NN}}}}\sum_{\left\langle {\left\langle {i,j} \right\rangle } \right\rangle } {{S}_i} \cdot {{S}_j} \nonumber\\ && - {{J}_{3NN}}\sum_{\left\langle {\left\langle {\left\langle {i,j} \right\rangle } \right\rangle } \right\rangle } {{S}_i} \cdot {{S}_j} - {{\mu }_i}\sum_i {{S}_i} \cdot B\end{eqnarray}$$$

Here S represents the normalized unit spin vector. The first term is single-ion anisotropy energy to describe the MAE and easy-axis magnetization. The second, third, and fourth terms represent the isotropic Heisenberg exchange couplings between two first, second, and third nearest-neighbored (NN, 2NN, and 3NN) localized spins, respectively, with details illustrated in inset of Figure 4d. The last term is the Zeeman energy describing the interaction of spins with an external magnetic field B (B = µ₀H is in units of Tesla). The derived parameters are A = 0.83, J_NN = 6.275, J_2NN = −1.981, and J_3NN = 2.975 meV (supplementary information F). Our first-principles calculations show that the magnetic moments of Fe sites are 3.3–3.6 µ_B, we thus treat them eligibly as classical spins and set µ_i = 4 µ_B for all spin sites. Then, the magnetic properties can be predicted by performing atomistic spin dynamics simulations based on the established spin Hamiltonian. The simulated temperature-dependent OP magnetization is shown in Figure 4d, from which T_C ≈35 K is predicted. The H_c can be predicted by tracking the magnetization dynamics of the spin system under cyclic OP magnetic field sweeping. Figure 4e shows the simulated hysteresis loops at different temperatures, indicating that H_c increases as the temperature decreases and reaches 7 T at 2 K. Our simulation results are qualitatively consistent with experimental observations, implying that the giant H_c at low temperature is an intrinsic property of the intercalated FePS₃, which is well captured by the established spin Hamiltonian. We note that the quantitative inconsistencies between experiment-measured and simulation-predicted Curie temperature and coercive field strength are expected considering the simplifications we used when constructing the material structure and spin Hamiltonian to ease the first-principles calculations and atomistic spin dynamics simulations, respectively.

Another potential explanation of the origin of strong FM is the formation of Fe²⁺ vacancies during the intercalation which can break the AFM ordering and induce FM.^[25] The employed intercalation method follows a direct ion exchange reaction, where Fe²⁺ ions are released directly into the solution. If the FM ordering indeed originates from a spin imbalance of up and down Fe²⁺ ions, the T_C should closely align with the Néel temperature of FePS₃. However, the induced T_C is significantly lower (≥35 K) than the Néel temperature values. Additionally, if the FM state originates from Fe²⁺ vacancies, it should persist even after the removal of PyH⁺. To investigate how PyH⁺ affects magnetism, intercalant PyH⁺ was intentionally removed (deintercalation) from the FePS₃ vdW gap by subjecting the sample to mild heating at 120 °C for 19 h in a nitrogen atmosphere. The XRD patterns before and after deintercalation in Figure 5a demonstrate the complete removal of PyH⁺ from P-phase. The broadening of the (0 0 1) peak of the deintercalated P-phase suggests less registered interlayer stacking, a common feature in deintercalated layered compounds. Both M–T and M–H curves in Figure 5c,d reveal the full recovery of the AFM state after deintercalation with T_N = 116 K. Notably, the Raman spectra (Figure 5b) remain unchanged despite the removal of PyH⁺ molecules, indicating that the Raman shift after intercalation is caused by the Fe²⁺ vacancies. As such, the disappearance of FM ordering upon deintercalation of PyH⁺ rules out the possibility of FM arising from Fe²⁺ vacancies.

[IMAGE OMITTED. SEE PDF]

Conclusion

In summary, we experimentally demonstrated that the strong out-of-plane ferromagnetism emerges from AFM FePS₃ single crystals after PyH⁺ intercalation. The studies have identified energetically stable and metastable configurations of PyH⁺ molecules within the FePS₃ vdW gap and established the interplay between those two configurations. Magnetic studies show that the resulting crystals exhibit almost hard-magnet like ferromagnetic ordering with T_C ≈72–87 K and the giant coercive field > 7 T. First-principles calculations and XPS analysis illustrate that the electron transfer from PyH⁺ to FePS₃ plays an important role in engineering the magnetic properties of the intercalated FePS₃. Atomistic spin dynamics simulations demonstrate the giant coercive field at low temperatures, indicating the hard magnetism an intrinsic property of the intercalated FePS₃ which is capable to be captured by our established effective spin Hamiltonian. Further deintercalation study rules out that the FM ordering originates from the imbalance of up and down spin caused by the Fe²⁺ vacancies.

Experimental Section

Crystal Growth

Single crystals of FePS₃ were synthesized using the chemical vapor transport (CVT) technique. No transport agent was used. Iron powder (Sigma–Aldrich, 99.99%), granular phosphorus (Sigma–Aldrich, 99.99%), and granular sulfur (Alfa Aesar, 99.9995%) were mixed in stoichiometric proportions in the nitrogen filled glovebox and sealed in an evacuated quartz ampoule at 10⁻⁵ Torr. The temperatures were set at 650 °C in the hot zone, and 550 °C in the cold zone. The ampoule was maintained at these temperatures for 7 days, after which it was allowed to cool to room temperature within the furnace. Following the cooling period, the crystal was removed from ampoule in the nitrogen filled glove box.

PyH⁺ Intercalation

In order to intercalate PyH⁺ into FePS₃, FePS₃ flakes (typical size: ≈3×3 mm²) were added into a solution made from ethanol (Sigma–Aldrich, 200 proof) dissolved PyHCl (Sigma–Aldrich, purum, ≥98.0% (AT)). The mixture was sealed in a vail and kept at 60 °C for certain days. The whole process was in the nitrogen filled glovebox. Then the flakes were removed from the solution, rinsed with ethanol three times, and dried naturally. As discussed in ref [31] 1 mg FePS₃ was kept in the solution of PyHCl in ethanol with the concentration of 1.5 g mL⁻¹ (0.5 g mL⁻¹) for 50 (60) hours, resulting in B-phase (P-phase).

X-Ray Diffraction

The X-ray diffraction (XRD) patterns were collected by a SmartLab high-resolution X-ray diffractometer (Rigaku, Japan) using Cu K_α radiation, λ = 1.5418 Å.

Raman Spectroscopy

Raman spectra were collected by confocal Raman microscope (Renishaw) using an excitation wavelength of 488 nm. The grading of 2400 in the reflection geometry. All measurements were collected with a 100x M APO chromatic objective with a NA of 0.24 for an acquisition time of 10 s and 5 accumulations.

Magnetic Properties Measurement

The magnetic properties were measured by a vibrating sample magnetometer (VSM) of a commercial physical property measurement system (PPMS, 2 K, 7 T, Quantum Design).

Scanning Electron Microscopy and Wavelength-Dispersive Spectrometry

The elemental analysis was conducted from JEOL JXA-8530F electron microprobe. Before the measurement was carried on samples, standard sample was used for calibration.

X-Ray Photoelectron Spectroscopy

The X-ray photoelectron spectra on samples were taken from Kratos Axis Supra+. The samples were exfoliated to reveal the fresh surface for analysis before loaded into vacuum. The pressure during the analysis was ≈10⁻⁹ Torr.

First-Principles Calculations

The first-principles calculations of FePS₃(PyH)_0.25 were performed based on density-functional theory as implemented in Vienna Ab initio Simulation Package (VASP).^[39] The projector augmented wave method^[40] was used to describe the interactions between core-valence electrons, and the local density approximation (LDA)^[41] was employed for the exchange-correlation functional. The energy cutoff for the plane wave expansion was set to 520 eV, and the Brillouin zone was sampled by the Γ-centered k-mesh with a grid density of 9 × 5 × 6, 9 × 5 × 3, and 5 × 5 × 6 for the unit cell, 1 × 1 × 2 and 2 × 1 × 1 supercell respectively. To describe the strongly correlated 3d electrons of Fe, the rotationally invariant DFT+U scheme proposed by Dudarev et al.^[42] was implemented, using an effective U value of 4 eV. The optB88-vdW method^[43] was adopted to account for the van der Waals interaction. Spin-orbit coupling effect was considered in all the calculations. Convergence criteria for energy and force were set to 10⁻⁶ eV and −0.01 eV Å⁻¹, respectively.

Atomistic Spin Dynamics Simulations

The atomistic spin dynamics simulations were performed by using the Spirit code^[44] based on the Landau-Lifshitz-Gilbert (LLG) equation.^[45] The Depondt's modified Heun method^[46] was used to numerically integrate the LLG equation over time, with the time step and Gilbert damping parameter set to 5 fs and 0.1, respectively. To simulate the out-of-plane magnetization M_Z as a function of temperature, a 50 × 30 supercell containing 6000 spin sites was adopted with periodic boundary conditions. The system was initially set to the paramagnetic state and was allowed to evolve for 10.35 ns at each considered temperature before the temperature changes. M_Z under each temperature was obtained via ensemble average over 700 spin configuration snapshots, which were collected after the system evolves for 10 ns and reaches equilibrium at that temperature. To simulate the temperature-dependent hysteresis loops, a smaller 24 × 14 supercell containing 1344 spin sites with periodic boundary conditions was adopted to save computational time. For each temperature, the simulation starts from a paramagnetic state under an out-of-plane magnetic field with the strength of 10 T. Then the strength of the magnetic field gradually decreases in steps of -1 T until it reaches −10 T. After that the strength of the magnetic field gradually increases in steps of 1 T until it reaches 9 T. This whole process simulates the magnetic field sweeping done by experiment. The system evolves for 10.35 ns before the magnetic field strength changes, and M_Z was obtained by ensemble averaging 700 spin configuration snapshots which were collected after the system evolves for 10 ns and reaches equilibrium under each magnetic field strength.

Acknowledgements

S.T. acknowledges direct support from DOE-SC0020653 (materials synthesis), NSF ECCS 2052527, DMR 2111812, and CMMI 2129412. The use of facilities within the Eyring Materials Center at Arizona State University was partly supported by NNCI-ECCS-1542160. X.L. and F.L. acknowledge financial support from the DOE-BES (No. DE-FG02-04ER46148). X.L. also acknowledges the National Science Foundation (No. 2326228). Computational resources for this work were supported by CHPC of the University of Utah and the DOE-NERSC.

Conflict of Interest

The authors declare no conflict of interest.

Author contributions

Y.O. and X.L. contributed equally to this work. All authors contributed substantially to the work. Y.O. and S.T. conceptualize and design the project. A.D. and G.G. grew the FePS₃ single crystal. J.K. and Y.O. collected the Raman spectra. A.D., M.S., and Y.O. collected the scanning electron microscopy and energy dispersive X-ray spectroscopy data. Y.O. conducted the intercalation of PyH⁺ into the FePS₃, X-ray diffraction, X-ray photoelectron spectroscopy and the magnetic measurements. X.L. and F.L. conducted the first-principles calculations and atomistic spin dynamics simulations. Y.O., X.L., F.L., and S.T. wrote the manuscript with the input from all authors.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, A. Kis, Nat. Nanotechnol. 2011, 6, 147.

O. Lopez‐Sanchez, D. Lembke, M. Kayci, A. Radenovic, A. Kis, Nat. Nanotechnol. 2013, 8, 497.

B. Huang, G. Clark, E. Navarro‐Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo‐Herrero, X. Xu, Nature 2017, 546, 270.

D. B. Trivedi, G. Turgut, Y. Qin, M. Y. Sayyad, D. Hajra, M. Howell, L. Liu, S. Yang, N. H. Patoary, H. Li, M. M. Petrić, M. Meyer, M. Kremser, M. Barbone, G. Soavi, A. V. Stier, K. Müller, S. Yang, I. S. Esqueda, H. Zhuang, J. J. Finley, S. Tongay, Adv. Mater. 2020, 32, [eLocator: 2006320].

D. Deng, K. S. Novoselov, Q. Fu, N. Zheng, Z. Tian, X. Bao, Nat. Nanotechnol. 2016, 11, 218.

M. Gibertini, M. Koperski, A. F. Morpurgo, K. S. Novoselov, Nat. Nanotechnol. 2019, 14, 408.

M. Wu, W. Zheng, X. Hu, F. Zhan, Q. He, H. Wang, Q. Zhang, L. Chen, Small 2022, 18, [eLocator: 2205101].

M. Z. Hasan, C. L. Kane, Rev. Mod. Phys. 2010, 82, 3045.

C.‐Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.‐L. Wang, Science 2013, 340, 167.

K.‐z. Du, X.‐z. Wang, Y. Liu, P. Hu, M. I. B. Utama, C. K. Gan, Q. Xiong, C. Kloc, ACS Nano 2016, 10, 1738.

J.‐U. Lee, S. Lee, J. H. Ryoo, S. Kang, T. Y. Kim, P. Kim, C.‐H. Park, J.‐G. Park, H. Cheong, Nano Lett. 2016, 16, 7433.

X. Wang, K. Du, Y. Y. Fredrik Liu, P. Hu, J. Zhang, Q. Zhang, M. H. S. Owen, X. Lu, C. K. Gan, P. Sengupta, C. Kloc, Q. Xiong, 2D Mater. 2016, 3, [eLocator: 031009].

C.‐T. Kuo, M. Neumann, K. Balamurugan, H. J. Park, S. Kang, H. W. Shiu, J. H. Kang, B. H. Hong, M. Han, T. W. Noh, J.‐G. Park, Sci. Rep. 2016, 6, [eLocator: 20904].

K. Kim, S. Y. Lim, J.‐U. Lee, S. Lee, T. Y. Kim, K. Park, G. S. Jeon, C.‐H. Park, J.‐G. Park, H. Cheong, Nat. Commun. 2019, 10, 345.

K. Kim, S. Y. Lim, J. Kim, J.‐U. Lee, S. Lee, P. Kim, K. Park, S. Son, C.‐H. Park, J.‐G. Park, H. Cheong, 2D Mater. 2019, 6, [eLocator: 041001].

M. J. Coak, D. M. Jarvis, H. Hamidov, C. R. S. Haines, P. L. Alireza, C. Liu, S. Son, I. Hwang, G. I. Lampronti, D. Daisenberger, P. Nahai‐Williamson, A. R. Wildes, S. S. Saxena, J.‐G. Park, J. Phys. Condens. Matter 2020, 32, [eLocator: 124003].

P. A. Joy, S. Vasudevan, Phys. Rev. B 1992, 46, 5425.

A. R. Wildes, V. Simonet, E. Ressouche, G. J. McIntyre, M. Avdeev, E. Suard, S. A. J. Kimber, D. Lançon, G. Pepe, B. Moubaraki, T. J. Hicks, Phys. Rev. B 2015, 92, [eLocator: 224408].

S. Jiang, J. Shan, K. F. Mak, Nat. Mater. 2018, 17, 406.

Y. Wang, C. Wang, S.‐J. Liang, Z. Ma, K. Xu, X. Liu, L. Zhang, A. S. Admasu, S.‐W. Cheong, L. Wang, M. Chen, Z. Liu, B. Cheng, W. Ji, F. Miao, Adv. Mater. 2020, 32, [eLocator: 2004533].

Y. Deng, Y. Yu, Y. Song, J. Zhang, N. Z. Wang, Z. Sun, Y. Yi, Y. Z. Wu, S. Wu, J. Zhu, J. Wang, X. H. Chen, Y. Zhang, Nature 2018, 563, 94.

N. Wang, H. Tang, M. Shi, H. Zhang, W. Zhuo, D. Liu, F. Meng, L. Ma, J. Ying, L. Zou, Z. Sun, X. Chen, J. Am. Chem. Soc. 2019, 141, [eLocator: 17166].

J. Zhou, Z. Lin, H. Ren, X. Duan, I. Shakir, Y. Huang, X. Duan, Adv. Mater. 2021, 33, [eLocator: 2004557].

A. Léaustic, J. P. Audière, D. Cointereau, R. Clément, L. Lomas, F. Varret, H. Constant‐Machado, Chem. Mater. 1996, 8, 1954.

R. Clment, A. Laustic, in Magnetism: Molecules to Materials, Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim 2003, pp. 397.

Y.‐D. Dai, Y. He, H.‐B. Huang, T. Shao, Y.‐F. Xia, Acta. Phys. Sin. 2003, 52, 3020.

L. Jun, D. Yao‐Dong, H. Yun, H. Hong‐Bo, X. Fan, X. Yuan‐Fu, Chin. Phys. Lett. 2004, 21, 386.

M. Amirabbasi, P. Kratzer, Phys. Rev. B 2023, 107, [eLocator: 024401].

Y. Gao, S. Lei, T. Kang, L. Fei, C.‐L. Mak, J. Yuan, M. Zhang, S. Li, Q. Bao, Z. Zeng, Z. Wang, H. Gu, K. Zhang, Nanotechnology 2018, 29, [eLocator: 244001].

D. Xu, Z. Guo, Y. Tu, X. Li, Y. Chen, Z. Chen, B. Tian, S. Chen, Y. Shi, Y. Li, C. Su, D. Fan, Nanophotonics 2020, 9, 4555.

A. A. El‐Meligi, Mater. Chem. Phys. 2005, 89, 253.

M. Cheng, Y.‐S. Lee, A. K. Iyer, D. G. Chica, E. K. Qian, M. A. Shehzad, R. dos Reis, M. G. Kanatzidis, V. P. Dravid, Inorg. Chem. 2021, 60, [eLocator: 17268].

K. Aramaki, M. Ohi, J. Uehara, J. Electrochem. Soc. 1992, 139, 1525.

A. McCreary, J. R. Simpson, T. T. Mai, R. D. McMichael, J. E. Douglas, N. Butch, C. Dennis, R. Valdés Aguilar, A. R. Hight Walker, Phys. Rev. B 2020, 101, [eLocator: 064416].

Y.‐J. Sun, J.‐M. Lai, S.‐M. Pang, X.‐L. Liu, P.‐H. Tan, J. Zhang, J. Phys. Chem. Lett. 2022, 13, 1533.

N. Preda, L. Mihut, M. Baibarac, I. Baltog, Acta Phys. Pol. A 2009, 116, 81.

W. Tang, E. Sanville, G. Henkelman, J. Phys. Condens. Matter 2009, 21, [eLocator: 084204].

B. L. Chittari, Y. Park, D. Lee, M. Han, A. H. MacDonald, E. Hwang, J. Jung, Phys. Rev. B 2016, 94, [eLocator: 184428].

G. Kresse, J. Furthmüller, Phys. Rev. B 1996, 54, [eLocator: 11169].

P. E. Blöchl, Phys. Rev. B 1994, 50, [eLocator: 17953].

J. P. Perdew, A. Zunger, Phys. Rev. B 1981, 23, 5048.

S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, A. P. Sutton, Phys. Rev. B 1998, 57, 1505.

J. Klimeš, D. R. Bowler, A. Michaelides, J. Phys. Condens. Matter 2010, 22, [eLocator: 022201].

G. P. Müller, M. Hoffmann, C. Dißelkamp, D. Schürhoff, S. Mavros, M. Sallermann, N. S. Kiselev, H. Jónsson, S. Blügel, Phys. Rev. B 2019, 99, [eLocator: 224414].

T. L. Gilbert, IEEE Trans. Magn. 2004, 40, 3443.

P. Depondt, F. G. Mertens, J. Phys. Condens. Matter 2009, 21, [eLocator: 336005].