Atomic-Scale Visualization of Multiferroicity in Monolayer NiI 2

Progress in layered van der Waals materials has resulted in the discovery of ferromagnetic and ferroelectric materials down to the monolayer limit. Recently, evidence of the ﬁrst purely 2D multiferroic material was reported in monolayer NiI 2 . However, probing multiferroicity with scattering-based and optical bulk techniques is challenging on 2D materials, and experiments on the atomic scale are needed to fully characterize the multiferroic order at the monolayer limit. Here, scanning tunneling microscopy (STM) supported by density functional theory (DFT) calculations is used to probe and characterize the multiferroic order in monolayer NiI 2 . It is demonstrated that the type-II multiferroic order displayed by NiI 2 , arising from the combination of a magnetic spin spiral order and a strong spin-orbit coupling, allows probing the multiferroic order in the STM experiments. Moreover, the magnetoelectric coupling of NiI 2 is directly probed by external electric ﬁeld manipulation of the multiferroic domains. The ﬁndings establish a novel point of view to analyze magnetoelectric eﬀects at the microscopic level, paving the way toward engineering new multiferroic orders in van der Waals materials and their heterostructures.


INTRODUCTION
0][11][12] Multiferroics can be classified based on the presence of weak or strong coupling between the orders as type I or type II, respectively. 7,8In the particular case of magnetic and ferroelectric orders, type-II multiferroics showing a strong magnetoelectric coupling 13,14 have a huge potential arXiv:2309.11217v1[cond-mat.mtrl-sci]6][17] Nevertheless, multiferroics displaying a strong magnetoelectric coupling at sufficiently high temperatures to build functional devices remain elusive.
The dawn of two-dimensional (2D) materials has enabled new strategies for the design of artificial multiferroics. 187][28][29][30][31] Recently, it was shown that the multiferroicity of the bulk van der Waals compound NiI 2 remains in the mechanically-exfoliated few-layer 32 and monolayer limits 3 .Theoretical analyses have addressed the origin of this type-II multiferroic order to the combination of the magnetic spin-spiral order of NiI 2 and the strong spin-orbit (SOC) coupling of the I atoms. 5To provide experimental evidence of 2D multiferroicity, circular dichroic Raman, birefringence, and second-harmonic-generation (SHG) measurements on both mono-and multi-layer NiI 2 were performed, establishing NiI 2 as the first purely 2D multiferroic with a transition temperature T C = 21 K. 3 However, optical measurements might not suffice to demonstrate the emergent ferroelectricity in the 2D limit. 4In addition, typical bulk techniques, such as neutron scattering, used to identify the magnetic spin spiral that gives rise to the multiferroic order, cannot be applied in the monolayer limit. 33These factors suggest that the origin of the multiferroic order in monolayer NiI 2 has not been experimentally established and its full characterization requires measurements at the microscopic level.
In this work, we present real-space visualization of the multiferroicity in monolayer NiI 2 .We show that the magnetoelectric coupling present in this kind of type-II multiferroics allows us to directly probe and characterize the multiferroic order by scanning tunneling microscopy (STM).
These results are supported by non-collinear ab initio calculations based on density functional theory (DFT) in monolayer NiI 2 .Moreover, we provide evidence of the magnetoelectric coupling of this 2D system by electrical manipulation of the multiferroic domains.Our findings demonstrate the visualization of multiferroic order in van der Waals materials with atomic-scale resolution, establishing a strategy that can be further applied to artificial multiferroic van der Waals heterostructures.

