Machine Learning Microscopic Form of Nematic Order in twisted double-bilayer graphene

TL;DR

This paper demonstrates how convolutional neural networks can analyze STM data to learn microscopic nematic order parameters in twisted double-bilayer graphene, distinguishing them from heterostrain effects.

Contribution

It introduces a CNN-based approach to extract microscopic nematic order from STM data in moiré superlattices, accounting for correlations and heterostrain.

Findings

CNN can learn microscopic nematic order parameters.

Including energy correlations improves model accuracy.

Neural networks distinguish nematic order from heterostrain.

Abstract

Modern scanning probe techniques, like scanning tunneling microscopy (STM), provide access to a large amount of data encoding the underlying physics of quantum matter. In this work, we analyze how convolutional neural networks (CNN) can be employed to learn effective theoretical models from STM data on correlated moir\'e superlattices. These engineered systems are particularly well suited for this task as their enhanced lattice constant provides unprecedented access to intra-unit-cell physics and their tunability allows for high-dimensional data sets within a single sample. Using electronic nematic order in twisted double-bilayer graphene (TDBG) as an example, we show that including correlations between the local density of states (LDOS) at different energies allows CNNs not only to learn the microscopic nematic order parameter, but also to distinguish it from heterostrain. These results demonstrate that neural networks constitute a powerful methodology for investigating the microscopic details of correlated phenomena in moir\'e systems and beyond.

Equations

$$\begin{array}{c}{\mathcal{H}}_{\mathbf{\Phi }}={\int }_{\mathbf{r}}{\int }_{\Delta \mathbf{r}}\mathbf{\Phi }\cdot {\mathbf{\varphi }}_{\sigma ,\ell ,s,\eta ;{\sigma }^{\prime },{\ell }^{\prime },s^{\prime },{\eta }^{\prime }}\left(\mathbf{r},\Delta \mathbf{r}\right)\\ \times c_{\sigma ,\ell ,s,\eta }^{†}\left(\mathbf{r}+\Delta \mathbf{r}\right)c_{{\sigma }^{\prime },{\ell }^{\prime },s^{\prime },{\eta }^{\prime }}^{}(\mathbf{r})+\text{H.c.},\end{array}$$

$$\begin{array}{c}{\mathcal{H}}_{\mathbf{\Phi }}^{\text{MN}}=\frac{1}{2}{\Phi }_{\text{MN}}{\int }_{\mathbf{r}}\sum _{m_1,m_2\in \mathbb{Z}}{\hat{\mathbf{\Phi }}}_{{\phi }_{\text{MN}}}\cdot {\mathbf{\varphi }}_{m_1,m_2}(\mathbf{r})\\ \times c_{\alpha }^{†}(\mathbf{r}+{\mathbf{R}}_{m_1,m_2})c_{\alpha }^{}(\mathbf{r})+\text{H.c.},\end{array}$$

$${\mathcal{H}}_{\mathbf{\Phi }}^{\mathrm{GN}}={\Phi }_{\text{GN}}{\int }_{\mathbf{r}}{\hat{\mathbf{\Phi }}}_{{\phi }_{\text{GN}}}\cdot {\mathbf{\varphi }}_{\ell ,s;{\ell }^{\prime },s^{\prime }}(\eta ;\mathbf{r})c_{\ell ,s}^{†}(\mathbf{r})c_{{\ell }^{\prime },s^{\prime }}^{}(\mathbf{r}).$$

$${\mathbf{\varphi }}_{\ell ,s;{\ell }^{\prime },s^{\prime }}(\eta ;\mathbf{r})={\delta }_{\ell ,{\ell }^{\prime }}{\psi }_{\ell }\left(\begin{array}{c}(e^{i{\alpha }_{\ell }\eta {\rho }_z}{\rho }_x{)}_{s{s}^{\prime }}\\ \eta (e^{i{\alpha }_{\ell }\eta {\rho }_z}{\rho }_y{)}_{s{s}^{\prime }}\end{array}\right),$$

$${\mathcal{E}}_2=-{\mathcal{E}}_1=\frac{1}{2}R({\theta }_{ϵ}{)}^{-1}\left(\begin{array}{cc}-ϵ& 0\\ 0& \nu ϵ\end{array}\right)R({\theta }_{ϵ}).$$

$$H_{AB-AB}=\left(\begin{array}{cccc}{H}_0\left({\mathbf{k}}_1\right)& S^{†}\left({\mathbf{k}}_1\right)& & \\ S\left({\mathbf{k}}_1\right)& H_0^{\prime }\left({\mathbf{k}}_1\right)& U^{†}& \\ & U& H_0\left({\mathbf{k}}_2\right)& S^{†}\left({\mathbf{k}}_2\right)\\ & & S\left({\mathbf{k}}_2\right)& H_0\left({\mathbf{k}}_2^{\prime }\right)\end{array}\right)+V$$

$$H_0\left(\mathbf{k}\right)=\left(\begin{array}{cc}0& -\hslash \nu k_{-}\\ -\hslash \nu k_{+}& d\end{array}\right)\text{and}{H}_0^{\prime }\left(\mathbf{k}\right)=\left(\begin{array}{cc}d& -\hslash \nu k_{-}\\ -\hslash \nu k_{+}& 0\end{array}\right).$$

$$S\left(\mathbf{k}\right)=\left(\begin{array}{cc}\hslash {\nu }_4{k}_{+}& {\gamma }_1\\ \hslash {\nu }_3{k}_{-}& \hslash {\nu }_4{k}_{+}\end{array}\right),$$

$$V=\text{diag}\left({\Delta }_1{\rho }_0,{\Delta }_2{\rho }_0,{\Delta }_3{\rho }_0,{\Delta }_4{\rho }_0\right)$$

$$U=\left(\begin{array}{cc}u& u^{\prime }\\ u^{\prime }& u\end{array}\right)+\left(\begin{array}{cc}u& u^{\prime }{\omega }^{-\eta }\\ u^{\prime }{\omega }^{\eta }& u\end{array}\right)e^{i\eta {\mathbf{G}}_1^M\cdot \mathbf{r}}+\left(\begin{array}{cc}u& u^{\prime }{\omega }^{\eta }\\ u^{\prime }{\omega }^{-\eta }& u\end{array}\right)e^{i\eta \left({\mathbf{G}}_1^M+{\mathbf{G}}_2^M\right)\cdot \mathbf{r}}$$

$${\Psi }_n\left(\mathbf{k}\right)={\left({\psi }_{n,\mathbf{k}}\left({\mathbf{G}}_{n_1,n_2}^1\right),\cdots ,{\psi }_{n,\mathbf{k}}\left({\mathbf{G}}_{n_1,n_2}^c\right)\right)}^T$$

$${\psi }_{n,\mathbf{k}}\left(\mathbf{G}\right)={\left(U_{n,\mathbf{k}}^{A_1}\left(\mathbf{G}\right),U_{n,\mathbf{k}}^{A_2}\left(\mathbf{G}\right),\cdots ,U_{n,\mathbf{k}}^{B_4}\left(\mathbf{G}\right)\right)}^T$$

$$\mathcal{D}\left(\mathbf{r},\omega \right)=\sum _{n,\mathbf{k}}\sum _{\mathbf{G},{\mathbf{G}}^{\prime }}{e}^{-i\left(\mathbf{G}-{\mathbf{G}}^{\prime }\right)\mathbf{r}}\delta \left(\omega -{\omega }_{n,\mathbf{k}}\right)\times \left({\left[U_{n,\mathbf{k}}^{A_4}\left({\mathbf{G}}^{\prime }\right)\right]}^{\ast }{U}_{n,\mathbf{k}}^{A_4}\left(\mathbf{G}\right)+{\left[U_{n,\mathbf{k}}^{B_4}\left({\mathbf{G}}^{\prime }\right)\right]}^{\ast }{U}_{n,\mathbf{k}}^{B_4}\left(\mathbf{G}\right)\right).$$

$$W\left(t,h\right)=\frac{1}{\sqrt{s}}\sum _{i=0}^{N_{\omega }-1}{\mathcal{D}}_{{\mathbf{r}}_0}({\omega }_i)\psi \left(\frac{i-t}{h}\right),\text{with}\psi \left(t\right)=e^{-t^2/2}\cos \left(5t\right),$$

