Introduction

Flexocoupling effects including flexoelectric and flexomagnetic coupling effects correspond to the response of microscopic polarized and magnetic orders to external gradient strain or stress1,2,3,4,5. The inhomogeneous deformation (i.e. gradient strain or stress) triggers the uneven distribution of bonding polarization and spin coupling, thereby altering the spin-orbital-lattice coupling and triggering exotic ferroelectric or magnetic phenomena such as effective Dzyaloshinskii-Moriya interactions (DMI)6,7. Since the prediction by Mashkevichand and Tolpygo8, flexocoupling has been experimentally confirmed1,9,10,11 and attracted wide interest in both experimental12,13 and theoretical fields14,15,16. The flexomagnetic phenomenon entails the interplay between the bending deformation and the magnetic excitations and renders additional dimensions to manipulate the spin orders.

In particular, spontaneous flexoelectric and flexomagnetic effects without the need of external electric/magnetic fields modify the magnitude of order parameters17,18,19. One effective method to assure an enhanced flexocoupling is the utilization of nanostructures like asymmetric thin films20,21,22 or strained interfaces23,24,25. The quantum confinement effect in nanostructures stabilizes the Skyrmions and is beneficial to room-temperature Skyrmion lattices without the need of an external magnetic field26,27. The presence of edges or surfaces in nanostructures, which breaks the periodicity, causes the displacement of the atoms from their equilibrium positions, thereby inducing the spatial gradient of order parameters28,29,30. Moreover, spatial discontinuity at edges facilitates the curvature-induced unconventional spin textures and anisotropy energy31,32,33,34. This has attracted great interest in generating curvilinear magnetism35,36,37,38,39. While tremendous phenomenological approaches40 have been utilized for bulk systems to underline the coupling of the intrinsic polar (polarization) and magnetic (magnetism) properties with extrinsic strain, such macroscopic methods are unable to sample the coupling near boundaries or surfaces due to the abrupt changes of electronic screening and correlation. Thus far, the effective coupling and interplay between spin re-ordering and graded deformation are unestablished.

In this study, working on the open boundary edge of emerging two-dimensional (2D) ferromagnetic CrI3, since the CrX3 family display intrinsic magnetic properties41,42 and is an ideal platform for studying 2D magnetism43. We show that bending-induced coexisting compressive and tensile strain (on the concave and convex surfaces) causes a transition from the collinear state to noncollinear state, which was proven to be unachievable by uniaxial or biaxial strain6. The dependence of exchange coupling on the structural parameters of nanoribbons of CrI3 is investigated with the maximum information coefficient statistical method. We reveal that the mechanism of superexchange in asymmetric structures adheres to the Goodenough-Kanamori rules.

Results

Pristine CrI3 monolayer and construction of CrI3 nanoribbons

Bulk CrI3 possesses an intralayer ferromagnetic coupling, and the spins of Cr atoms lie inside the a–b plane. In contrast, monolayer CrI3 entails the magnetic easy axis (EA) along the c-axis (Fig. 1a) with the hybrid magnetic interaction of short-range Heisenberg-like superexchange and the direct exchange interaction among atomic neighbors carrying the spins. Each Cr magnetic center has a highly localized magnetic moment equal to 3 µB. The calculated lattice parameters are a = b = 7.08 Å, which are mostly consistent with the previous DFT results44,45 but slightly overestimated compared with the experimental results46 and the calculation results without GGA + U47 because the Coulomb interaction U can enlarge the lattice parameters.

Fig. 1: Structure and bending deformation of a CrI3 nanoribbon.
figure 1

a Top and side views of a CrI3 monolayer. The unit cell is indicated by the black dashed quadrilateral. b Schematic pathways for spin coupling between neighboring Cr atoms. c Construction of CrI3 zigzag nanoribbons in flat and bent forms. d Geometric bending diagram, where d0 and d are the straight distances between Cr1 and Cr9 before and after bending. C corresponds to the curvature of the nanoribbon induced by bending. The Cr and I atoms are represented in blue and violet, respectively.

Full size image