MULTIFERROIC ORDER IN NiI 2
In bulk, NiI 2 displays a magnetic spin spiral at low temperatures. 33DFT calculations have determined that the q vector characterizing the spin spiral order in bulk is a consequence of the competition between magnetic exchange interactions between Ni atoms 5,34,35 .In particular, ferromagnetic intralayer first-neighbor J 1 , antiferromagnetic intralayer third-neighbor J 3 , and antiferromagnetic interlayer second-neighbor J2 magnetic exchange interactions are the most relevant.
These interactions lead to a q vector in bulk whose in-plane component lies on the ΓM segment in reciprocal space, and whose value can be modified by external parameters, such as pressure, that change the ratio between the different magnetic exchanges. 34,36In the monolayer limit, there are no interlayer interactions, hence the q vector is determined by the competition between J 1 and J 3 (Fig. 1b).DFT calculations have predicted that in this scenario the q vector lies on the ΓK segment in reciprocal space, corresponding to the third-neighbor bond direction. 37The classical Heisenberg model solution predicts that q = (q, q, 0), with q in units of the reciprocal lattice of the structural single-unit cell with lattice vector a, is determined by the J 1 /J 3 ratio. 38,39art from the magnetic order, theoretical analyses have shown that monolayer NiI 2 develops a ferroelectric polarization due to the combination of this spin spiral order and the strong spin-orbit coupling of the I atoms. 5These results are based on the theoretical derivation for the emergence of ferroelectricity in spiral magnets. 40This approach predicts that a net ferroelectric polarization P, which is proportional to the SOC strength λ SOC , arises perpendicular to the q and e rotation vector that characterizes the spin spiral (as shown in Fig. 1a).At the more microscopic level, this relationship shows that given a spin spiral with magnetization M = M(− sin(qr), cos(qr), 0) in cartesian coordinates (Fig. 1c) the emergent electric polarization is given by where Λ is proportional to λ SOC and some physical constants. 41Equation (1) leads to P = Λq(− sin(2qr), sin 2 (qr), 0) for the given spin spiral M. The spatial average of this equation leads to a net polarization P with non-zero components only in the direction perpendicular to the q and e vectors. 40,41Equation (1) shows that the emergent electric polarization has a real space modulation whose periodicity is half of the periodicity of the spin spiral (Fig. 1d).
The modulation of the induced ferroelectric order, which is produced by the coupling between spin and charge degrees of freedom, is the property that can be used to characterize the multiferroic order of monolayer NiI 2 with an STM.The emergence of ferroelectricity is accompanied by ferroelectric displacements of the atomic positions in the directions of the electric polarizations and results in the emergence of a modulated electrostatic potential.This electrostatic potential displays the same periodicity as the ferroelectric polarization and modulates the energy bands of NiI 2 .The modulation results in an observable contrast in STM images, thus providing a signature of the multiferroic order allowing its characterization (Fig. 1e).The magnetoelectric coupling occurring in this multiferroic allows the visualization of the magnetic spin spiral periodicity without the need to use more complex spin-polarized scanning microscopy techniques.In the real space images, the light blue rectangle depicts the unit cell of the commensurate spin spiral shown in Fig. 1a.In contrast to the theory, the experimental stripes are not aligned precisely with the atomic lattice.In the FFT images, the red-circled peaks correspond to the single unit cell giving rise to the depicted Brillouin zone, while the green-circled peaks are associated with the ferroelectric modulation.

