Electrically tunable correlated domain wall network in twisted bilayer graphene

We investigate the domain wall network in twisted bilayer graphene (TBG) under the influence of interlayer bias and screening effect from the layered structure. Starting from the continuum model, we analyze the low-energy domain wall modes within the moiré bilayer structure and obtain an analytic form representing charge density distributions of the two-dimensional structure. By computing the screened electron–electron interaction strengths both within and between the domain walls, we develop a bosonized model that describes the correlated domain wall network. We demonstrate that these interaction strengths can be modified through an applied interlayer bias, screening length and dielectric materials, and show how the model can be employed to investigate various properties of the domain wall network and its stability. We compute correlation functions both without and with phonons. Including electron–phonon coupling in the network, we establish phase diagrams from these correlation functions. These diagrams illustrate electrical tunability of the network between various phases, such as density wave states and superconductivity. Our findings reveal the domain wall network as a promising platform for the experimental manipulation of electron–electron interactions in low dimensions and the study of strongly correlated matter. We point out that our investigation not only enhances the understanding of domain wall modes in TBG but also has broader implications for the development of moiré devices.

By identifying relevant degrees of freedom at low energy, the models provide a bosonic description for the strongly correlated twisted bilayer systems [21, 66-68, 70-72, 74], offering an alternative perspective to other theoretical works on these systems [31,34,.Notably, this bosonic description integrates with the renormalization-group (RG) analysis, serving as an effective tool for exploring correlated phenomena, including the correlated insulating phase and superconductivity (SC) observed in earlier studies [98,99] and various unconventional states of matter or features in subsequent works .While these observations are typically associated with the quasiflat bands in devices near a certain twist angle (referred to as the magic angle), here we explore a regime where the system exhibits correlation effects even when away from the magic angle.
Crucially, the moiré systems provide an opportunity to systematically determine the interaction strengths in devices consisting of one-dimensional channels.This is a notable contrast to earlier studies on sTLL and csTLL in non-moiré systems [5][6][7][8], which often relied on ad hoc assumptions regarding the specific forms of the interaction strength.Our investigation is also inspired by recent studies [49,50,73,113] that have shown the tunability of interaction strength in graphene-based devices.In [113], a control of electron-electron interactions was achieved by altering the distance between TBG and a metallic screening layer, observing a suppression of correlated insulating phases in devices with a smaller screening distance.Additionally, domain wall modes in Bernal-stacked bilayer graphene can be manipulated through voltage gating [49], which subsequently alters the properties of TLL [50].Similar domain wall modes appearing in large-angle TBG can also be electrically controlled [73].These insights further substantiate the systematic controllability of the interaction strength in the domain wall network.
In this article, we investigate the domain wall network in TBG under the influence of interlayer bias.We compute the energy dispersion and density profile of the low-energy domain wall modes within the moiré bilayer structure.By globally fitting the charge density in the two-dimensional structure, we derive an analytical expression representing charge density distributions.This allows for efficient calculations of electron-electron interaction strengths within and between domain walls.Our investigation extends to how these interaction strengths can be modulated by externally controllable parameters, highlighting the potential for experimental manipulation through device design and preparation.Building on this model, we explore the spectroscopic and transport properties of the network, uncovering potential Anderson localization in moderately disordered devices.We also examine the instability of the network towards various density wave (DW) states and SC in two scenarios, without and with phonons.In purely electronic systems, repulsive electronelectron interactions favor the formation of charge density waves (CDW) and spin density waves (SDW).However, longitudinal acoustic phonons can enhance pairing instability within the domain wall network, potentially driving these systems towards SC.In the presence of very strong electron-phonon coupling, we identify the Wentzel-Bardeen (WB) singularity that eventually destabilize the network.Our findings establish the TBG network as an electrically tunable, low-dimensional platform for the systematic investigation of strongly correlated phenomena.
The article is structured as follows.In section 2, we introduce the continuum model in the presence of an applied interlayer bias.In section 3, we calculate the spatial charge density distribution from the continuum model and derive an analytic form via global fitting, examining the influence of external parameters on this distribution.Section 4 focuses on computing the effective interaction strength for both intrawire and interwire terms, including their dependence on interlayer bias, screening length, and dielectric materials.In section 5, drawing insights from the single-particle model, we construct an interacting model for the correlated domain wall network.Here, we investigate its spectroscopic and transport properties, along with quantifying the localization length and temperature for the Anderson localization of the network.In section 6, we examine the instability of the network towards DW and SC, either in purely electronic systems in section 6.1 or in the presence of phonons in section 6.2.Our findings and their broader implications are discussed in section 7.
Appendix A provides an analysis on the influence of the chemical potential variation on the domain wall network.Additional details on the effective action and correlation functions are provided in appendix B for purely electronic systems and in appendix C in the presence of phonons, respectively.

Continuum model of TBG in the presence of interlayer bias
The continuum model of TBG, originally developed in [26], has been widely employed to study the band structure of the system at the single-particle level [29,31,32,34,62].As a starting point, we incorporate the interlayer bias in this model [30], with H 0 describing the top and bottom graphene layers with a twist angle θ and H hyb the hybridization The vectors q j link the Dirac cones between layers.The vectors q M1 = q 0 − q 1 and q M2 = q 1 − q 2 are the reciprocal lattice vectors of TBG with magnitude qM = √ 3k θ and k θ = 4π/(3λM).
between the two layers.For given valley and spin, this single-particle Hamiltonian can be expressed as ⊺ with the fermion field c η αγσ , the transpose operator ⊺, and the indices corresponding to the layer η ∈ {1, 2}, sublattice α ∈ {A, B}, valley γ ∈ {K ≡ +1, K ′ ≡ −1}, and spin σ ∈ {↑, ↓}.In the above, the diagonal blocks incorporate the interlayer bias 2V d and contain the following term representing the Dirac cone at the γ valley, where we have the Fermi velocity v F for monolayer graphene, the momentum operator p = p x , p y , and the location K γ of the Dirac point in momentum space.We denote τ γ = (γτ x , τ y ) with the component µ ∈ {x, y, z} of the Pauli matrix τ µ acting on the sublattice index.The off-diagonal blocks in equation (1b) describe the spatially dependent interlayer hybridization and act as a periodic moiré potential with moiré wavelength λ M defined in the caption of figure 1(a).
Setting the origin at one of the AA-stacking region centers, we have T γ hyb = 2 j =0 e i γq j •r T γ j , with the twoby-two matrix )] . (3) Here, we have q j = R 2π j/3 (0, −k θ ), with the rotation operator R ϕ and k θ defined in figure 1 caption.We account for the effect of lattice relaxation by allowing for different amplitudes for hoppings between identical (w AA ) and distinct (w AB ) sublattices between the layers [29,31,32,62,127].Following [32], we adopt an value for the ratio w AA /w AB between 0.2 and 0.5, assuming moderate relaxation effect.As the formulations in [26,33,34], we define the dimensionless parameters α AA = w AA /(hv F k θ ) and α AB = w AB /(hv F k θ ).We will term the latter as the effective hybridization parameter, which plays a key role in the band structure.
As in [26,32,33,128], we proceed by expressing the single-particle Hamiltonian in equation ( 1) in a plane wave basis, implementing a momentum cutoff at the scale determined by the larger of α AB q M and V d q M /(hv F k θ ) with q M specified in the caption of figure 1.The resulting expression allows us to perform numerical diagonalization.In addition to the band structure, the model introduced here allows us to compute the spatial distribution of the charge density, which we investigate next.

