Multiferroic kinks and spin-flop transition in Ni 2 InSbO 6 from first principles

,


I. INTRODUCTION
In magnetoelectric multiferroics, a magnetic order coexists and interacts with a ferroelectric one.Several microscopic scenarios of why such coexistence may occur and how the magnetic order can affect the electric polarization have been established [1][2][3][4][5] and the work is rapidly progressing in this direction.Understanding such interactions between the magnetic and electric degrees of freedom is of great importance from both the fundamental and practical points of view.A special attention is paid to the mutual control of the magnetic structure and the electric polarization by applying the magnetic or electric field.For instance, the external magnetic field can control the magnetic structure, while also changing the ferroelectric polarization.Conversely, the external electric field can be used to control the magnetic structure [6].
There are two main types of multiferroic (MF) materials [4].In type-I multiferroics the crystal structure itself is ferroelectric, irrespectively of the magnetism.However, the electric polarization can still be controlled by changing the magnetic structure.In type-II MF the crystal structure is centrosymmetric, but the inversion symmetry can still be broken by a magnetic order, which leads to the ferroelectric polarization.An interesting aspect of the type-I materials is that many of them develop chiral magnetic structure, driven by antisymmetric Dzyaloshinskii-Moriya (DM) interactions in the non-centrosymmetric crystal structure.This chirality can be controlled by the magnetic field, presenting another interesting avenue for magnetoelectric control in type-I materials.For instance, a very special type of the chiral magnetic structure is the skyrmion lattice, which has been intensively studied in the context of MF applications in Cu 2 OSeO 3 and GaV 4 S 8 [7,8].
Ni 2 InSbO 6 (NISO) is one of such chiral MFs.Its low temperature structure has a polar (non-centrosymmetric) rhombohedral R3 space group [9,10] similar to well known MF corundum derivative Ni 3 TeO 6 [11,12].A number of polar corundum derivatives have recently been synthesized and materials with above room temperature magnetism have been found [13,14].These are promising candidates for magnetoelectric applications, in some of which polarization switching has been predicted [15].Previous studies [16] revealed an incommensurate antiferromagnetic proper-screw spiral (helix) within each Ni layer, with a long periodicity of 30 unit cells.The polarization along the threefold rotation axis changes quadratically with the magnetic field due to variations in the spiral order induced by the field [16].In addition, recent experiments have revealed a spin-flop (SF) transition upon applying the magnetic field along the threefold rotation axis [17,18].However, the detailed microscopic analysis of these observations is lacking.The phenomenological mechanisms of the magneto-electric coupling considered so far [19][20][21] are not universal and are influenced by system-dependent factors.Therefore, it is important to construct realistic models of these materials with interacting spins and magnetically-induced polarization, starting from the modern theory of polarization in terms of Berry phases and Wannier centers [22][23][24] and using model parameters obtained from first-principles calculations.According to such calculations in the generalized gradient approximation (GGA), NISO can be regarded as an S = 1 material.In the corundum structure, each Ni 2+ ion is located in a distorted octahedron of O 2− ions.Therefore, the Ni 3d states split into triply degenerate t 2g and double degenerate e g manifolds.The t 2g states are fully occupied and do not significantly contribute to magnetism.On the other hand, the e g states are half-filled and form a group of narrow bands near the Fermi level, which are mainly responsible for the S = 1 physics.Thus, we can greatly simplify our analysis by constructing the 2-orbital model for these bands and extracting all parameters of the model from the first principles calculations in the Wannier basis.
Here we model the exchange interactions and magnetically-induced polarization emerging from such realistic 2-orbital e g model at the half-filling.After extracting the parameters of electronic Hubbard-like model from the first principles calculations, we employ the superexchange theory, which in our case is formulated as a firstorder perturbation theory for the Wannier functions with respect to the hopping parameters.For the exchange interactions, the treatment is equivalent to the standard second-order perturbation theory for the magnetic energy, which in the Wannier basis results in the expression for the spin dependent electric polarization.Our models highlight the emergence of intriguing cross-coupling phenomena in NISO.Specifically, we explore the SF transitions and a crossover to the multiferroic kink array state, both induced by the external magnetic field along the threefold axis.The kinks contribute ferroelectric polarization opposite to that of the collinear state and their energetics can be rationalized in terms of repulsion through the Yukawa-like potential and the competition between the DM and magnetic field setting their chemical potential.Additionally, using a continuous theory, we explore the possibility of cross-control of magnetic (electric) order by an electric (magnetic) field.