STM CHARACTERIZATION OF MONOLAYER NiI 2
Figure 2a shows a large area STM scan of our sample (see SI for details).It reveals the presence of a well-formed monolayer NiI 2 on HOPG.The dI/dV spectrum of monolayer NiI 2 is shown in Fig. 2b.The valence and conduction bands can be distinguished at −1.9 V and 0.4 V respectively, indicating that NiI 2 is an insulator with a gap of 2.3 eV.Figures 2c-d show small-area STM images at 0.3 V (within the NiI 2 bandgap) and at 0.9 V (within the NiI 2 conduction band).A systematic analysis with scans at different energies can be found in the SI.Inside the gap, the contrast is not affected by the polarization of the NiI 2 , and we observe a hexagonal moiré modulation arising from the lattice mismatch between HOPG and NiI 2 (Figure 2c).However, as the conduction band is modulated by the ferroelectric polarization, the stripe-like pattern expected for this kind of type-II multiferroic appears at biases corresponding to the NiI 2 conduction band (Fig. 2d).
Two distinct regions with different stripe orientations can be observed.These correspond to two different multiferroic domains.Finally, there are different types of atomic scale defects that are analysed in detail in the SI.
Figure 2e shows an atomic-resolution scan in a single domain region, the FFT of this scan is shown in Fig. 2f.The spin-spiral q vector can be determined from the relationship between the single unit cell a × a and the stripe modulation.From our experimental results, we obtain a = 3.85 Å for monolayer NiI 2 .The periodicity of the stripes is L S = 17.8 Å corresponding to about 4.6a.According to Eq. (1) the modulation of the stripes caused by the emergent electric polarization displays at half the periodicity of the spin spiral (Figs.1c-d), leading to a spin spiral periodicity of 9.2a in the experiment.The spin-spiral q vector can be directly extracted from the peaks associated with the stripes in the FFT (green circles in Fig. 2f) as half of the peak value.We have obtained q = (0.069, 0.041, 0) in units of the reciprocal lattice vectors, showing a small deviation from ΓK segment expected for the J 1 − J 3 spin model.The projection in the ΓK segment is 0.173ΓK, corresponding to q = (0.057, 0.057, 0).It can be observed that this small variation in the experimental q vector from the ΓK segment leads to a small tilt (∼ 4.6 • ) of the stripe modulation (Fig. 2f).This effect could be a consequence of local defects in the sample, the proximity of other multiferroic domains, or the influence of other parameters such as intralayer second-neighbor magnetic exchange, Kitaev, or biquadratic interactions not considered in the J 1 − J 3 spin model. 39,44Finally, the diffuse peaks around the M points are associated to the defect seen in the scan.The analysis of the filtered FFT images is shown in the SI.
Extracting the information on the spin spiral from the STM measurements allows us to model monolayer NiI 2 using DFT.Neglecting the small tilt in the ferroelectric modulation, the unit cell describing the multiferroic order in monolayer NiI 2 can be well approximated by the commensurate 9a × √ 3a supercell shown in Fig. 1a, corresponding to q = (0.064, 0.064, 0) and √ 3/9ΓK ∼ 0.192ΓK.From the experimental q = (0.057, 0.057, 0) we can estimate J 3 /J 1 = −0.263. 38,39This result shows that monolayer NiI 2 is on the verge of a spin-spiral to ferromagnetic transition (J 3 /J 1 = −0.25), 39which might explain why in some studies the presence of ferroelectric polarization was only reported down to the bilayer limit but not in the monolayer. 32It also points out that DFT tends to overestimate the |J 3 /J 1 | ratio. 36,37This can be rationalized by the typical limitations of the different DFT functionals to deal with electronic correlations and electronic localization, thus introducing an overestimation of the long-range magnetic exchanges.This could explain why ab initio calculations predict a spin spiral for NiCl 2 while the observed ground state is ferromagnetic. 37 order to shed light on the microscopic origin of the STM modulation observed in the experiments, we have performed non-collinear DFT calculations in the commensurate 9a × √ 3a supercell and the q = (0.064, 0.064, 0) spin spiral (see SI for computational details).When SOC is included in the atomic relaxations, we observe the emergence of ferroelectric forces that drive the atoms to new positions following the same periodicity as the electric polarization (Fig. 1d).Due to this symmetry breaking in real space, the computed STM image (Fig. 2g) and its FFT (Fig. 2h) display a modulation with the periodicity of the electric polarization, i.e. half of the magnetic spiral, in a good agreement with the experimental observations.
Having characterized the spiral order of monolayer NiI 2 , we now proceed to determine the strength of the ferroelectric polarization.6][47] The emergence of an inhomogeneous electrical polarization in NiI 2 is expected to produce similar band bendings, but in this case within the same multiferroic domain and with the periodicity of the ferroelectric modulation (Fig. 1d).
This behavior is experimentally demonstrated and theoretically reproduced in Fig. 3.
Figure 3b shows dI/dV spectra recorded along a line in a single domain region (green arrow shown in Fig. 3a) parallel to the q vector.Figure 3d shows the DFT simulated line spectra for the commensurate 9a × √ 3a supercell and the q = (0.064, 0.064, 0) spin spiral.Both plots show the conduction bands.A clear underestimation of the band gap is observed for the DFT calcula-tions, which is typical for the LDA approximation.However, both experiment and theory show a band bending modulation of the conduction band that becomes more apparent when analyzing the dI/dV spectra at a maximum and a minimum of the stripe modulation in Fig. 3c (spectra taken at the positions indicated by the blue and red dots in Fig. 3a, respectively).A clear shift in energies between these two positions occurs as a consequence of the ferroelectric modulation introducing a band bending.The shift in energies is also seen in the DFT calculations (Fig. 3f).
The band bending can be analyzed by tracking the shift in the differential conductance maxima as a function of the position across the modulation.This is shown in Figs.3d and 3g for the experimental and theoretical data, respectively.The observed oscillatory behavior can be fitted to E = E 0 + E P sin(2πx/L S + ϕ), where E 0 corresponds to an average of the differential conductance maximum at the bottom of the conduction band, E P is the band bending modulation caused by the inhomogeneous emergent polarization, L S is the stripe periodicity that has already been analyzed and, ϕ is a trivial phase related to the starting point considered for the fitting.
We obtain E P = 16.8 mV from the experiment and E P = 5.0 meV from the DFT calculation.
This underestimation of the ab initio calculations could be a consequence of the limitations of the LDA approximation for the atomic relaxations.However, apart from this small difference, the qualitative agreement between theory and experiment is remarkable.For phonon-driven ferroelectrics such as SnSe and SnTe, STM experiments have shown the band bending is on the order of several hundreds of meV [45][46][47][48] and an associated electric polarization P ∼ 10 -10 C/m has been estimated. 46,48Therefore, considering that in monolayer NiI 2 we observe a band bending of 1 or 2 orders of magnitude smaller, a P ∼ 10 -12 C/m can be estimated for monolayer NiI 2 , consistent with previous theoretical estimations. 5