Charge density distribution under interlayer bias
Motivated by domain wall network discussed in [30], we explore the gapless modes in the lowest two bands near the charge neutrality point around the K and K ′ points in the original first Brillouin zone of monolayer graphene.More precisely, for a given valley (for instance, K), we focus on the low-energy modes around K and K ′ points in the mBZ, as indicated in figure 1(b).Aiming at exploring domain wall modes, we set V d to nonzero values and focus on the parameter range where V d /hv F k θ ≳ α AB > 1 for the rest of the article.By numerically diagonalizing the single-particle Hamiltonian, we compute the wave functions and subsequently obtain the spatial distribution of the charge density.An example is illustrated in figure 2, where we display the computed charge density alongside the corresponding state in the mBZ.It can be observed that the charge density is predominantly concentrated along the domain walls in the direction of propagation.Furthermore, within these domain walls, the density exhibits higher values around the AA-stacking regions.This latter feature aligns with findings in previous studies [31,32], but not well described by the simple approximations in [30].We note that a smaller twist angle leads to a shift of density from the AA-stacking regions to the domain wall segments, resulting in a more uniform distribution along the entire domain wall [43].However, this smaller twist angle also increases α AB , thereby requiring a higher momentum cutoff in our numerical diagonalization.For illustrative purposes, we mainly use a larger twist angle θ = 0.5 • , unless specified otherwise.As depicted in the main panel of figure 2, in the momentum space, for each valley and spin, there are two energy band branches hosting these lowenergy modes.Reversing the valley inverts the velocity of the modes, indicating electrons propagating in the opposite direction.This feature is consistent with findings in TBG [30] and also noted in nontwisted samples [52][53][54].The charge density profiles presented in the insets of figure 2 correspond to µ = 0. We have verified that the results remain largely unchanged when varying the chemical potential, provided it stays within the narrow bands depicted in the main panel.A more detailed analysis on the effects of chemical potential on the density distributions of the domain wall modes is provided in appendix A. It is important to note that our findings apply to systems with a rather small twist angle, which contrasts with systems near the magic angle, where the chemical potential can induce nonnegligible effects [85,[129][130][131][132][133][134][135][136].
To resolve the domain wall modes more clearly, we focus on one of the domain walls along the y direction.In figure 3(a), we show the spatial dependence of the charge density perpendicular to the domain wall for several V d values.In addition to confirming the confinement at the domain wall, we observe that the density profile can be electrically modified.Namely, as the interlayer bias increases, the charge density becomes more confined within the domain wall, as a result of the local spectral gap induced by the bias at Bernal-stacking domains [30,137,138].This feature can be examined more quantitatively.In figure 3(b), we plot the density at the center of the AB-stacking region (labeled by r AB ) as a function of the interlayer bias, observing a decrease in density at r AB due to the increasing local gap with the interlayer bias.To illustrate the density shift with respect to the interlayer bias, we also show the ratio of the density at the center of the domain wall segment (labeled by r dw ) to that at r AB .As the local gap within the domain regions is enhanced by increasing interlayer bias, the reduction in the density at r AB leads to increase in the density at r dw .Since the interactions within the domain wall network are governed by the screened Coulomb potential, which involves essentially short-range densitydensity interactions, this feature suggests the tunability of these interactions through external electrical control.
In addition to locating the low-energy modes in real space and determining their spatial profile, we are also able to extract the bandwidth ∆ a and the Fermi velocity v dw of the domain wall modes through the dispersions of domain wall modes in the spectrum.We summarize the results in figure 4, where we show how the two quantities influenced by the interlayer bias and effective hybridization parameter (quantified by α AB ).Again, the dependence on the interlayer bias demonstrates the electrical tunability of the domain wall network.As V d increases, the bands around the charge neutrality point flatten, and the bandwidth ∆ a decreases, a feature also observed in [45].Crucial for our subsequent analysis, this leads to a decreased v dw value.On the other hand, the dependence on the effective hybridization parameter indicates variations in these quantities across different samples, though with a weaker dependence in the studied range.Specifically, since α AB is influenced by the hybridization strength w AB and the moiré wavelength λ M , the properties of the network can be altered through device preparation, either by adjusting w AB through applied pressure (experimentally achieved in [106]) or by varying λ M through the twist angle.Having established the presence of domain wall modes, we introduce a fitting function in order to express the density using an analytic form, which will be useful when we estimate the interaction strength below.Specifically, we take the following functional form, with the normalization constant C 0 , domain wall length L ∥ , the extent L z of the wave function perpendicular to the layered structure, the Heaviside function Θ(x) and the index m labeling the domain wall for a given domain wall direction indicated by j.Here, we define the index 1 δ ∈ {1, 2} to denote the branches of the gapless modes in the spectrum, which are represented by the blue and yellow curves in figure 2.
In the above, we adopt a functional form that allows for the separation of the two orthogonal directions locally defined through the domain wall orientation.The parameter κ ⊥ is introduced to characterize the confinement perpendicular to the domain 1 However, our numerical analysis indicates no practical differences in the spatial profile of ρ dw δ,m between the branches δ = 1 and δ = 2. Consequently, we refrain from introducing redundant labels for the parameters κ ⊥ , κ ∥ , and c 0∥ .Therefore, the primary distinction between the two branches lies in their wave functions in momentum space, which will be distinguished using different Fermi momenta in the subsequent sections.
wall.The form of ρ ⊥ is inspired by the Jackiw-Rebbi solution as described in [30].Partially motivated by the overall C 3z symmetry of TBG [33], we adopt an empirical formula ρ ∥ of a form similar to the transverse component to describe the density profile along the domain wall, with a higher density near the AAstacking regions characterized by the parameter κ ∥ .In the above expression, κ ⊥ , κ ∥ and c 0∥ serve as the fitting parameters, which can be determined by globally fitting2 the two-dimensional density obtained from the continuum model to ρ 2D (r).
An example of such fitting is illustrated in figure 5, where we show the computed density profile (blue dots) and the fitting results (red curves) along two directions, one transverse to the wall (in the x direction), and the other along the domain wall (in the y direction).We see that the expression ρ ∥ along the domain wall fits well, including the peaks around the AA-stacking regions.Perpendicular to the domain wall, the fit is reasonably good, although some discrepancies are observed.It is important to note that these discrepancies are of minor significance in the two-dimensional global fitting due to the smaller overall magnitude compared to the AA-stacking centers.
Through our fitting procedure, we establish how the system parameters influence the fitted values, as illustrated in figure 6.The local spectral gap in the Bernal-stacking regions, which increases with the interlayer bias V d [137,138], leads to a corresponding increase in both κ ⊥ and κ ∥ (approximately linearly).This observation further confirms the enhanced confinement within the domain wall and AA-stacking regions.Consequently, these findings support the  idea that it is possible to manipulate charge properties and interaction strength within the domain walls using externally controllable parameters.

Interaction strength in the domain wall network
With the analytic form of the spatial density profile and its parameters deduced from the fitting procedure, we now compute the effective interaction strength for both intrawire (within a domain wall) and interwire (between parallel domain walls) contributions.The empirical formula in equation ( 4) allows for a straightforward separation of variables parallel and perpendicular to the domain walls.To proceed, we analyze a device configuration where the TBG is positioned on top of a dielectric layer with thickness d (setting z = 0 at the interface).This layered structure is then situated above a metallic back gate, located in the region z ⩽ −d, as depicted in figure 7. The presence of the metallic gate results in a screened Coulomb interaction, making it effectively short-ranged.Consequently, the thickness of the dielectric layer serves as the screening length, quantifying the degree of the screening effect.
In general, screening effects in such a three-layer structure can be formally incorporated through an infinite sum of image charges [140,141].Here, since in our setup the electron wave function in the TBG has a narrow distribution in the direction perpendicular to the layered structure (more quantitatively, L z ≪ d), we simplify our problem by considering an image charge distribution (caused by the charges in TBG) in −L z ⩽ z ⩽ 0 with a screening factor of −a scn [142] and another distribution (caused by the TBG and dielectric layer) in the metallic layer at  6), with the curves serving as a visual guide.Unless otherwise specified, we adopt L ∥ = 0.5 µm, Lz = 6.68 Å (twice the interlayer separation of graphite), and ϵr = 6.93, corresponding to a hexagonal Boron Nitride (hBN) dielectric layer [139] for the rest of this article.The adopted values of the other parameters are given in the caption of figure 2.
Here, we introduce the parameter, a scn = (ϵ r − 1)/(ϵ r + 1), with the relative dielectric constant ϵ r .With this approximation, for the distribution ρ dw δ,m of the branch δ in a given domain wall labeled by m, we have the following expressions for the image charge density, The intrabranch and interbranch electrostatic energy between the mth domain wall and (m + n)th domain wall can be expressed as where M ∈ {dw, diel, met} denotes different contributions and δ denotes the branch opposite to δ.
As mentioned earlier, as we do not find numerical difference between ρ dw 1,m and ρ dw 2,m , the above integrals are practically the same for any δ.To proceed, we reformulate the denominator in the integrand of equation ( 6) using a Gaussian integral [143][144][145].Using this approach, the integral over z becomes straightforward owing to the uniform distribution in this direction.In contrast, the integrals over the in-plane spatial coordinates will be evaluated numerically.
The main results of this numerical procedure are summarized in figures 8 and 9, for the intrawire (U ee,n=0 ) and interwire (U ee,n̸ =0 ) terms, respectively.For the intrawire interaction strength in figure 8, we observe an increase with a larger interlayer bias due to the stronger confinement of the electron wave functions in the domain walls, consistent with the computed density profile in section 3.In addition, there is a decrease in U ee,0 with a larger dielectric constant and a shorter distance from the metallic layer owing to the stronger screening effect.Notably, the interaction strength can be enhanced by more than seven fold in the investigated parameter regime, demonstrating great electrical tunability.
As to the interaction strength between electrons in parallel domain walls, we have calculated the interwire terms up to ninth neighboring domain walls, as shown in figure 9.The dependence of U ee,n for a given n on interlayer bias and screening length follows the same trend as the intrawire term.For different n, while the interwire interaction is more severely screened by a shorter d for domain walls that are farther apart (see figure 9(a)), the V d dependence behaves similarly (see figure 9(b)).The quantitative differences observed between the two panels are expected.Specifically, the interaction between the nth nearest neighbor domain walls is significantly reduced when the screening length is comparable to their separation distance (that is, d ∼ √ 3nλ M /2).This reduction in interaction strength is less significant when d is much larger.On the other hand, the dependence of the interaction strength on the bias voltage originates from the confinement of electron wave functions induced by V d .Consequently, this control parameter influences the strength uniformly across different values of n, resulting in a minor dependence on the latter.Since the strength of the interwire interaction decays with the separation between domain walls, as expected, it allows for the exclusion of interwire interaction terms at large n in the construction of the interacting network model.The rate of decay changes with the distance d from the metallic gate, with a smaller d value leading to a faster decay.
Before concluding this section, we note previous experimental studies that have demonstrated tuning the interaction strength in Bernal-stacked bilayer graphene through V d [50], as well as in unbiased TBG by varying d [113].With the demonstrated tunability here, we expect future efforts towards systematic investigations in biased TBG, either by utilizing the aforementioned control methods or by exploring the use of various dielectric materials.Beginning with the continuum model, we have analyzed the charge density and the effective interaction strengths within the domain wall network.This analysis provides the foundation for constructing our model describing interacting electrons in the network, which we introduce in the following section.