$$W_{\pm }(\mathbf{r})\propto \exp \left(\mp c_{1}y\left(y^2-3x^2\right)-c_2{\left(x^2+y^2\right)}^2\right),\mathbf{r}=(x,y{)}^T,$$

$${\mathcal{H}}_{\text{tb}}=\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}{h}_{\mathbf{k}}^0{c}_{\mathbf{k}}^{},h_{\mathbf{k}}^0=\left(\begin{array}{cc}{h}_{\mathbf{k}}^{W_{+}}({\alpha }_2,\Delta )& h_{\mathbf{k}}^C\\ h_{\mathbf{k}}^{C†}& h_{\mathbf{k}}^{W_{-}}(-{\alpha }_2,-\Delta )\end{array}\right),$$

$$h_{\mathbf{k}}^W({\alpha }_2,\Delta )=\left(\begin{array}{cc}\Delta +f\left({\alpha }_2,t_2\right)& g\left(t_1,t_3\right)\\ g^{†}\left(t_1,t_3\right)& -\Delta +f\left(-{\alpha }_2,t_2\right)\end{array}\right),$$

$$h_{\mathbf{k}}^C=\left(\begin{array}{cc}{e}^{i\Theta \omega }{\omega }_1& {\omega }_2\left(1+2e^{i\frac{\sqrt{3}}{2}{k}_x}\cos (k_y/2)\right)\\ {\omega }_2\left(1+2e^{-i\frac{\sqrt{3}}{2}{k}_x}\cos (k_y/2)\right)& e^{i\Theta \omega }{\omega }_1\end{array}\right),$$

$$f\left({\alpha }_2,t_2\right)=t_2\left(\cos (\mathbf{k}\cdot {\mathbf{a}}_2-{\alpha }_2)+\cos (\mathbf{k}\cdot ({\mathbf{a}}_1-{\mathbf{a}}_2)-{\alpha }_2)+\cos (\mathbf{k}\cdot {\mathbf{a}}_1+{\alpha }_2)\right)$$

$$g\left(t_1,t_3\right)=t_1\left(1+e^{i\mathbf{k}\cdot {\mathbf{a}}_1}+e^{i\mathbf{k}\cdot {\mathbf{a}}_2}\right)+t_3{e}^{i{\alpha }_3}\left(e^{-i\mathbf{k}\cdot \left({\mathbf{a}}_2-{\mathbf{a}}_1\right)}+e^{-i\mathbf{k}\cdot \left(-{\mathbf{a}}_2-{\mathbf{a}}_1\right)}+e^{-i\mathbf{k}\cdot \left(-{\mathbf{a}}_2+{\mathbf{a}}_1\right)}\right),$$

$${\mathbf{g}}_{\mathbf{k}}={\alpha }_0{\rho }_0{\sigma }_0\left(\begin{array}{c}{X}_{\mathbf{k}}\\ Y_{\mathbf{k}}\end{array}\right)+{\alpha }_1{\rho }_0{\sigma }_x\left(\begin{array}{c}{X}_{\mathbf{k}}\\ Y_{\mathbf{k}}\end{array}\right)+{\alpha }_2{\rho }_0{\sigma }_y\left(\begin{array}{c}-Y_{\mathbf{k}}\\ X_{\mathbf{k}}\end{array}\right)+{\alpha }_3{\rho }_z{\sigma }_z\left(\begin{array}{c}-Y_{\mathbf{k}}\\ X_{\mathbf{k}}\end{array}\right),$$

$$(X_{\mathbf{k}},Y_{\mathbf{k}})=\frac{8}{3}\left(\cos k_y-\cos \frac{\sqrt{3}{k}_x}{2}\cos \frac{k_y}{2},\sqrt{3}\sin \frac{\sqrt{3}{k}_x}{2}\sin \frac{k_y}{2}\right).$$

$$\frac{\text{dI}}{\text{dV}}\left(\mathbf{r},\omega \right)\propto \text{Im}\left[\sum _{j,k,\alpha ,\beta }{W}_{{\mathbf{R}}_{j\alpha }}\left(\mathbf{r}\right)G_{\alpha \beta }^R\left({\mathbf{R}}_j-{\mathbf{R}}_k,\omega \right)W_{{\mathbf{R}}_{k\beta }}^{\ast }\left(\mathbf{r}\right)\right]=\mathcal{D}\left(\mathbf{r},\omega \right)$$

$$G_{\alpha \beta }^R=\left(\mathbf{R}-{\mathbf{R}}^{\prime },\omega \right)=\frac{1}{V}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\left(\mathbf{R}-{\mathbf{R}}^{\prime }\right)}\underset{\eta \to 0^{+}}{\lim }{\left(\frac{1}{\omega -h_{\mathbf{k}}+i\eta }\right)}_{\alpha ,\beta }.$$

Full text

I Introduction

Driven by the impressive improvements in machine learning (ML) in the last couple of years, exploring its potential for quantum many-body physics has recently become subject of intense research Carleo et al. (2019); Dawid et al. (2022). For instance, ML provides powerful tools to solve inverse problems that occur frequently in physics Lee et al. (2023); Dubois et al. (2022); Berthusen et al. (2021); Chertkov and Clark (2018): given a model, it is often straightforward with conventional many-body techniques to compute observables that can be measured experimentally, whereas the often needed inverse problem of extracting the model and underlying microscopic physics from observations is much more challenging and typically even formally ill-defined. A second example of a large class of applications of ML in physics is ML-assisted analysis of experiments, in particular of those yielding image-like data like scanning tunneling microscopy (STM) Choudhary et al. (2021); Joucken et al. (2022); Wang et al. (2020); Zhang et al. (2019), photoemission Liu et al. (2022a), and others Khan et al. (2023); Chen et al. (2022); Ede (2021); Iwasawa et al. (2022); Burzawa et al. (2019); Basak et al. (2022).

In the context of applying ML algorithms to data from imaging techniques like STM, van der Waals moiré superlattices Balents et al. (2020); Andrei and MacDonald (2020) are particularly promising for three reasons: (i) they display a huge variety of correlated quantum-many-body phenomena, such as interaction-induced insulating phases Cao et al. (2018a), magnetism Sharpe et al. (2019), superconductivity Cao et al. (2018b), electronic nematic order Kerelsky et al. (2019); Jiang et al. (2019); Choi et al. (2019); Cao et al. (2021), which can also coexist microscopically Cao et al. (2021); Lin et al. (2022). Despite intense research on these phenomena over several decades, e.g., in the pnictides or cuprates, their origin and relations are still subject of ongoing debates. However, compared to these microscopic crystalline quantum materials, moiré superlattices are (ii) highly tunable; for instance, the density of carriers can be varied within a single sample just by applying a gate voltage (as opposed to chemical doping) and even the interactions can be tuned Liu et al. (2021). This allows producing large data sets of measurements on a single sample, containing a lot of information on the microscopic physics. This aspect, which is crucial for data-driven approaches, is further enhanced by (iii) the large moiré unit cells of these systems compared to that of microscopic crystals, increasing the relative spatial resolution of scanning-probe techniques significantly. This enables experiments to probe the structure of the wave functions within the unit cell, and thus provides unprecedented access to the microscopic physics compared to conventional quantum materials. For instance, in the extreme limit of only one degree of freedom (Wannier state or pixel) per unit cell, the broken rotational symmetry of the electron liquid—the defining property of electronic nematic order Fradkin et al. (2010); Fernandes et al. (2014)—is not visible as a consequence of translational symmetry and thus requires a careful analysis of the behavior around impurities Goetz et al. (2020).

