Kekul´e spirals and charge transfer cascades in twisted symmetric trilayer graphene

.

Introduction.-Thediscovery of superconductivity proximate to correlated insulating behaviour [1][2][3][4][5][6][7][8][9] in twisted bilayer graphene (TBG) in the 'magic angle' regime [10,11] has stimulated intense investigation of a host of other multilayer graphene systems.Among these, alternating-twist multilayers exhibit identical moiré patterns between each pair of adjacent layers and have welldefined magic angles in the chiral limit [12].The simplest member in the group, twisted symmetric trilayer graphene (TSTG), which is composed of three layers of graphene stacked together with the inner layer rotated, has been experimentally studied [13][14][15][16][17][18][19][20][21][22][23][24][25], especially in the context of its superconducting properties, and has received extensive theoretical investigation .In the absence of external fields, the single-particle Hilbert space of TSTG can be decomposed into a direct sum of monolayer graphene and TBG parts [12], allowing many phenomena in TBG to be reproduced in TSTG.If one applies a perpendicular electric field, the TBG and graphene sectors are hybridized, allowing for further tunability of the system.Among the various correlated states in TBG, a novel form of translation-symmetry breaking order dubbed the incommensurate Kekulé spiral (IKS) has been proposed by some of us to characterize the many-body ground state at all integer fillings except charge neutrality [53], as well as a range of non-integer fillings [54].The signature of the proposed state, which is stabilized by experimentally realistic uniaxial strains, consists of a graphene scale Kekulé-like pattern with complex spatial dependence on the moiré scale, and has recently been observed with high-resolution scanning tunnelling microscopy (STM) [55].A similar pattern has also been observed in TSTG [25], but Kekulé spiral order in TSTG is theoretically unexplored, motivating the present work.
Here, we perform extensive Hartree-Fock (HF) simulations of TSTG that incorporate the effect of strain and allow for translation symmetry breaking.As in the case of TBG, we find IKS to be ubiquitous for non-zero integer and non-integer fillings of TSTG at finite strain.Furthermore, we find that commensurate Kekulé spiral (CKS) order in unstrained TSTG can be accessed even at zero strain by applying a large interlayer potential -a possibility absent in TBG.These results should be contrasted with previous studies of TSTG [36,40,44,50,52] at zero strain.We also explore the process of 'charge transfer' between the TBG and graphene sectors of TSTG at zero interlayer potential, and elucidate subtleties in defining commensurate fillings in light of this mechanism.
Model and Symmetries.-TSTG is made by stacking arXiv:2310.16094v2[cond-mat.str-el]15 Nov 2023 three sheets of graphene with the middle layer twisted by an angle θ relative to the top and bottom layers [Fig.1(a)].The moiré patterns generated within the top and bottom pairs of layers are identical, preserving the overall moiré periodicity of the system.The single-particle Hamiltonian of TSTG can be decomposed [12] into a mirror-symmetry even sector, which is effectively a renormalized version of TBG, and a mirror-symmetry odd sector, which resembles monolayer graphene [Fig.1(b)].The magic angles of TSTG and TBG are related via θ TSTG = √ 2θ TBG ≈ 1.56 • .In contrast to TBG, an electric displacement field normal to the layers has a significant effect: it produces an interlayer potential that breaks the mirror symmetry and mixes the odd and even sectors [Fig.1(c)].This modifies the dispersion, expanding the parameter space for new phases.
We also incorporate the effect of strain in the system.Homostrain (where all layers experience the same strain) has only a small effect on TBG [56], and we expect it to play a similarly minor role in TSTG.Therefore, we restrict our attention to heterostrain, produced by unequal strain in the different layers.In principle, strain could break the mirror symmetry of the system and introduce incommensurate moiré patterns for the top and bottom pairs of layers, potentially leading to a supermoiré pattern.However, we focus here on the case where top and bottom layers experience identical strains, thus preserving the moiré periodicity of the system.This choice can be justified as STM experiments observe well defined moiré patterns in large regions of the TSTG samples [16,17,25].We also ignore the small-angle Paulimatrix rotations in the Dirac kinetic terms in each layer; then, absent strain or interlayer potential, the system has exact particle-hole symmetry, which is lost when both perturbations are included.
Phase Diagram.-Weperform band-projected selfconsistent HF calculations on the model discussed above at the magic angle θ = 1.56 • .We use realistic hopping parameters w AA = 75 meV and w AB = 110 meV, and include strain ϵ and interlayer potential ∆V as free parameters of the model.We model interactions via the dual-gate screened Coulomb potential V (q) = (e 2 /2ϵ 0 ϵ r q) tanh qd, with a screening length d = 25 nm and relative permittivity ϵ r = 10.To avoid doublecounting interactions, we use the "average" subtraction scheme [57].Despite the lack of exact particle-hole symmetry for ∆V, ϵ ̸ = 0, we find that the phase diagrams at ν and −ν are qualitatively similar.Hence we only present results for |ν| ≥ 0 in the main text, and refer readers to Ref. [57] for ν < 0 results.
Unless otherwise noted, we assume collinear spin configurations, but we have also checked that, for the parameter range in Fig. 2, relaxing this condition does not lead to different ground states up to SU(2) K ×SU(2) K ′ valleydependent spin rotations.In order to capture the incommensurate Kekulé spirals (IKS) proposed in Ref. [53] or its commensurate counterparts, we allow for coherence between states of momentum k − τ q/2 in valley τ with states of momentum k − τ ′ q/2 in valley τ ′ (here τ and τ ′ take on values of ±1) for some Kekulé spiral vector q in our HF calculations, and then optimize over all q.One could further consider the slightly more general situation where the two spin sectors have different q.As such, the one-particle density matrix takes the form where s is spin index, a and b are single particle band indices, and q s is the Kekulé spiral vector with possible spin-dependence.In the majority of cases, this spin-dependence is not relevant and we simply scan over q ↑ = q ↓ ≡ q.Spin-dependence of Kekulé spiral vectors becomes relevant when both spin sectors have intervalley coherence and the two sectors are nonequivalent, as is the case for ν = 1.As such, we have allowed q ↑ ̸ = q ↓ when computing the phase diagram at ν = 1.For the range of interlayer potentials considered in the main text, we found that allowing for more generic translation symmetry breaking beyond Eq. 1 does not lead to different results.At higher interlayer potentials, we observe a complex pattern of translational symmetry breaking with nontrivial spin structure [57].Fig. 2 shows the phase diagrams for ν = 0, 1, 2, 3. Kekulé states involve non-zero Kekulé vector q and finite intervalley coherence.To diagnose the presence of these states and demonstrate their robustness, we show energy differences per unit cell between the lowest-energy q = 0 solution and the energy achieved by minimizing over all q: non-zero values indicate Kekulé spiral ground states.
At zero interlayer potential, the single-particle Hamiltonian decomposes into mirror-symmetry even (TBG) and odd (graphene) sectors.For all fillings and strains examined here, we observe no spontaneous breaking of this symmetry, and the two sectors are only coupled by Hartree terms.Compared to pure TBG, the inclusion of the graphene sector can lead to 'charge-transfer cascades' from the TBG sector, discussed below, but otherwise the phase diagram is consistent with the corresponding one for TBG [53].We find that IKS order appears at ν = 2, 3 at finite strain (Fig. 2).While for the range of strains shown in Fig. 2 we find Kramers intervalley-coherent (K-IVC) state at ν = 0, 1, for larger strains the system transitions into IKS at ν = 1 and a symmetric state at ν = 0, again consistent with the situation in TBG.
A non-zero interlayer potential mixes the two symmetry sectors, and the simple picture of Hartree-coupled TBG and graphene sectors is no longer applicable.The mixing of relatively flat TBG bands with graphene bands leads to more dispersive central bands.This favors the formation of IKS by reducing the minimal strain required for its stabilization, except at charge neutrality, where no Kekulé spiral is observed for either TBG or TSTG.At high interlayer potential of about 200 meV, a common feature except at ν = 3, is the breaking of C 2 T symmetry and the opening of a charge gap.While strong C 2 T -breaking has not been observed in monolayer graphene under ordinary conditions, as graphene orbitals are hybridized with TBG orbitals, the Dirac point becomes less dispersive and interaction-induced spontaneous symmetry breaking becomes possible.For ν = 3, C 2 T symmetry only spontaneously breaks in one of the two spin sectors, leaving the system gapless.We note that Ref. [21] finds a charge gap at ν = 2, in qualitative agreement with our result.
The phase diagram is most complex at ν = 1, where the two spin sectors have densities of approximately 0 and 1 electron per unit cell relative to their respective charge neutrality points (finite but small charge transfers between the two sectors could occur).The two sectors are only coupled with Hartree interaction, and otherwise exhibit their own physics.To capture the physics correctly, we allow for different Kekulé vectors for the two sectors [58].In particular, we find ground states with K-IVC (q ↑ = 0) in the spin sector at charge neutrality and IKS (q ↓ ̸ = 0) in the other spin sector.
Commensurate Kekulé Spirals.-Anew phase that we identify in TSTG is the commensurate Kekulé spiral (CKS*) [59] at high interlayer potential and zero strain.In CKS*, the spiral wavevector relates the K M -point of valley K (Fig. 1(c)) where TBG and graphene hybridization is most prominent, to the Γ-point of valley K ′ .As such, the resulting spiral order is commensurate with moiré periodicity.We note that CKS* further differs from IKS in that it breaks C 2 T symmetry.This is also true of other 'starred' states, including the KS* at ν = 1, 2, 3.At zero strain all KS* states have commensurate q (except ν = 3 and large ∆V ), but this commensurability is lost once any finite strain is intro-duced (though the graphene scale lock-in mechanism introduced in Ref. [25] may stabilize commensuration at finite strain).While CKS* is possible at all of ν = 1, 2, 3, it is only robust (i.e.significant energy advantage over competing states) at ν = 2, where it is the ground state up to about ∆V = 240 meV.
Charge-Transfer Cascades.-Atzero interlayer potential, the TBG and graphene sectors are coupled only by Hartree terms as the mirror symmetry is unbroken.While the total electron density is fixed, electronic charge could transfer non-trivially between the two sectors, as first noted in Ref. [36] at ν = 3 and large chiral ratio.We now track this charge-transfer mechanism across general fillings, and with finite strain.We use a HF interpolation scheme [36] to mitigate finite-size effects on the dispersion and hence more accurately compute the charge transfer.
Fig. 3 shows the filling of the TBG sector ν TBG as a function of total filling ν.We observe plateaus at integer values of ν TBG near ν = 2 and ν = 3.This can be rationalized in terms of the formation of a correlated insulator with a finite charge gap in the TBG sector, such that injection of a small amount of additional charge only changes the filling in the graphene sector.The precise width of the plateau depends on the size of the TBG charge gap and the density of states of graphene at the chemical potential.We also find that the sum of unsigned Fermi volumes of all bands -a rough proxy for how 'metallic' the system is -can reach a local minimum at non-integer total fillings, meaning that the most insulator-like filling does not necessarily occur at integer ν due to the charge transfer.
In generating Fig. 3, we have determined the chargetransfer from self-consistent HF calculations of the full system.However, if we treat the Hartree interaction between the two sectors as a function of only the densities of the respective sectors (which can be justified by assum-FIG.3. 'Charge transfer plateaus' in TSTG.We plot the fillings of TBG and graphene sectors from 10 × 10 self-consistent HF, interpolated to 60 × 60, as a function of total fillings.We notice plateaus at integer TBG fillings, which are due to a correlated charge gap in the TBG sector.In the inset, we additionally show the total unsigned Fermi volumes of all bands ('Fermi Volume') as a measure of the metallicity of the system.We notice that, due to non-trivial charge transfer between the two sectors, minima in 'Fermi Volume' need not occur exactly at total integer fillings.
ing the graphene sector charge density is approximately uniform), we can write where n denotes the total electron density and U 0 = V (q = 0).We can re-write this as where ẼTBG (n TBG ) = E TBG (n TBG ) − U 0 n 2 TBG /2 is the energy of TBG sector with the electrostatic energy of uniform electron gas subtracted away, and similarly for Ẽgraphene .We have also dropped terms that depend only on total electron density.Minimization of the total energy corresponds to dE/dn TBG = 0, giving μTBG = μgraphene . (4) The graphene sector can be modelled as a non-interacting Dirac cone, so that we can estimate the charge transfer from the interacting physics of TBG sector alone.The result is consistent with full TSTG HF calculations [57].
Comparison with Experiments.-Ref.[25] observes Kekulé patterns using STM in a TSTG device with heterostrain of −(0.12 ± 0.04)% and finds doping-dependent q Kekulé vectors, which are equivalent to our q vectors up to choice of reference point, at fillings from −2 to FIG. 4. (a) Doping-dependent q Kekulé vectors, found from 10 × 10 HF, of a system with ϵ = −0.12%and φ = 87 • , consistent with the experimental parameters in Ref. [25].We set ∆V = 0 as displacement field in the single-gate geometry is weak, and use single-gate screened Coulomb interaction V (q) = (e 2 /2ϵ0ϵrq)(1 − e −2qd ) with screening length d = 25 nm and relative permittivity ϵr = 20.We take the middle-layer graphene to be rotated counter-clockwise, as in experiment but opposite to our usual convention.(b,c) q Kekulé -dependence of the energies of IKS states for ν = −2 and ν = −2.5 respectively, where ∆E = E(q Kekulé ) − Emin.
−2.5.From our HF calculations with relative permittivity ϵ r = 20 (elsewhere in this paper, ϵ r = 10), at the experimental angle and magnitude of strain, we find that the ground states of TSTG for experimental parameters have IKS order, and we show the corresponding q Kekulé in Figure .4. We observe doping-dependent q Kekulé , in qualitative agreement with experiment.However, the values of q Kekulé are highly sensitive to modelling details, and so we do not expect precise quantitative agreement.
Concluding Remarks.-Motivated by the experimental discovery of Kekulé spiral order in TSTG [25], we have systematically studied its ground state phase diagram under strain and interlayer potential.At zero interlayer potential, since the system decomposes into TBG and graphene sectors, the physics is similar to that of TBG, including the emergence of IKS order at finite strain.At finite interlayer potential, we find new effects such as the breaking of C 2 T symmetry and the formation of a commensurate Kekulé spiral state.Our prediction that the Kekulé spiral state at ν = 2 breaks C 2 T symmetry and opens up a charge gap is in agreement with experimental results [21].The commensurate spiral state at zero strain is an addition to the family of Kekulé spiral states introduced by the discovery of the incommensurate Kekulé spiral, and hints at a broader relevance of KS order beyond strained samples.While the high interlayer potential required to stabilize the phase in TSTG necessitates a dual-gate device, precluding the direct detection of its Kekulé pattern using STM, the mechanism behind the formation of CKS order is more general, and could emerge in other moiré systems, especially in cases where C 2 symmetry is broken explicitly [60].We have further investigated the normal (nonsuperconducting) state of TSTG at non-integer fillings.We establish a charge-transfer cascade mechanism between the TBG and graphene sectors in the absence of interlayer potential, and show that this is consistent with a picture where the two sectors are coupled by Hartree interactions and the TBG sector density of states tracks that of standalone TBG.Notably, we find that integer fillings of the TBG sector and 'maximally insulating' densities need not coincide with overall integer filling, which has implications for experimentally mapping the phase diagram of TSTG.
Since the BM model band structures have already taken into account of some of the interaction, it is necessary to subtract away such effect.In practice, this amounts to replacing the density matrix P with P − P 0 in computing the interaction term, where P 0 is the reference matrix.A few different methods exist, but in this work we have employed the average scheme: two central bands per spin/valley have occupation 1/2, while bands above (below) the central bands are empty (filled).This is similar to the scheme employed in Ref. [36].

