Tunable quantum criticalities in an isospin extended Hubbard model simulator

Device fabrication

We fabricated the twisted double bilayer graphene by a modified polymer-based cut-pick technique. We mechanically exfoliated bilayer graphene, graphite and hexagonal boron nitride (h-BN) onto a highly doped Si wafer covered by a 300-nm-thick SiO₂ layer, with the thickness of these flakes identified with optical contrast and a Bruker MultiMode 8 atomic force microscope. The bilayer graphene was cut into two halves by using the atomic force microscope tip in contact mode. We then used stamps consisting of poly(bisphenol A carbonate) film and polydimethylsiloxane to pick up top graphite and h-BN at around 80 °C. After that, we picked up one half of the bilayer graphene and then the other half with a rotation of 180°+0.8° at room temperature. This aimed rotation angle is slightly larger than the twisted angle of our measured device, which is attributed to the inevitable angle relaxation in the fabrication process. Finally, we picked up bottom h-BN at 80 °C and the five-layer heterostructures  were transferred onto bottom graphite at around 135 °C. The heterostructures were then shaped into Hall bar geometry by using standard electron-beam lithography and dry etching in an inductively coupled plasma system. The edge contact electrodes (10 nm Cr/10 nm Pd/30 nm Au) to cTDBG and top electrodes (5 nm Ti/40 nm Au) to the top and bottom graphite gates were deposited by standard electron-beam evaporation. The temperature was always kept below 150 °C during stacking and fabrication processes.

Electrical measurements

The devices were measured in an Oxford cryostat with a base temperature of about 1.5 K. The resistance signals were collected using a low-frequency (17.7 Hz) SR830 Lock-in Amplifier (Stanford Research). Gate voltages were applied with Keithley 2400 Source Meters. The differential resistance signals were collected by applying an a.c. excitation current of 1 nA added on top of a variable d.c. bias current up to 1 μA. We apply d.c. currents and record amplified voltages by adopting a data acquisition system (National Instruments 6251).

Continuum model of twisted double bilayer graphene

We consider AB-BA-stacked twisted double bilayer graphene: an AB-stacked bilayer graphene is placed on top of another BA-stacked bilayer graphene and they are twisted with respect to each other by an angle θ, forming a moiré superlattice. The moiré lattice constant L_s = a/2sin(θ/2), in which a = 0.246 nm, is the atomic lattice constant of monolayer graphene. We further consider atomic corrugation in the system, that is, the interlayer distance d(r) between the two twisted layers at the interface (between the AB and BA bilayers) varies as a function of in-plane position r in the moiré supercell, which can be approximately modelled as \(d\left({\bf{r}}\right)={d}_{0}\left({\bf{r}}\right)+2{d}_{1}\mathop{\sum }\limits_{i=1}^{3}\cos {{\bf{g}}}_{i}\cdot {\bf{r}}\), in which g₁, g₂ and g₃ = g₁ + g₂ are the three reciprocal lattice vectors of the moiré supercell. We take d₀ = 0.3433 nm and d₁ = 0.00278 nm to reproduce the interlayer distances of the AA-stacked and AB-stacked bilayer graphene at the ABBA point and the ABAC point in the moiré supercell, respectively.

The low-energy electronic structure of cTDBG can be well described by the Bistritzer–MacDonald continuum model⁴¹, in which the low-energy states from the two atomic valleys K and K′ are assumed to be completely decoupled from each other. To be specific, the continuum model for the K valley of the cTDBG system is expressed as

in which \({{\rm{H}}}_{{\rm{AB}}}^{K}\) and \({{\rm{H}}}_{{\rm{BA}}}^{K}\) are the Hamiltonians of the untwisted bilayers:

and

in which K₁ (K₂) denotes the Dirac point of the AB-stacked (BA-stacked) bilayer. The Pauli matrices σ = (−σ_x, σ_y) are defined in the space of the A, B sublattices of graphene. h₊ denotes the interlayer hopping matrix between the untwisted bilayer,

in which \({t}_{\perp }=0.48\,{\rm{eV}}\), t₂ = 0.21 eV and t₃ = 0.05 eV are the nearest-neighbour, second-neighbour and third-neighbour interlayer hopping amplitudes, respectively, and \({h}_{-}={h}_{+}^{\dagger }\).