The geometric information of each Cr-Cr pair and the mediated I atoms are labeled in Fig. 1b. The bond length of Cr-I is 2.77 Å, and the direct distance between the nearest neighboring Cr atoms is 4.09 Å. Two symmetrical Cr-I-Cr angles measured in the top and bottom I sub-lattices are equal to 95.03°. Due to the octahedral coordination, the 3d orbitals of the Cr atom are split into the lower triply degenerate t2g orbitals and upper doubly degenerate eg orbitals. Direct exchange interaction between neighboring Cr species leads to a weak antiferromagnetic (AFM) coupling because of distant separation and suppressed parallel spin hopping of Cr3+ (t2g3eg0) under the restriction of the Pauli principle. In addition, there exists another coupling mechanism associated with the Cr-I-Cr superexchange interaction (Fig. 1b). According to the Goodenough-Kanamori rules48,49, a near-90° bridging configuration between the two 3d orbitals of Cr3+ ions separated by I- ions favors the ferromagnetic (FM) superexchange coupling. The existence of the two competing interactions, together with the flexible Cr-I-Cr hinge structure, renders flexomagnetic response in CrI3 as shown hereafter.

To demonstrate the vast spin ordering possibilities under the bending manipulation, we construct one-dimensional nanoribbons of CrI3 where the edge discontinuity facilitates the decoupling of the spin strips from local aperiodic deformation. As shown in Fig. 1c, the CrI3 monolayer is constructed into nanoribbons with exposed zigzag edges. Flexo-bending deformation is applied along the width (a||) direction of the nanoribbon, while its axis (a⊥) direction is free of bending, making its structure one-dimensional. Six differently deformed nanoribbons associated with consecutive bending angles, named BX (X = 0, 20, 50, 80, 150 and 250), are examined, where B stands for bending deformation and the number X is the force constant K0 (kcal·mol−1Å−2) used in the Harmonic potential \(E = K_0\left( {\theta - \theta _0} \right)^2/2\) with θ being the target angle of each bent nanoribbon (Fig. 1c) and θ0 the equilibrium angle. The “B0” stands for the original flat nanoribbons without bending.

The deformation is quantitatively described by the curvature of the nanoribbon induced by bending as illustrated in Fig. 1d. The evolution of the geometric characteristic with bending is listed in Table 1. The distance between two Cr atoms (Cr1 and Cr9) shortens by ~12.2% as the curvature increases to 0.062 Å−1. The bending triggers inhomogeneous deformation across the nanoribbon. The average in-plain strain of the top and bottom I-layers are asymmetric at the concave/convex, defined as \(\varepsilon _{{{{\mathrm{top}}}}/{{{\mathrm{bottom}}}}} = \left( {\overline {d_{{{{\mathrm{I}}}} - {{{\mathrm{I}}}}}} - \overline {d_{0{{{\mathrm{I}}}} - {{{\mathrm{I}}}}}} } \right)/\overline {d_{0{{{\mathrm{I}}}} - {{{\mathrm{I}}}}}}\) through calculating the ratio of the change of the average I-I distance on the top and bottom I-layers over the I-I length of unbent structure.

Table 1 Structural information of CrI3 nanoribbons under different degrees of bending.
Full size table

Table 1 shows εtop ranging from −0.5% to −4.4%, and εbottom from 0.7% to 3.2%. The maximum local strain (−4.4% and 3.2%) occurs at a curvature of 0.062 Å−1. A previous study showed that a co-planar deformation under such an amount of homogeneous strain is not large enough to change either the magnetic ground state or the spin configuration47. In contrast, hereafter we show that the bending and asymmetric deformation with an equivalent degree of strain triggers re-ordering of spins.

Flexomagnetic effect

Figure 2a shows the bending energy density, which is calculated as the difference in the total energy between the bent and flat nanoribbons divided by the a–b plane area of the flat nanoribbon. We find that when the curvature is <~0.05 Å−1, the bending energy density of 2D CrI3 is smaller than the van der Waals interlayer binding energy density of bulk CrI3 (14.3~17.9 meV·Å−2)50. Such a small energy density associated with the bending deformation implies that bending, rippling, and folding of 2D CrI3 sheets can easily occur. Recently, bending of 2D materials has been realized by compressing or stretching their supporting piezoelectric/dielectric substrate materials through the application of an electric field for a changeable deformation51. Other methodologies such as using tips to cause out-of-plane deformation during sliding on a graphene sheet52 or draping graphene over atomically sharp steps of hexagonal boron nitride53 have also been tested through experiments.