Symmetry Properties of States Found in Phase Diagram
In Table S1, we have summarized the symmetry properties of various states found in the main phase diagram.
. The symmetries of various states that appear in the phase diagram.Since all states considered here can be reduced to colinear spin configurations with appropriate SU(2)K × SU(2) K ′ rotations, we denote the filling of the spin sector where the state could occur as νs.For states that break UV (1) symmetry, an appropriate UV (1) rotation is allowed before C2T symmetry is considered.Asterisk(*) labels variants of states that break C2T symmetry (e.g.VP preserves C2T -symmetry but VP* breaks C2T -symmetry).Note about notation: Kekulé spiral states could either break or conserve C2T -symmetry.The C2T -breaking version is termed KS*, which could either have commensurate or incommensurate q.The C2T -conserving version is called IKS, as there is no commensurate version of it in the phase diagrams.†: At some regions, KS* states could have finite valley polarization and break time-reversal symmetry along the way.Refer to Figure 2 of the main text and the next section on phase boundaries for more discussion.

Competition between K-IVC and T-IVC at |ν| = 2
Ref. [64] predicts that the ground states of TBG at zero strain for fillings |ν| = 2 are K-IVC states, and Ref. [53,61] finds that this remains true at small enough strains.In this paper, we also find K-IVC to be the ground state for TSTG at small interlayer potential and strain at |ν| = 2.However, for ultra-low strained TBG samples at fillings near ν = −2, Ref. [65] found a Kekulé density pattern, despite the fact that such a pattern is expected to vanish in the K-IVC from symmetry reasons [66,67].One candidate strong-coupling state that could manifest a Kekulé pattern is the T-IVC state, which was previously predicted to be energetically unfavorable [64].To resolve this contradiction, some of us proposed [68] that electron-phonon coupling could provide a mechanism that favors T-IVC over K-IVC at |ν| = 2 based on self-consistent HF calculations [69].Ref. [70] proposes an alternative mechanism based on renormalization group arguments.
While the aforementioned experimental measurements were made in a TBG system, TSTG with small interlayer potential resembles strongly a TBG system.As such, it is possible that the ground states of TSTG at |ν| = 2 and at small interlayer potential and strain are actually T-IVC.
Phase Boundaries at ν = 2 At ν = 2, in the main phase diagram, we show the transition from K-IVC to KS* as we increase the interlayer potential at zero or small strain.In fact, we also observe intermediate phases of K-IVC* (which is the C 2 T -symmetry breaking variant of K-IVC) and KS*-VP (which is a KS* state with finite valley polarization) in our HF numerics, as shown in Figure S2.We believe that the order of transitions depends on the details of modelling, and in an actual experiment the intermediate phases may not be present.

