Twirling and Spontaneous Symmetry Breaking of Domain Wall Networks in Lattice-Reconstructed Heterostructures of Two-Dimensional Materials

Lattice relaxation in twistronic bilayers with close lattice parameters and almost perfect crystallographic alignment of the layers results in the transformation of the moiré pattern into a sequence of preferential stacking domains and domain wall networks. Here, we show that reconstructed moiré superlattices of the perfectly aligned heterobilayers of same chalcogen transition metal dichalcogenides have broken-symmetry structures featuring twisted nodes (“twirls”) of domain wall networks. The analysis of twist-angle dependence of strain characteristics for the broken-symmetry structures shows that the formation of twirl reduces the amount of hydrostatic strain around the nodes, potentially weakening their influence on the band edge energies of electrons and holes.

This study addresses a detailed analysis of domain wall networks (DWN) which form in long-period moiré patterns characteristic for highly aligned heterostructures of same-chalcogen transition metal dichalcogenides MoX 2 /WX 2 (TMDs with X=S or Se).As same chalcogen TMDs have very close lattice constants, moiré patterns at their interface have long periods offering a sufficient space for creating preferential stacking areas (domains).That is, the energy gain due to better adhesion can surmount the cost of intralayer strain in each of the constituent crystals.The reconstruction of small-angle twisted bilayers into an array of domains [1] has been observed [2][3][4][5][6] both in MoS 2 /WS 2 and MoSe 2 /WSe 2 heterostructures.The observed [2][3][4][5][6] and theoretically modelled [1, 7,8] structures feature hexagonal for anti-parallel (AP) orientation of unit cell and triangular for parallel (P) orientation of unit cells DWN.
Here, we show that lattice relaxation in P/AP-MoX 2 /WX 2 (X=S,Se) bilayers in the limit of θ = 0 • twist angles undergoes symmetry breaking of the domain wall patterns resulting in formation of twirled structures of DWN nodes (see Fig. 1 for P-MoX 2 /WX 2 bilayers).The emergence of twirled DWN nodes is specific for heterobilayers, in contrast to previously studied marginally twisted homobilayers (MX 2 /MX 2 ) [1, 7,8], where only symmetric star-like nodes have been found.This difference stems from the hydrostatic strain component determined by a small lattice mismatch of the two constituent 2D crystals.By adjusting lattice constants of MoX 2 and WX 2 inside the domains, each monolayer compression/expansion is inflicted onto domain walls creating hot spots of hydrostatic strain in the DWN nodes with XX (chalcogen over chalcogen) stacking [9].An excessive energy costs of a large hydrostatic strain around XX nodes (as compared to shear deformations dominating DWN structure in homobilayers) is, then, negotiated by left/right-handed twirling of the domain walls, leading to a broken-symmetry configuration sketched in Fig.