In this work, we explore these advantages of moiré superlattices for extracting or ‘learning’ effective field-theoretical descriptions of their correlated many-body physics from STM data. This can be viewed as an inverse problem and is also conceptually related to the goal of ‘Hamiltonian learning’ in quantum simulation Granade et al. (2012); Wiebe et al. (2014); Wang et al. (2017); Valenti et al. (2019); Kokail et al. (2021); Yu et al. (2022), albeit in rather different regimes and based on different measurement schemes. As a concrete example, we use electronic nematic order in twisted double-bilayer graphene (TDBG) Cao et al. (2020); Liu et al. (2020); Shen et al. (2020); Rubio-Verdú et al. (2022); He et al. (2021); Kuiri et al. (2022); Su et al. (2022). This moiré system consists of two AB-stacked bilayers of graphene that are twisted against each other; as one can see in Fig. 1(a), it exhibits the point group $D_{3}$, generated by three-fold rotation $C_{3}$ along the out-of-plane $z$ axis and two-fold rotation $C_{2x}$ along the in-plane $x$ axis. Evidence of electronic nematic order has been observed in previous STM experiments Rubio-Verdú et al. (2022); Samajdar et al. (2021) which clearly exhibit stripe-like features breaking the $C_{3}$ symmetry spontaneously for certain electron concentrations. While simple limiting cases have been compared with the data in Ref. Samajdar et al., 2021, there is no systematic analysis of the microscopic form of nematicity in the system. To fill this gap, we consider the more general case in which all leading terms on the graphene and moiré scale describing nematic order in a continuum-model description of TDBG are included. In addition, as it is common in graphene moiré systems Huder et al. (2018); Kerelsky et al. (2019); Jiang et al. (2019); Choi et al. (2019); Rubio-Verdú et al. (2022), we also allow for finite strain. The Hamiltonian defining the changes in TDBG resulting from nematic order and strain depends on a set of parameters $\beta$, which we reconstruct from STM data using convolutional neural networks (CNN) in a supervised learning procedure. As such, our study differs significantly from recent works, which focused on detecting the presence or absence of nematic order Goetz et al. (2020) or performed a phenomenological data analysis of STM measurements Taranto et al. (2022) with ML, rather than extracting the underlying microscopic physics as we do here.

II Results

II.1 Nematic order in TDBG

The non-interacting band structure of TDBG features two moiré minibands per spin and valley close to charge neutrality, where a variety of correlation-driven phenomena can emerge Cao et al. (2020); Liu et al. (2020); Shen et al. (2020); Rubio-Verdú et al. (2022); He et al. (2021); Kuiri et al. (2022); Su et al. (2022). In Fig. 1(b), these minibands are denoted as valence (VFB) and conduction flat bands (CFB). The band structure shown is obtained from continuum model calculations close to half filling of the CFB (band filling $\nu=0.475$), where electronic nematic order was observed to be the strongest Rubio-Verdú et al. (2022), see Appendix A for more details. STM experiments probe the band structure and wave functions of a system by providing direct access to the spatial and energy dependence of the LDOS. Most commonly, the LDOS is studied either for a fixed position $\bm{r}_{0}$ over a range of different energies, $\mathcal{D}_{\bm{r}_{0}}(\omega)$, or for a fixed energy $\omega_{0}$ covering a spatial region of the system, $\mathcal{D}_{\omega_{0}}(\mathbf{r})$. The behavior of $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ and $\mathcal{D}_{\bm{r}_{0}}(\omega)$ following from the continuum model for TDBG for three different energies and high-symmetry positions in the moiré unit cell is shown in Fig. 1(c). The $C_{3}$ rotational and translational symmetry of the moiré lattice can be clearly seen in $\mathcal{D}_{\omega_{0}}(\mathbf{r})$. Meanwhile, $C_{2x}$ is broken, albeit weakly, as a consequence of the electric field required to control the electron filling to be close to the middle of the CFB in an open-faced STM sample geometry Rubio-Verdú et al. (2022).

In graphene moiré systems, there are two fundamentally distinct sources of $C_{3}$ symmetry breaking—strain and electronic nematic order. Postponing the discussion of the former to below, electronic nematic order Fradkin et al. (2010); Fernandes et al. (2014) refers to the spontaneous rotational symmetry breaking as a result of electronic correlations. While recent works also indicate the possibility of nematic charge-density wave states in TDBG Wilhelm et al. (2022); He et al. (2021), where moiré translational symmetry is simultaneously broken, we here focus on translationally symmetric nematic order since the STM data of Ref. Rubio-Verdú et al., 2022 preserves moiré translations. The underlying nematic order parameter we study is a time-reversal- and moiré-translation-invariant vector $\mathbf{\Phi}=\Phi\hat{\mathbf{\Phi}}_{\varphi}$, $\hat{\mathbf{\Phi}}_{\varphi}=\left(\cos 2\varphi,\sin 2\varphi\right)$, transforming under the irreducible representation $E$ of $D_{3}$ (or of $C_{3}$, taking into account the weak $C_{2x}$ breaking); $\Phi$ and $\varphi$ stand for the intensity and orientation of the nematic director, respectively. The microscopic form of nematicity can be modeled by a coupling of $\mathbf{\Phi}$ to a fermionic bilinear and reads in its most general form in a continuum-model description as Samajdar et al. (2021)

[TABLE]

where $c^{\dagger}$ and $c$ are the electronic creation and annihilation operators. This general form encompasses couplings between the two sublattices $s=A,B$ of the microscopic graphene sheets, the four graphene layers $\ell=1,\dots,4$, the valley $\eta=\pm$ and spin $\sigma=\>\uparrow,\downarrow$ degrees of freedom in the tensorial form factor $\bm{\phi}_{\sigma,\ell,s,\eta;\>\sigma^{\prime},\ell^{\prime},s^{\prime},\eta^{\prime}}(\bm{r},\Delta\mathbf{r})$; its two components are required to transform in the same way as $\bm{\Phi}$ under all symmetries of the system. In the following, we will take $\bm{\phi}$ to be trivial in the spin and diagonal in the valley indices, $\bm{\phi}_{\sigma,\ell,s,\eta;\>\sigma^{\prime},\ell^{\prime},s^{\prime},\eta^{\prime}}=\delta_{\sigma,\sigma^{\prime}}\delta_{\eta,\eta^{\prime}}\bm{\phi}_{\ell,s;\>\ell^{\prime},s^{\prime}}(\eta)$. This is motivated by the weak spin-orbit coupling in graphene Kane and Mele (2005); Min et al. (2006) and the lack of indications of interaction-induced spin-orbit coupling, which is also strongly constrained Kiselev et al. (2017). Furthermore, the intervalley-coherent nematicity is known to lead to stronger effects on the remote bands Samajdar et al. (2021) that were not observed experimentally Rubio-Verdú et al. (2022).

Since we are working with a continuum theory, the space of possible couplings $\bm{\phi}$ in Eq. (1) is technically infinite dimensional. As such, a complete reconstruction of $\bm{\phi}$ from experimental data is impossible given the finite resolution and energy range of the available data. On top of this, it is not required either as we are primarily interested in understanding the low-energy behavior of the system. In the spirit of gradient expansions commonly used in continuum low-energy field theories, we will therefore only keep the leading terms in $\bm{\Phi}$. There is, however, a subtlety associated with the presence of an additional moiré length scale. We will therefore have to consider two basic classes of nematic orders, referred to as graphene (GN) and moiré (MN) nematicity Rubio-Verdú et al. (2022); Samajdar et al. (2021).

In the case of MN, nematic order is associated with the moiré scale, i.e., we choose $\Delta\bm{r}=\mathbf{R}_{m_{1},m_{2}}=m_{1}\mathbf{L}_{1}^{M}+m_{2}\mathbf{L}_{2}^{M}$ in Eq. (1), $m_{j}\in\mathbb{Z}$, with moiré lattice vectors $\mathbf{L}_{j}^{M}$, to represent the non-trivial transformation behavior of $\bm{\phi}$ under $C_{3}$. We can thus take it to be diagonal in the remaining internal indices, yielding

[TABLE]

with multi-index $\alpha=\left(\sigma,\ell,s,\eta\right)$. We further focus on the lowest moiré-lattice harmonic by setting $\phi_{m_{1},m_{2}}(\mathbf{r})=\phi_{m_{1},m_{2}}$ and only keeping the terms with the shortest possible $\mathbf{R}_{m_{1},m_{2}}$. Intuitively, MN order can be thought of as a distortion of the effective inter-moiré-unit-cell hopping matrix elements, as illustrated schematically in the lower right panel of Fig. 1(d).

Conversely, GN acts as a local order parameter, $\Delta\bm{r}=0$ in Eq. (1), without any explicit reference to the moiré scale,

[TABLE]

Here, the correct transformation properties of $\bm{\phi}$ result from its structure in the internal indices. Focusing on the local intra-layer contributions and the leading (constant) basis function, the most general form reads as