Momentum Space Properties of Commensurate Kekulé Spiral
In Figure . S3, we show the strength of intervalley coherence (IVC) in the moiré BZ for the CKS* state.We note that the regions where IVC strength is significantly reduced correspond to the K M -point of K-valley or K ′ M -point of K ′ -valley.

Behaviours with Large Strains
Here we consider the behaviour with strains larger than those considered in the main text, up to 0.6%, at zero interlayer potential.As shown in Figure S4, at ν = 0, there is a transition from K-IVC to symmetric semi-metal when going from the unstrained system to a system with large strain of about 0.5%.The Kekulé vector q remains zero wherever IVC is non-vanishing.At ν = 1, there is onset of IKS from 0.4% strain onwards.FIG.S4.Systems with fillings ν = 0 and ν = 1 with strain up to 0.6% and zero interlayer potential.The results are obtained from 12 × 12 self-consistent HF.On the left, "IVC" is defined as ||P+−|| 2 /(N1N2), the Frobenius norm squared of the intervalley part of the density matrix divided by the system size.On the right, ∆E is the energy difference between the q = 0 state with the overall energy minimum.For the plot of ν = 0, we choose to retain, instead of all bands closest to charge neutrality, 2 TBG bands and 2 graphene bands respectively closest to charge neutrality per spin and valley.

