Ferroelectric Switching at Symmetry‐Broken Interfaces by Local Control of Dislocations Networks

Semiconducting ferroelectric materials with low energy polarization switching offer a platform for next‐generation electronics such as ferroelectric field‐effect transistors. Recently discovered interfacial ferroelectricity in bilayers of transition metal dichalcogenide films provides an opportunity to combine the potential of semiconducting ferroelectrics with the design flexibility of 2D material devices. Here, local control of ferroelectric domains in a marginally twisted WS2 bilayer is demonstrated with a scanning tunneling microscope at room temperature, and their observed reversible evolution is understood using a string‐like model of the domain wall network (DWN). Two characteristic regimes of DWN evolution are identified: (i) elastic bending of partial screw dislocations separating smaller domains with twin stackings due to mutual sliding of monolayers at domain boundaries and (ii) merging of primary domain walls into perfect screw dislocations, which become the seeds for the recovery of the initial domain structure upon reversing electric field. These results open the possibility to achieve full control over atomically thin semiconducting ferroelectric domains using local electric fields, which is a critical step towards their technological use.


Introduction
Overcoming challenges in modern computer engineering and telecommunication technologies requires compact multifunctional components.Atomically thin films of semiconducting transition metal dichalcogenides (TMDs) offer multiple advantages for such a development.They retain robust semiconducting properties down to sub-nanometer thickness, which allows for an efficient transistor operation, and they exhibit strong light-matter coupling.[3][4][5][6][7][8] Semiconducting ferroelectrics with low energy barrier polarization switching are of interest, as they enable ferroelectric field-effect transistors, devices that combine memory storage and logic processing into a single device. [9]The possibility of assembling vertical stacks of designer sequences from atomically thin planes of layered van der Waal materials allows for unprecedented flexibility in engineering novel electronic devices. [10,11]14] Lack of inversion symmetry among individual monolayers of TMDs permits control over the symmetry of interfaces formed in layer-by-layer assembled heterostructures.For example, when assembling two layers with inverted orientation of unit cells, one obtains bilayers with a local inversion-symmetric structure, which reconstructs [12,14,15] (at small twist angles) into a honeycomb set of 2H stacking domains, [14,16,17] separated by domain walls that have the form of intra-layer shear solitons. [18][4][5]19] The bilayers with parallel orientation of monolayer unit cells exhibit two energetically equivalent favorable stacking orders, which are mirror images of each other (illustrated in Figure 1a).In one type of stacking (XM'), the metal site M' of the bottom layer appears In area B, the triangular moiré superstructure is homogeneous, with equilateral triangles, and has periodicity given by the translation vectors ' t 1,2 .In area A, the moiré pattern is skewed due to a small strain in one of the layers (whose effect is magnified by the moiré superlattice effect [23] ).d) Spatial displacement field map for the inter-layer deformations, obtained using Eq. ( 1) from the variation of moiré superlattice pattern in c).
under the chalcogen site X of the top layer, whereas the chalcogen/metal sites (X'/M) in the bottom/top layers sit under/over the middle of the honeycomb of the other layer's lattice.In the second type of stacking (MX'), the metal site M of the top layer appears over the chalcogen site X' of the bottom layer, whereas the chalcogen/metal sites (X/M') in the top/bottom layers sit over/under the middle of the honeycomb of the other layer's lattice.The distribution of transferred charge density between the layers is shown in in Figure 1a for XM' stacked WS 2 bilayers (see methods for details on the calculation).This demonstrates formation of an electrical dipole layer between the inner chalcogen sublayers, that is characterized by electric polarization areal density P = 3.26 × 10 −3 e/nm, suggesting that this system features an interfacial ferroelectric effect.For small-angle twisted P-bilayers, in-plane relaxation [12,14,16,17] leads to the alternating twinning of XM' and MX' domains carrying opposite electric polarization. [14,16,17,19,20]These domains are separated by domain walls that are similar to partial dislocations in the 3R polytype of bulk TMDs.In the context of the interfacial ferroelectric effect, by coupling the domain polarization to an out-of-plane electric field, it should be possible to change the polarization of the sample by domain switching, which involves structural changes.In the context of the interfacial ferroelectric effect in twistronic structures this appears as a lateral interlayer shift of the lattices transforming MX' into XM' stacking (or vice versa, depending on the direction of the out-of-plane electric field).This repolarization via structural modification happens primarily at the boundaries of domains with the opposite polarization, leading to the shift in the domain wall position (as it also happens in bulk ferroelectrics).

Results and Discussion
Below, we study regimes and evolutions of interfacialferroelectric domain structures using P-aligned WS 2 samples.Those were assembled as schematically indicated in Figure 1b (optical micrographs and fabrication details are available in Figure S1 in Supporting Information).A thin graphite film was placed under the TMDs to improve electrical contact.Scanning tunneling microscopy and spectroscopy (STM/STS) was used to image and control the reconstruction domains (XM' and MX') formed in the twisted WS 2 samples.As indicated in Figure 1b, such an experiment implies the presence of a perpendicular electric field between the tip and the sample, which can couple to the local polarization of the XM' and MX' stacking areas.A typical scanning tunneling topography for our samples (Figure 1c) shows large triangular domains of XM' and MX' stacking (Figure 1a), separated by a network of domain walls with clearly identifiable nodes (to which one can attribute [14] energetically unfavorable XX' stacking characterized by vertical alignment of chalcogens in the top and bottom layers).According to the earlier observations [16,17,19] and theories, [14] such domains form in P-aligned bilayers of TMDs for the twist angles below 2.5°.We note that, as such images are acquired in constant current mode, the large apparent height misrepresents the actual topography of the sample surface (in particular, reflecting the effect of electric polarization on the value of tunneling current from the tip into the substrate).
Within the triangular domains of the network in Figure 1c, we first focus on the area B containing a network of approximately equilateral triangular domains.These equilateral domains have periodicity of ≈78 nm corresponding [16,17,19,21] to a twist angle 0.23°.This area will be used as a reference for the analysis of a larger-scale aperiodic domain pattern (domains of various sizes and with pronounced anisotropy).Such a distortion of the domain network is caused by a small accidental strain in either of the assembled WS 2 monolayers, inflicted in the fabrication process. [22]Theoretically, it has been demonstrated [23] that a moiré pattern acts as a magnifying glass for small deformations in the constituent monolayers of the bilayer.Strain-induced shifts of the atomic registry between the layer -even as small as a fraction of the lattice constant -lead to a shift of the equivalent stacking areas in the moiré pattern by distances comparable to the moiré superlattice period.Here, we use this generic property of moiré patterns to map the displacement field distribution (due to the accidental strain) and analyze the position of domain wall network nodes, which feature XX' stacking.For this, we successively enumerate all nodes by integers (n 1 ,n 2 ), starting from the referenced point O = (0, 0) inside area B in Figure 1c with approximately equilateral domains.For this area, we determine a twist angle  = Then, we take into account that, without strain, the positions, ⃗ r n 1 ,n 2 , of the XX' network nodes in the rest of the sample would correspond to the mutual shift of atomic registry such that whereas the difference between the actual, ⃗ R n 1 ,n 2 , and expected, ⃗ r n 1 ,n 2 , positions of the network sites enables us to determine the strain-induced interlayer shift (superimposed onto the local lattice reconstruction into domains) of the monolayer WS 2 lattices as, [23] ⃗ u The resulting deformations map is shown in Figure 1d, indicating a vortex in top left corner and discernible deformations from the sides of the analyzed area.To further elaborate the character of the moiré periodicity distortion caused by the deformations we analyzed distribution of two independent, hydrostatic and shear, strain components calculated in local principal axes.We found that the deformations produce both the hydrostatic and shear strains (Figure S7 in Supporting Information) that, as expected from the deformation map, attain the highest values in the vicinity of top left corner in area A in Figure 1c.
The effect of the reconstruction domains on the local electronic properties, as measured using STM/STS, is discussed in Figure 2. The STM topographic map in Figure 2a shows a typical region with triangular domains.By observing the domain wall bending (domain expansion/shrinking) in response to the applied biased voltage, the bright and dark domains were identified from the ferroelectric polarization as MX' and XM' stacking regions, respectively. [24]The contrast in the STM topography of the two stackings is consistent with previously reported STM measurements of twisted TMD homobilayers with similar scanning parameters [25] .The tunneling spectrum in Figure 2b, averaged over all domains, shows the conduction and valence band edges at ≈0.6 and −1.5 eV respectively, indicating an electronic band gap consistent with previous measurements of this system. [26]The differential tunneling conductance (dI/dV) maps acquired at the indicated voltages in Figure 2c, demonstrate the dependence of the local density of states on the stacking order inside the domains and at the domain walls.
Although a scanning tunneling microscope uses an atomically sharp tip end as a local probe, the overall shape of the tip and its holder make the electric field between the tip and the sample approximately uniform over a typical domain size (Figure 1b).To verify this, we compared forward and reverse topographic scans (Figure S2 in Supporting Information) that are found to be identical, supporting the scenario of an approximately uniform out-ofplane electric field.Such a perpendicular electric field can be controlled by varying the bias voltage between the tip and the sample.A finite element simulation of the electric field under the STM tip is presented in Supporting Information.In this case, an out-ofplane electric field lifts the energetic equivalence of the XM' and MX' domains, expanding the areas of those that have favored direction of polarization, thus, curving the domain walls. [12]n Figure 3a we experimentally demonstrate that for a sufficiently large bias voltage, reversing the direction of the perpendicular electric field resulted in significant shape changes of the domain walls, while the nodes of the domains remained in the same position.This is due to the topological nature of domain walls, which precludes nucleation of new domains, with the node density determined by the local twist angle.All changes in the domain pattern are produced by the displacements of existing domain walls, [19] whose ends remain pinned to the nodes of the initial domain wall network.In the data presented in Figure 3a we find that the brighter domains tend to shrink for the negative bias voltage at the expense of the darker domains; in contrast, for a positive bias voltage, the darker domains shrink at the expense of the brighter domains, who tend to enlarge.We quantitatively analyze the shapes of such equilateral triangular domains using a previously developed model [7] for the calibration of the electric displacement field between tip and sample to the applied bias voltage.We will subsequently use the resulting calibration to capture the behavior of domain redistribution by vertical electric field across a complex network of elongated domains (e.g., in area A of Figure 1c).
For the equilateral domains, experimentally accessible electric fields lead only to a small bending of the domain walls (partial dislocations) with ends clamped at the domain wall network nodes.To analyze this behavior, we used a previously developed analytical model [5,7] predicting shape of equilateral domain walls for a given field value based on an energy functional capturing the trade-off between energy costs on domain wall bending and energy gain, arising due to area expansion of preferentially polarized domains.The latter is provided by coupling of the ferroelectric polarization, P =  0 Δ (parametrised using the potential drop across the resulting double-charge layer) with electric field (parametrised using the displacement field controlled by external charges), − 2PE = − 2DΔ.This considers gain in energy due to growth of domain area with preferential orientation of of out-of-plane D (see equation (S1) in Supporting Information).The resulting shape of each domain wall can be described by its transverse deflection, y (x dx ′ , from a straight line, 0 ≤ x ≤ ℓ, connecting a pair of consecutive nodes.In this expression, the function f(x') is a root of a cubic polynomial, f 3 − Af − B (x −  2 ) = 0 where A = w w + 2 ≈ 3.52, with w = 1.05 eV nm −1 characterizing minimal single domain wall energy per its unit length (realized for armchair direction in WS 2 ) and w = 0.69 eV nm −1 taking into account orientation-dependent parts (energy of partial dislocations is maximal for zigzag direction).The value of electric displacement field, D, is incorporated into another parameter in the solution, B = 2DP  0 w , which determine coupling of the field to areal density of ferroelectric polarization P in XM (MX) domains shown in Figure 1a.
Using this model, we analyzed STM topographic maps (e.g., such as in area B in Figure 1b) that showed almost equilateral triangular domains at different voltages and at different length scales.By adjusting the value of the electric displacement field D for each analyzed domain wall, we obtained the curves represented by white lines in Figure 3a and b.We note that this model adequately captures the shape of domain bending at different electric fields (positive and negative) as well as at different lengths between the nodes of the network.Figure 3c summarizes the values obtained for the electric displacement field as a function of bias voltage by performing such fits of domain wall shape for five areas on two different samples and with two different tips.The relationship (see Figure 3c) between the bias voltage V B and the electric displacement field, D =  0 V B d * , gives us an effective distance, d* = 1.0 ± 0.7 nm, between the tip and the sample, which is close to a typical tip-sample separation in the reported STM experiment.
We note that the experimentally observed trend of the domain curvature increasing with the domain size ℓ is also a way to rule out electrostriction as a possible mechanism for the observed phenomena.In the scenario of electrostriction, the bending of the domains would either not be dependent on the direction of the electric field (if each layer couples individually to the electric field), or it would be smaller for smaller domains (if the bilayer couples to the electric field), as we detail in Supporting Information.
As we have pointed out earlier, even a small strain in one of the layers can distort the equilateral domains into elongated ones, such as those indicated in region A area of Figure 1c.Examples of such measured domains are highlighted in Figure 4b, where -for a strong enough electric field D>D * -the long stretches of the curved domain walls (which structurally are  , 0.050 e/nm 2 .We note the presence of defects (indicated by arrows), which served as a marker for the areas of interest. [27,28]The domain wall dynamics was realized in an imperfect crystal, where the domain walls were not pinned by defects, and where defects did not act as nucleation sites.In addition, the defects were found to be stable the same a partial screw dislocations) merge together into a "full" perfect screw dislocation (PSD), highlighted by red lines in Figure 4b.These PSDs separate pairs of identically polarized MX' domains.Hence, we extend the domain network modelling to deformed triangular domains, where two of the inter-node distances, ℓ 1,2 , are longer than the third one, ℓ 3 .For the nonequilateral domains, the anisotropy of DW energy is described by w + wsin 2 ( 1,2,3 + arctan(y ′ 1,2,3 )), where angles  1,2, 3 characterize deviations of DWs from the closest armchair direction at D = 0, (see SI), and arctan(y′ i ) appears due to a local transversal deflection, y i induced by electric field.The DW shapes are found by minimizing the energy functional elaborated in detail in Supporting Information with the following boundary conditions: y 1,2 ( cos( 1,2 )) = sin( 1,2 ) and y 1,2, 3 (ℓ 1,2, 3 ) = y 3 (0) = 0.The first of those accounts for the emergence of a PSD streak of length  (see Supporting Information), where we also define angles  1,2 .To find these PSD parameters, we consider that the transversal shifts y 1 and y 2 bring the DW together at the PSD end, x 1,2 =  cos( 1,2 ).For  1,2 < <1, justified for the strongly elongated domains in Figure 4a, the condition can be approximated by, where A 1,2 = A − 3sin 2  1,2 .Equations (2) enable us to describe all the details of the field-induced changes of the domain wall network, as marked on the maps in Figure 4.For a quantitative comparison, in Figure 4c we compare the theoretically computed lengths  of the PSD streaks with the measured values L exp for different domains and electric fields.

Conclusion
To summarize, we have demonstrated that for each marginally twisted sample the analysis of the ferroelectric domain network can give detailed information about the deformations in the assembled layers.These deformations lead to the variation of domain sizes across the network, producing a variety of switching conditions for individual domains.We demonstrate that one can achieve full local control over the domain network, by squeezing out the unfavorable polarization from the larger and elongated domains.This happens through merging pairs of partial dislocations into a perfect dislocation.Note that, despite an apparent disappearance of the unfavorably polarized domains, the surface polarization remains reversible, as, upon the reversed bias, perfect "full" dislocations split up into pairs of partial dislocation, nucleating domains of opposite stacking order.We also note that usually switching of ferroelectric polarization across a large volume/area of ferroelectric material is associated with hysteresis.In the context of interfacial ferroelectricity, such hysteresis has been observed (both in SEM imaging and transistor operation) in marginally twisted MoS 2 bilayers with extremely large domain sizes. [5]In that case, domain wall pinning occurred due to defects encountered by the domain wall propagation to long distances.In 0.2°− 1.5°twisted structures with relatively small domains, local deformation of domain walls, which ends are pinned to the XX sites of the domain wall network, explore local clean areas of high-quality material without pinning centers.As a result, no hysteresis appears, and we achieve reversible local control of ferroelectric dipole moments in the studied transition metal dichalcogenide bilayer.

Experimental Section
We computed the charge density in Figure 1a using density functional theory implemented in Quantum ESPRESSO, [29,30] constructing a supercell with a pair of mirror-reflected XM and MX bilayers separated by a large vacuum interval, as in Ref. [4] We exploited full-relativistic projector augmented wave pseudo-potentials with Perdew-Burke-Ernzerhof exchange correlation functional, kinetic energy and charge density cutoffs -90 Ry and 1000 Ry, respectively, and 23 × 23 × 1 reciprocal space grid.Parameters of WS 2 bilayer crystal structure were taken from. [31]he sample was assembled using the tear-and-stack technique. [32,33]S 2 crystals were purchased from HQ Graphene and were mechanically exfoliated onto Si/SiO 2 substrates.Monolayer WS 2 was identified optically, and the thickness was confirmed with atomic force microscopy (AFM).Suitable graphite flakes were identified and were picked up with a polydimethylsiloxane (PDMS) stamp covered with Polypropylene carbonate (PPC) film on a glass slide.The graphite on the polymer stamp was then used to pick up monolayer WS 2 .The two layers of WS 2 correspond to two halves of a large flake, allowing for precise rotational control.The PPC film was then cut and the 2D material stack was transferred onto a Si/SiO 2 substrate coated with Ti/Au (5 nm/150 nm).The sample was then annealed under ultra-high vacuum at 300 °C for several hours.The surface was cleaned using an AFM tip. [34,35]STM measurements were performed at pressures below 10 −9 Torr, at room temperature.For STS we used bias modulation of 5.0 mV, 1.325 kHz.

Figure 1 .
Figure 1.a) Top and side views of XM' and MX' stackings, with arrows indicating the direction of electric polarization.Inset shows in-plane-averaged charge density transferred between the layers, which is equal to antisymmetric (with respect to the middle plane between the layers) part of total charge density for XM-stacked WS 2 bilayer.b) Schematic of the STM experiment, indicating the electric field produced by applying a bias voltage between the tip and the sample.c) Large STM topographic map (V B = −1.3V, I T = 50pA), showing two distinct regions in the sample indicated by A and B.In area B, the triangular moiré superstructure is homogeneous, with equilateral triangles, and has periodicity given by the translation vectors ' t 1,2 .In area A, the moiré pattern is skewed due to a small strain in one of the layers (whose effect is magnified by the moiré superlattice effect[23] ).d) Spatial displacement field map for the inter-layer deformations, obtained using Eq.(1) from the variation of moiré superlattice pattern in c).

Figure 2 .
Figure 2. a) STM topographic map of a typical region with equilateral triangular stacking order domains (I t = 100 pA, V B = −1.3V). b) STS averaged over bright and dark domains (set point: I T = 60 pA and V B = −1.9V).The arrows represent energies at which the data in (c) was collected.c) dI/dV maps at the indicated bias voltages (I t = 100pA).

Figure 3 .
Figure 3. a) STM topographic maps acquired at I T = 100 pA and V B = ±1.3V.The white lines represent fits with string-like model using as parameter the electric displacement field values D = −0.075e/nm2 , 0.050 e/nm 2 .We note the presence of defects (indicated by arrows), which served as a marker for the areas of interest.[27,28]The domain wall dynamics was realized in an imperfect crystal, where the domain walls were not pinned by defects, and where defects did not act as nucleation sites.In addition, the defects were found to be stable, even when crossed by a moving domain wall.b) STM topographic maps of different areas with domains of indicated lengths (I T = 100 pA and V B = −1.4V).The white lines represent fits with the equilateral triangle domain wall model (from left to right: D = −0.047,−0.070, −0.053, −0.08 e nm 2 ).c) The electric displacement field as a function of bias voltage for five areas.Colors indicate different areas, squares correspond to Sample 1 Tip 1, circles correspond to Sample 2 Tip 2 (Supporting Information).
Figure 3. a) STM topographic maps acquired at I T = 100 pA and V B = ±1.3V.The white lines represent fits with string-like model using as parameter the electric displacement field values D = −0.075e/nm2 , 0.050 e/nm 2 .We note the presence of defects (indicated by arrows), which served as a marker for the areas of interest.[27,28]The domain wall dynamics was realized in an imperfect crystal, where the domain walls were not pinned by defects, and where defects did not act as nucleation sites.In addition, the defects were found to be stable, even when crossed by a moving domain wall.b) STM topographic maps of different areas with domains of indicated lengths (I T = 100 pA and V B = −1.4V).The white lines represent fits with the equilateral triangle domain wall model (from left to right: D = −0.047,−0.070, −0.053, −0.08 e nm 2 ).c) The electric displacement field as a function of bias voltage for five areas.Colors indicate different areas, squares correspond to Sample 1 Tip 1, circles correspond to Sample 2 Tip 2 (Supporting Information).

Figure 4 .
Figure 4. a) STM topography with observed angle of armchair direction indicated.Data acquired at V B = −1.3V, I T = 50pA.b) STM topography of a region in top-left of panel a).Images taken at I T = 100 pA and V B = ± 1.7 V. Single domain wall fits from string-like models plotted in white, predicted double domain walls plotted in red with A = A() where  is the angle of the triangle side relative to the nearest armchair direction.D values of −0.067and 0.070 e/nm 2 are used for the two images.As we already determined the effective distance that characterized the equilateral triangles (Figure3), in the case of the elongated triangles, drawing the shapes of the domain walls was done free of tunable parameters: the only input is the position of the network nodes.c) Predicted PSD length  compared to observed length for elongated domains (see Supporting Information for details on estimating L exp ) The error bars reflect the pixel size of our images. 2