[TABLE]

where Pauli matrices in sublattice space are represented by $\rho_{j}$; $\alpha_{l}$ and $\psi_{l}$ are real-valued parameters. As shown schematically in the upper left panel of Fig. 1(d), one can think of GN as the nematic distortion of the bonds of the individual graphene layers in a way that preserves the graphene translational symmetry.

We emphasize that GN and MN should not be viewed as distinct phases; they break the same symmetries and as such in general mix. We thus take $\mathcal{H}_{\bm{\Phi}}^{\text{MN}}+\mathcal{H}_{\bm{\Phi}}^{\mathrm{GN}}$ to describe nematicity in TDBG in the following, which depends on the set of parameters $\beta=\{\alpha_{\ell},\psi_{\ell},\Phi_{\text{MN}},\Phi_{\text{GN}},\varphi_{\text{MN}},\varphi_{\text{GN}}\}$. The computation of the LDOS for a specific set of parameters can be done straightforwardly from the continuum model. The resulting spatial dependence of the LDOS, $\mathcal{D}_{\omega_{0}}(\mathbf{r})$, is also shown in Fig. 1(d) for two different values of $\beta$. As opposed to the plots without nematic order, $C_{3}$ is now broken, leading to stripes in the VFB, while translational symmetry is still preserved. The inverse problem—inferring the value of the parameters $\beta$ from a given LDOS pattern—is a much more challenging task. Our goal in the following sections will be to use ML, in particular, CNNs to learn the set $\beta$ directly from LDOS images.

II.2 ML architecture

Using CNNs to solve this inverse problem can be interpreted as a supervised learning task Dawid et al. (2022), i.e., a regression-like procedure using synthetic LDOS data labeled by their respective value of nematicity parameters $\beta$. More specifically, our CNNs take as inputs $65\times 65$ pixels of LDOS images and apply consecutive transformations (represented by a set of weights between each layer) in order to extract meaningful correlations that represent the set $\beta$. One example of the CNN image inputs is shown in Fig. 2(a). The complete dataset consists of 12000 images which are divided into training ($60\%$), validation $(20\%)$ and test $(20\%)$ subgroups. Each image is generated for a randomly sampled set of nematic parameters $\beta$ and the intensities in the LDOS are modified with the addition of Gaussian noise (see Appendix A). The motivation for noise is twofold: to avoid overfitting Goodfellow et al. (2015) and to test the stability against and performance of the procedure with noise, which is inevitably present in experimental data.

The ML architecture chosen for this task is represented in Fig. 2(a) and its implementation was done with the TensorFlow library Abadi et al. (2015). Each convolutional layer is followed by batch normalization and max pooling layers (Conv-Batch-MaxPool). The batch normalization layers normalize the input weights in each stage, and also reduce the number of epochs necessary for convergence Ioffe and Szegedy (2015). This process is repeated four times, with the convolutional layers having a kernel size of $3\times 3$ and strides set to $1$. The filters follow a sequence of $16-32-32-16$ with rectified linear unit (ReLU) activation functions Fukushima (1975). Padding is set to zero such that the reduction of dimensionality is performed only by the MaxPool layers. In turn, these have both strides and pool sizes set to $2\times 2$. After a Flatten stage, dense layers lead to a dropout before the final layer with filters equal to the number of parameters in $\beta$. The Flatten layer transforms the data to a one-dimensional shape, and the Dropout reduces overfitting by setting a percentage of 20% adjusted weights to zero Srivastava et al. (2014). Tests on variations of this architecture are briefly described in Appendix B.

The learning procedure is then defined by the minimization of the loss function with respect to the CNN’s weights in a backward propagation procedure Rumelhart et al. (1986). The loss function can be represented as the mean squared error (MSE), which is defined as the difference between the true and expected set of parameters $\beta$ in $\text{MSE}=\sum_{j}^{N}\left(\beta_{j}^{\text{true}}-\beta_{j}^{\text{predicted}}\right)^{2}/N$, with $N$ representing the number of samples in the test dataset. Finally, we consider the adaptive moment estimation (ADAM) for the minimization of the loss function, with a learning rate of 0.001 and batch size equal to 64 Kingma and Ba (2017). After the completion of the training stage, the algorithm is ready to be deployed to previously unseen data, returning as outputs the parameters $\beta^{\text{predicted}}$.

II.3 Orientation of nematic director

As a first investigation, we consider the task of predicting the orientation $\varphi$ of the nematic director from $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ images at a single energy in the VFB [$\omega_{0}=-15$ meV, see Fig. 1(b)]. For this, we consider a dataset with randomly generated MN and GN intensities $\Phi^{\text{MN}},\Phi^{\text{GN}}\in[0.001,0.1]$ eV, and $\varphi^{\text{MN}}=\varphi^{\text{GN}}=\varphi\in[0,\pi]$. Furthermore, $\psi_{l}=1$ and $\alpha_{l}=0$ for all layers. The relation between the shape of the LDOS at a single energy $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ and $\varphi$ is highly non-trivial for two reasons: even for a given form of nematicity, changing $\varphi$ generically not just merely rotates the LDOS pattern, due to the lattice, but leads to complex distortions of its structure. Additionally, by sampling $\mathcal{H}_{\bm{\Phi}}^{\text{MN}}+\mathcal{H}_{\bm{\Phi}}^{\mathrm{GN}}$, even if the same bond direction is favored over the $C_{3}$-related ones in the LDOS pattern of two samples, the underlying $\varphi$ can be rather different. As can be seen in the three sample LDOS plots in Fig. 2(b) with different values of $\varphi$, the correspondence between $\varphi$ and $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ is complex and not apparent to the human eye.

Using the angles $\varphi$ as labels to the data is the most straightforward choice, but leads to inaccurate predictions around [math] and $\pi$ due to the periodicity in the definition of the nematic order parameter, $\hat{\bm{\Phi}}_{\varphi}=(\cos 2\varphi,\sin 2\varphi)=\hat{\bm{\Phi}}_{\varphi+\pi}$. To circumvent this feature, we use the two-component label $\hat{\bm{\Phi}}_{\varphi}$ instead of $\varphi$ in the training process and then fold the network’s prediction back to $\varphi$ with the $arctan2$ function Fischer et al. (2015). The results, shown in Fig. 2(b), are consistent with the true labels, including at the boundaries of $\varphi$’s domain. This shows that even when the precise nature of nematicity (predominantly MN or GN or an admixture of the two) is not known, the director orientation $\varphi$ can be accurately predicted with our CNN setup from $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ at a single energy. We have checked that the few outliers in Fig. 2(b) are directly related to small nematic intensities, where $\varphi$ has virtually no impact on the LDOS and is, thus, impossible to predict.

II.4 Form of nematicity

After successfully learning the director orientation $\varphi$ in the presence of different nematicities, we proceed into investigating finer details of these couplings by learning the parameters $\beta=\{\Phi^{\text{MN}},\Phi^{\text{GN}},\alpha_{l}\}$ defined in Eqs. (2-4). To this end, we consider $\psi_{l}=1$ and $\alpha_{l}=\alpha$ for all layers. For concreteness, we set $\varphi^{\text{MN}}=\varphi^{\text{GN}}=\varphi=2\pi/3$, which is one of the possible discrete orientations ($\varphi^{\text{MN}}=\varphi^{\text{GN}}=2\pi/3,\pi/6$ and symmetry related) of the nematic director in presence of $C_{2x}$. The dataset now consists of randomly generated MN and GN intensities $\Phi^{\text{MN}},\Phi^{\text{GN}}\in[0.001,0.1]$ eV, and $\alpha\in[0,\pi]$. The intensity values are chosen such that the stripes in the VFB resemble the experimental results Rubio-Verdú et al. (2022). As with $\varphi$, instead of learning the angular variable $\alpha$ directly, the $arctan2$ mapping from Section II.3 is applied.

Using only the LDOS at a single energy (i.e. one $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ channel) in the ML architecture for this task does not produce accurate predictions. Additionally, both hyperparameter optimization and architecture modifications did not lead to any significant improvement, implying that nematic order impacts the electronic structure in complex ways that cascade across energy scales. In fact, this is also intuitively clear since, for example, the samples marked by a star and pentagon in Fig. 3(a) have fundamentally different nematic couplings and yet exhibit visually similar $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ images at the VFB energy.