Comparison of Two Methods to Compute Charge Transfer
As explained in the main text, the charge transfer between the TBG sector and graphene sector of a TSTG system can be understood in terms of matching the chemical potential of the two sectors.This allows the determination of the charge transfer with only Hartree-Fock calculations of TBG, where in addition the graphene sector can be modelled as simply non-interacting.The results of this method, labelled as 'TBG HF', are compared with full TSTG HF calculation, shown in Figure S5, and we obtain good agreement.The remaining discrepancy is likely the result of nonequivalent band cut-off implemented for TBG and TSTG calculations.

Particle-Hole Symmetry Breaking and Negative Fillings
As was shown in Ref. [31], TSTG with interlayer potential but without strain has exact particle hole symmetry, which is given by a combined operator m z C 2x P , where P is the usual particle-hole symmetry, m z is mirror reflection in the z-plane, and C 2x is π-rotation with respect to the x-axis.Now, applying strain would break C 2x symmetry for Hamiltonian of the TBG sector ( ĤTBG ), and break the combined symmetry for ĤTBG (strained graphene sector ĤD also breaks m z C 2x P symmetry, though the reasoning is less trivial as it starts without C 2x symmetry in the unstrained case).As such, there is no particle-hole symmetry for TSTG with both strain and interlayer potential [71].
It is easily seen that under C 2x symmetry, for the TBG sector, strain is mapped from (ϵ, φ) to (−ϵ, −φ), where ϵ is the magnitude of strain and φ is the direction of strain measured from the x-axis.It can be shown that under the combined m z C 2x P the strain of the graphene sector transforms an the identical way, establishing a correspondence between system at (ν, ϵ, φ) and system at (−ν, −ϵ, −φ), where ν is the filling factor as usual.
Lacking exact particle-hole symmetry, we computed the phase diagrams for the negative fillings as well, as shown in Figure S6.We see that there is no qualitative difference to that of the positive fillings.From the last paragraph, we understand that the results could also be applied to systems with positive fillings but negative strains.