Correlated domain wall network
With the interaction strengths within the network determined, we are set to establish an interacting model for the domain wall modes.This model distinctively considers two branches of gapless modes for each direction of motion and therefore differs from the previous bosonization models for moiré networks [21,67,68,72], which neglected these additional degrees of freedom.With our description, each domain wall is reminiscent of the bosonization model for metallic carbon nanotubes [146][147][148], where the doubling in degrees of freedom is due to the two Dirac cones in the spectrum.Here, the situation is further complicated by the formation of the moiré network, leading our model to effectively describe interacting electrons in coupled nanotubes forming a mesoscopic network.
To be explicit, we describe electrons as fermion fields ψ ( j) σ,m in each domain wall in figure 10 and bosonize them as follows, with the array index j ∈ {0, 1, 2} corresponding to the three q j directions introduced in figure 1 with the index ν ∈ {c, s} labeling the charge/spin sector and P ∈ {s, a} for the symmetric/antisymmetric combinaiton of the two branches.The above relation implies that ϕ j νP,m and ∂ x θ j νP,m are conjugate to each other.
In terms of the boson fields, our model for the correlated domain wall network takes the form, where ϕcP,n and U ( j) θcP,n quantify the strengths of the density-density and current-current interaction terms between the nth neighbor domain walls, respectively.Their specific forms are given by with ee,n and V ee,n computed from equation ( 6).Our numerical analysis shows U ee,n within the relevant parameter regime, suggesting that the charge antisymmetric sector (labeled by 'ca') is essentially noninteracting.However, we retain its notation for generality.Finally, for the spin sector, we have identical velocity u s = v dw /K s and interaction parameter K s for all the domain walls.
While our description above preserves the C 3z rotational symmetry, leading to an independence from the j index, we retain this label for generality.This notation can be useful in scenarios where spatial inhomogeneity or anisotropy arises, such as from disorder or spontaneously broken symmetry.On the other hand, without imposing the translational symmetry in the direction perpendicular to the domain walls, there is a visible dependence on m, which we will demonstrate in the following sections.The application of bosonization in our model enables the computation of physical quantities nonperturbatively in the interaction strength.
To proceed, we introduce orthogonal matrices M ( j) ϕcP and M ( j) θcP that diagonalize the charge sector in equation (9b), leading to where the interaction strength and fields in the new basis read and similarly for U  8), provided that we have M θcP .For the sake of notational clarity, we will continue to use the distinct notations in the following discussion.It should be noted that the stability of the model presented in equation ( 9) is ensured by the positive definiteness of the interaction matrix, which can be inferred from equation (9b).Additionally, this stability can be confirmed by verifying that all values of U ( j) ϕcP,m in equation (12a) are positive.Indeed, we have numerically verified that these criteria are met across the entire range of parameters we have examined in this work.In section 6.2, however, it is shown that including phonons in the system can lead to a scenario where strong electron-phonon coupling destabilizes the model, resulting in the WB singularity.
Having established the bosonic model for our correlated domain wall network, we now explore its properties.To this end, we investigate physical quantities, including the local density of states (DOS), impurity-induced conductance correction, and tunneling current between parallel domain walls or nonparallel domain walls at intersections.In addition, we look into the Anderson localization induced by potential disorder in these domain walls.