Fig. 2: Energetics and spin configurations of CrI3 with a series of curvatures.
figure 2

a Bending energy density as a function of the curvature. The light green area is the range of the van der Waals interlayer binding energy (IBE) of bulk CrI3. b Energy differences between AFM and FM configurations. The total energy of the FM state at zero curvature, E0, is used as the energy reference. c Schematics of four main spin states where the cones are the classical precession trajectories of thermally excited spins (arrows). d Energy of different spin states as a function of curvature relative to the energy of the parallel state at zero curvature E0.

Full size image

The energetics of AFM and FM spin configurations in the CrI3 nanoribbons are examined. For the FM state, all magnetic moments point to the positive z direction (all↑), while for the AFM state, the neighboring spins have opposite orientations; i.e., one is in the positive z direction, and the neighbors are in the negative z direction (2↓1↑4↓3↑6↓5↑8↓7↑9↓10↑, with the numbers and arrows standing for different Cr sites in Fig. 1c and their associated directions of spins, respectively). As shown in Fig. 2b, the FM state is more energetically favorable than the AFM state for all bending conditions, implying a robust FM under bending. Traditional co-planar biaxial strain engineering shows that the energy difference between FM and AFM ordering decreases with increasing compressive strain44,47 resulting from a shortened Cr-Cr pair and accordingly enhanced direct exchange. In contrast, here in our work, the energy difference between the FM and AFM states is nearly constant (25 ± 1 meV per Cr) for all the bending situations. The coaxial bending barely influences the Cr-Cr distance; for instance, the most bent structure B250 has an average Cr-Cr distance of 4.051 Å, which is almost the same as that of B0 structure (4.055 Å). Therefore, the flexo-bending deformation on CrI3 is largely accompanied by varying the hinge angle of Cr-I-Cr to regulate the superexchange. This structural flexibility facilitates the manipulation of orientationally rearranged spin patterns while its intrinsic FM order is maintained.

We next reveal the detailed spatial rearrangement of the spin configurations after the flexo-bending. Four typical spin states are constructed (Fig. 2c): spins aligned along the c direction (parallel state, Ec||), the radial direction (ERadial), the in-plane (EIn-plane) direction, and the out-of-plane direction (EOut-of-plane). The in-plane magnetization keeps the magnetic moment of each Cr atom aligned co-planarly within the triangular Crmn-Cr-Crmn plane formed by the most neighboring (mn) Cr atoms, and the out-of-plane magnetic moment is set normal to this basal plane. The radial state here is in analogy to the collinear state, and the parallel state is noncollinear and similar to a bunch of plumb lines hanging on each Cr atom. More details about the definition of these spin states can be found in Supplementary Note 1.

Figure 2d shows the energy of these spin states as a function of the curvature. When the nanoribbon is flat (B0), we recover the magnetic ground state of the CrI3 monolayer which has an out-of-plane easy axis i.e., EOut-of-plane < EIn-plane. In this situation, Ec|| is the same as EOut-of-plane and ERadial due to a coincidence of the structures. This EA anisotropy forces the magnetic moment to align in the out-of-plane crystallographic direction, counteracting the thermal agitation and stabilizing magnetic ordering in 2D materials43,54,55. With an increasing curvature, the Crmn-Cr-Crmn plane starts to tilt inward and develops consecutive bending angles, thereby lifting the degeneracy of different spin configurations. The flexo-deformation induces a deviation of the plane formed by the ends of spin vectors from the crystallographic atomic lattice plane. The energy Ec|| now becomes lower than EOut-of-plane and ERadial, leading to a cross point (Fig. 2d) beyond the curvature of 0.042 (cross point under different computational conditions are shown in Supplementary Note 4). The stabilization of such a noncollinear state Ec|| can be attributed to the spin-orbital interactions on the curvature-induced anisotropy6. In the unconstrained calculation, the spin direction of the ground states is more divergent than the initial radial direction after bending (Supplementary Fig. 4), creating a plumb-line-shape-like response of spin configuration to the flexo-bending deformation. This tendency is more prominent on the edge atoms and enhanced with the increase in curvature, thus inducing the transition from a pristine collinear state (radial state) to a noncollinear state (parallel state).