Translation Symmetry Breaking at High Interlayer Potential
We perform self-consistent Hartree-Fock calculations allowing for all translation and spin symmetry breaking, and we observe the ground state to have complex spin textures at ν = 1 with large interlayer potential (> 200 meV) and finite strain.The system has significant spatial fluctuations of spin densities, with much smaller total charge density modulation, which may be an artefact of our numerical methods.The spin patterns partially resemble spin spirals, but the spin rotations are not co-planar, and the results show strong finite size dependence.As this is not the main purpose of our work, the result remains inconclusive, but this could be the subject of future studies with careful control of finite size effects.

FIG. 1 .
FIG. 1. Twisted symmetric trilayer graphene.(a) The large hexagons denote graphene Brillouin zones of individual layers and the small hexagon denotes the moiré BZ for valley K. (b) Without interlayer potential, TSTG decomposes into a TBG sector (red) and a graphene sector (blue) with opposite mirror eigenvalues.(c) A finite interlayer potential (100 meV) breaks mirror symmetry and hybridizes the two sectors.The color represents the relative weight of the wavefunction on the zero-field TBG and graphene sectors.(d,e) TBG sector conduction band at zero interlayer potential without and with strain (ϵ = 0.1%) respectively.Strain breaks C3-symmetry, unpinning the Dirac points, and increases the bandwidth.The corresponding band structures in valley K ′ can be found using time-reversal symmetry.