Local DOS
The local DOS at the position x in the mth domain wall of the jth array can be computed by generalizing the calculation for one-dimensional systems [149,150].Keeping the forward scattering contributions, which give the signature universal scaling behavior, we have where we define ⟨• • • ⟩ ee with respect to the effective action in equation (B.1).With the bosonization formula in equation ( 7), we compute the local DOS, which can be expressed as a function of energy E and temperature T, with the Boltzmann constant k B and the interactiondependent parameters defined as We see that, due to the coupling between the domain walls, the local DOS of a given domain wall depends on the intrawire and interwire interactions of the entire network.The local DOS in equation ( 14) follows a universal scaling curve and serves as spectroscopic features that can be verified through scanning tunneling spectroscopy as illustrated in figure 11(c).It can be checked that, at low T, the formula reduces to a power law ρ m .It is noteworthy that in the noninteracting limit, both ∆ ( j) ϕcP,m and ∆ ( j) θcP,m approach unity, while β ( j) m tends towards zero.These deviations thus serve as an experimentally accessible metric for quantifying the correlation in the network.
To understand how system parameters influence the above exponents, we calculate their dependence on the screening length and interlayer bias for a domain wall at the center of a mesoscopic device, and present the results in figure 11.Our numerical analysis reveals significant deviations from unity for both ∆ ( j) ϕcP,m and ∆ ( j) θcP,m in most regions of the plot, highlighting the presence of substantial correlation effects that can be observable through spectroscopic probes.As anticipated, the deviations become more pronounced with increased screening length and interlayer bias.This is particularly evident in the dependence of interlayer bias, which we attribute to the substantial reduction in v dw with increasing interlayer bias (see figure 4).
Besides the domain wall located near the center of a device, we extend our analysis to various domain walls denoted by different m.We focus on the interlayer bias dependence, as depicted in figure 12. Notably, since our approach does not assume translational invariance in the direction perpendicular to the domain walls, it leads to weak but visible quantitative differences between domain walls located at the edges and those in the interior of the twodimensional system.We observe that both ∆   potential remains within the low-energy bands under investigation.We conclude this subsection by noting that a systematic investigation of the local DOS at domain wall segments (away from the AA-stacking regions) using scanning tunneling spectroscopy could be highly useful in quantifying correlations within the domain wall network.

Charge transport features at microscopic scales
Here we discuss the charge transport features from three microscopic processes, which we illustrate in figure 10.To begin with, we consider an isolated impurity in the mth domain wall of the jth array, illustrated in figure 13(a).It can be modeled as a delta function V imp (x) = v b δ(x), which induces the backscattering of the electrons within the domain wall, ] .
(16) In the above, we have intrabranch terms involving ϕ j ca,m field and interbranch terms involving θ j ca,m field.While their scaling behaviors are identical in our case with a noninteracting charge antisymmetric sector, in a more general situation the two terms might scale differently.We therefore focus on the scaling of the intrabranch terms, as these can be the most RG relevant term in either scenario.From equation (16), one can obtain the RG flow equation for the backscattering strength [4,151,152], with v b = v b /∆ a and the dimensionless length scale l.It leads to a correction in the (differential) conductance, with an interaction-dependent exponent.
Next, assuming weak coupling between nonparallel domain walls, the intersections at the AA-stacking regions allow for tunneling process between electrons in two of the crossing domain walls, say, the mth domain wall of the array j and the m ′ th domain wall of the array j ′ ; see illustration in figure 13(b).To the leading order, the tunnel current between these domain walls can be obtained by generalizing the formula in one-dimensional systems [153,154], with the field operators defined at the coordinates of the intersection.Here, the fields ψ σ,m and ψ σ,m ′ belong to the two domain walls involving the tunneling process, the tunnel amplitude t ( j,j ′ ) ×,mm ′ , and the voltage difference V between the two domain walls.
Generalizing the algebra in [155] for networks, we get the current-voltage curve for general V and temperature T, which involves the local DOS of the two crossing domain walls.Alternatively to the current-bias curves, one can derive the RG flow equation for the tunnel amplitude, where t ×,mm ′ /∆ a is the dimensionless coupling.From the RG flow equation, we obtain the differential tunneling conductance in the high-bias or high-temperature regime, consistent with equation ( 20) in these limits.Finally, we can also discuss the tunneling between two parallel wires within a given array j, as plotted in figure 13(c).In the high-V or high-T regime, we get the following power laws for the differential tunneling conductance, where we incorporate the tunneling taking place at any spatial points along the parallel domain walls [63,157].
While a comprehensive analysis of the transport properties of the entire TBG device requires further exploration at the mesoscopic level, similar to previous studies in the single-particle picture [158][159][160], we can examine the exponents characterizing three distinct microscopic processes relevant to transport.To this end, we demonstrate how α

and α ( j)
⊥,mm ′ vary with experimentally controllable parameters, as shown in figures 13(d)-(f), considering a noninteracting spin and charge antisymmetric sectors (that is, ∆ θca,m = K ss = K sa = 1 for any j and m).Once again, we observe that all the exponents are electrically tunable through the interlayer bias and screening length, with a more pronounced dependence on the former.This is consistent with the behaviors of ∆ ( j) ϕcs,m and ∆ ( j) θcs,m discussed in figure 11.Specifically, a larger interlayer bias leads to stronger interactions, resulting in a larger deviation of the exponents from their noninteracting values.Among the three exponents, the behavior of α As demonstrated in this work, electric transport at microscopic scales can be influenced by external parameters, which should consequently affect the transport properties of mesoscopic devices, in which the three explored microscopic processes take place.A detailed analysis of these mesoscopic networks should include the TLL-induced correlation effects discussed here, along with the magnitudes of the bare couplings, which depend on material or device specifics and are not included above.While such detailed analysis exceeds the scope of the present work, it is important to highlight that, as devices are scaled up, the localization effect induced by potential disorder may play a more pronounced role.This effect is particularly notable in one-dimensional interacting channels [4,151,152].We will therefore explore this phenomenon in the following section.