Exchange interaction in asymmetrical systems

The outmost edged Cr atoms (Cr2 and Cr9) are partially saturated with two Cr-I bonds. The local magnetic moments of Cr2 and Cr9 are around 4.10 µB and 3.20 µB, respectively, while the local magnetic moments of the interior Cr atoms (from Cr3 to Cr8) are equal to (3.40 ± 0.01) µB, which is similar to that of the 2D monolayer. It is worth noting that the position of the edged Cr atoms (Cr1, Cr2, Cr9, and Cr10) is fixed during the structural relaxation. The distance between the Cr1 and Cr2 atoms remains around 4.05 ± 0.03 Å, which results in strong ferromagnetic coupling for the Cr1-Cr2 pair (J12 = 4.35 ± 0.5 meV). Such extraordinary interactions are less likely due to the bending effect. Thus, to suppress the effect of the edge states, in the following discussion, we focus on the interior Cr pairs (pairs belong to from Cr3 to Cr8). The calculated Curie temperature Tc slightly increases for bending structures with relatively small curvatures; for instance, Tc is calculated to be 36 and 45 K for the B0 and B150 bent nanoribbons, respectively (Supplementary Fig. 5a). The small changes in Tc can be understood from the fact that the spin couplings J of two neighboring Cr atoms for all the bent nanoribbons fluctuate around 2.25 meV (2.37 meV for the 2D monolayer). Nevertheless, the magnitude of J shows a significant site-dependence on the location across the width direction, and this spatial variation becomes more significant with an increasing curvature (Supplementary Fig. 5b).

To further explain this variation with curvature, we investigate the dependence of the exchange coupling on structural parameters under bending. The maximal information coefficient (MIC)56 method (Calculation details refer to Supplementary Note 3) is adopted to analyze the correlation between the structure parameters of Cr-Cr pairs and the exchange coupling. A larger value of the MIC (ranging from 0 to 1) associated with one specific structural parameter reflects a larger amount of information about the exchange coupling. The maximum value of the MIC corresponds to the Cr-I-Cr angle on the top I layer (θtop), which is much higher than the second one (Fig. 3a). This result is consistent with the Goodenough-Kanamori rules. However, by comparing the MIC of θtop (Cr-I-Cr angle measured at the concave) and θbottom (Cr-I-Cr angle measured at the convex), it is found that for the asymmetric bent structures in which the top and bottom ligand layers have different metal-ligand-metal angles, the effects of these two differently distorted structures on the exchange coupling dramatically differ from each other. Our calculation shows that after bending, the superexchange path (dtop) at the concave accompanied by larger metal-ligand-metal angles is always shorter than the path at the convex (dbottom) (Fig. 3b). Similar structure was also found in the calculation of the Janus structure Cr(X,Y)3 by Albaridy et al.44. Compared to pristine CrBr3, lighter ligands such as Cl atoms result in shorter Cr-X bonds and a larger Cr-X-Cr angle, similarly to the top I layer of bent nanoribbons in our work. From the previous work33, the exchange energy values for CrBr3, Cr(Cl,Br)3 and Cr(I,Br)3 are 32, 27 and 32 meV, respectively. The angle dependence of the exchange energy is consistent with our results about the angle sensitivity of the superexchange interaction in asymmetrical bent CrI3 with the distorted Cr-I-Cr angles. Therefore, a more correlated relationship between θtop and J is identified as being associated with a shorter exchange path on the top iodine atomic layer which has larger Cr-I-Cr angles.

Fig. 3: Variation of exchange interaction with distorted bonds under bending.
figure 3