1.
Methods.To study lattice relaxation in MoX 2 /WX 2 bilayers we use multiscale modelling approach developed in Refs.[1,10].This combines mesoscale elasticity with analytically interpolated description of adhesion energy of two layers, computed using density functional theory (DFT).The latter is represented below upon taken into account interlayer distance relaxation (see SM): Here, r 0 is a locally defined in-plane vector determining stacking arrangement between layers (r 0 = 0 for XX stacking, r 0 = (0, −a/ √ 3) for XW/2H stacking and r 0 = (0, a/ √ 3) for MoX/MoW stacking in P/APheterostructures), phases γ P = π/2, γ AP = 0 account r 0 → −r 0 symmetry of the adhesion energy in P-and 1,2,3 are three sets of the shortest reciprocal lattice vectors (|G where δ ≈ 0.2% for MoS 2 /WS 2 bilayers and δ ≈ 0.4% for MoSe 2 /WSe 2 ; u Mo (r) and u W (r) are in-plane displacement fields describing lattice relaxation in MoX 2 and WX 2 layers, respectively.Then, we minimise total energy of moiré superlattice, with respect to the displacement fields u W and u Mo .
Here, λ Mo,W , µ Mo,W are elastic moduli of MoX 2 and WX 2 monolayers (see Table I), and u l ij = (∂ j u l i + ∂ i u l j )/2 are components of strain tensors.
Table I.Adhesion [1] and elastic [11,12] energy parameters (in eV/nm 2 ) used in Eq. ( 3).X w To find optimal distributions of displacement fields, u W and u Mo , we solve a system of Lagrange-Euler equations, with periodic boundary conditions, implemented via the finite difference method.
In particular, we chose rectangular supercell with √ 3/2 and sides ( = a/ √ δ 2 + θ 2 is the moiré superlattice period, a is averaged monolayer lattice parameter), setting the following boundary conditions: Here, prime superscripts indicate that x Oy reference frame was rotated by an angle φ = π/6 + arctan ( √ 3δ − θ)/( √ 3θ + δ) with respect to the fixed frame for which Ox and Oy axes are along zigzag and armchair directions in the crystal, respectively.Such a rotation allows us to fix boundary conditions (note that the Lagrange-Euler equations themselves are rotational-invariant).
In the numerical analysis we use a sufficiently dense grid (one point per nm) to provide convergence of the computed displacement fields.This requires solving ∼ 30000 coupled nonlinear equations.This solution is ob-tained in an annealing-type computational scheme: first, we choose the starting point as u Mo/W = 0 and scale down the adhesion energy parameters w (s,a) 1,2,3 by a small factor η = 10 −3 , finding slightly relaxed structure.Then, we gradually enlarge this factor, up to η = 1 (at the final step), using solutions obtained at the previous iterations as starting points at each next step.
Twirls in P-MoX 2 /WX 2 heterostructures.Lattice relaxation in P-MoX 2 /WX 2 heterostructures results in formation of DWN separating triangular domains with MoX and XW stackings, similar to those found in 3R-TMD polytypes [2,3] (see Fig. 2).Minimising functional (3) we find the following three competing distributions of strain.One of them (symmetric), shown on the middle panels, is characterized by straight edge dislocation lines linking XX stacking nodes of DWN.The other two (left and right panels) are twirled (broken-symmetry) structures of DWN nodes with a left/right-handed twist.The computed total energies per supercell, gathered in Table II, show that the twirled DWN structure has lower energy than the symmetric one.Note that left/right-handed twirls have the same energy suggesting that DWN in perfectly aligned heterostructure (θ = 0 • ) undergoes a spontaneous symmetry breaking into twirled state.
Values of adhesive and elastic energies compared in Table II also suggest that the dominant energy gain for the twirled structures comes from lowering elastic energy.The latter is determined by hydrostatic strain, divu, and shear deformations, characterized by a vector, A = (u xx − u yy , −2u xy ).Their distribution across DWN is shown Fig. 2, in the form of color maps for divu and |A|.Here, the difference between symmetric and twirled structures is such that hydrostatic strain component, concentrated around nodes, is higher for symmetric one.Since energy costs of the hydrostatic strain (∝ λ+µ) are higher than those of shear strain (∝ µ), the twirled structures emerge as energetically favourable.
In Fig. 3 we analyse dependence of divu and A on the twist angle.The data shown in Fig. 3(i  that the maximal value of divu at the hot spot of strain decreases with the twist, which also explains why twirling is weaker for larger θ's.Another deformation field characteristic which is interesting to consider is

Mo(W) yy
).This characteristic determines the size of piezoelectric charges, generated by inhomogeneous strain in each layer, and pseudomagnetic field that would be experienced by charge carriers at the K-valley band edge in the heterostructure [1,10].Note that in P-heterostructures piezocharges (∝ e 11 B * ) have opposite signs in MoX Twirls in AP-MoX 2 /WX 2 heterostructures.Lattice relaxation in AP-MoX 2 /WX 2 heterostructures leads to the formation of hexagonal domains of 2H-like stacking (simultaneous metal-on-chalcogen and chalcogen-on-metal like in bulk 2H TMD crystals).Those domains are separated by a hexagonal DWN with two inequivalent nodes: one with XX stacking (chalcogen-on-chalcogen) and the other with two metallic site on the top of each other, MoW.Upon minimising energy in Eq. ( 3) for crystallographically aligned bilayers, we identify two candidates for the lowest energy configuration of DWN: symmetric and left/right-handed twirled structures, Fig. Energies of symmetric and twirled structures compared in Table III show that the broken-symmetry structure is energetically favourable.In Fig. 5 we show twist angle dependences of strain field characteristics around a single twirl.The data for divu in Fig. 5(i) indicate decrease of its magnitude (in each of the layers) with the twist angle.This trend corresponds to the decay of twirling with growth of layers' misorientation (Fig. 5(a-d)) as already discussed about P-heterostructures in the previous section.We also point out that, in contrast to P-bilayers, for APheterostructures shear strain results in the same signs of piezocharge densities (∝ e ).Note that DWN nodes host area of non-zero B * , which reverses sign at θ ≈ δ, with magnitude at MoW nodes substantially higher than that of twirled nodes (see Fig. 5(j)).
Conclusion.The main finding of this study consists in the prediction of a spontaneous symmetry breaking in the form of domain wall network in MoX 2 /WX 2 heterostructures.For almost perfectly aligned MoX 2 /WX 2 bilayers, this symmetry breaking consists in the formation of twirl-shaped nodes, accompanied by substantial reduction of a hydrostatic strain component in the vicinity of those nodes -as compared to earlier-studied symmetric (non-twirled) DWN structures [9].The reduction of hydrostatic strain component by twirling could have a pronounced effect in the energetics of the band edge states of the electrons and holes, making those less confined at the twirled nodes, as compared to symmetric nodes.As a result, the earlier studies [9] overestimated electron, hole and exciton bindings by the XX DWN nodes, so that the analysis of electron, hole and interlayer excitons localisation at the twirls requires farther studies.Graphene Flagship Project, EC-FET Quantum Flagship Project 2D-SIPC, EPSRC grants EP/S030719/1 and EP/V007033/1, and the Lloyd Register Foundation Nanotechnology Grant.Author Contributions.MAK, VVE, and VIF contributed equally to this paper.MAK and VVE developed formalism and software and carried out calculaions.VIF conceived the paper.All authors contributed to the data analysis and manuscript writing.

Figure 2 .
Figure 2. Maps of hydrostatic (a-c) and shear (d-f) components of strain tensor in principal axes for symmetric (b) and twirled (a,c) structures of reconstructed moiré superlattice in aligned P-MoX2/WX2 bilayers.Sketch on the top shows energy gain from formation of twirled structures.Scale bar is 30 nm and 45 nm for X=Se and X=S, respectively.Dashed rectangles on (c) and (f) show area displayed in Fig. 3.
2 and WX 2 layers, as B * Mo ≈ −B * W , and piezocoefficients have the same signs, e Mo 11 ≈ e W 11 .For the displayed range of angles, 0 • ≤ θ < ∼ 0.4 • , distributions of B * shows change of sign and zero value in the middle of the domain wall.For larger twist angles θ δ piezocharge and pseudomagnetic field distributions take the form of those established earlier for homobilayers [1].
Mo 11 B * Mo ≈ e W 11 B * W ) in both layers (opposite signs of piezocoefficients, e Mo 11 ≈ −e W 11 , due to anti-alignment of layers, are compensated by signs inversion for B * Mo ≈ −B * W

Figure 4 .
Figure 4. Maps of hydrostatic (a-c) and shear (d-f) components of strain tensor in principal axes for symmetric (b) and twirled (a,c) structures of reconstructed moiré superlattice in aligned AP-MoX2/WX2 bilayers.Sketch on the top shows energy gain from formation of twirled structures.Scale bar is 30 nm and 45 nm for X=Se and X=S, respectively.Dashed rectangles on (c) and (f) show area displayed in Fig. 5.

Table III .
Values of adhesion, elastic, total energies (in eV/supercell) computed for symmetric (sym) and twirled (twirl) structures and total energy differences, ∆, for aligned AP-MoX2/WX2 bilayers.Twirled structures form around XX nodes, whereas DWN around MoW nodes is almost unchanged (only rotated as a whole by ≈ 5 • ), which is because those nodes host metastable MoW stacking with much weaker hydrostatic strain component as compared to XX nodes).