Anderson localization of the network
Apart from an isolated impurity, one can also consider potential disorder in the mth domain wall of the jth array as a random potential V ( j) dis,m , and investigate the possibility of Anderson localization in the network.We further assume a random potential of the form, with the disorder average ⟨⟨• • • ⟩⟩ and the disorder strength D ( j) b,m = h2 v 2 F /(2π λ mfp ) related to the mean free path λ mfp of the sample.For estimation purposes, we employ the experimentally extracted mean free path from typical two-dimensional devices [156], since a corresponding estimate for the one-dimensional domain wall mode is not available.The actual disorder strength may thus be further decreased due to reduced intervalley scattering, similar to the helical channels in quantum spin Hall insulators [161][162][163].As a side remark, while we examine a particular type of potential disorder locally coupled to the domain wall modes, other forms of disorder, such as local twist angle variation and distortion [114,164,165], might also influence the transport properties of the system.
The random potential considered here leads to elastic backscattering of electrons within the domain wall, which can be described as Lδ ′ σ,m (x) + H.c.. (25) Upon bosonization and performing disorder average [4], it leads to a contribution to the action, as given in equation (B.3) in appendix B, which can be quantified with the effective backscattering strength, Similar to single wires [4,151,152], one can derive the following RG flow equation for the effective backscattering strength, The flow equation indicates the presence of a gapped phase when ∆ ϕca,m < 6 − K ss − K sa = 4, indicating a possibility for the Anderson localization in the network for a system even with (weakly) attractive interactions.
Given that the renormalization of other parameters, such as the interaction strength and the velocity, contributes at higher orders [4], we neglect these subleading contributions.This allows us to compute the localization length and localization temperature directly from the RG flow equation in equation (27) and get and respectively.
In figures 13(g)-(h), we illustrate the influence of screening length and interlayer bias on the physical quantities.Remarkably, the results are a culmination of the analyses we have discussed so far.Specifically, they are influenced by multiple factors, including the velocity v dw of the domain wall modes presented in figure 4, as well as the interaction strength U ee,n summarized in figures 8 and 9.For the localization length, we observe a more pronounced dependence on the interlayer bias compared to the screening length.As the interlayer bias increases, it not only strengthens the interaction but also suppresses v dw .These two effects collectively contribute to the localization length, as evidenced in equation (28a).The localization temperature, on the other hand, exhibits a quantitatively different behavior.Within the explored parameter range, we note only a mild variation in the localization temperature.This is attributed to a balanced competition between the effects of interaction strength and v dw , as indicated by equation (28b), resulting in a region with weaker dependence on the external parameters.
Our findings indicate that, under typical conditions, the domain wall network is likely to transition into an Anderson insulator for domain walls longer than the order of 0.1 µm at temperatures of several kelvins or below.It is important to note that these results are contingent upon the actual disorder strength, which may vary from sample to sample.In our estimation, we have assumed a mean free path λ mfp = 1 µm to represent a sample with moderate mobility [156].Considering cleaner samples or reduced intervalley backscattering strength, we anticipate a longer localization length and a correspondingly lower localization temperature for the domain wall network.Moreover, while the disorder strength influences these quantities, the RG relevance is entirely determined by equation ( 27) and thus the electron-electron interaction encoded by the parameters ∆ ( j) ϕcs,m and ∆ ( j) ϕca,m .Consequently, as one attempts to scale up the device, the system can ultimately reach the localization regime at sufficiently low temperatures, provided that the disorder strength remains finite.

Correlation functions in the domain wall network
As one of the key features of TLL systems, various correlation functions exhibit power-law decay [1][2][3][4].These functions are defined by exponents highly dependent on the interaction strength and indicate the instability of the TLL towards various phases.In this section, we examine the instabilities of the correlated domain wall network towards the DW and SC phases.We analyze two scenarios, one without phonons in section 6.1 and the other with phonons in section 6.2.This investigation not only addresses the tendency to these instabilities but also illustrates the utility of the developed model in such analyses.

Correlation functions of purely electronic systems
As in previous works on TLL-like model [4], we start with DW correlation functions.To this end, we introduce the following two operators, for the charge and spin (number) density, respectively.In the above, we introduce the identity matrix σ 0 = 1, the Pauli matrix σ µ of the component µ ∈ {x, y, z}, the coordinate x = (x, ut) with the velocity u ≈ v dw .
The backscattering terms of the above can be defined as the CDW and SDW operators.Owing to the existence of two branches for each moving direction and spin, we obtain both intrabranch and interbranch terms, akin to those in carbon nanotubes [147,148,166].We distinguish these two types of contributions and apply bosonization to the corresponding operators using equation (7).The explicit expressions of the operators are presented in equation (B.4) in appendix B. With these expressions, we evaluate the intrabranch and interbranch CDW and SDW correlation functions, where the exponents corresponding to the CDW and µ component of the SDW are given in table 1, with the subscripts oe ö and ò denoting intrabranch and interbranch terms, respectively.
Similar to the DW correlation functions, we can also examine the pairing correlation functions.To this end, we introduce the following pairing operators for singlet superconductivity (SSC) and triplet superconductivity (TSC), which describe pairings with zero total momentum.Since interbranch pairing would involve Cooper pairs carrying nonzero momentum, we do not include it in the present analysis.With the explicit expressions given in equation (B.5) in appendix B, we obtain the pairing correlation functions, where the exponents corresponding to SSC and the µ component of the TSC are given in table 1.
Similar to the single-wire case [4], each of these correlation functions is characterized by a powerlaw behavior.Owing to the correlations between the domain walls of the network, the quantities for a given domain wall are also influenced by others, thereby reflecting the correlated characteristics of the entire network.From the above correlation functions, we get the instability conditions listed in table 1.
To proceed, we determine the dominant instability tendencies toward various phases by identifying the slowest decaying correlation function(s) from equations ( 30) and (32) in the parameter regime where the corresponding condition(s) in table 1 Table 1.Various correlation types in the domain wall network, the exponents of the corresponding correlation functions, their explicit forms, and the instability conditions.The correlation types sharing identical exponents are grouped in the same row.The correlation functions are given in equations ( 30), ( 32), (35), and (38).The interaction-dependent parameters ∆ ϕcP,m , and Γ ( j) θcP,m are given in equations (15b), (15c), (36), and (39) θcs,m and ∆ ( j) ϕcs,m < 1 for the parameter range under investigation, we expect the dominant instability of the correlated domain wall network is CDW and SDW along the domain walls, compared to SSC and TSC.Our model does not incorporate interactions or perturbations that break the spin rotational symmetry; hence, with K ss = K sa = 1, both the charge and spin sectors exhibit the same power-law behavior.We additionally examined the interwire correlation functions for both CDW and SDW (see equation (B.7)) and SC (see equation (B.9)).We found that these functions decay even more rapidly.Consequently, the explicit formulas are presented in appendix B, and we do not discuss them further in the main text.
To conclude, the correlation functions can also be controlled through their dependence on the interaction-dependent exponents.While in the pure electronic system the instability towards DW phases always prevails that towards SC phases, below we demonstrate the appearance of the phonon-induced SC in the correlated domain wall network.