In experiments, one can typically obtain single point spectra [$\mathcal{D}_{\bm{r}_{0}}(\omega)$] and real space LDOS images at fixed energies [$\mathcal{D}_{\omega_{0}}(\mathbf{r})$]. We can therefore include additional input channels corresponding to $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ and $\mathcal{D}_{\bm{r}_{0}}(\omega)$ for different energies $\omega_{0}$ and points $\bm{r}_{0}$, respectively. In the second case, the individual point spectra are transformed to scaleogram images for consistency with the input data for CNNs Berthusen et al. (2021); Mallat (1999), see upper left inset in Fig. 3(a) and Appendix A. The new architecture is then formed by four channels with $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ inputs at fixed energies $\omega_{0}=\left(-35,-15,1,23\right)$ meV within the flat and remote bands, such that they resemble visually the corresponding ones in the experimental data of Ref. Rubio-Verdú et al., 2022, and three channels for $\mathcal{D}_{\bm{r}_{0}}(\omega)$ scaleogram inputs at stacking positions $\bm{r}_{0}=\{\text{ABAB, BAAC, ABCA}\}$, cf. Fig. 1(c). Each channel is passed through parallel Conv-Batch-MaxPool layers as in Fig. 2(a), but instead of flattening each channel separately, they are concatenated to a Dense-Dropout stage before the last layer [Fig. 3(a)].

In Fig. 3(b-d), predictions on the test data set are represented for (b) $\alpha$, and (c) the moiré and (d) graphene nematic intensities; as can be seen, very good agreement is found between the reconstructed and true parameters. The outliers in $\alpha$ are related to small $\Phi_{\text{GN}}$ (brighter colors). From Eqs. (3) and (4), it is clear that for small $\Phi_{\text{GN}}$, minimal changes will be induced in the LDOS, irrespective of the true value of the phase governed by $\alpha$. This is a similar behavior to what was observed for outliers in the nematic director prediction. If $\alpha$ is maintained fixed, we observed (not shown) that predictions for $\Phi_{\text{GN}}$ and $\Phi_{\text{MN}}$ get more accurate. The results of Fig. 3 demonstrate that the microscopic form of nematicity can be extracted from the LDOS if significant energy dependence is included in the input data set.

II.5 Including strain

As already alluded to above, another possible source of $C_{3}$ breaking is strain Nguyen and Dollfus (2015); Yan et al. (2013); Huder et al. (2018); Bi et al. (2019), which is believed to be a ubiquitous property of graphene moiré superlattices at small twist angles. Breaking the same symmetries as nematic order, strain can obscure the experimental identification of nematic order and their precise interplay is still under debate Kerelsky et al. (2019); Jiang et al. (2019); Choi et al. (2019); Scheurer (2019). Experiments indicate Huder et al. (2018); Kerelsky et al. (2019); Jiang et al. (2019); Choi et al. (2019); Rubio-Verdú et al. (2022) that the most relevant form of strain in graphene superlattices such as twisted bilayer graphene or TDBG is uniaxial heterostrain. In this case, the matrices $\mathcal{E}_{j}$ describing the in-plane metric deformation of the coordinates in the $j$th rotated bernal bilayer of TDBG are of the form

[TABLE]

Here $\nu=0.16$ is the Poisson ratio for graphene and $R(\theta_{\epsilon})$ is the $2\times 2$ matrix describing rotations of 2D vectors by angle $\theta_{\epsilon}$. We see that uniaxial heterostrain is characterized by two variables, the strain intensity $\epsilon$ and the direction of strain, parameterized by the angle $\theta_{\epsilon}$.

In the following, we allow for the simultaneous presence of uniaxial heterostrain and nematic order, leading to two additional parameters, $\epsilon$ and $\theta_{\epsilon}$, in $\beta$. We will study whether our ML approach is still able to extract the microscopic form of nematicity and also learn the relative strength and direction of strain. Note that the form of nematicity is still given by Eqs. (2-4), with the only difference that we replace $\mathbf{L}_{j}^{M}$ in the definition of $\mathbf{R}_{m_{1},m_{2}}$ by the strained moiré lattice vectors. The data set for this task is built with nematic intensities $\Phi^{\text{MN}},\Phi^{\text{GN}}\in[0.001,0.1]$ eV, with the addition of strain parameters $\epsilon\in[0,0.8]\>\%$ and $\theta_{\epsilon}\in[0,\pi/3]$. Here, $\alpha_{l}=0$, $\psi_{l}=1$ and $\varphi=\varphi_{\text{MN}}=\varphi_{\text{GN}}=2\pi/3$. The domain for the strain intensities is chosen based on typical values observed in TBG Kerelsky et al. (2019), and for $\theta_{\epsilon}$ on the periodicity of the unstrained system as $\theta_{\epsilon}\rightarrow\theta_{\epsilon}+\pi/3$ Bi et al. (2019). The ML architecture employed in this section is the same as in the previous investigation [Fig. 3(a)].

In Fig. 4(a-d), predictions on the test data set are shown for $\epsilon$ (a), $\theta_{\epsilon}$ (b), and the nematic intensities (c-d). At first sight, the result for the strain angle in Fig. 4(b) looks as if the procedure ceased to work, since there are many data points where the true and predicted value of $\theta_{\epsilon}$ differ significantly. However, when indicating the true strain intensity label $\epsilon$ for each prediction, it becomes clear that the outliers are related to small values of $\epsilon$ (brighter colors). As such, this behavior is not a shortcoming of the learning procedure but actually a feature of strain: for small enough $\epsilon$ in Eq. (5), the angle $\theta_{\epsilon}$ has no meaning. We have checked that removing the samples with small strain $\epsilon$ from the training and test data set will lead to accurate predictions of $\theta_{\epsilon}$. The stability that we find for our learning procedure in the presence of virtually vanishing $\epsilon$ is, however, important when applying it to experimental data, where the strength of strain is unknown.

Most importantly, we see in Fig. 4(c-d) that the nematic couplings can still be accurately predicted when varying strain is present. The MAE is equally distributed in these cases, in contrast to the strain intensity prediction. This shows that not only nematic order can be identified when strain is present, but also its internal structure and the strength of strain that is present at the same time can be resolved when using different channels consisting of both $\mathcal{D}_{\bm{r}_{0}}(\omega)$ and $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ as inputs. This allows the networks to take into account correlations between different energies in the STM data, which in turn conveys the crucial microscopic physics, enabling the model to disambiguate between lattice and electronic effects.

II.6 Experimental data

After demonstrating the effectiveness of CNNs on learning microscopic parameters $\{\beta_{i}\}$ from a synthetic (theoretical) data set $D_{\text{th}}\left(\beta_{1},\cdots,\beta_{N_{\text{th}}}\right)$ with $N_{\text{th}}$ samples, we now proceed into applying the trained ML architecture for predictions of the a priori unknown sets of parameters $\{\beta^{\prime}_{i}\}$ in an experimental data set $D_{\text{exp}}\left(\beta_{1}^{\prime},\cdots,\beta_{N_{\text{exp}}}^{\prime}\right)$. For concreteness, we use the same synthetic training data set as in Appendix B, where only the nematic and strain intensities are predicted, i.e., $\beta=\{\Phi_{\text{MN}},\Phi_{\text{GN}},\epsilon\}$. The data set $D_{\text{exp}}$ is constituted of both scaleograms $\mathcal{D}_{\mathbf{r}_{0}}(\omega)$ and $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ maps for different fillings of the CFB ($n_{s}$). More details about the preprocessing of the experimental data $D_{\text{exp}}$ can be found in Appendix C.

In Fig. 5, predictions of the trained CNN for the set $\{\beta_{i}^{\prime}\}$ show non-zero values of nematicity (a) and strain (b) for all fillings of the CFB. For $n_{s}\geq 0.47$ (gray region), the experimental data shows the most pronounced signatures of broken rotational symmetry to the human eye, which was previously interpreted as electronic nematic order Samajdar et al. (2021); Rubio-Verdú et al. (2022). Here the CNN predicts MN to dominate over GN, although both are finite (as expected by symmetry). As can be seen in Fig. 5(c), the parameters predicted by the CNN nicely reproduce the key features in the experimental data, including the strong stripes in the VFB and the much weaker, albeit finite, signatures of nematicity in the other bands.