We consider only the nearest-neighbour interlayer coupling between the two twisted bilayers, that is, the coupling between the top layer of the AB bilayer and the bottom layer of the BA bilayer, which is expressed as

in which ΔK = K₂ − K₁ = (0, 4π/3L_s) is the shift between the two Dirac points of the two sets of twisted bilayers⁴². U is the interlayer hopping term between the two Dirac states of the twisted layers⁴³,

in which \({{\bf{r}}}_{{\rm{AB}}}{\boldsymbol{=}}\left(\sqrt{3}{L}_{{\rm{s}}}/3,\,0\right)\), \(g\left({\bf{r}}\right)={\sum }_{j=1}^{3}{{\rm{e}}}^{{{\rm{i}}{\bf{q}}}_{j}\cdot {\bf{r}}}\), with q₁ = (0, 4π/3L_s), \({{\bf{q}}}_{2}=\left(-2\pi /\sqrt{3}{L}_{{\rm{s}}},-2\pi /3{L}_{{\rm{s}}}\right)\) and \({{\bf{q}}}_{3}=\left(2\pi /\sqrt{3}{L}_{{\rm{s}}},-2\pi /3{L}_{{\rm{s}}}\right)\) (ref. ⁴⁴). In cTDBG, the corrugation effects result in \({u}_{0} < {u}_{0}^{{\prime} }\), in which u₀ ≈ 0.078 eV and \({u}_{0}^{{\prime} }\approx 0.097\,{\rm{eV}}\) are the intra-sublattice and inter-sublattice interlayer coupling terms, respectively⁴³. We note that all the parameters of the continuum model, as given in equations (1)–(6), are derived from a realistic Slater–Koster tight-binding model introduced in ref. ⁴⁵. We have also compared the results calculated from the continuum model with those calculated using the atomic Slater–Koster tight-binding model and find excellent agreement with each other.

The displacement field D is introduced by applying a homogeneous vertical electrostatic potential drop across the four layers

in which U_d = eD·d_t/ϵ is the electrostatic potential energy difference between the top and bottom layers, D is the displacement field, d_t ≈ 1.005 nm is the thickness of the entire system and ϵ ≈ 4 is the dielectric constant of the BN substrate. We note that several different values for the dielectric constant of h-BN has been widely used, which can lead to markedly different conclusions in the theoretical calculations. In our experiment, the adoption of ϵ ≈ 3.8 gives us a consistent result, so we choose the value of 4 for our theoretical calculations.

We note that the Coulomb potentials from the fully occupied bands below the target flat band may strongly enhance the bandwidth of the flat band and renormalize the flat-band wavefunction. This question has recently attracted a lot of attention. For example, Vafek and Kang recently performed comprehensive renormalization group studies on magic-angle twisted bilayer graphene⁴⁶ and found that both the kinetic energy and the Coulomb interaction energy increase as the Fermi level approaches the charge neutrality point. In the end, when the flat bands are partially occupied, the ratio between the interaction energy projected to the (renormalized) flat-band subspace and the renormalized bandwidth is roughly unchanged, and the system is still in the strong-coupling regime. A similar argument could also be applied to the twisted double bilayer graphene system.

To calculate the numerical prefactor C of the transport ‘scattering rate’, we have used the m_e at the valence band edge (in which the effective mass can be well defined) and the carrier density corresponding to the 8n₀ filling. Because the 7+2/3 filling in our device is in very close proximity to the 8n₀ filling, it is therefore reasonable that the C value calculated for the 8n₀ filling applies to that at the 7+2/3 filling.

Criterion for the Wigner crystal state

Here we introduce a dimensionless quantity r_s = V/W, the ratio between the Coulomb interaction energy (V) and kinetic energy (W) of the electrons, as the key criteria for the formation of Wigner crystal state⁴⁷. The average distance between the electrons in the moiré bands is denoted as r_e, which satisfies \(\pi {r}_{{\rm{e}}}^{2}=1/{n}_{{\rm{e}}}\), in which n_e = v/Ω is the carrier density, with Ω the area of the moiré primitive cell and ν the filling factor. Then the characteristic Coulomb interaction energy is V = e²/4πϵ₀ϵr_e, in which e is the elementary electronic charge, ϵ ≈ 4 is the dimensionless background dielectric constant and ϵ₀ is the vacuum permittivity. The characteristic kinetic energy of the system can be effectively described by the free electron kinetic energy \(W={\hbar }^{2}{k}_{{\rm{F}}}^{2}/2{m}_{{\rm{e}}}^{* }\), in which \({m}_{{\rm{e}}}^{* }\) is the effective mass of the carrier at the Fermi level and k_F is the Fermi wave vector. Notably, \({m}_{{\rm{e}}}^{* }\) is calculated to be about 0.2m_e, based on the formula \(\frac{1}{{m}^{* }}=\frac{1}{{\hbar }^{2}}{\nabla }_{{\bf{k}}}^{2}\,E({\bf{k}})\). The carrier density is related to the Fermi wave vector:

in which the factor of 4 is from the fourfold valley–spin degeneracy. Then we obtain:

in which a_B = 4πϵ₀ħ²/m₀e² is the Bohr radius and m₀ is the bare electron mass. According to previous studies^48,49, electrons crystallize for r_s ≥ 31 in the one-valley two-dimensional electron gas (2DEG) system and for r_s ≥ 30 for the two-valley 2DEG system⁴⁹.

We note that the small bandwidth of the second conduction flat band (around 6–8 meV) gives rise to r_s of about 35 (see Extended Data Fig. 6a), which exceeds the value for Wigner crystallization in 2DEGs⁴⁹. The zero Chern number indicates a trivial band topology different from the topological charge density wave and fractional Chern insulator residing in the first moiré band of other moiré graphene systems^50,51,52. Notably, the Chern number is non-zero in twisted double bilayer graphene of the stacking of AB-AB at the same fillings and D field, indicating that the stacking order plays a crucial role in forming the Wigner crystal. In addition, the Wigner crystal state only emerges at a large displacement field, at which there is only one Fermi surface with large r_s centred at the Γ_s point for each valley (see Extended Data Fig. 6b). By contrast, as |D| is reduced to zero, other Fermi surfaces with much smaller r_s emerge (see Extended Data Fig. 6c,d), making the system still behave as a Fermi liquid.

Scaling analysis

The procedure of scaling analysis is detailed described next. (1) We first measure the resistance R at different temperatures T and displacement fields D, obtaining several R–T curves at different D, with typical results shown in Fig. 3a. (2) We choose a trial critical field D* near the phase boundaries at the base temperature in the measurements, which has R–T curve denoted by R*(T), and normalize the R–T curves at different D in the Wigner crystal regime by R*(T). We then try to make the normalized R/R^*(T) curves at different D collapse onto a single branch after scaling the temperatures by D-dependent T₀. If this fails, then find a nearby point to repeat this process and find a point that has the optimal collapse of R–T curves. The critical D_c for the Wigner crystal regime is determined to be −0.325 V nm⁻¹ in our device in zero magnetic field. A similar method has been used for other scaling analyses in this work.

Determination of T
_n

Now we explain in detail the algorithm for how we determine the phase boundary T_n between the strange metal and normal metal phases. In the normal metal region (T < T_n), the R–T relation follows R(T) = AT² + R₀, whereas in the strange metal region (T > T_n), R is proportional to T. It follows that the temperature boundary T_n between strange metal and normal metal can be determined by fitting the functional relationship of the resistance R and temperature T. We first attempt to fit R to the form R(T) = AT² + R₀ from the base temperature to a trial maximum temperature (T_trial) and extract the r_p² representing the correlation coefficient of such parabola fitting. Meanwhile, we fit R to the form R(T) = AT + R₀ for T > T_trial to extract another r_l² representing the correlation coefficient of this linear fitting. Then we try several different T_trial for the fitting to search for the temperature simultaneously having high r_p² and r_l² values, which is then determined as T_n. Extended Data Fig. 7 shows a typical example of fitting at \(n=7\frac{2}{3}{n}_{0}\) and D = −0.21 V nm⁻¹. We find that the resistance is fit well (r² = 0.98) by a T² form up to a temperature of 14.7 K and is also fit well (r² = 0.996) by a T-linear form between 14.7 and 30 K. Here we can also determine the error of the fitting. If there exist a few trial temperatures that all give good fitting results, then we can define this range of T as the error of fitting. If there is only one trial temperature that can give good fitting, then we define the error by the step length of the temperature in our measurements. Such error analysis gave the error bars in Figs. 3c and 4b.