Correlation functions in the presence of phonons
In this section we consider the presence of phonons, which create a displacement field and thus coupled to the low-energy electronic modes (that is, the domain wall modes in our system).The Hamiltonian now contains three terms, H ee + H ph + H ep , with the electron part in equation ( 9), the phonon part H ph , and the coupling H ep between the two.
While the electron-phonon coupling in twisted structures is currently an important topic of research [77,79,[167][168][169][170][171][172], a detailed microscopic model for this coupling in domain wall networks is not fully established.We therefore consider an effective model describing longitudinal acoustic phonons coupling to gapless modes that can move in onedimensional channels [173,174], with the effective mass density ρ a distributed along the domain wall, the phonon velocity c ph , the conjugate field Π j ph,m of the displacement field d j ph,m generated by the phonons.The displacement field couples to the charge density in the domain walls, which in the bosonic language can be expressed as with the effective coupling g ep between the domain wall mode and the phonon.Here, since we assume the same coupling strengths for both branches of a given domain wall, the effective coupling only acts on the charge symmetric sector.In a more general case, however, phonons can also couple to the ϕ j ca,m field and thus affect the exponent in the charge antisymmetric sector.
In our analysis, we include the nonperturbative effects of electron-phonon coupling, thereby extending beyond the perturbative regime for g ep , which we treat as a free parameter.It is essential to recognize that electron-phonon coupling may be significantly enhanced in the moiré structures, as highlighted by [79,167,170,172].Clearly, the incorporation of the phonon contributions leads to a breakdown in the duality between ϕ j cs,m and θ j cs,m typically existed in TLL [152].This indicates a change in the properties of the domain wall modes, potentially leading to different electronic behavior of the network.
With the phonon-induced contribution, we derive the effective action in equations (C.1)-(C.3),as well as the correlators modified by phonons in equation (C.4).To investigate the correlation functions for various DW phases, we focus exclusively on intrabranch correlations here, motivated by the fact that interbranch correlations are, at best, as RG relevant as intrabranch correlations.With the details presented in appendix C, we compute the following correlation functions, with ⟨• • • ⟩ ee+ph with respect to the action in equation (C.1).In the above, the exponents are given in table 1, with the phonon-modified parameters defined as Here, the phonon-induced modification stems from the factor η j ϕ,m , which contains the following parameters obtained in the course of deriving the effective action in appendix C, In the absence of electron-phonon coupling (that is, v j ep,m → 0), the velocities u j ±,m correspond to the velocity u j m of the charge symmetric mode in the mth domain wall of the jth array and that of the phonon.A nonzero electron-phonon coupling hybridizes the two, resulting in the renormalization of both velocities.This, combined with the characteristic features of TLL, can influence physical quantities via the exponents of various correlation functions discussed below.Similar to one-dimensional systems [173][174][175][176][177], in the limit of very strong electron-phonon coupling, there is a WB singularity where the velocity u −,m approaches zero, beyond which our model becomes unstable.
In the presence of phonons, we are particularly interested in the potential formation of SC phases within the domain wall network.This can be explored through the following SSC and TSC correlation functions, with the exponents given in table 1 and the modified parameters, The above results show that the pairing correlation in the network can be substantially influenced by the phonons.By comparing equations (36a) and (39a) with equations (15b) and (15c), it is evident that the phonon-induced effect is captured by the factor η j ϕ,m for the DW correlations and η j θ,m for the SC correlation.We therefore present their values in figure 14(a) and analyze how the effective electronphonon coupling influences these factors, which are unity in the absence of phonons.A key observation is their opposite trend (that is, η j ϕ,m ⩾ 1 and η j θ,m ⩽ 1), implying that the phonons effectively induce attractive interaction between electrons in the domain wall network.The stronger the electron-phonon coupling is, the more pronounced the deviation of η  of the modes reaches this divergence.Concerning the instability of the network, the η's significantly influence the parameters Γ The correlation functions of the network obtained here share qualitative similarities with those in strictly one-dimensional channels; see, for instance, [4].This similarity is evident in the power-law scaling, indicative of TLL characteristics, and in the general trend with the enhancement of DW and the suppression of SC due to repulsive electron-electron interactions.These findings align with the qualitative similarities between the one-dimensional systems and sTLL or csTLL [5][6][7][8].Nonetheless, the inclusion of interwire correlations modifies the network instability on a quantitative level.For a quantitative analysis, we compute the exponents of the correlation functions and derive phase diagrams, the latter determined by the dominant instability exponents among those listed in table 1.It is important to note that the scope of analysis can be extended to include various phases induced by backscattering operators, such as the moiré correlated states discussed in [21].Additionally, a unique aspect of the network here is the tunability of electron-electron interaction strength, adjustable through sample preparation and interlayer bias.
With the goal of providing a concrete estimation and inspired by the experimental setup described in [113], we consider a device with a short screening length, corresponding to a relatively strong screening effect.For our analysis below, the screened interaction strength can be subsequently modified by an interlayer bias.In addition, we examine three regimes characterized by the relative magnitude of phonon velocity to the Fermi velocity in a graphene monolayer.We use the latter quantity because the velocity of domain wall modes, being dependent on external parameters, does not serve well as an overall magnitude.
In figures 14(b)-(d), we present three types of phase diagrams, resulting from different values of the phonon velocity.As in figure 14(a), the WB singularity is indicated, which marks the strong electronphonon coupling limit.Beyond this limit, the model in equation ( 9) becomes unstable.For velocities comparable to the graphene Fermi velocity, figure 14(b) shows a phase diagram featuring a DW phase and a correlated domain wall network phase.The latter phase corresponds to the parameter regime where none of the instability conditions listed in table 1 is fulfilled, and the network remains gapless, with its properties discussed in section 5.In [173,174], the corresponding phase is identified as a metallic phase.It is also noteworthy that the phase diagram shown in figure 14(b) closely resembles the phase diagrams of systems that do not include phonons.
Interestingly, when the phonon velocity is sufficiently high, a new region emerges, as shown in figure 14(c).Specifically, with a sufficiently strong electron-phonon coupling strength, the dominant instability is SC at low V d values.At slightly higher V d values, a gapless network phase is observed, transitioning to a DW phase as electron-electron interactions increase at even higher V d values.Finally, an extremely high (and unrealistic) phonon velocity effectively suppresses the gapless network regime, as displayed in figure 14(d).Consequently, the region for SC becomes enlarged, as phonons mediate an instantaneous interaction between electrons.
In addition to the notable characteristics of electrically tunable network, which allows for in-situ driving of the system through phase transitions, we have several remarks.First, as has been identified in strictly one-dimensional interacting electronic systems [173][174][175][176][177], the WB singularity represents a phonon-driven instability in low-dimensional interacting electronic systems.Here, the location of this singularity in the parameter space, as defined in equation (37), is influenced by electron interaction and can thus be also electrically tuned, a feature absent in the previous literature [173][174][175][176][177]. Second, the phase diagrams are based on the exponent of a single domain wall in the interior of the device.As a result of the absence of translational invariance here, the exponents generally differ across the network, and the tendency towards SC or DW phases can vary.Consequently, the domain wall modes, including the two branches in each wall, progressively become unstable towards these phases in a nonuniform manner.We therefore expect a crossover, rather than a sharp transition in the network.Third, since the spin sectors are noninteracting (K ss = K sa = 1), the phase diagrams do not distinguish between spin-singlet and triplet phases.Introducing spin-dependent interactions, such as through proximity-induced spin-orbit coupling, could deviate from this specific case and further enrich these phase diagrams.
Before concluding this section, we note that the coupling strength between phonons and domain wall modes is a free parameter in our analysis, with its value not specified here.However, we observe that for a sufficiently strong coupling, the electronic properties can be adjusted to reach SC.As shown in figure 14(c), the range of SC broadens with decreasing electronic interaction strength.Additionally, as stated in section 3, our fitting procedure for the charge density is less accurate at very low V d values (V d ≲ 1.2hk θ v F ) due to weaker confinement of domain wall modes in that regime.Nonetheless, based on the trends we have observed, it is reasonable to expect an extended SC phase at even lower V d values.

Discussion
In the present analysis, we show how the physical properties of the domain wall network in TBG can be tuned electrically.Motivated by the network observation in minimally TBG [38], in most parts of this work we consider a relatively small twist angle, specifically θ = 0.5 • , compared to the magic angle.This value corresponds to a larger effective hybridization parameter3 α AB ≈ 1.4.However, our findings are not limited to this particular value.As demonstrated in figures 4 and 6, we have explored a range of α AB ∈ [1,2], corresponding to 0.35 • ≲ θ ≲ 0.7 • .Therefore, our analysis remains applicable across various parameter sets that result in the appearance of gapless modes within the domain walls.
We note that an even smaller twist can lead to a larger separation between the parallel domain walls, thereby reducing the interwire interaction strength.Conversely, a larger twist angle may bring the domain walls closer together, resulting in illdefined one-dimensional channels.To ensure welldefined domain wall channels for the TLL description, the criterion ρ2D(rAB) ρ2D(r dw ) ≪ 1 should be fulfilled.From the fitting function in equation ( 4), this criterion simplifies to e −2κ ⊥ ≪ 1, a condition that is fulfilled across the entire range of investigation.To deduce the twist angle range satisfying the criterion, we examine a simplified expression from [30].In our notation, it reads e −wABλM/(2πhvF) ≪ 1, indicating a twist angle on the order of O(1 • ) or smaller.It is noteworthy to highlight that the analysis in [30] did not include relaxation effects.Given that stronger confinement of the domain wall modes is expected in the presence of relaxation, larger twist angles are possible for the feasible range of TLL description.Furthermore, while we focus on TBG in the present work, the procedure outlined above can be extended to other systems [61][62][63][64] in which similar onedimensional channels have been identified.
From a broader perspective, we establish a mesoscopic network model from a single-particle model that describes interlayer hybridization at a more microscopic scale.This approach not only simplifies the subsequent analysis of correlated moiré systems with enlarged unit cells [31,83], but also provides a broader implication of the network description in various settings.Namely, we present a systematic method for determining the interaction strength in networks formed by domain wall modes.Moreover, extending beyond the previous works of the moiré network [21,67,68,72], here we explicitly incorporate the two branches for each moving direction and spin in a given domain wall.This consideration effectively establishes a mesoscopic network of coupled carbon nanotube TLLs.
Our approach further deviates from existing literature on sTLL or csTLL [5][6][7][8], which often rely on certain assumptions to build the bosonized model [21].Specifically, these studies usually adopt a predefined form of the interaction.In contrast, we avoid presuming the forms of density-density and currentcurrent interactions.Instead, we start directly with equation ( 9) and incorporate inputs from the outcomes derived from the single-particle continuum model in equation (1a).Additionally, we do not presume translational invariance perpendicular to the domain walls.Such invariance is not very realistic for a mesoscopic network with a limited number of parallel domain walls (fewer than 20 here), as opposed to the sTLL or csTLL description for CuO chains in typical cuprate samples.
In terms of observable features, our predictions include scaling behaviors in spectroscopic and transport measurements at microscopic scales.These behaviors are notably distinct from the topological edge modes of the quantum anomalous Hall states, where the scaling is determined by universal fractions [21].In contrast to the topological nature of the edge modes, the scaling exponents here are directly influenced by the interaction strength, making them adjustable using the experimental parameters outlined above.This tunability consequently affects physical properties of mesoscopic devices, such as localization length and temperatures.
Furthermore, we examine the stability of this correlated domain wall network.We demonstrate that, while the screened Coulomb interaction in purely electronic systems favors the formation of CDW and SDW, the presence of longitudinal acoustic phonons can still induce pairing instability, driving the system towards SC instability and even the WB singularity.Our focus here has been on the properties of the network model as described by the forward scattering terms of the screened electron-electron interactions.Therefore, we have not covered those driven by backscattering operators, which can be analyzed through the (perturbative) RG procedure, as shown in [21,67,68,72,178].Similarly, the 2k F phonons [179][180][181] not considered here can also induce instability by coupling to the 2k F component of the charge density.The analysis presented here can provide a more realistic view of the phase diagrams in those previous studies, where interaction strengths were introduced as free parameters.
Finally, our analysis indicates that twisted structures can manifest correlated phenomena by hosting spatially confined modes, without necessarily being tuned to magic angles.These modes are characterized by enhanced correlations and electrically tunable properties, forming the basis for an expanding theoretical literature employing bosonization [21,[66][67][68][70][71][72].In addition to TBG [30,32,35,36,38,40,45,73], the explored one-dimensional modes also appear in a wider range of nanoscale systems, including twisted bilayer WTe 2 [63,64,74] and twisted trilayer graphene [61,62].Our work suggests that the exploration of domain wall network in twisted structures can open avenues for understanding and manipulating correlated phenomena in nanoscale systems.In the presence of potential disorder, we obtain a perturbation term in equation (25), which contains the intrabranch and interbranch contributions.In general, the two terms have different scalings.Since they contain mutually conjugate fields, they cannot be ordered simultaneously.We therefore focus on the intrabranch contribution, motivated by its RG relevance, in the main text.To proceed, we apply the replica method, perform the disorder average [4,151,152], and arrive at the following contribution to the action, where we have introduced a replica index r, r ′ and the dimensionless coupling parameter D ( j) b,m (see equation ( 26)).The above formula is used to derive the RG flow equation in equation (27) in the main text.
Here we give the expressions for the intrabranch ( oe ö ) and interbranch (ò) terms of the CDW and SDW operators, ) ) These interwire correlation functions decay faster than the intrawire correlations in equations ( 30) and (32); we therefore do not discuss them in the main text.