For smaller fillings, $n_{s}<0.47$, the experimental data still exhibit distortions that break $C_{3}$, see Appendix C, but no clear stripe-like features appear. The CNN tries to assign different anisotropy sources to these distorted regions, but the agreement between theoretical prediction and experiment is less accurate than for larger $n_{s}$. It is clearly possible that, indeed, a crossover from primarily MN to GN occurs when lowering $n_{s}$, as predicted by the neural network, see Fig. 5(a), in particular, since nematic order is also a plausible instability in non-twisted bilayer graphene Liu et al. (2021); Cvetkovic et al. (2012). However, we believe that additional experimental data and refined theoretical models are required to conclude whether this is really the case.

In contrast to this interplay between the nematic couplings, strain remains relatively constant for all $n_{s}$, and slightly decreases in Fig. 5(b) for $n_{s}\geq 0.47$ as it approaches the same order of magnitude of $\epsilon\in[0.003-0.1\%]$ that is expected for the experimental samples in $D_{\text{exp}}$ Rubio-Verdú et al. (2022). We note that at low fillings the value of strain that is predicted by the neural network is nevertheless significantly greater than the value extracted from experimental topography. This is likely a consequence of subtle differences between the continuum model calculations and the experimental spectroscopy, which the network attempts to accommodate by including finite strain.

III Discussion

We constructed and demonstrated a ML procedure that can extract the form of the nematic order parameter in TDBG from LDOS data. The key ingredient was the use of several channels that capture the correlations among different energies. Our work has several important implications. First, it shows that the presence and even the strength and internal structure of nematic order can be extracted when the sample exhibits significant heterostrain; this is a crucial aspect for moiré systems where the issue of distinguishing between nematicity and strain has been subject of debate. Second, our analysis also shows which type of STM data is needed and most useful to extract information about nematicity: as we have seen, the LDOS maps at a single energy, $\mathcal{D}_{\omega_{0}}(\bm{r})$, are not enough to deduce the form of the nematic order parameter and—contrary to what one might have expected—point spectra, i.e., $\mathcal{D}_{\bm{r}_{0}}(\omega)$, contain a lot of helpful complementary information for that task (see also the second model discussed in Appendix D). We emphasize that this form of ‘solid-state Hamiltonian learning’, i.e., of parameterizing the leading terms of a set of microscopic order parameters (like nematic order) or perturbations (such as strain) and extracting their form using multi-channel CNNs can be more broadly applied—to other systems, see Appendix D where we discuss a toy model for twisted-bilayer graphene, and other forms of instabilities. As such, this could open up completely new ways of revealing the form and role of nematic order and other phases for the physics of quantum materials.

Acknowledgements.

J.A.S. and M.S.S. acknowledge funding by the European Union (ERC-2021-STG, Project 101040651—SuperCorr). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Salary support is also provided by the National Science Foundation via grant DMR-2004691 (S.T.) and by the Office of Basic Energy Sciences, Materials Sciences and Engineering Division, U.S. Department of Energy under Contract No. DE-SC0012704 (A.N.P.). J.A.S. is grateful for discussions with J. P. Valeriano, Sayan Banerjee, Patrick Wilhelm, Igor Reis and Pedro H. P. Cintra. M.S.S. also thanks R. Samajdar, R. Fernandes, and J. Venderbos on a previous collaboration on nematic order in TDBG Samajdar et al. (2021).

Appendix A Continuum model and LDOS maps

Twisted double-bilayer graphene (TDBG) consists of two bilayer graphene stacks with primitive lattice vectors $\bm{a}_{1}=a\left(1,0\right)$ and $\bm{a}_{2}=a\left(1,\sqrt{3}\right)/2$, where $a\simeq 0.246$ nm is the lattice constant of graphene, and corresponding reciprocal lattice vectors $\bm{b}_{j}$, following from $\bm{a}_{i}\cdot\bm{b}_{j}=2\pi\delta_{ij}$. After applying a twist angle $\theta$ between the Bernal stacks, these vectors are modified by the two-dimensional rotation matrix $R(\theta)$ as $\left(\bm{a}^{l}_{i},\bm{b}^{l}_{i}\right)=R\left(\mp\theta/2\right)\left(\bm{a}_{i},\bm{b}_{i}\right)$ with $-,+$ for each stack $l=1,2$. The corresponding Dirac cones in each valley $\eta=\pm 1$ of the individual graphene layers are located at $\bm{K}^{l}_{\eta}=-\eta\left(2\bm{b}_{1}^{l}+\bm{b}_{2}^{l}\right)/3$. The emerging moiré pattern [with triangular Bravais lattice shown as green domains in Fig. 1(a)] is represented by the difference of the new lattice vectors from each bilayer stack in reciprocal space as $\bm{G}_{i}^{M}=\bm{b}_{i}^{1}-\bm{b}_{i}^{2}\>\left(i=1,2\right)$, with corresponding primitive lattice vectors $\bm{L}_{j}^{M}$ obtained from the relation $\bm{G}_{i}^{M}\cdot\bm{L}_{j}^{M}=2\pi\delta_{ij}$.

We consider a description of the low-energy physics for TDBG via the continuum Hamiltonian of Ref. Koshino, 2019. In the Bloch basis given by carbon’s pz orbitals $\left(A_{1},B_{1},\dots,A_{4},B_{4}\right)$, with sublattices $s=\{A_{\ell},B_{\ell}\}$ and layers $\ell=1,2,3,4$, the continuum Hamiltonian in valley $\eta$ for small twist angles ($\theta\ll 1$) in AB-AB double bilayer graphene can be written as

[TABLE]

with Bloch wavevectors generated by $\bm{k}_{j}=R(\mp\theta/2)\left(\bm{k}-\bm{K}_{\eta}^{j}\right)\left(j=1,2\right)$ and single-layer graphene Hamiltonians with $k_{\pm}=\eta k_{x}\pm ik_{y}$ as

[TABLE]

The explicit matrix structure in Eq. (7) refers to sublattice space (associated with Pauli matrices $\rho_{j}$ in the main text). In turn, $H_{0}$ and $H_{0}^{\prime}$ are coupled via

[TABLE]

with parameters $\{d,\hbar\nu/a,\gamma_{1},\nu_{3},\nu_{4}\}=\{0.050,2.776,0.4,0.32,0.044\}\>\text{eV}$ and $\nu_{i}=\left(\sqrt{3}/2\right)\gamma_{i}a/\hbar\left(i=3,4\right)$. For more details about the physical significance of each term, see Ref. Koshino (2019). Considering a self-consistently calculated screened electric field, the interlayer potential matrix reads as

[TABLE]

where $\rho_{0}$ is the unit matrix in sublattice space. For a filling fraction of $\nu=0.475$, on-site potentials representing the electrostatic energy between adjacent layers are given by $\Delta=(\Delta_{1},\Delta_{2},\Delta_{3},\Delta_{4})=(4.079,1.021,-1.537,-3.563)$ meV Samajdar et al. (2021), which we use in our numerical calculations.

Finally, the moiré interlayer coupling defined between the twisted layers $\ell=2-3$ is given by

[TABLE]

with $\omega=\exp 2\pi i/3$, $u=0.0797$ eV and $u^{\prime}=0.0975$ eV. These parameters were chosen in accordance with Ref. Samajdar et al. (2021); Koshino (2019).

From the continuum Hamiltonian in Eq. (6), a numerical diagonalization in momentum space is performed by selecting a finite number of $\bm{q}$ wavevectors in a cutoff circle $|\bm{q}-\bm{q}_{0}|<q_{c}$, with radius $q_{c}=4|\bm{G}_{i}^{M}|$ around the midpoint $\bm{q_{0}}=\left(\bm{K}_{\eta}^{1}+\bm{K}_{\eta}^{2}\right)/2$ between Dirac cones $\bm{K}_{\eta}^{j}$. Here, the Bloch vector $\mathbf{k}$ in the moiré Brillouin zone is hybridized with the graphene eigenstates at $\bm{q}=\bm{k}+\mathbf{G}_{n_{1},n_{2}}$ due to the coupling between Bernal bilayers via Eq. (10), with $\mathbf{G}_{n_{1},n_{2}}=n_{1}\mathbf{G}_{1}^{M}+n_{2}\mathbf{G}_{2}^{M}\left(n_{1},n_{2}\in\mathbb{Z}\right)$. Since we do not consider the intervalley graphene nematicity Samajdar et al. (2021), the calculations are performed with a fixed valley index, e.g., $\eta=+1$. The corresponding band structure for $\eta=-1$ can be obtained by a time-reversal symmetry transformation. For a certain band $n$, the wave functions, truncated up to a wavevector $\mathbf{G}_{n_{1},n_{2}}^{c}$ in the reciprocal lattice, are represented as