II. RESULTS
A. Basic electronic structure, electronic and spin model for NISO Results of electronic structure calculations using the experimental crystal structure, shown in Fig. 1a [10], in GGA [25,26] with the spin-orbit coupling (SOC) for a nonmagnetic state are illustrated in Fig. 1e.These calculations clearly reveal two groups of the Ni 3d bands: six t 2g bands (per two Ni sites) around −1 eV and four e g bands around the Fermi level.In NISO, the t 2g bands are fully occupied and nonmagnetic, while the magnetic properties are mainly associated with the e g bands.Therefore, we pick this group of states to construct a realistic Hubbard-type model, which would capture the magnetic behavior of NISO.Such Hubbard model has the form: where a and b (= 1 or 2) label the e g orbitals, and H on−site stands for the on-site Coulomb interactions, which are specified by the intra-orbital interaction U , Hund's coupling J H and inter-orbital interaction U ′ = U − 2J H [27].
ĉ † iaσ (ĉ iaσ ) in Eq. ( 1) stands for the creation (annihilation) of an electron on the Wannier orbital a of the Ni site i with the spin σ.The parameters of the one-electron part, t abσσ ′ ij , are defined as the matrix elements of GGA Hamiltonian in the Wannier basis.Since the basis is complete for the e g bands, these parameters perfectly reproduce the original band structure in GGA (see Fig. 1e).The main sources of the SOC in NISO are the 5p states of the heavy In/Sb atoms.Therefore, it is important to include the SOC before the wannierization, at the level of regular GGA calculations.Then, although the Wannier functions are formally associated with the Ni e g states, the SOC of the heavy In/Sb atoms will still contribute to the matrix elements t abσσ ′ ij , which can be diagonal as well as off-diagonal with respect to the spin indices.
Next, we map this low energy electronic model onto the spin model employing for these purposes the superexchange theory.In the atomic limit, the ground state corresponding to the high-spin S = 1 state at half-filling is described by the single Slater determinant.The same holds for the one-electron and one-hole e g states emerging in the superexchange theory in the process of virtual excitations.Therefore, here we essentially deal with the one-electron theory.This results in the spin Hamiltoninan: where ⃗ e i stands for the classical spin vector at site i, the first term is the isotropic interaction, the second term is the anisotropic Dzyaloshinskii-Moriya interaction and the last term is the (traceless) symmetric exchange anisotropy [28].
Similarly, the expression for the magnetically-induced polarization is given by where the first term describes isotropic exchange striction, while the following terms originate form the antisymmetric and the (traceless) symmetric anisotropy [29][30][31].This model is derived in the framework of the modern theory of polarization in solids [22][23][24], using perturbation expansion of the Wannier functions.
Since the R3 group has only one three-fold rotation axis (along z in the chosen coordinates), all bonds with the same distance should be transformed into each other by the Ĉz 3 rotation.To illustrate this, consider the bonds surrounding a Ni ion labeled as 0 in Fig. 1c.We index its Ni neighbors in the layer above as j = 1, 2, 3, which are classified as bond type α = 1.Conversely, the Ni ions in the layer below, indexed as j = 4, 5, 6 and categorized as bond type α = 2, are at a slightly longer distance.The vectors connecting ion 0 with ions j = 1, 2, 3 are expressed as ⃗ ϵ 0j = (ϵ , where ϵ ∥ 0j and ϵ ⊥ 0j are the lengths of the vector components parallel and perpendicular to the Ni layer, respectively (see Fig. 1d).For j = 4, 5, 6, similar formula applies, with the same ϵ ∥ 0j , but a slightly different ϵ ⊥ 0j and with the arguments of sin and cos incremented by π.Then, the DM vectors are as follows: where θ ′ j,α = 2πj/3 + θ α , d ⊥ α and d ∥ α are bond-dependent parameters.The DM vectors are antisymmetric: Similarly, contributions from the isotropic (Heisenberg) exchange to the magnetically-induced polarization (first term of Eq. ( 3)) are as follows: with symmetric ⃗ P α j0 : ⃗ P α j0 = ⃗ P α 0j .An additional single-ion term is allowed in systems S > 1/2.However, in the subsequent discussion, we neglect this term, expecting its effect to be minor since the orbital angular momentum matrix elements vanish between e g orbitals, making SOC inactive.Indeed, in our 2-orbital model, the energy splitting due to SOC inside S = 1 triplet results in very tiny ∆E ≈ 7 µeV (see Supplementary Note 2).Turning to the magnetically-induced polarization, we calculate the single-ion term by the formula given in Ref. [32].However, we also neglect this term as it is spin-independent in this 2-orbital model.
We derive analytical formulas for the symmetric and antisymmetric exchange constants in (2) and magneticallyinduced electronic component of the polarization in (3) from the 2-orbital model following the strategy used in Refs.[30,31] (see Methods for the derivation).For exchange interactions we have: while the parameters for magnetically-induced electronic component of the polarization have the form: where ∆E = U + J H , t αβ ij and ⃗ r αβ ij are the spin dependent matrix elements of the hopping and position operator, correspondingly, expanded in Pauli matrices as tij = t0 ij ⊗ σ0 + x,y,z γ is a 3 × 3 identity matrix.J ij and ⃗ P ij are the symmetric isotropic (Heisenberg) exchange (scalar) and exchange striction (vector), while the consecutive terms describe antisymmetric (Dzyaloshinskii-Moriya) exchange and symmetric anisotropic (Ising-like) interactions, respectively.The exchange parameters obtained with this method are summarized in Tables I and II.We find that the significant exchange interactions in NISO predominantly originate from the bonds between Ni2 and its Ni neighbors in the layer above (ions 1,2,3, exchange constant J 1 ) and below (ions 4,5,6, exchange constant J 2 ), as shown in Fig. 1c.The obtained parameters obey J 2 < J 1 , which can be understood considering geometric Ni-O-Ni angles, ∠(Ni-O-Ni, J 1 ) = 129.38 • and ∠(Ni-O-Ni, J 2 ) = 136.75• .Typically, the half-filled 2-orbital model predominantly yields AFM interactions, with FM contributions manifesting as effective suppressions in the AFM iteraction constants.According to the Goodenough-Kanamori rule [33], the bond angle ∠(Ni-O-Ni, J 1 ) being close to 90 • compared to ∠(Ni-O-Ni, J 2 ), leads to J 1 receiving compensations by the ferromagnetic contributions.Consequently, this results in a smaller magnitude of J 1 compared to J 2 .
We note that the 2-orbital model parameters also provide insight into importance of different terms in Eq. ( 2) and Eq. ( 3).As for the magnetic energy, isotropic and antisymmetric exchange contributions are found to be nonnegligible (see Supplementary Note 1 for the values of the symmetric anisotropy tensor components).Actually, the relative magnitudes of isotropic and DM interactions indicate that the symmetric anisotropic exchange is inherently small.This fact aligns with Moriya's paper [28] and our analytical formulas (6).These state that the DM interaction is first-order in SOC, whereas the symmetric anisotropy is second-order in SOC.Among the contributions to the polarization, only the isotropic term is found non-negligible.Therefore, in the following we focus mainly on these terms.The actual values of the anisotropic terms of the magnetically-induced polarization and their small effects are detailed in Supplementary Note 3. Consequently, in the following we primarily focus on the isotropic exchange, the DM exchange and the isotropic term of the magnetically-induced polarization.
We note that the antiferromagnetic contribution overestimates the exchange parameters as the 2-orbital model does not take into account ferromagnetic contributions from t 2g orbitals and Hund's couplings on non-magnetic ions accurately.Additionally, the relatively small value of the denominator in superexchange theory, specifically ∆E = 2.9 eV, further contributes to potential overestimations in the overall exchange parameters.Using these parameters in the Monte-Carlo simulations in a 30 × 30 × 1 simulation box with periodic boundary conditions, we obtain the transition temperature, overestimating the experimental one by approximately a factor of two.The transition temperature was identified by the peak in the heat capacity in Monte-Carlo simulations using the calculated parameters (Supplementary Note 5).Nevertheless, we emphasize that the main physics discussed in the following is purely originating from the relative strength between isotropic and anisotropic exchanges.We now discuss the magnetic ground state obtained from the single-q spiral Ansatz analysis in NISO.The 2orbial parameters show that the J 1 and J 2 are very strong and AFM.Thus, we can expect spins in the neighboring layers to be opposite.Since there are four Ni layers in the hexagonal cell, Fig. 1a, the period of the AFM order coincides with the length of the hexagonal c-axis, thus, we have ⃗ q C = (0, 0, 0).DM interactions modify q-vector from a commensurate phase to an incommensurate one by δ⃗ q.The 2-orbital model parameters show that those bonds fulfill a relation d α .This stabilizes spin-spiral states with the propagation vector within the layer.Namely, the q-vector is modified as ⃗ q IC = ⃗ q C + δ⃗ q = (δq x , δq y , 0).Then, we can consider the situation given in Fig. 2a, where the spin spiral propagates within a Ni layer.
Magnetic energy contributions from isotropic (E iso ) and DM (E DM ) interactions are Then, we take derivative of the energy with respect to δq x and δq y , and find the energy minimum at δ⃗ q, The wavevector components are plotted in Fig. 2a using the parameters from the 2-orbital model.We see that, approximately, (cos ϕ, sin ϕ, 0) ⊥ n⊥ , thus, a cycloidal spiral state (Fig. 2b), while the experiment reported a proper-screw spiral state.However, the wave vector and the spiral period, δq Ans ≈ 0.034 (29 unit cells) are in a good agreement with reported experimental values δq ≈ 0.029 (30 unit cells) [9].Additionally, the symmetric anisotropic interactions tilt the spin-spiral plane.The calculated parameter actually give a small rotation ≈ π/8 (see Supplementary Note 1).Another parameter set, calculated from Green's function method [34], gives a properscrew type spiral with a very long wave length of 142 unit cells (δq GF ≈ 0.007; see Supplementary Note 4).Since the spin spiral type is not important for the following discussion of magnetic kink generation and magneticallyinduced polarization, we use a cycloidal spiral as the ground state of NISO.We note that in this mean-field analysis, the energy does not depend on the rotation of the spiral plane ϕ.This can be straightforwardly confirmed by substituting the analytical formula for the spiral wave vector Eq. ( 9) into the energy Eq. ( 8).The resulting formula yields In order to understand the magnetic ground state and the response to the external magnetic fields in NISO, we also perform classical Monte-Carlo simulations (MCS) of the spin model Eq. ( 2) with the 2-orbital parameters.Incorporating the Zeeman term − i H z e z i , simulations are carried out on a 30 × 30 × 1 supercell with periodic boundary conditions.
Our MCS reveal the emergence of an anticipated spiral ground state, see Fig. 3a.This result validates the applicability of our single-q spiral Ansatz.Upon the application of an external magnetic field perpendicular to the layer, a SF transition is found.This supports the observations reported in Refs.[16][17][18].Remarkably, before the SF transition, we observe the solitons, favored by DM interaction, shrinking, while the regions between them, with spins along the magnetic field, experiencing a notable increase (Fig. 3b).In comparison to the experimentally determined transition field of approximately 20 T, our calculations yield an overestimated field of approximately 100 T. This discrepancy is attributed to the overestimated exchange parameters and finite size effects in the MCS.Despite the success of the MCS, finite size effects significantly influence important spiral deformations by forcing the spiral to be commensurate with the simulation box.Therefore, in the following we derive and employ a continuous Ginzburg-Landau-type theory.

C. 1D continuous model
Application of the external magnetic field can modify the spiral period.Such a phenomenon, often elusive in MCS due to finite size limitations, can result in overlooking of interesting physics like spiral period changes.For instance, De Gennes investigated modifications of the solitonic lattice period by an external magnetic field in liquid crystals [35] and concluded that the spiral period follows the analytical expression: where L(0) is the original period of the perfect spiral at zero field, and K(k), E(k) are elliptic integrals of the first and second kind, respectively.Here, we formulate a continuous theory to capture such period changes and resulting magnetically-induced polarization without any restriction of the periodicity.As we see from MCS, the obtained magnetic structures are quasi-1D, therefore we can effectively model NISO as a 1D-chain.This is due to the following reasons.First, the spiral lies in the xy-plane.Second, important exchanges are only J 1 and J 2 bonds.They involve the group of bonds connected by the Ĉz 3 operation.Those give rise to the energy expression identical to that of a 1D chain.Fig. 4 shows the spin-spiral ground state and the states it deforms into under an external magnetic field perpendicular to the layer, as obtained by minimizing 1D AFM continuous chain energy in Eq. (37).First, application of an external magnetic field generates a kink array state (Fig. 4 c, g, k): instead of rotating in space with the constant wave vector, antiferromagnetic order parameter now rotates in a non-uniform fashion.Indeed, when ⃗ A is pointing perpendicular to the magnetic field ⃗ H, sublattices cant along the field, which leads to a gain of Zeeman energy, linear in the canting angle (at the expense of the quadratic loss in the antiferromagnetic exchange energy between the sublattices).Thus, the regions with ⃗ A ⊥ ⃗ H expand, forming plateaus in A x (Fig. 4 c), which increases the gain of the Zeeman energy on canting.In contrast, between the plateaus, kinks in A x (or solitons) occur where ⃗ A quickly rotates through the field direction.At these points, one sublattice aligns with the field and anotheropposite to it, and therefore the gain of Zeeman energy on the canting is not possible.The shape of these kinks is analogous to the solitons in nonlinear dynamics because they are solutions of the Sine-Gordon equation.The continuous model reveals that with further increase of the field strength, the kinks are pushed apart.At a higher field, the flat spiral turns into a conical one and its plane flops perpendicular to the field (Fig. 4 b, f, j), which enables a higher gain of Zeeman energy from spin canting along the field.The kinks, present both in the flat and conical spiral phases, can be viewed as particles interacting with each other via exchange of magnons.Further increase of the magnetic field strength drives the transition into SF phase.These transitions are very similar to the one found in our MC simulations (Fig. 3b).However, in this continuous theory, the transition is much smoother due to unrestricted simulation box size.In this case, the interactions between kinks are found by solving Euler-Lagrange equation for the magnetic structure deformation, induced by the soliton, and result in the following interaction between kinks, where β = H z /J, and r is the distance between the kinks.This Yukawa-like interaction can be understood as a long-range repulsion caused by exchanging virtual massive magnons (similar to how Coulomb interaction is caused by exchanges of virtual massless photons).High magnetic field results in the enlargement of the area with ⃗ A ⊥ ⃗ H because of the exponential dependence on H z in Eq. (11).This is reminiscent of physics found in multiferroic TbFeO 3 [36].Fig. 5 illustrates the relationship between spiral period and the polarization change.Our continuous theory predicts a period change that closely matches De Gennes's analytical formula Eq. (10).Then, the difference in the polarization (computed as a dipole moment of the cluster in Fig. 1c divided by the volume) between spiral and SF states is given as: The model parameters result in a very small change due to the isotropic exchange, ∆P iso z = 0.923 µC/m 2 .Such a minor polarization change is likely to be imperceptible in experiments.In the first step (flat spiral phase), a change in polarization is induced by a change in the spiral period (kink distance).In the second step (conical phase, above 16 T, indicated by a dotted line in Fig. 5), the curve shows a sharp increase of polarization resulting from a ferromagnetic component developing in xy-plane due to the SF contribution.
Fig. 6 illustrates the phase diagram in the (H z , E z ) plane.Here we account for the electric field perturbatively by introducing the lowest order energy density correction −E z P z .We find that a small electric field along the ferroelectric polarization stabilizes SF state while an opposite field favors the flat spiral.Asymmetry with respect to the sign of E z is caused by the non-centrosymmetric structure of NISO, with the pyroelectric polarization set by a specific cation ordering.The P z contributions from the spiral and FM states are opposite, therefore they are stabilized by opposite electric fields.We note that the large magnetic field continuously deforms the SF state into the FM state.

D. Polarization change during spin flop-to-FM transition
Our MCS results indicate the transition from SF phase to the FM state as we apply an external magnetic field (Fig. 7a).Notably, the experimental data reveals a much greater polarization change during this transition when compared to the change between the spiral and the SF phases.During this process, we can express the isotropic exchange contribution to the polarization as follows: where θ is the angle between the spin and the xy-plane.In our MCS this angle exhibits a linear dependence on the external magnetic field strength, as depicted in Fig. 7b.This linear dependence is also found in the experiment [16].Consequently, we anticipate that the magnetically-induced polarization will resemble the behavior shown in Fig. 7c.

Conical SF
Importantly, the magnitude of this polarization change far exceeds that observed during the spiral → SF phase transition.The polarization change across the sequence of phase transitions, spiral → conical spiral → SF → FM, reproduces well the experimental data [16].

III. DISCUSSION
Our study of the complex interplay of magnetic and electric fields, magnetic exchanges and DM interactions in corundum nickelate derivative NISO has uncovered intriguing properties of multiferroic kinks and their implications for the field of MFs.The methodological advancements allowed a quantitative modelling of key magnetoelectric phenomena in NISO.
We have derived analytical formulas for magnetic exchanges and magnetically-induced polarization, starting from a minimal 2-orbital model, suitable for computing these parameters from first principles.These formulas can be applied to other S = 1 MFs (e.g.Ni 3 TeO 6 , Ni 2 ScSbO 6 ).Additionally, we have developed a continuous theory to address complex changes in the spin-spiral structures, which extends our understanding of spiral MFs.This theory should be applicable to a wide class of materials where spiral state is stabilized by Dzyaloshinskii-Moriya interactions, for example BiFeO 3 .Our analysis starts from the ground state with spiral ordering within the layers.We find that a spiral structure deforms into an array of multiferroic kinks upon the application of an external magnetic field in the plane of the spiral.These kinks play a crucial role in determining the period of the spiral, leading to a small polarization that is opposite to the FM contribution.Importantly, external magnetic fields can directly control the distances between these kinks.When a greater magnetic field is applied parallel to the three-fold axis, it induces the SF transition, consistent with experimental observations.Using analytical and numerical models with parameters calculated from first principles, we have identified a three-step phase transition (spiral → conical spiral → SF → FM) under an applied magnetic field.Our ab-initio model shows that the most important mechanism of the magnetically-induced polarization in NISO is the isotropic exchange striction.During the first phase transition, the polarization undergoes a small change, closely related to the magnetic kink distance (spiral period).The kink density produces a contribution to the polarization that opposes the ferromagnetic contribution, resulting in polarization change similar to the spiral period change.In the second step, the polarization curve shows a sharp increase related to the FM component development in the xy-plane due to SF contribution.In the third step, the polarization changes significantly, simultaneously with the linear change of the canting angle θ with the magnetic field.These polarization changes imply the possibility of switching between these magnetic structures by an external electric field, enabling the electric control of magnetism.The proposed scenario is supported by the excellent agreement with the experimental data on the field dependence of the polarization.
In summary, the work advances the understanding of magnetoelectric effect in spiral multiferroics.NISO shows a rich phase diagram with flat spiral ground state, magnetic kink arrays, conical kink phase and SF phase.We reveal their connection to the magnetically-induced polarization, and the possibility of electric switching between magnetic phases via an external electric field.Given the importance of the manipulation of the spiral MFs by external fields, we expect that our findings will motivate new experiments and facilitate the applications of spiral MFs in the next-generation memory and spintronic devices.

A. First principles calculations
Density functional theory calculations have been performed for the experimental crystal structure [10] (see Fig. 1a for the hexagonal cell) using a rhombohedral cell.These calculations employed norm-conserving pseudopotentials within the Quantum ESPRESSO package [37].The plane wave cutoff was set to 1088 eV, the Brillouin zone was sampled by a 7 × 7 × 7 Monkhorst-Pack k-point mesh.
The screened Coulomb and exchange interactions of Eq. ( 1) are calculated using constrained RPA technique [42], which yields U = 2.3 eV and J H = 0.6 eV.
The crystal structures were visualised with VESTA [43].

B. Analytical formula for the parameters of spin models
The second order perturbation energy with respect to electron hopping in the 2-orbital model is formulated as: where |α + (−), i⟩ indicates occupied (unoccupied) spin-orbital at site i and ∆E = U + J H . Now, we assume Hamiltonian matrix elements tij are ordered as |1⟩ |2⟩ (pairs of Kramers' states).Consequently, the corresponding orbital ket vectors are straightforwardly defined as: |1⟩ = 1 0 and |2⟩ = 0 1 .Furthermore, when considering the occupied and unoccupied spin states, we assume them to have general spinor functions: and These states have the maximal spin projection in the direction of an associated classical spin ⃗ e = ⟨+| ⃗ σ |+⟩ = (sin θ cos ϕ, sin θ sin ϕ, cos θ).
Next, we decompose general 4 × 4 hopping matrices (with the proper phase) into their orbital and spin components, as described by the following equations: with the orbital part where σ0 is the 2 × 2 identity matrix and σγ ; γ = x, y, z are the spin Pauli matrices.Here, we choose the phases such that non-SOC related terms are represented solely by pure real coefficients, while SOC terms are represented solely by pure imaginary coefficients.Then, taking energy difference for different spin orientations and mapping, it is straightforward to find analytical formual Eq. ( 6).
The spin model can also be derived with respect to spin operators acting on S = 1 multiplets.The procedure is simply finding the correspondence between a general spin Hamiltonian: and the second order perturbation energy: Here, . This approach results in the following expressions: where the difference in signs between the second-order energies arises from particle exchanges.This tensor can generally be decomposed into isotropic interaction: , and symmetric anisotropy: However, this approach results in the exactly same result as the classical vector approach Eq. ( 6).The equivalence of the two approaches can be easily demonstrated by considering the relationship between the hopping integrals For example, Substituting this definition into Eq.(9-17) yields the same formula as Eq. ( 6).Here, we provide an explicit notation for the diagonal elements as an example: Therefore, the isotropic interaction is given as All the other expression can be restored following the same procedure.
A similar theory can be developed for magnetically induced electronic polarization.The main idea is to expand the Wannier functions with respect to first-order of the hopping and apply the general theory of electronic polarization in solids.This expansion corresponds to spin-dependent modulations of the Wannier density, resulting in changes in the Wannier centers.Starting point of the theory is the general theory for polarization in solids [23]: where ⟨w i |r|w i ⟩ represents the position operator's diagonal elements in the Wannier basis, and V is the cell volume.
Next, the first-order perturbation expansion of the Wannier function with respect to electron hopping is: By substituting Eq. ( 31) into Eq.( 30), we obtain pair and single-ion terms as ⃗ P = i ⃗ P i + <ij> ⃗ P ij .The pair interaction terms are as follows: Similar to the hopping matrix, the position matrix can be decomposed into a linear combination of the Pauli matrices as: Then, it is straightforward to find the pair interaction formula Eq. ( 3) and the corresponding analytical formula Eq. (7). Figure 2: a Definition of the spin-spiral parameters.b The components of δq minimizing the energy as a function of the spin rotation plane orientation (given by a polar angle ϕ, where the rotation plane normal is n⊥ = (− sin ϕ, cos ϕ, 0)).The energy is minimized by δq ∼ (cos ϕ, sin ϕ, 0) which is perpendicular to n⊥ , giving rise to a cycloidal spiral.The horizontal axis indicates the rotation of the spin-spiral plane.c The cycloidal spiral ground state for the Hamiltonian, with magnetic exchange constants derived from our superexchange theory.Black lines connecting the sites and spins of neighbors are drawn as a guide.Figure 4: a-d Continuous spin, e-h real space continuous spin textures and i-l sphere area covered by each spin array as calculated using the continuous model Eq.( 37) at several external magnetic field strengths H z .Panels a, e, i: SF state; panels b, f, j: conical spiral; panels c, g, k: kink state; panels d, h, l: flat spiral state.
Figure 5: The response to an external magnetic field H z , applied along the three-fold axis in NISO.(Upper panel) Spin spiral period calculated from the continuous model (red line) and using the analytical formula in Ref. [35] (blue line).(Lower panel) magnetically-induced polarization resulting from the isotropic exchange striction.The dotted line indicates phase transition between the flat spiral state and the conical state shown in Fig. 6.
Figure 6: Phase diagram of NISO as obtained by minimization of the continuous model, Eq.( 37), with a small electric field E z .Yellow, green, and red areas are indicating flat spiral phase, conical phase, and SF phase, respectively.The tendency towards FM state is indicated with white in color gradients.To confirm the magnitude of the single-ion anisotropy (SIA) in our model, we perform a diagonalization of the on-site Hamiltonian, given as H on−site = H CF + H SOC + H U , where H CF , H SOC and H U indicate crystal-field, SOC, and on-site Coulomb interaction, respectively.Obtained energy spectrum is depicted in Fig. 2.After the diagonalization, the Sz values are reduced to Sz eff ≈ 0.7 due to SOC.The SIA is determined by the energy difference within the S = 1 multiplet.It is important to note that, due to the influence of SIA, these states do not strictly form a triplet.The energy difference between doubly degenerate ±Sz eff states and non-degenerate Sz = 0 state is E ±Szeff − E Sz=0 ≈ 7 µeV.Indicating very tiny easy-plane anisotropy, that is consistent with the expectation from the experiment [1].Consequently, in our 2-orbital model, the SIA can be considered negligible.

SUPPLEMENTARY NOTE 3: ANISOTROPIC MECHANISM FOR MAGNETICALLY INDUCED POLARIZATION
Here, we discuss the effect of the anisotropic mechanism for the magnetically-induced polarization (the second and the third term in Eq. (3) of the main text) in NISO.We summerize the parameters for the antisymmetric mechanism of the magnetically-induced polarization in Table II.In all bonds, the z-component (the third row of the tensor) of this term results in small compared to the in-plane component (the first and the second row of the tensor).When we take only z-component and write it as ⃗ P α 0j,z for each bond type α, this vector transform between the same α exactly same as the DM vector (without unphysical numerical error due to position operator).Namely, using the same definition of DM vector (Eq.( 4) of the main text), one can write it as ⃗ P α 0j,z = p ∥ α,z cos θ ′ j,α , p ∥ α,z sin θ ′ j,α , p ⊥ α,z .Now, we take into account the long-wave length spin-spiral ground state propagate in the Ni plane as obtained from our analysis and the prior experimental observation.Then, neighboring sites for the same α can be given as Since the spiral rotates in a plane formed by (0, 0, 1) and an arbitrary vector in xy-plane (cos ϕ, sin ϕ, 0), the total polarization along z-direction will vanish.On the other hand, the xy components of the polarization can be finite.Using the values of antisymmetric mechanism given in Table II, regardless of the spiral type (cycloidal or proper-screw), P x and P y are evaluated as P x ≈ −16 µC/m 2 and P y ≈ −9 µC/m 2 , respectively.Therefore, polarization from the antisymmetric term can be neglected in the main discussion.
The parameters for the symmetric anisotropic term of the magnetically-induced polarization are shown in Table III.The symmetric anisotropic term for the magnetically-induced polarization results in very tiny.This is expected from the formula Eq. ( 7) of the main text.The symmetric aniostropy term of the polarization originates from spin off-diagonal term of hoppings and position operators.As one can expect from the relation between isotropic and antisymmetric parameters, we obtain few µC/m 2 order for the symmetric anisotropic term.Thus, we can safely ignore this term in the main discussions.By employing Lloyd's formula, it is possible to represent changes in the electron density of states resulting from perturbations to the Hamiltonian and Green's function.Then, the idea is to establish a connection between this perturbation energy and the energy change in the spin Hamiltonian due to infinitesimal rotations of spins.As a result, for example, isotropic exchange is evaluated as: where the trace runs over the magnetic quantum number m, ∆ i represents the magnetic splitting at site i, and G ij denotes the Green's function.The expansion of this method to include SOC can be derived using the Dyson equation of the Green's function, as discussed in Ref. [3].We use this Green's fucntion method as it is implemented in the TB2J package [4].The calculated exchange parameters from a layer AFM state are given in Tab.IV.The anisotropic exchange parameters yield magnitudes quite different from those obtained by our analytical theory.This suppression of anisotropic exchanges leads to a long-wavelength solution of the spin-spiral Ansatz.By employing the Ansatz equations Eq. ( 9) of the main text, we find the ground state to be approximately properscrew state (Fig. 3).The wavelength results in a very long δq GF ≈ 0.007 (equivalent to 142 unit cells) due to the above-mentioned suppression of the anisotropic exchange.

SUPPLEMENTARY NOTE 5: CRITICAL TEMPERATURE BY MONTE-CARLO SIMULATION
To investigate the transition temperature, we calculate the heat capacity.The heat capacity calculated from the Monte-Carlo simulations is shown in Fig. 4. The result indicates the phase transition around 150 K.This is exaggerated compared to the experimental value 77 K [5].This inconsistency arises from finite size dependence and the overestimated exchange parameters in the 2-orbital model.Given the relationship between the calculated transition temperature and the experimental value, the exchange parameters appear to be overestimated by approximately a factor of two.
Table II: Values of isotropic term for polarization and corresponding parameters in NISO [µC/m 2 ] In the kink array phase, fast change of the angle θ near the kink at position X results in a delta function-like change of the gradient ∂θ ∂x ∝ δ(x − X).

Figure 1 :
Figure 1: a Hexagonal cell of NISO.b Ni ions in the hexagonal cell.Only the closest neighbors coupled with exchange constants J 1 and J 2 are shown with yellow and gray lines, respectively.c J 1 and J 2 bonds from side view and top view around a Ni ion.d Definition of DM vector parameters.This figure explicitly depicts parameters for α = 1 bond type.e Electronic structure of a rhombohedral unit cell of NISO around Fermi level calculated within GGA (solid black line) and 2-orbital model constructed by MLWF method (cyan dashed line).The inset shows a schematic of the k-path in the Brillouin zone for the rhombohedral unit cell.Figure2: a Definition of the spin-spiral parameters.b The components of δq minimizing the energy as a function of the spin rotation plane orientation (given by a polar angle ϕ, where the rotation plane normal is n⊥ = (− sin ϕ, cos ϕ, 0)).The energy is minimized by δq ∼ (cos ϕ, sin ϕ, 0) which is perpendicular to n⊥ , giving rise to a cycloidal spiral.The horizontal axis indicates the rotation of the spin-spiral plane.c The cycloidal spiral ground state for the Hamiltonian, with magnetic exchange constants derived from our superexchange theory.Black lines connecting the sites and spins of neighbors are drawn as a guide.Figure 3: a Magnetic ground state of NISO as obtained from classical Monte-Carlo simulations.The red rectangle indicates the area, for which we show the spin texture in Panel b. b Magnetic structure inside the area, indicated by the red box in Panel a for different values of an external magnetic field, perpendicular to the spiral plane (three-fold rotation axis).Figure4: a-d Continuous spin, e-h real space continuous spin textures and i-l sphere area covered by each spin array as calculated using the continuous model Eq.(37) at several external magnetic field strengths H z .Panels a, e, i: SF state; panels b, f, j: conical spiral; panels c, g, k: kink state; panels d, h, l: flat spiral state.Figure5: The response to an external magnetic field H z , applied along the three-fold axis in NISO.(Upper panel) Spin spiral period calculated from the continuous model (red line) and using the analytical formula in Ref.[35] (blue line).(Lower panel) magnetically-induced polarization resulting from the isotropic exchange striction.The dotted line indicates phase transition between the flat spiral state and the conical state shown in Fig.6.Figure6: Phase diagram of NISO as obtained by minimization of the continuous model, Eq.(37), with a small electric field E z .Yellow, green, and red areas are indicating flat spiral phase, conical phase, and SF phase, respectively.The tendency towards FM state is indicated with white in color gradients.Figure 7: a Schematic representation of the SF phase.b Evolution of the angle θ after SF transition under the external magnetic field as obtained from the MCS.c Polarization change during the transition from SF phase to FM state.

Figure 3 :
Figure 1: a Hexagonal cell of NISO.b Ni ions in the hexagonal cell.Only the closest neighbors coupled with exchange constants J 1 and J 2 are shown with yellow and gray lines, respectively.c J 1 and J 2 bonds from side view and top view around a Ni ion.d Definition of DM vector parameters.This figure explicitly depicts parameters for α = 1 bond type.e Electronic structure of a rhombohedral unit cell of NISO around Fermi level calculated within GGA (solid black line) and 2-orbital model constructed by MLWF method (cyan dashed line).The inset shows a schematic of the k-path in the Brillouin zone for the rhombohedral unit cell.Figure2: a Definition of the spin-spiral parameters.b The components of δq minimizing the energy as a function of the spin rotation plane orientation (given by a polar angle ϕ, where the rotation plane normal is n⊥ = (− sin ϕ, cos ϕ, 0)).The energy is minimized by δq ∼ (cos ϕ, sin ϕ, 0) which is perpendicular to n⊥ , giving rise to a cycloidal spiral.The horizontal axis indicates the rotation of the spin-spiral plane.c The cycloidal spiral ground state for the Hamiltonian, with magnetic exchange constants derived from our superexchange theory.Black lines connecting the sites and spins of neighbors are drawn as a guide.Figure 3: a Magnetic ground state of NISO as obtained from classical Monte-Carlo simulations.The red rectangle indicates the area, for which we show the spin texture in Panel b. b Magnetic structure inside the area, indicated by the red box in Panel a for different values of an external magnetic field, perpendicular to the spiral plane (three-fold rotation axis).Figure4: a-d Continuous spin, e-h real space continuous spin textures and i-l sphere area covered by each spin array as calculated using the continuous model Eq.(37) at several external magnetic field strengths H z .Panels a, e, i: SF state; panels b, f, j: conical spiral; panels c, g, k: kink state; panels d, h, l: flat spiral state.Figure5: The response to an external magnetic field H z , applied along the three-fold axis in NISO.(Upper panel) Spin spiral period calculated from the continuous model (red line) and using the analytical formula in Ref.[35] (blue line).(Lower panel) magnetically-induced polarization resulting from the isotropic exchange striction.The dotted line indicates phase transition between the flat spiral state and the conical state shown in Fig.6.Figure6: Phase diagram of NISO as obtained by minimization of the continuous model, Eq.(37), with a small electric field E z .Yellow, green, and red areas are indicating flat spiral phase, conical phase, and SF phase, respectively.The tendency towards FM state is indicated with white in color gradients.Figure 7: a Schematic representation of the SF phase.b Evolution of the angle θ after SF transition under the external magnetic field as obtained from the MCS.c Polarization change during the transition from SF phase to FM state.

Figure 7 :
Figure 1: a Hexagonal cell of NISO.b Ni ions in the hexagonal cell.Only the closest neighbors coupled with exchange constants J 1 and J 2 are shown with yellow and gray lines, respectively.c J 1 and J 2 bonds from side view and top view around a Ni ion.d Definition of DM vector parameters.This figure explicitly depicts parameters for α = 1 bond type.e Electronic structure of a rhombohedral unit cell of NISO around Fermi level calculated within GGA (solid black line) and 2-orbital model constructed by MLWF method (cyan dashed line).The inset shows a schematic of the k-path in the Brillouin zone for the rhombohedral unit cell.Figure2: a Definition of the spin-spiral parameters.b The components of δq minimizing the energy as a function of the spin rotation plane orientation (given by a polar angle ϕ, where the rotation plane normal is n⊥ = (− sin ϕ, cos ϕ, 0)).The energy is minimized by δq ∼ (cos ϕ, sin ϕ, 0) which is perpendicular to n⊥ , giving rise to a cycloidal spiral.The horizontal axis indicates the rotation of the spin-spiral plane.c The cycloidal spiral ground state for the Hamiltonian, with magnetic exchange constants derived from our superexchange theory.Black lines connecting the sites and spins of neighbors are drawn as a guide.Figure 3: a Magnetic ground state of NISO as obtained from classical Monte-Carlo simulations.The red rectangle indicates the area, for which we show the spin texture in Panel b. b Magnetic structure inside the area, indicated by the red box in Panel a for different values of an external magnetic field, perpendicular to the spiral plane (three-fold rotation axis).Figure4: a-d Continuous spin, e-h real space continuous spin textures and i-l sphere area covered by each spin array as calculated using the continuous model Eq.(37) at several external magnetic field strengths H z .Panels a, e, i: SF state; panels b, f, j: conical spiral; panels c, g, k: kink state; panels d, h, l: flat spiral state.Figure5: The response to an external magnetic field H z , applied along the three-fold axis in NISO.(Upper panel) Spin spiral period calculated from the continuous model (red line) and using the analytical formula in Ref.[35] (blue line).(Lower panel) magnetically-induced polarization resulting from the isotropic exchange striction.The dotted line indicates phase transition between the flat spiral state and the conical state shown in Fig.6.Figure6: Phase diagram of NISO as obtained by minimization of the continuous model, Eq.(37), with a small electric field E z .Yellow, green, and red areas are indicating flat spiral phase, conical phase, and SF phase, respectively.The tendency towards FM state is indicated with white in color gradients.Figure 7: a Schematic representation of the SF phase.b Evolution of the angle θ after SF transition under the external magnetic field as obtained from the MCS.c Polarization change during the transition from SF phase to FM state.

FigFigure 1 :
Figure 1: Spin-spiral plane rotation induced by the symmetric anisotropic exchange interactions.a The energy change by the spiral plane rotation solely due to the symetric anisotropic exchange interactions.b Total energy change by the spiral plane trotation.c A schematic for the spin-spiral state given as the minima of b.

Figure 2 :
Figure 2: Energies of two-electron state on a Ni 2-orbital model as obtained from diagonalization of the on-site Hamiltonian.

Figure 3 :
Figure3: Single-q Ansatz result for the exchange parameters calculated using Green's function method.

Figure 4 :
Figure 4: Heat capacity as obtained from Monte-Carlo simulation

Table I :
Values of isotropic exchanges and DM parameters in NISO [meV]

Table II :
Values of antisymmetric mechanism of the magnetically-induced polarization for bonds given in Fig.1cin the main text [µC/m 2 ]

Table III :
Values of the symmetric anisotropic mechanism of the magnetically-induced polarization for bonds given in Fig. 1c [µC/m 2 ]

Table IV :
Values of magnetic exchanges in NISO [meV ] and corresponding parameters calculated by Green's function method