Appendix C. Details about the calculation in the presence of phonons
In this section, we present the details about the calculation in the presence of phonons.Starting from equations ( 9), (33), and (34), we can perform straightforward algebra similar to [173,174,182]  which includes electronic degrees of freedom, phonon fields, and their couplings.In the bosonic language, all these terms remain at most bilinear in the fields, still allowing for direct diagonalization.By integrating out the phonon fields, one can derive the effective action that incorporates phonon-mediated interaction in the following form,  36) and (39), respectively.In addition, the interwire CDW and SDW correlations in the presence of phonons take the form of equation (B.7) upon the replacement of ∆

Figure 1 .
Figure 1.(a) TBG real space lattice.The moiré lattice vectors λ Mj have a magnitude of moiré wavelength λM = a0/(2 sin(θ/2)) with the graphene lattice constant a0.The domain wall network is marked in green, with local AB-or BA-stacking configurations indicated.(b) TBG reciprocal space.The original first Brillouin zones of the two monolayers are shown in red and orange.The inset displays the moiré Brillouin Zone (mBZ) of TBG, highlighting new symmetry points (denoted with bars).The vectors q j link the Dirac cones between layers.The vectors q M1 = q 0 − q 1 and q M2 = q 1 − q 2 are the reciprocal lattice vectors of TBG with magnitude qM = √ 3k θ and k θ = 4π/(3λM).

Figure 2 .
Figure 2. Spatial profiles of charge density for representative states within the mBZ, marked by a hexagon, computed from equation (1).The high symmetry points are also indicated.The blue and yellow curves represent the lowest energy bands intersecting at the chemical potential µ = 0 along three inequivalent K-K ′ lines.The inset provides a detailed view of the computed spatial charge density profiles (with domain walls marked by dashed lines) for six states (indicated by pink dots) at the charge neutrality point.Unless otherwise specified, we adopt V d /hvFk θ = 1.5 and αAB = 1.4 (corresponding to a twist angle of θ = 0.5 • ) throughout the article.