a MICs between the exchange interaction J and the structural parameters of each Cr-Cr pair, dtop = dtop_L + dtop_R, dbottom = dbottom_L + dbottom_R and dR-L = dtop_R + dbottom_R − dtop_L − dbottom_L. b Sketch of a bending deformation of a Cr-Cr pair with all the structural parameters listed in a. c Contour plot for exchange interaction as the function of Cr-I-Cr angles. The PDOSs of d the Cr7-I24-Cr8 pair in the B0 nanoribbon and e the Cr5-I18-Cr6 pair in the B250 nanoribbon corresponds to points A and B in c, respectively. The inset in e is the sketch of the FM superexchange interaction through t2g-p-eg orbitals and the AFM direct exchange between two Cr half-occupied t2g orbitals.

Full size image

The underlying mechanism can be explained through the θ-J mapping diagram in Fig. 3c showing the response of J to θtop and θbottom. The points A, B and C in Fig. 3c occur at the Cr7-Cr8 pair in the B0 nanoribbon, the Cr5-Cr6 and Cr8-Cr5 pairs in the B250 nanoribbon, respectively. At point A, the angles are the same as the CrI3 monolayer which is equal to 95.03°, and the exchange interaction J (2.21 meV) is comparable with the 2D monolayer (2.37 meV). Points B and C correspond to the maximum (2.89 meV) and minimum J (1.32 meV), respectively. One claim of the Goodenough-Kanamori rules is that when the angle of cation-ligand-cation is >90°, the coupling of both d orbitals of the two cations (Cr) to the same p orbitals of the ligand (I) becomes possible according to the Slater-Koster rules, thus initiating the AFM superexchange which competes with the FM superexchange mediated by the Coulomb exchange on the ligand. This effect is more prominent with a further increase in the angle from 90° to 180°, finally shifting the superexchange coupling from FM to AFM. The change from point B to point C obeys this rule. This signifies our finding that θtop has a stronger effect than θbottom. Figure 3e, f show the projected density of states (PDOSs) of the Cr7-I24-Cr8 pair in the B0 nanoribbon and the Cr5-I18-Cr6 pair in the B250 nanoribbon corresponds to points A and B, respectively. Compared to that of point A, the eg orbitals in point B is spatially close to the t2g orbitals. The enhanced hopping between the half-occupied t2g orbitals and the empty eg orbitals of Cr atoms in point B leads to a stronger FM superexchange.

Discussion

In summary, we have revealed the noncollinear response of the spin configurations toward the flexo-deformation of the CrI3 nanoribbons. Our results suggest a divergent magnetic ground state at small curvature, and it develops into a parallel noncollinear state when the curvature is above 0.06 Å−1. Most importantly, we reveal a gradient dependence of the exchange coupling on the bending angles, and the dominant role of the the superexchange governed by the asymmetric Cr-I-Cr angle at the concave/convex region. Such findings allow the manipulation of spins via intentional control of unique Cr-I-Cr hinges through introducing a proper degree of surface curvature for magnetic sensors and spintronics.

Methods

Construction of bent structure

In this work, we add an angle restraint to the Cr(out-most)-Cr(center)-Cr(out-most) angle (see Fig. 1c). A harmonic type of functional form with different force constants K0 is applied to generate initial structures with different magnitudes of bending curvature. Then, a geometry optimization task that considers the applied restraints is conducted to obtain the preliminary curved structures for the density functional theory (DFT) calculations.

DFT calculation

All DFT calculations are performed using the projector augmented wave (PAW) method as implemented in the Vienna Ab initio Simulation Package (VASP)57. Perdew-Burke-Ernzerhof (PBE) as a particular parametrization of generalized gradient approximation (GGA) is chosen for the exchange-correlation functional58. The onsite Coulomb interaction U determined by the linear response approach59 is set at 2.65 eV for Cr atoms. The energy cutoff is set at 400 eV. The 8 × 8 × 1 and 1 × 8 × 1 k-point grids are utilized to sample the Brillouin zone for the unit cell and nanoribbons, respectively. Periodic boundary conditions are applied along the axial direction, while a 20 Å vacuum space is inserted in the in-plane and out-of-plane directions. All structures are fully relaxed until the force was <0.01 eV Å−1. Specifically, during the structural relaxation of bent nanoribbons, the outermost edge Cr atoms are fixed to maintain the bending condition, and the lattice parameters are only optimized along the axial direction. Spin orbital coupling is taken into account in dealing with noncollinear magnetic state.