[TABLE]

with each term

[TABLE]

containing elements in layer $\ell$ and sublattice $s$ spaces. From these wave functions, the LDOS mappings $\mathcal{D}_{\bm{r}_{0}}(\omega)$ and $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ are computed from

[TABLE]

where $\omega_{n,\bm{k}}$ is the corresponding eigenvalue to the wave function $\Psi_{n}\left(\bm{k}\right)$. This is already projected onto the topmost graphene layer $\ell=4$, where tunneling of electrons from the STM tip are expected to occur in the experimental setup.

To transform the one-dimensional map of LDOS as a funtion of energies $\mathcal{D}_{\bm{r}_{0}}(\omega)=\{\mathcal{D}(\omega_{0}),\cdots,\mathcal{D}(\omega_{N-1})\}_{\bm{r}_{0}}$ into image inputs for the CNN, we use continuous wavelet transforms (CWT) Berthusen et al. (2021); Mallat (1999). These are defined as

[TABLE]

representing the mother wavelet function in a real Morlet form. Here, the transformation is linearly spaced, with $h=1,2,\dots 65$ being equivalent to the spacing in energy taken in the maps $\mathcal{D}_{\bm{r}_{0}}(\omega)$. This scale factor is analogous to frequency in Fourier transforms. Besides $h$, there is a time scale $t$ which is also taken as $t=1,2,\dots 65$, such that $W\left(t,h\right)$ produces $65\times 65$ pixel images. An example of such a “scaleogram” is shown in Fig. 6.

Finally, the LDOS pixel intensities in both maps $\mathcal{D}_{\bm{r}_{0}}(\omega)$ and $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ are modified by the addition of Gaussian random noise via $p_{g}\left(z\right)=\exp(-z^{2}/2\sigma^{2})/\sqrt{2\pi\sigma^{2}}$ with $\sigma=0.31$ Liu et al. (2022b). For $\mathcal{D}_{\bm{r}_{0}}(\omega)$ images, the noise must be added before the CWT for physical consistency.

Appendix B Including strain with fixed $\theta_{\epsilon}$, and variations of the ML architecture

In this appendix, we discuss the changes in the performance of the ML procedure when the training data set or the ML architecture are modified. First, we tested in all cases (Sections II.3-II.5) the performance for the predictions of intensities for fixed angles ($\varphi$, $\alpha$ and $\theta_{\epsilon}$, respectively). For concreteness, we here focus on predicting the microscopic nematic form in the presence of strain (Sec. II.5). The data set for this task is, as before, built by randomly sampling nematic and strain intensities $\Phi^{\text{MN}},\Phi^{\text{GN}}\in[0.001,0.1]$ eV, and $\epsilon\in[0,0.8]\>\%$. Here, $\theta_{\epsilon}=0$, $\alpha=0$, $\psi_{l}=1$ and $\varphi_{MN}=\varphi_{GN}=\varphi=\pi/3$. In this case, all intensities are easily distinguishable and with high accuracy, see Fig. 7. An identical behavior was observed for investigations with fixed $\varphi$ and $\alpha$ for the predictions of the microscopic form of nematicity (Sec. II.4), showing that this is a general feature of the considered CNN architecture. Naturally, in the absence of outliers, the precision of the CNN can be further increased (i.e., by reaching a lower MAE for the predictions) with increasing size and variability of the data set Chollet (2021).

We next address whether the complex architecture in Fig. 2(a) for each channel is really necessary to solve this inverse problem, or in other words, whether simpler architectures could have the same performance and whether modifications of it could produce significant changes in the predictions. For this, we compared the results from the main text with two other architectures in the case of learning the microscopic form of nematicity (Sec. II.4): (i) a very simple sequential neural network that takes the images as inputs, followed by a flatten and dense layer which predicts the parameters $\beta=\{\alpha,\Phi_{\text{GN}},\Phi_{\text{MN}}\}$; (ii) the architecture in Ref. Berthusen et al., 2021, which has a similar structure, i.e., Conv-Batch-MaxPool channel followed by dense layers, but with different number of filters in each layer; we refer to Ref. Berthusen et al., 2021 for details of the architecture. We have found that, even if there is some clear correlations between the true and the predicted values of $\alpha$ in simple architectures, such as (i), it fails completely on predicting the nematic intensities. Additionally, (ii) does not lead to any significant improvement in the predictions.

We also investigated the performance of the CNN with respect to hyperparameter optimization for both Sections II.4 and II.5. These included using different activation functions (SELU, ELU, LeakyReLU, PReLU, ReLU and Sigmoid) Abadi et al. (2015), batch sizes (9, 16, 32, 48 and 96), different number of filters and convolution layers in the Conv-Batch-MaxPool channels, different learning rates ($10^{-1}$, $10^{-2}$, $10^{-3}$ and $10^{-4}$) and optimizers (RMSprop, SGD and ADAM). The architecture described in Sec. II.2 and its variation in Fig. 3(a) already correspond to the optimal configuration. Even though these investigations are not an exhaustive treatment with respect to all possible parameters and correspondent combinations, it shows that certain elements play a major role in the CNN’s performance, such as choosing ReLU as activation functions, setting padding to zero in the convolution layers and using a learning rate of $10^{-4}$. Finally, even though the four Conv-Batch-MaxPool channels in the main architecture may not be necessary for simpler cases (e.g. predicting the nematic director with fixed nematic intensities), it is essential for increasing parameters and complexity, such as learning strain and the internal structure of nematicity simultaneously.

Appendix C Preprocessing of the experimental data and further implications

The experimental data set $D_{\text{exp}}(\beta_{1}^{\prime},\cdots,\beta_{N_{\text{exp}}}^{\prime})$ consists of $N_{\text{exp}}=8$ samples for fillings of the CFB equal to $n_{s}=\{-0.58,-0.45,-0.32,0,0.34,0.47,0.61,0.67\}$. Each sample has $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ images in an energy interval of $\omega\in[-100,100]$ meV with resolution of $2$ meV. From these, the $\mathcal{D}_{\mathbf{r}_{0}}(\omega)$ channels can be calculated by taking an average of intensities at the corresponding BAAC, ABCA and ABAB sites, see Fig. 8(a).

In order to obtain consistent results, the experimental data set $D_{\text{exp}}$ needs to be fed into the trained CNN as similar as possible to the training data in $D_{\text{th}}$. For this, the preprocessing of $D_{\text{exp}}$ consists of:

• (1)

Transforming the experimental plots $\mathcal{D}_{\mathbf{r}_{0}}(\omega)$ into scaleograms as described in Fig. 6. Here, these plots are considered for $\omega\in[-70,60]\,\text{meV}$. This is necessary in order to have scaleograms of $65\times 65$ pixels. In this energy range, these channels contain information about CFB, VFB, RV1 and RC1.

• (2)

Normalizing each image to have the distribution of pixel intensities with same mean $\mu$ and standard deviation $\sigma$ in both $D_{\text{exp}}$ and $D_{\text{th}}$. Here, we have chosen $\mu=0$ and $\sigma=1$ Chollet (2021). This step is essential to produce meaningful predictions on $D_{\text{exp}}$, since the trained CNN have weights associated to the scale of $D_{\text{th}}$.

• (3)

Cropping the images $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ from $D_{\text{exp}}$ such that they show roughly the same number of moiré unit cells as in the corresponding ones in $D_{\text{th}}$. Additionally, the orientations of each of these images in both data setss also need to be consistent pair-wisely, see Fig. 8(b-d) and Fig. 5(c).