Figure 3 .
Figure 3. (a) Density profile of gapless domain wall modes perpendicular to the domain wall in the y direction for various interlayer biases (V d ).(b) Interlayer bias (V d ) dependence of the density (blue) at the AB-stacking region center (rAB) and the ratio (purple) of the density at the domain wall segment center (r dw ) to that at rAB.The three V d values corresponding to Panel (a) are marked with dashed lines.The adopted values of the other parameters are given in the caption of figure 2.

Figure 4 .
Figure 4. (a) Interlayer bias (V d ) and effective hybridization parameter (αAB) dependence of the bandwidth (∆a) of the domain wall modes.(b) A similar plot for the velocity (v dw ).The values of these quantities are extracted from the energy spectrum from equation (1).

Figure 5 .
Figure 5. Density profile (in arbitrary units) of the gapless modes in the direction (a) perpendicular and (b) parallel to a domain wall along y direction.The blue dots are computed from equation (1) and the red curves are fitted with equation (4).The adopted values of the other parameters are given in the caption of figure 2.

Figure 6 .
Figure 6.Interlayer bias (V d ) and effective hybridization parameter (αAB) dependence of fitted parameter values from equation (4).For the parameter range of interest, the simple ansatz presented in equation (4) provides a fairly good fit for the computed density.The white regions are excluded due to large fitting errors.

Figure 7 .
Figure 7.A schematic of the device is depicted, including the TBG stacked with a dielectric layer (blue) of thickness d and a metallic gate (light blue).For a given direction of parallel domain wall modes (navy blue) within the TBG, image charges induced in both the dielectric (light purple) and metallic (orange) layers are indicated.

Figure 8 .
Figure 8.(a) Screening length d and interlayer bias V d dependence of the intrawire interaction strength (Uee,0).(b) Interlayer bias V d and (c) screening length d dependence of the intrawire interaction strength (Uee,0) for three dielectric materials.The dots indicate values calculated from equation(6), with the curves serving as a visual guide.Unless otherwise specified, we adopt L ∥ = 0.5 µm, Lz = 6.68 Å (twice the interlayer separation of graphite), and ϵr = 6.93, corresponding to a hexagonal Boron Nitride (hBN) dielectric layer[139] for the rest of this article.The adopted values of the other parameters are given in the caption of figure 2.

Figure 9 .
Figure 9. Interwire interaction strength calculated from equation (6), for various (a) dielectric layer thickness, with V d /hvFk θ = 1.5, and (b) interlayer bias, with d = 20 nm.The adopted values of the other parameters are given in the captions of figures 2 and 8.

Figure 10 .
Figure 10.Illustration of various microscopic processes in the domain wall network discussed in section 5.The inset shows the linearized energy spectrum near the Fermi level in the mth domain wall of the jth array.Blue and yellow indicate the two branches shown in figure 2, with solid and dashed lines representing opposite spin states.The operator ψ ( j) ℓδσ,m represents a fermion field.For ℓ = R (ℓ = L), it corresponds to the right-(left-) moving mode.The outer (inner) branch is denoted by δ = 1 (δ = 2), and σ indicates the spin state.
caption, the domain wall index m ∈ [1, N ⊥ ], the movingdirection index ℓ ∈ {R ≡ +, L ≡ −} along a given domain wall, the spin index σ ∈ {↑ ≡ +, ↓ ≡ −}, the local coordinate x, the Fermi wave vector k ( j) Fδ,m , and δ ∈ {1 ≡ +, 2 ≡ −} labeling the outer/inner branch of the domain wall spectrum.In the bosonic description, we introduce the Klein factor U j ℓδσ,m , the shortdistance cutoff set by a = hv dw /∆ a , and the boson fields satisfying ( j) θcP,m and θ j cP,m .It is straightforward to check that in the new basis the fields ϕ j cP,m and θ j cP,m follow the same commutation relation as equation ( θcP,m show increased deviations from unity near the sample edges.Additionally, we explore the impact of the chemical potential µ on ∆ ( j) ϕcP,m and ∆ ( j) θcP,m , finding a relatively weak dependence (see figure A1 for m = 9).The results indicate the robustness of the exponents (including those discussed below) with respect to the doping level, provided that the chemical

Figure 11 .
Figure 11.(a), (b) Dependence of ∆ ( j) ϕcs,9 and ∆ ( j) θcs,9 on screening length (d) and interlayer bias (V d ), calculated from equation (15) for a domain wall labeled by m = 9 in a device comprising 18 domain walls per array (N ⊥ = 18).(c) Illustration of the scanning tunneling spectroscopy measurement.(d) A similar plot for the exponent (β ( j) 9 ) of the local DOS calculated from equation (15), with ∆ ( j) ϕca,m = ∆ ( j) θca,m = Kss = Ksa = 1.The adopted values of the other parameters are given in the captions of figures 2 and 8.

Figure 12 .
Figure 12.Spatial dependence of the exponents (a) ∆ ( j) ϕcs,m and (b) ∆ ( j) θcs,m in equation (15) under various interlayer biases (V d ).Here we show the results for domain walls labeled by m in a device with d = 60 nm and N ⊥ = 18.The adopted values of the other parameters are given in the captions of figures 2 and 8.

Figure 13 .
Figure 13.Illustrations of the (a) intrawire backscattering induced by an isolated impurity (gray dot), (b) tunneling process between two nonparallel domain walls at their intersection and (c) tunneling process between two parallel domain walls (which can take place at any point along the domain wall direction).(d)-(h) Screening length (d) and interlayer bias (V d ) dependence of the exponents α ( j) imp,m in equation (18b), α ( j,j ′ ) ×,mm ′ in equation (20b), α ( j) ⊥,mm ′ in equation (23b), the localization length (ξ ( j) loc,m ) in equation (28a), and the localization temperature (T ( j) loc,m ) in equation (28b).The numerical results for (d)-(h) are obtained for domain wall(s) representing the interior of a device with N ⊥ = 18 and λ mfp = 1 µm [156].The adopted values of the other parameters are given in the captions of figures 2, 8, and 11.
is rather simple, as it exhibits a constant shift from ∆ ( j) ϕcs,m and thus follows an identical trend.Given that both exponents α( j,j ′ ) ×,mm ′ and α ( j) ⊥,mm ′ are linear functions of the sum ∆ ( j) ϕcs,m + ∆ ( j)θcs,m , they behave similarly.As previously mentioned, in an isotropic and (relatively) homogeneous network like the one studied here, both ∆ ( j) ϕcs,m and ∆ ( j) θcs,m show only a weak dependence on m.
j ϕ ,m and η j θ,m is from unity.Beyond a critical value, a divergence occurs in the parameters, with η j ϕ ,m → ∞ and η j θ,m → 0. Consequently, we identify the WB singularity at the coupling strength where the first
crucial in determining various correlation functions, as demonstrated in the inset.With a strong electron-phonon coupling, a lower η j θ,m value reduces the exponent in the pairing correlation function.Therefore, one can in general expect a slower decay of pairing correlations in the presence of phonons, indicating a stronger tendency towards SC.