MANIPULATION OF MULTIFERROIC DOMAIN BOUNDARIES
We have also directly probed the magnetoelectric coupling in monolayer NiI  been achieved before in ferroelectric SnTe. 47Now, we show here that multiferroic domains in this kind of magnetic spin-spiral multiferroics are also tunable by external electric fields, thus showing direct evidence of the magnetoelectric coupling in monolayer NiI 2 .

CONCLUSIONS
In conclusion, we have probed and characterized the multiferroic order in monolayer NiI 2 down to the atomic scale.Using a combination of STM experiments and non-collinear ab initio calculations, we show that in this spin-spiral SOC-driven multiferroic, the emergent ferroelectric polarization induces a modulation of the electrostatic potential.This is directly reflected in the local density of states allowing its direct visualization by STM.This modulation has half of the spin-spiral periodicity, allowing us to characterize the spin-spiral vector as q = (0.069, 0.041, 0) and the parameter relation FIG. 1. a: Unit cell of monolayer NiI 2 .9a × √ 3a supercell is required to describe the magnetic spin-spiral order with propagation vector q in the x-direction and rotation vector e in the z-direction.The presence of spin-orbit coupling λ SOC induces a net electric polarization P in the y-direction.b: Ferromagnetic first neighbor J 1 and antiferromagnetic third neighbor J 3 magnetic exchange interactions between Ni atoms giving rise to the spin spiral order in the monolayer.c: Magnetization components of the spin spiral showing a periodicity of 9a in the x-direction.d: Electric polarization components associated with the magnetic spin spiral in the presence of SOC showing a periodicity of 9a/2 in the x-direction.e: Modulated electrostatic potential that can be probed with an STM.It shows a 9a/2 periodicity in the x-direction as a consequence of the multiferroic order.

FIG. 3 .
FIG. 3. a: Small area STM scan showing a multiferroic domain where a line spectrum was performed along the green arrow (V = 1 V, I = 100 pA).b: dI/dV line spectrum capturing the conduction band of monolayer NiI 2 along the green arrow in panel (a) parallel to the q vector.c:dI/dV spectra showing the shift in energy between a maximum and a minimum of the stripe modulation (blue and red dots in panel (a))d: dI/dV intensity maximum of the bottom of the conduction band along the green line showing a modulated band bending as a consequence of the ferroelectric modulation.e-g: Corresponding DFT calculations for panels (b-d) computed in the 9a × √ 3a supercell and the q = (0.064, 0.064, 0) spin spiral.

2
by external electric field manipulation of the multiferroic domains.Figs.4a-4dshow the movement of the multiferroic domain wall induced by voltage pulses from the STM tip (the bias is swept from 1 V to 4 V with a closed feedback loop and a current set point of 100 pA).Starting in Fig.4a, the multiferroic domain wall (black-dashed line) and different kinds of defects (highlighted with circles, detailed analysis in the SI) can be identified.The neutral defects (pink circle) can be used as a spatial reference for the domain manipulation.A first voltage pulse is performed in the position indicated by the sky blue dot (in the left domain) and it leads to the configuration shown in Fig.4b.We

FIG. 4 .
FIG. 4. a-d: Manipulation of the multiferroic domains by sweeping the bias from 1 V to 4 V with closed feedback loop (image size 30 × 27 nm 2 , V = 1 V, I = 20 pA).Sky blue dots mark the position of the tip during the process and the numbers on them track the order of events.The black-dashed lines highlight the domain boundary, the pink circle shows a neutral defect on the surface used as a spacial reference, and the dashed red-and-green circle corresponds to a type of mobile charge defect.
J 3 /J 1 = −0.263 of the associated spin model.The observed band bending allows us to estimate the polarization P ∼ 10 -12 C/m for monolayer NiI 2 .Finally, we have probed the magnetoelectric coupling of NiI 2 by manipulating the multiferroic domains by local electric fields induced by the STM tip.Our results firmly establish the atomic scale origin of multiferroicity in NiI 2 and pave the way for future studies on multiferroics and 2D van der Waals materials at the microscopic level.