We have also investigated the influence of additional preprocessing of the $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ channels of the data set $D_{\text{exp}}$ on the predictions. We introduce more contrast to the images [Fig. 8(c)] and reduce noise by smoothing the pixel distribution with a multidimensional Gaussian filter [Fig. 8(d)]. In Fig. 8(e-f), the resulting predictions for the nematicities and strain using these augmented $D_{\text{exp}}$ are shown. Here, every dot represents an average over the predictions on 10 variations of $D_{\text{exp}}$ with Gaussian filter with standard deviations of the Gaussian kernel in $\sigma_{GF}=\{0,1,2,5,10\}$ with and without higher contrast. The overall behavior described in Sec. II.6 is unaffected by these modifications, but the predictions with the lowest strain intensity in the gray region of Fig. 5(a-b) were found for the raw data inputs [Fig. 8(b)] - see Fig. 9.

We emphasize the importance of the multi-channel CNN architecture in Fig. 3(a) for the predictions in Fig. 5(a-b). While using only $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ at the flat bands already seems to capture the interplay between MN and GN as a function of fillings of the CFB, the addition of channels for the remote bands and scaleograms is crucial to discern the influence of strain and nematicity in the experimental samples. This can be intuitively understood by the relative stability of the remote bands with respect to heterostrain over a wide range of fillings $n_{s}$ in $D_{\text{exp}}$ Rubio-Verdú et al. (2022). These results indicate that the inclusion of more channels from even more remote bands could potentially produce more accurate predictions for the strain intensity; we leave this for future work.

Appendix D Applicability of the CNN to different models and moiré systems

To demonstrate that our ML approach of extracting microscopic parameters based on CNNs with multiple channels works more generally, we here apply it to a different moiré system. To further increase the variability of models studied in this work, we do not use a continuum model but, instead, consider a tight-binding model on the moiré scale that captures the symmetries and topological features of the twisted bilayer graphene (TBG).

As is well known Kang and Vafek (2018); Koshino et al. (2018); Po et al. (2018), the representations of the flat bands of TBG at high-symmetry momenta requires taking a model on the honeycomb lattice. To be able to study the valleys separately, we take the valley quantum number to be conserved such that the associated (fragile) topological obstructions necessitates taking at least four bands Po et al. (2019). We, therefore, place two Wannier orbitals $W_{\pm}(\bm{r})$ at every site of the honeycomb lattice. We choose them to be invariant under $C_{3}$ (three-fold rotation perpendicular to the graphene layers) and transform into one-another under $C_{2x}$ (two-fold rotation along $x$) and $\Theta C_{2}$ (the product of time-reversal and two-fold rotation perpendicular to the layers); this specifies the behavior of the Wannier orbitals under all symmetries of TBG that act within a given valley. Here, our goal is not to provide a quantitatively accurate description of the LDOS of TBG but rather to demonstrate our ML procedure. It is therefore sufficient to take the simple, phenomenological forms of the Wannier states given by

[TABLE]

which obey the required symmetry constraints, as shown in Fig. 10(a). In Eq. (15), $c_{1}$ and $c_{2}$ are real-valued constants that we set to $\{c_{1},c_{2}\}=\{1.5,0.7\}$ for concreteness.

Including symmetry-allowed intra-orbital (inter-orbital) hopping processes up to third-nearest (nearest) neighbor leads to a tight-binding model with momentum-space form (for valley $\eta=+$ and a given spin flavor)

[TABLE]

with orbital Hamiltonians

[TABLE]

coupled via

[TABLE]

with

[TABLE]

and

[TABLE]

where $c^{\dagger}_{\bm{k}}$ are four-component electronic creation operators with the first two (second two) indices referring to the two sublattices of Wannier orbitals $W_{+}$ ($W_{-}$). Here, the two primitive lattice vectors from monolayer graphene are given by $\mathbf{a}_{1}/a=\left(\sqrt{3},1\right)/2$ and $\mathbf{a}_{2}/a=\left(\sqrt{3},-1\right)/2$. In Fig. 10(b,c), we show the band structure for $\{t_{2}/t_{1},\Delta/t_{1},\omega_{1}/t_{1},\omega_{2}/t_{1},\alpha_{2},t_{3}/t_{1},\alpha_{3},\Theta\omega\}=\{0.6,0,0.6,-0.5,0.3\pi,0.1,0.6\pi,0\}$, which are also used in the ML calculations presented below. We will take the lower two, isolated bands [indicated in black in Fig. 10(b)] as a phenomenological description of the quasi-flat bands of TBG; they exhibit Dirac cones at K and K’, which can be shown to have the same chirality—exactly as in TBG.

We next add different forms of nematicity to $h^{0}_{\bm{k}}$, i.e., $h^{0}_{\bm{k}}\rightarrow h_{\bm{k}}=h^{0}_{\bm{k}}+\Delta h_{\bm{k}}$, which have very different structure compared to those discussed in the main text since the model is very different. The nematic order parameter $\bm{\phi}=\left(\phi_{1},\phi_{2}\right)^{T}=\left(\cos 2\varphi,\sin 2\varphi\right)^{T}\in\mathbbm{R}^{2}$ will couple as $\Delta h_{\mathbf{k}}=\bm{\phi}\cdot g_{\mathbf{k}}=\phi_{1}g_{1,\mathbf{k}}+\phi_{2}g_{2,\mathbf{k}}$. Here, the matrix-valued functions $g_{l,\mathbf{k}}$ play a similar role as the tensorial form factor $\bm{\phi}_{\sigma,\ell,s,\eta;\>\sigma^{\prime},\ell^{\prime},s^{\prime},\eta^{\prime}}(\bm{r},\Delta\mathbf{r})$ in the continuum model nematic coupling in Eq. (1). Denoting by $X_{\bm{k}}$ and $Y_{\bm{k}}$ Brillouin-zone-periodic, real-valued functions that transform as $x$ and $y$ under TBG’s point group $D_{3}$, we can write

[TABLE]

where $\alpha_{j}\in\mathbbm{R}$ are parameters and $\sigma_{j}$ ($\rho_{j}$) are Pauli matrices in Wannier (sublattice) space. Technically, the explicit form of $X_{\bm{k}}$ and $Y_{\bm{k}}$ in each of the four terms in Eq. (21) can be different. However, the functional space of possible $X_{\bm{k}}$ and $Y_{\bm{k}}$ is technically infinite dimensional and we will focus only on the leading contribution which then also becomes identical for all four terms in Eq. (21) and reads as

[TABLE]

Consequently, there are four parameters, $\beta=\{\alpha_{0},\alpha_{1},\alpha_{2},\alpha_{3}\}$, describing the microscopic form of nematicity in our model. Our goal will be to reconstruct their values from LDOS images, which we compute via

[TABLE]

with

[TABLE]

The indices $\alpha$ and $\beta$ of $W_{\mathbf{R}_{j\alpha}}$ in Eq. (23) correspond to four different realizations of the Wannier functions in each unit cell $\bm{R}_{j}$, as each of the two orbitals $W_{\pm}$ can be placed on each of the two sublattices.

To reconstruct $\beta=\{\alpha_{0},\alpha_{1},\alpha_{2},\alpha_{3}\}$, we consider a variation of the ML architecture in Fig. 3(a) with four channels for $\mathcal{D}_{\omega_{0}}(\mathbf{r})$ with $\omega_{0}/t_{1}=\{-2,-1,1,2\}$, and one for the scaleograms from $\mathcal{D}_{\mathbf{r}_{0}}(\omega)$. The complete data set consists of 12000 images which are divided into training (78.5%), validation (15%) and test (6.5%) subgroups. These are generated with randomly sampled $\alpha_{j}\in[0.01,0.1]$ for $j=0,\cdots,3$ and for a fixed nematic director $\varphi=5\pi/6$. All the images in the data set were modified by the addition of Gaussian noise with a standard deviation of $\sigma=0.05$. In Fig. 11(a-d), one can see that all four parameters can be accurately predicted. Additionally, we have also observed that even training the CNN with only a single $\mathcal{D}_{\mathbf{r}_{0}}\left(\omega\right)$ channel in this case is sufficient to also yield very good predictions, evidencing the fundamental role of point spectra as an additional source of information. These results indicate that the framework proposed in this work could be successfully applied to a plethora of correlated phenomena and different moiré systems.