FIG. 2 .
FIG.2.Phase diagrams as a function of strain and interlayer potential at integer fillings ν.We have performed 10 × 10 selfconsistent Hartree-Fock calculations, scanning over all possible Kekulé spiral vectors q. ∆E is the energy difference between the q = 0 state with the overall energy minimum.The acronyms are: K-IVC for Kramers-intervalley coherence, VP for valley polarization, (I)KS for (incommensurate) Kekulé spiral, SM for (compensated) semi-metal.Asterisk (*) denotes states that break C2T -symmetry.The KS*(VP) region in ν = 3 consists of Kekulé spiral states with possible finite valley polarization.The red lines on the zero-strain axis indicate commensurate Kekulé spirals.The shaded regions have a finite charge gap.For ν = 1, we have shown the phases and charge gaps of the two spin sectors separately, where νs denotes the filling of each spin sector.As ∆E cannot be separately defined for each spin sector, the data are simply duplicated in the two plots of ν = 1.
FIG. S2.Properties and ground states at ν = 2 without strain from 10 × 10 HF results.The phases in each parameter range is labelled above the graph.
FIG. S3.Momentum resolved intervalley coherence strength of CKS* from 24 × 24 self-consistent HF.The blue dots mark KM /ΓM position (for K and K ′ valleys), and the red dot marks ΓM /K ′ M position.
FIG. S5.Comparison of fillings of TBG and graphene sectors obtained from two different methods.The results are from 10×10 HF calculations with strain of 0.3% and zero interlayer potential.
FIG.S6.The phase diagram at negative fillings.The acronyms are the same as that used in the main text, with non-IVC* denotes a C2T -breaking state that preserves other symmetries.Despite lack of exact particle-hole symmetry, there is little qualitative difference to the phase diagrams of positive fillings.Due to the rather coarse grid used in parameter space, the phase boundaries should be interpreted as a rough guide.