Quantum Hall ferroelectrics and nematics in multivalley systems

TL;DR

This paper investigates how anisotropic Fermi surfaces in multivalley quantum Hall systems lead to ferroelectric and nematic broken symmetry states, revealing new insights into their electronic properties and quasiparticle excitations.

Contribution

It demonstrates the formation of quantum Hall ferroelectric states linked to Fermi surface geometry and identifies robust skyrmion quasiparticles in nematic states.

Findings

Ferroelectric states are driven by Fermi surface anisotropy.

Dipole moments relate to the parent metal's Fermi surface geometry.

Skyrmion quasiparticles are notably robust in nematic states.

Abstract

We study broken symmetry states at integer Landau level fillings in multivalley quantum Hall systems whose low energy dispersions are anisotropic. When the Fermi surface of individual pockets lacks twofold rotational symmetry, like in Bismuth (111) and in Sn$_{1-x}$Pb$_x$Se (001) surfaces, interactions tend to drive the formation of quantum Hall ferroelectric states. We demonstrate that the dipole moment in these states has an intimate relation to the Fermi surface geometry of the parent metal. In quantum Hall nematic states, like those arising in AlAs quantum wells, we demonstrate the existence of unusually robust skyrmion quasiparticles.

Figures

Equations

$$v({\mathbf{r}}_1-{\mathbf{r}}_2)=\frac{e^2}{ϵ|{\mathbf{r}}_1-{\mathbf{r}}_2|}.$$

$$P_{LL}v({\mathbf{r}}_1-{\mathbf{r}}_2)P_{LL}=\frac{1}{A}\sum _{\mathbf{q}}{v}_{\mathbf{q}}{F}_I(\mathbf{q})F_J(-\mathbf{q})e^{i\mathbf{q}\cdot ({\mathbf{R}}_1-{\mathbf{R}}_2)},$$

$$|\Psi ⟩=\prod _{a=1}^{\nu }\prod _{k=1}^{N_{\varphi }}\left(\sum _{I=1}^M⟨{\chi }_a|I⟩c_{kI}^{†}\right)|\mathrm{O}⟩,$$

$$\begin{array}{rl}& E[P]=-N_{\varphi }\sum _{I,J}{X}_{IJ}⟨I|P|J⟩⟨J|P|I⟩,\\ & P=\sum _{a=1}^{\nu }|{\chi }_a⟩⟨{\chi }_a|,X_{IJ}=\frac{1}{2}\int \frac{d^2\mathbf{q}}{(2\pi {)}^2}{v}_{\mathbf{q}}{F}_I(\mathbf{q})F_J^{\ast }(\mathbf{q}),\end{array}$$

$$⟨F_I,F_J⟩\equiv X_{IJ}=\frac{1}{2}\int \frac{d^2\mathbf{q}}{(2\pi {)}^2}{v}_{\mathbf{q}}{F}_I(\mathbf{q})F_J^{\ast }(\mathbf{q}).$$

$$E=-N_{\varphi }\sum _{IJ}{p}_I{X}_{IJ}{p}_J,p_I=|⟨I|\chi ⟩{|}^2,$$

$${\mathbf{p}}_1\cdot {\mathbf{p}}_2\equiv \sum _{IJ}{p}_{1I}{X}_{IJ}{p}_{2J}.$$

$$\begin{array}{rl}||\mathbf{p}(t)|{|}^2& \equiv \mathbf{p}(t)\cdot \mathbf{p}(t)=\\ & ||{\mathbf{p}}_1-{\mathbf{p}}_2|{|}^2{t}^2+2t({\mathbf{p}}_1\cdot {\mathbf{p}}_2-||{\mathbf{p}}_2|{|}^2)+||{\mathbf{p}}_2|{|}^2.\end{array}$$

$$\mathbf{D}=-|e|(⟨k,J|\mathbf{r}|k,J⟩-⟨k,J|\mathbf{r}|k,J{⟩}_0).$$

$$\mathbf{r}\equiv \mathbf{R}-\hat{\mathbf{z}}\times \mathbf{p}.$$

$$\mathbf{D}=-|e|l^2\hat{\mathbf{z}}\times ⟨J|\mathbf{p}|J⟩.$$

$${\psi }_{k_y}(x,y)=\frac{e^{i{k}_{y}y}}{\sqrt{L_y}}{u}_{k_y}(x).$$

$$I_y=-|e|\frac{\partial H_{k_y}}{\partial k_y},$$

$$D_y=-i|e|\frac{l^2}{L_x}{\int }_0^{L_x/l^2}d{k}_y⟨u_{k_y}|{\partial }_{k_y}|u_{k_y}⟩.$$

$$\begin{array}{rl}{D}_y& =-i|e|⟨u_{k_y}|{\partial }_{k_y}|u_{k_y}⟩{|}_{k_y=0}\\ & =i|e|l^2\int dx{u}_0^{\ast }(x){\partial }_x{u}_0(x).\end{array}$$

$$\begin{array}{rl}{D}_y& =i|e|l^2\int d^{2}r{\psi }_{k_y}^{\ast }(x,y){\partial }_x{\psi }_{k_y}(x,y)=|e|l^2⟨{\psi }_{k_y}|p_x|{\psi }_{k_y}⟩\\ & =-|e|⟨{\psi }_{k_y}|(y-R_y)|{\psi }_{k_y}⟩.\end{array}$$

$$H\left(\frac{\hslash }{i}\frac{d}{dx},k_y-eBx\right)u_{k_y}(x)={ϵ}_n{u}_{k_y}(x).$$

$$u_0(x)=\exp \left(\frac{i}{\hslash }{\varphi }_0(x)+{\varphi }_1(x)+\mathcal{O}(\hslash )\right),$$

$$u_0(x)\approx \sum _s\frac{1}{\sqrt{v_{x,s}^{\mathrm{cl}}(x,ϵ)}}\exp \left(\frac{i}{\hslash }{\int }^{x}d{x}^{\prime }{p}_{x,s}^{\mathrm{cl}}(x^{\prime },ϵ)\right),$$

$$v_{x,s}^{\mathrm{cl}}(x,ϵ)=\frac{\partial H(p_x,-eBx)}{\partial p_x}{|}_{p_x\to p_{x,s}^{\mathrm{cl}}(x,ϵ)},$$

$$\sum _s\frac{s}{eB}{\int }_{p_y^{-}}^{p_y^{+}}d{p}_y{p}_{x,s}^{\mathrm{cl}}(p_y,ϵ)=l^{2}A({ϵ}_n)=2\pi (n+1/2).$$

$$\begin{array}{rl}& \int dx\mathrm{Im}(u_0^{\ast }(x){\partial }_x{u}_0(x))\approx \\ & -\sum _s\frac{s}{eB}{\int }_{p_y^{-}}^{p_y^{+}}d{p}_y\frac{1}{|v_{x,s}^{\mathrm{cl}}(p_y,ϵ)|}{p}_{x,s}^{\mathrm{cl}}(p_y,ϵ)\\ & =-\oint dt{p}_x(t,ϵ).\end{array}$$

$$\int dx{u}^{\ast }(x)u(x)\approx \sum _s\frac{1}{eB}{\int }_{p_y^{-}}^{p_y^{+}}d{p}_y\frac{1}{|v_{x,s}^{\mathrm{cl}}(p_y,ϵ)|}=\oint dt.$$

$$\mathbf{D}\approx -|e|⟨\mathbf{r}(t,{ϵ}_n){⟩}_t\equiv -|e|\frac{\oint dt\mathbf{r}(t,{ϵ}_n)}{\oint dt},$$

$$H_0=\pm v_x{\sigma }_x{p}_x+v_y{\sigma }_y{p}_x\pm \delta v_x{p}_x,$$

$$E/A=-\beta n_z^2,$$

$$\begin{array}{rl}& \beta =\frac{{\tau }^2}{4\pi }\int \frac{d^{2}q}{(2\pi {)}^2}{v}_q{e}^{-|\bar{q}{|}^2/2}{\left(3-\frac{|\bar{q}{|}^2}{2}\right)}^2{\bar{q}}_y^2,\end{array}$$

$$\begin{array}{rl}& D_y\approx \pm |el|\mathrm{sign}(B)\frac{3\delta v_x}{\sqrt{2}{v}_y},D_x=0.\end{array}$$

$${\rho }_{qh}(\mathbf{r})=\frac{1}{2\pi l^2}-|{\psi }_{k,J}(\mathbf{r}){|}^2,$$

$$E/A=-\alpha n_z^2.$$

Full text

I Introduction

Quantum Hall liquids with nearly degenerate internal degrees of freedom have long been the source of a rich variety of phenomena. Aside from multilayer systems, there exist today an array of multi-valley two-dimensional electron systems exhibiting quantum Hall effect, including quantum wells of AlAs AlAsrev , monolayer and bilayer graphene Kim ; Novoselov , Silicon surfaces Si111 , PbTe (111) surfaces PbTe , surface states of the topological crystalline insulator Sn1-xPbxSe Madhavan , and more recently the (111) surface of Bismuth Yazdani . These systems all have multiple valleys related to each other by discrete crystal symmetries.

Among these multivalley systems, AlAs heterostructures, Si surfaces, PbTe (111) quantum wells and Bi(111) surfaces all have highly anisotropic pockets, whose orientations are valley-dependent. The shape of the Fermi surface is crucial in determining the pattern of valley symmetry breaking. For AlAs, and PbTe, the pockets are centered at time-reversal-invariant momenta, and therefore posses an elliptical shape with two-fold symmetry. This Fermi surface anisotropy favors valley-polarized states where a subset of valley degenerate Landau levels are fully occupied  Abanin ; AllanSnTe . This state spontaneously breaks the larger crystal rotational symmetry down to twofold rotations, and therefore the resulting state is a nematic quantum Hall state.

However, the recently studied Bi(111) surface and Sn1-xPbxSe (001) surface, brings in a novel ingredient into this problem. As shown in (Fig. 1), Bi (111) surface has six tadpole-shaped hole pockets Bi2001 ; Bi111 ; Sn1-xPbxSe (001) surface has four crescent-shaped pockets PbSnSe . These pockets come in time-reversed pairs located at opposite momenta away from time-reversal-invariant points in the Brillouin zone. Importantly, the Fermi surface of each pocket does not have a twofold rotational symmetry. Therefore, we expect Landau orbitals associated with each valley generically to carry an in-plane dipolar component. This implies that a valley-polarized state will be an insulator that spontaneously breaks inversion symmetry. We will refer to this phase as a quantum Hall ferroelectric state. The giant fermi surface anisotropy of Bi(111) and Sn1-xPbxSe (001) surface states makes them promising candidate systems to study this novel ferroelectricity in quantum Hall states (Fig. 1).

In this letter we provide a unified description of these ferroelectric and nematic states in multi-valley integer quantum Hall systems. We establish that due to the Fermi surface anisotropy, the long-range Coulomb interaction generically favors states with full valley polarization. We compute the electric polarization of quantum Hall ferroelectrics using a quantum mechanical approach based on Berry phase and also using a semiclassical approach that directly relates the dipole moment to the underlying Fermi surface geometry. For quantum Hall nematic states, we study the energetics of Skyrmion-type charged excitations Sondhi using the density-matrix renormalization group method, and find they are suprisingly robust against mass anisotropy of the valleys, shedding new insights on the experiments on AlAs quantum wells AlAsskyrm .

II General setup

We start from symmetry considerations on two-dimensional (2D) multivalley systems. We consider systems with three, four, or six fold rotational symmetry, which rules out any principal axis within the 2D plane. In order for each valley to have an anisotropic dispersion, we require that the symmetry group that leaves each valley invariant, or the “little group”, is only a subgroup of $C_{2v}$, or equivalently, contains no more than a twofold rotation $(x,y)\rightarrow(-x,-y)$ and at most two mirror planes that are orthogonal to each other, $(x,y)\rightarrow(-x,y)$ and $(x,y)\rightarrow(x,-y)$.

Within this wide class of systems it is convenient to distinguish two types two-dimensional systems with multiple anisotropic valleys: type-I are those whose little group contains at most a single mirror plane; type-II are those with larger symmetry. For systems of type-I, each valley is of such low symmetry that the electron dispersion at zero field has no inversion center, i.e., $\epsilon({\bf k})\neq\epsilon(-{\bf k})$ where ${\bf k}$ is the “small” momentum within the valley. As we will show, in the quantum Hall regime and at odd-integer fillings, Coulomb interaction tend to induce a valley-polarized state which breaks all rotational symmetry of the crystal and is therefore a ferroelectric. For systems of the type-II, the electron dispersions at each valley have twofold symmetry $\epsilon({\bf k})=\epsilon(-{\bf k})$, and the resulting valley polarized state will retain $\pi$ rotations. In this case, valley-polarized quantum Hall states exhibit nematic instead of ferroelectric orders. The Bi(111) surface Yazdani is a representative of type-I systems and the AlAs quantum well AlAsrev a representative of type-II.

To study interacting electrons in multi-valley quantum Hall systems we assume that the magnetic field is large enough so that the Hamiltonian can be projected into a set of $M$ degenerate Landau levels associated with different valleys. The system is at some integer filling $\nu\in\mathbb{Z}$ ($\nu\leq M$). The interaction between an electron with position ${\bf r}_{1}$ from valley $i$ and an electron with position ${\bf r}_{2}$ from valley $j$ is dominated by the long-range part of the Coulomb interaction 111We neglect short distance valley dependent corrections. Note also that $v$ is independent of the relative orientation of the two electrons, because systems considered here all have isotropic dielectric properties.:

[TABLE]

Despite the bare Coulomb interaction being valley independent, it effectively becomes valley dependent after projection into the Landau levels of interest:

[TABLE]

where ${\bf R}\equiv{\bf r}-l^{2}\hat{\bf z}\times{\bf p}$ is the intra-Landau level guiding center operator, ${\bf p}={\bf\nabla}/i-e{\bf A}$ is the mechanical momentum, and $F_{J}({\bf q})\equiv\langle J|e^{-il^{2}\hat{\bf z}\cdot{\bf q}\times{\bf p}}|J\rangle$ is the form factor determined by the wavefunction, $|J\rangle$, of the Landau level of interest associated with valley $J$ ($\hbar\equiv 1$, $l^{2}\equiv 1/eB$). These form factors are crucial for the energetics of valley symmetry breaking quantum Hall states and are described for specific cases of type-I and type-II multi-valley systems in Appendices A, B and C.

We would like to note that our approach rests on the assumption that the long range Coulomb force projected into a set of valley-degenerate Landau levels is the dominant term that is able to break the many-body degeneracy and select the ground states. This assumption is justified in systems with anisotropic Fermi surfaces (which are the subject of this work), provided that short-range interactions are small. For systems with nearly isotropic Fermi surfaces such as graphene, long range Coulomb interaction leaves large remnant continuous symmetries, hence short distance corrections are needed to select the ground state (see e.g. Refs. Kharitonov, ; Murthy, ). Additionally, for massless Dirac fermions systems the Landau level projection is only perturbatively enforced by the smallness of the effective fine structure constant and hence corrections beyond naive degenerate perturbation theory that include the role of the negative energy sea might be important when such parameter is not small.

III Hartree-Fock Theory

We study the Hamiltonian (2) using both the Hartree-Fock approximation and the DMRG numerical method. Within Hartree-Fock, we consider translationally invariant trial Slater determinant states that are arbitrary coherent superpositions of the $M$ valleys. These states can be parametrized as follows:

[TABLE]

where $|{\rm O}\rangle$ is the reference vacuum in which the Landau levels of interest are empty, $k$ is an intra Landau level guiding center label, $I$ labels the $M$ valleys, $|\chi_{a}\rangle$ are $N$ orthonormal vectors describing the occupied coherent combinations of such valleys. The expectation value of the energy (with the Hartree part removed by including a neutralizing background) is:

[TABLE]

where $I$ labels the valleys and $|\chi_{a}\rangle$ are $\nu$ orthonormal vectors describing the occupied coherent combinations of the valleys. Trial states are found by minimizing $E[P]$ as a function of $P$.

It is well known that generically at integer fillings the ground states of multicomponent systems with repulsive interactions are spontaneous integer quantum Hall states with coherence between different components known as quantum Hall ferromagnets. We have found a generalized and stronger version of quantum Hall ferromagnetism at $\nu=1$ for the kind of multicomponent systems with anisotropic fermi surfaces that we are considering. In particular, we will show that the state that minimizes the Hartree-Fock energy is maximally polarized into a single valley and those states with coherence between different valleys are energetically disfavored by the exchange energy. In fact, we will prove the following theorem:

Full valley polarization theorem— If the set of form factors $F_{I}(\bf q)$, $I=\{1,...,\nu\}$, are linearly independent functions of ${\bf q}$ and the interaction is strictly repulsive, $v_{\bf q}>0$ for all ${\bf q}$, the minimum of Hartree-Fock energy $E[P]$ at $\nu=1$ for any multi-valley system is a state maximally polarized into a single valley.

Proof—We begin by noting that provided the interactions are repulsive, $v_{\bf q}>0$, the exchange integral defines an inner product for form factors as functions of ${\bf q}$:

[TABLE]

$X_{IJ}$ is a real symmetric matrix whose entries are the inner products, known in linear algebra as a Gram matrix. A classic theorem establishes that if the vectors used in the Gram matrix are linearly independent, the matrix is strictly positive definite, namely, all its eigenvalues are real and positive. Therefore, we conclude that the linear independence of $F_{I}({\bf q})$ implies that $X_{IJ}$ is a strictly positive matrix.

On the other hand, at $\nu=1$, the expression for the exchange energy reduces to:

[TABLE]

where $|\chi\rangle$ is the single trial coherent combination of the valleys that is occupied at $\nu=1$. The task is to find the minimum of this energy as a function of probabilities $p_{I}$ that define a point in $\mathbb{R}^{M}$ that is further restricted to a region defined by the constraints $\sum_{I}p_{I}=1$, $p_{I}\in[0,1]$. Let us denote this region by $\mathbf{P}$, which is illustrated in Fig. 2. Now, since $X_{IJ}$ is strictly positive definite, it defines itself an inner product in ${\mathbb{R}}^{M}$:

[TABLE]

Now consider a line in ${\mathbb{R}}^{M}$ parametrized by $t\in[0,1]$, ${\bf p}(t)=t{\bf p}_{1}+(1-t){\bf p}_{2}$, connecting two distinct points, ${\bf p}_{1}\neq{\bf p}_{2}$. It is easy to show that if ${\bf p}_{1}\in{\bf P}$ and ${\bf p}_{2}\in{\bf P}$, then ${\bf p}(t)\in{\bf P}$, $\forall t\in[0,1]$. Notice that:

[TABLE]

Since the coefficient of $t^{2}$ is strictly positive, then, the maximum of $||{\bf p}(t)||^{2}$ is always achieved at the end points of the line and never strictly within its interior. Since the energy is $E=-N_{\phi}||{\bf p}||^{2}$, it follows that any point that can be viewed as being strictly in the interior of a line contained in ${\bf P}$ cannot be a minimum of the energy. Therefore the minima of the energy must necessarily be located at the corners of $\mathbf{P}$, because these are the only points which cannot be viewed as being strictly in the interior of a line contained within ${\bf P}$. The $J$-th corner is defined by $p_{J}=1$ and $p_{I}=0$ for $I\neq J$. The state with the minimum energy is therefore fully valley polarized and the valley with minimum energy will be the one with the largest exchange integral $X_{JJ}$. In the case in which the valleys are related by a discrete crystal symmetry the exchange integral $X_{JJ}$ will be the same for every $J$ and we get a discrete degeneracy of all the fully polarized states. And thus we conclude our proof that the system spontaneously fully polarizes into a single valley.

A simple way to test that form factors are linearly independent is to check that $\det(X_{IJ})\neq 0$. It is easy to imagine that form factors that arise from valleys with anisotropic Fermi surfaces are generically linearly independent. This behavior is the generic case for valleys with anisotropic Fermi surfaces that are distinct. We use the word distinct here in the strong sense that the fermi surfaces are not identical upon translation in momentum space. For example, two elliptical Fermi surfaces that are rotated with respect to each other, as in AlAs, are distinct. Similarly, all the Fermi surfaces at the Bi(111) surface are distinct (Fig. 1). Therefore, we expect generically that the ground states of type-I and type-II systems are fully valley polarized ferroelectric or nematic states at $\nu=1$.

It is possible, however, that two or more Fermi surfaces for distinct valleys are, in a small momentum approximation, identical (as in the case of graphene within the leading linear $k$ “dot” $p$ Dirac approximation), hence the form factors are linearly dependent. These system in the presence of long-range Coulomb interaction have an emergent SU($N$) symmetry, between the $N$ valleys with identical form factors. Also, an important condition for the preceding analysis to apply is that all electrons reside in the same layer so that the interactions before projecting to the Landau level are identical. Otherwise, non-trivial Hartree terms may favor superpositions of components in different layers, as it happens for the exciton condensate in quantum Hall bilayers Moon .

In fact, the fully valley polarized states are exact eigenstates of the many-body Hamiltonian projected into the Landau levels of interest. This follows from the constraint of valley number conservation in our model implies that, after projection into the multivalley Landau level of interest, they have no other states in the Hilbert space to mix with. It is more subtle to demonstrate that they are the lowest energy states. As we have seen, however, this is the expectation based on Hartree-Fock theory. We have confirmed via density matrix renormalization group numerical simulations on the torus geometry Shibata2001 ; Feiguin2008 ; Xiang1996 ; White1992 ; Zhao2011 , that they are indeed the exact ground states, as discussed in Section VII.

Although the valley polarization theorem predicts fully valley polarized states in a general class of systems one should be mindful of potential subtleties in specific systems that would require going beyond the simple model we have described. One consideration is the importance of short-distance corrections to the Coulomb interaction, which typically are reduced from the energy scale of the Coulomb interaction by a factor of the order of $a/l$, where $a$ is a lattice length scale. A simple criterion for when these terms are important is to compare their size with the energy difference per particle between a valley polarized state and valley coherent state within the description we have provided. When such difference is small the tendency to select a unique ground state by the Coulomb interaction alone would be weak. In such cases one should consider the role of these short interactions in selecting ground states on an equal footing to the asymmetry of the form factors we have described, in an analogous fashion to how it has been done in Graphene Kharitonov .

IV Quantum Hall ferroelectric dipole moment

Following our full valley polarization theorem, we expect Type-I systems to exhibit valley polarized states at filling $\nu=1$ that fully breaks rotational symmetry. The presence of a macroscopic dipole moment is a subtle issue. Unlike, conventional ferroelectrics, quantum Hall ferroelectrics are always accompanied by a chiral metallic edge. Therefore one expects that any polarization charge that accumulates at the boundary of the sample to be screened by the metallic edge 222We thank A.H. MacDonald for pointing out to us this subtlety on the measurement of the Dipole moment for this state.. In the forthcoming discussion we will ignore this subtlety in the measurement of the dipole moment. We will discuss in Section VI.1 how the dipole moment that we will compute can be experimentally measured from the charge distribution of quasi-particles.

IV.1 General expressions for the dipole moment

The dipole of the insulating state at $\nu=1$ can be defined as the change of the average position of the electrons relative to a reference state with inversion symmetry. Since the state at $\nu=1$ is a Slater determinant the dipole moment per electron can be computed for each single orbital as a single particle property. We will illustrate two different ways of computing the dipole moment which result in identical final results but are useful as they illustrate different points of view on the quantum mechanical theory of polarization.

The first approach starts from the expression that the dipole moment associated with an orbital $k$ in a valley $J$ is given by the change of the expectation value of the position operator:

[TABLE]

Here $k$ is an intra-Landau level guiding center label and $J$ labels the valley that is spontaneously chosen in the ground state. For purposes of defining the dipole moment we consider a hypothetical reference hamiltonian that has inversion symmetry, whose eigenfunctions are $|k,J\rangle_{0}$, and imagine that this Hamiltonian can be adiabatically deformed into the Hamiltonian of interest describing the dispersion of valley $J$ which breaks inversion symmetry. The dipole is obtained as the change in the position relative to the reference inversion symmetric state. The position operator can be decomposed into cyclotron and guiding center components as follows:

[TABLE]

where ${\bf p}=\nabla/i-{\bf A}({\bf r})$ is the gauge invariant mechanical momentum. The key observation is that the guiding center variables are adiabatic invariants because guiding center operators commute with the mechanical momentum, $[R_{a},p_{b}]=0$, and the Hamiltonian describing each valley is a function of ${\bf p}$ throughout the entire adiabatic path because the path is assumed not to break translational symmetry. Thus, by choosing the origin of coordinates to be the initial position of the orbital, so that $\langle k,J|{\bf r}|k,J\rangle_{0}=0$, we arrive at the following expression for the dipole moment:

[TABLE]

Here we have omitted explicit reference to the guiding center label of the eigenstate in question as the formula makes manifest that the dipole moment is independent from it. The formula just derived offers great versatility for computations of the dipole moment for specific systems when the cyclotron eigenstates are known explicitly in terms of the canonical cyclotron raising and lowering operators, because ${\bf p}$ is a simple linear function of these operators. We will take advantage of this to perform a swift computation of the dipole moment in tilted Dirac cones that are relevant for the surface of topological crystalline insulators described in Section V.

The second approach follows from adapting the formalism of the modern theory of polarization Thouless ; Vanderbilt to our current problem. To compute the dipole moment along a given direction within this approach, we chose a gauge so that there is translational invariance along such direction. For example, to compute the dipole along $y$, we choose $A_{x}=0$ and $A_{y}=Bx$. In this case the single particle states within the Landau level of interest can be labeled by a momentum $k_{y}$ and have a real space form:

[TABLE]

Periodicity along the $y$ direction imposes a discretization of $k_{y}$ ($\Delta k_{y}=2\pi/L_{y}$), and the finite size and the periodicity of the system along the $x$ direction determines the size of the effective the Brillouin zone to be: $k_{y}\in[0,L_{x}/l^{2}]$. To determine the polarization we compute the integral over time of the current operator:

[TABLE]

in response to an adiabatically changing Dirac cone tilt $\delta v_{x}(t)$ over a period of time $T$, so that $\delta v_{x}(t=0)=0$ and $\delta v_{x}(t=T)=\delta v_{x}$. Following the arguments of Refs. Thouless, ; Vanderbilt, , one finds the following expression for the change in the dipole moment per particle:

[TABLE]

In the present case, the real space wave-functions at different $k_{y}$ differ by a translation: $u_{k_{y}}(x)=u_{0}(x-k_{y}l^{2})$. As a consequence the formula reduces to:

[TABLE]

To make manifest the equivalence of the two approaches to compute the dipole moment, we recast Eq. (15) as follows:

[TABLE]

The last expression makes manifest that the dipole moment computed in this way coincides with the average position of the single particle orbital measured with respect to the guiding center $R_{y}$.

IV.2 WKB approximation for the dipole moment

We now introduce a semiclassical approach to establish a direct connection of the dipole moment with the underlying Fermi surface geometry at zero field. Consider a valley described by a single-band Hamiltonian $H(p_{x},p_{y})$ of arbitrary form. For simplicity we will imagine that the eigenstates near this valley have negligible Berry phase, or equivalently that the Hamiltonian has only differential operators but no pseudo-spin matrix structure. This is typically a good approximation near the bottom of a band for a single-band system provided other bands are sufficiently far in energy. As in the previous subsection, by choose the Landau gauge $A_{x}=0$ and $A_{y}=Bx$, we view the eigenvalue problem as effectively one-dimensional. The eigenvalue problem reads as:

[TABLE]

without loss of generality we can set $k_{y}=0$, since other solutions are obtained by a global translation. We can search for an approximate WKB solution to this differential equation which formally is an expansion in $\hbar$:

[TABLE]

substituting this expression into Eq. (17) and keeping corrections up to linear order in $\hbar$ leads to the WKB approximation of the problem Falkovs :

[TABLE]

where $p^{\rm cl}_{x,s}(x^{\prime},\epsilon)$ is the $s$ root of the algebraic equation defined by the classical dispersion relation: $H(p^{\rm cl}_{x,s},-eBx)=\epsilon$, and

[TABLE]

is the group velocity evaluated on the classical trajectory. Typically we have two roots $s={+,-}$ and two turning points $x_{+,-}$ that separate the classically allowed from the classically forbidden regions. In the classically forbidden regions $p^{\rm cl}_{x,s}$ is imaginary and the wavefunction decays exponentially. As is well known, the quantization condition for the WKB solution is obtained by matching boundary conditions between classically allowed and forbidden regions and leads to the Onsager quantization condition, which, in the absence of Berry phase, reduces to:

[TABLE]

where we made the change of variables $p_{y}=-eBx$, $A(\epsilon_{n})$ is the area in momentum space inside the iso-energetic contour $H(k_{x},k_{y})=\epsilon_{n}$ and $n$ is the Landau level index. From Eq. (15) we write the Dipole as $D_{y}=-el^{2}\int dx\ {\rm Im}(u^{*}_{0}(x)\partial_{x}u_{0}(x))$. There are two types of contributions to the dipole moment, namely the sum of the contributions coming from each classical root separately and the contribution corresponding to the interference-like crossed terms between the two classical roots. The semiclassical approximation is expected to work when the oscillations of the wavefunction occur over a length much smaller over that in which the confining potential changes, therefore we restrict to states that have a high energy and hence a large Landau level index. For these states the crossed terms have rapid oscillatory factors of the form $\exp\left(i/\hbar\int^{x}dx^{\prime}(p^{\rm cl}_{x,1}-p^{\rm cl}_{x,2})\right)$ and can be neglected. In addition in this case the integrals can be approximated to be taken over the classically allowed region, because the states with rapid oscillations in the classically allowed region correspondingly have fast decay in the classically forbidden region. Therefore the terms in which the oscillatory part cancels within the classically allowed region have the dominant contribution. As a consquence we obtain:

[TABLE]

where in the last line we have established the connection to the classical cyclotron motion described in the previous section, in which $d{\bf r}(t)/dt={\bf v}=l^{2}d{\bf k}(t)/dt\times\hat{z}$, from which we have that $v_{x}=l^{2}dk_{y}/dt$. Additionally we must normalize the WKB solution. Using the same approximations employed to compute the dipole one obtains that the approximate normalization of the WKB solution is:

[TABLE]

Thus the normalization is simply the total time that takes to complete the classical cyclotron orbit. A similar analysis can be followed to find the dipole along the $x$ direction $D_{x}$. Therefore within the approximations here outlined, the WKB quantized wave-function predicts that the dipole moment is simply given by time averaged position of the electron over the semiclassical orbit:

[TABLE]

where the integral is perfomed over the semiclassical cyclotron orbit ${\bf r}(t)=l^{2}{\bf k}(t)\times\hat{z}$, with ${\bf k}(t)$ tracing a constant energy contour $H({\bf k}(t))=\epsilon_{n}$ at a speed $v=|\partial\epsilon/\partial{\bf k}|$ Kittel . This picture predicts intuitively and generically that the dipole is orthogonal to the direction of the distortion of the Fermi surface since the real space orbit is a rotated version of the Fermi surface.

V Quantum Hall Ferroelectrics at the surface of Topological crystalline insulators

There are various material platforms for multi-valley systems of the type-I, where anisotropic valleys are located at non-time-reversal-invariant momenta and each valley lacks twofold symmetry. One candidate platform are multi-valley systems on the (111) surface of bismuth Bi2001 ; Bi111 . Another interesting platform, on which we will focus in this section, is the (001) surface of topological crystalline insulators (TCI) SnTe, SnxPb1-xSe and SnxPb1-xTe LiuSnTe . In both Bi(111) and in TCI’s, surface state Landau levels have been observed by means of scanning tunneling microscopy Yazdani ; Madhavan .

To study ferroelectricity in at the surface of TCI’s, we consider a model with two valleys ($M=2$) of Dirac fermions located at opposite momenta, at filling $\nu=1$. This model describes the low-energy dispersion of the [001] surface of SnTe and PbSnTe topological crystalline insulators in the low temperature phase Serbyn . Here, the Dirac cones are generically “tilted” Dipole and described by the following Hamiltonian:

[TABLE]

The only discrete symmetry that leaves each valley invariant is a single mirror plane along the line that connects the Dirac cones Serbyn ; Dipole . We consider the $n=1$ Landau level 333The zeroth Landau level does not carry dipole even in the presence of the tilt because the Hamiltonian retains an intra-valley particle-hole times inversion symmetry which forbids a dipole moment for the zeroth level. However, all other Landau levels will carry a dipole.. The wavefunctions associated with each valley break rotation and inversion symmetry as illustrated in Fig. 3. Using the form factors described in Appendix C and applying Eq. (4) to the case of $M=1$ at $\nu=1$, the Hartree-Fock energy can be found to be (up to a global constant):

[TABLE]

where $n_{z}={\rm tr}(\sigma_{z}P)$ and $\sigma_{z}$ is a Pauli matrix in the valley indices, and $A$ is the system area. $\beta$ is a positive constant, in agreement with the theorem of full valley polarization, and it is explicitly given by:

[TABLE]

where $\tau=\delta v_{x}/(\sqrt{2}v_{y})$, and $(\bar{q}_{x},\bar{q}_{y})=(\sqrt{v_{y}/v_{x}}q_{y},\sqrt{v_{x}/v_{y}}q_{x})$. Therefore, the ground state is an Ising type ($n_{z}=\pm 1$) ferroelectric. The dipole moment, to leading order in the Dirac cone tilt, can be computed using the form of the Landau levels in terms of cyclotron raising and lowering operators (Appendix C), and using the Eq. (11), and it is found to be:

[TABLE]

The dipole of each valley is orthogonal to the direction of the Dirac cone tilt and reverses with the direction of the perpendicular field. This is allowed since the magnetic field breaks the mirror symmetry present at zero field. In our discussion we have neglected the Zeeman coupling. This approximation is justified when the Zeeman energy is much smaller than the cyclotron spacing of the Dirac Landau levels in question, but including the Zeeman term should not change qualitatively the results in these systems.

VI Experimental manifestations of quantum Hall ferroelectrics

We now discuss potential experimental manifestations of the quantum Hall ferroelectrics. In the case of Bi(111) surface Yazdani , we expect these states to appear at odd integer fillings. The degree of broken inversion asymmetry in this system is expected to increase with Landau level index because the lowest Landau levels might be well described by an elliptical dispersion near the bottom of the hole bands, for which no dipole moment is expected, but only nematic states reported experimentally Yazdani . In this system, because the system has coexisting ferroelectric and nematic character, one signature is the appearance of additional ferroelectric domain walls in addition to those associated with nematicity seen at even fillings. One interesting possibility is the existence of gapless edge states at the ferroelectric domain walls. If one neglects inter-valley scattering, the valley pseudospin is conserved and one expects gapless charge carrying counter-propagating edge modes that arise at the domain walls. Due to their ferroelectric nature, these domains may be manipulated by STM bias voltage or in-plane external electric field. Finally, we note that Bi (111) hole pockets located at opposite momenta carry opposite in-plane spin-polarizations, which cancel when the two are equally occupied. Therefore, valley-polarized states at odd-integer filling also carry an in-plane spin polarization that could be manipulated with in-plane magnetic fields. However, the broken inversion symmetry of this state would manifest directly in the charge distribution of quasiparticles which would carry a dipole moment given by the formulae described in section IV. We note in passing that other types of quantum Hall ferroelectrics arising in wide quantum wells widewells and in bilayer graphene Barlas have been previously proposed, but to our knowledge these proposals have not been experimentally realized so far.

VI.1 Measurement and physical meaning of the Dipole moment

The dipole moment we have computed is defined as the change of the average position of every electron in a Landau level along a path that connects the Hamiltonian of interest with a reference initial Hamiltonian that has inversion symmetry. The argument presented in Section IV.1 demonstrates that, as long as there is translational invariance, such a change is independent of the guiding center variables and hence it is the same for any orbital constructed within a given Landau level.

As we noted in section IV, the presence of the metallic edge of quantum Hall states prevents the build-up of charge at the boundaries of quantum Hall ferroelectrics. This may invalidate conventional approaches to measure the polarization used in ordinary ferroelectrics that rely on the measurement of macroscopic charge accumulation at the boundaries. Experimental manifestations of the formation of this inversion symmetry breaking quantum Hall state, therefore require looking at other kinds of observables. One instance of the manifestation of the broken inversion symmetry will be the charge distribution of quasiparticles. One consequence of inversion symmetry is that the charge distribution of elementary quasiparticles will be inversion symmetric. However, this will no longer be true when the system breaks inversion symmetry. Consider, the elementary quasi-hole of the system which is created by removing the electron from a single orbital $\psi_{k,J}({\bf r})$, with an intra-Landau label $k$, in the Slater determinant associated with completely filling a Landau level that we label by $J$. The particle density in the quasi-hole state is:

[TABLE]

where $1/{2\pi l^{2}}$ is the background density of the $\nu=1$ state. Now, from this equation we can see that the change of the average position of the quasi-hole state computed along an adiabatic path that connects the Hamiltonian of interest with a reference inversion symmetric Hamiltonian will coincide with the change of position of the single particle state $\psi_{k,J}({\bf r})$ with an overall opposite sign. Therefore the change of the dipole moment of the quasi-hole is exactly minus the dipole moment that we have described in Section IV.

Broadly speaking, the quasi-hole can be considered as a limiting version of a hole inside the sample that separates the quantum Hall ferroelectric and vacuum. When the size of this hole is very large there are gapless excitations that appear accompanying the formation of the chiral edge that will screen the build up of charge. However, when the size of the hole is small or comparable to the magnetic length these excitations will acquire a finite size gap and we expect that dipole moments will build up in such limit. Therefore, a dipole moment should build up for cavities or holes inside the sample that are on the order of the magnetic length, in resemblance to the cavity electric fields that are used as a classic method in elementary discussions of dielectrics as a way to measure the electric dipole moment Zangwill . Strictly speaking a cavity is different from a free quasi-hole because there is an explicit potential that expels electrons from its interior and breaks translational symmetry, which was one of our assumptions on the derivation of the dipole moment. Studying in detail the dipole moment build up in such cavities is an interesting open question that we hope future work will address.

VII Quantum Hall nematics

In this section we consider multi-valley systems of type-II that develop nematicity but which carry no dipole moment. We focus on the case with two valleys whose anisotropic mass tensors are rotated by a $\pi/2$ angle relative to each other, as realized in AlAs AlAsrev ; Abanin . Our main objective is to study non-trivial charged excitations, but we briefly review ground state properties. At filling $\nu=1$ one finds that the energy is (up to a global constant) Abanin :

[TABLE]

where $n_{z}={\rm tr}(\sigma_{z}P)$, and $\alpha>0$. The system has Ising character in agreement with the theorem of full valley polarization. To further asses the validity of the Hartree-Fock approach we have performed DMRG simulations that are described in Appendix D. It is in fact easy to covince one-self that the Ising nematic state is an exact eigenstate of the Hamiltonian projected into the valley degenerate Landau level. We confirmed that it is the ground state via DMRG employing a torus with a square aspect ratio (Appendix D) for as many as $24$ electrons and for mass anisotropies as large as $m_{x}/m_{y}=16$.

VII.1 Valley skyrmions in quantum Hall nematics

Having established that the ground state is a nematic state, we now turn to study its charged excitations. In the SU(2) invariant limit, i.e. when the mass tensors are identical for both valleys, we expect the lowest energy excitation to be infinite sized skyrmions Sondhi . When the mass tensors are slightly different, we expect that the skyrmions will have a finite size dictated by the competition between the Ising anisotropy which wants to shrink them and the Coulomb energy which wants to expand them to smear the charge over large distances 444The stiffness cost is scale invariant and hence indifferent to the size of the skyrmion.. To be able to study accurately the properties of skyrmions we resort to DMRG 555The energetics can be captured by a non-linear sigma model Sondhi ; Moon , but using long-wavelength estimates for non-linear sigma model parameters can lead to a substantial underestimation the size of the skrymions. This is a known pitfall of this type of estimates for small skyrmions Abolfath .. The skyrmion quasi-electron can be obtained as the ground state of the Hamiltonian when $N=N_{\phi}+1$ 666Quasi-hole and quasi-electron excitations are related by particle hole symmetry.. Figure 4 shows that non-trivial skyrmions (those involving at least one spin flip) survive up to large mass ratio $m_{x}/m_{y}\approx 3.8$. AlAs has mass anisotropy of about $m_{x}/m_{y}\approx 5$ AlAsrev . The experiment of Ref. AlAsskyrm found a non-linear dependence of the charge gap on strain, which was interpreted as evidence for skyrmions. Our findings suggest that this is scenario is not unlikely since the critical mass ratio to observe skyrmions is ultimately dependent on the details of interactions and can easily change by effects beyond our model (e.g. finite well widths and Landau level mixing).

We have also computed the charge gaps for these quasiparticles in Fig. 4 . The charge gap of a system can be defined as the jump in chemical potentials. The chemical potential is the energy change for adding a particle. Thus the chemical potential to add a quasi-electron on top of the ground state is $\mu_{+}=E_{N=N_{\phi}+1}-E_{N=N_{\phi}}$, similarly the chemical potential for quasiholes is $\mu_{-}=E_{N=N_{\phi}}-E_{N=N_{\phi}-1}$. With this the charge gap is found to be:

[TABLE]

A detailed comparison of the charge gaps of the non-trivial skyrmions we have found with those expected within Hartree-Fock theory is presented in Appendix E. Finally, although we focused on the skyrmions on the nematic case, skyrmions could also be present in the ferroelectric states. These Skyrmions will carry also a dipole moment although its magnitude cannot be inferred from the simple formulae we have developed. The precise pattern of the skyrmion texture in this case might be very complex due to the interplay of breaking of inversion symmetry and the long range nature of Coulomb interactions. We hope future work addresses the nature of these interesting quasiparticles.

VIII Summary

In summary, we have shown that multivalley systems with anisotropic dispersions lead generically to ferroelectric or nematic quantum Hall states at odd integer Landau level fillings. The ferroelectric states arise when the parent Fermi surface of a single valley lacks an inversion center, such as the case of Bi$(111)$ Yazdani . We have shown that the resulting dipole moment has an intimate relation with the underlying fermi surface geometry of the parent metal. We also demonstrated the existence of non-trivial charged skyrmion excitations in the nematic states with an unexpectedly large stability to the Ising symmetry breaking terms, shedding light into the question of the presence of these excitations in AlAs AlAsskyrm .

Acknowledgements.

We would like to thank Vidya Madhavan, Benjamin Feldman, Malika Randeria, Ali Yazdani, Donna Sheng, Heun Mo Yoo, Joonho Jang, Ray Ashoori, Fengcheng Wu, and Allan H. MacDonald for stimulating discussions. We specially thank Pok Man Tam for discussions and for catching a mistake in the numerical prefactor of Eqs. (28) and (49) in an earlier unpublished version of the manuscript. While at MIT I.S. was supported by the Pappalardo Fellowship. Z.Z. is supported by the David and Lucile Packard foundation. L.F. is supported by the DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award No. DE-SC0010526.

Appendix A Landau levels and density form factors for anisotropic parabolic dispersions

Consider the problem of Landau levels (LL) with an anisotropic mass tensor:

[TABLE]

Where $p=\nabla/i-eA$, and $g$ is a tensor that can be diagonalized as $g=Q^{T}S^{2}Q$, where $Q\in SO(2)$ and $S$ is diagonal with positive eigenvalues and $\det S=1$. Explicitly $S={\rm diag}\{(m_{x}/m_{y})^{1/4},(m_{y}/m_{x})^{1/4}\}$ and $m^{*}=(m_{x}m_{y})^{1/2}$. We can define rescaled momenta along the principal axes of the tensor to be: $\pi_{a}=(SQ)_{ab}p_{b}$, these satisfy the same commutation relations as the original ones:

[TABLE]

Which allows to solve the LL problem by defining the LL raising operators:

[TABLE]

satisfying $[a,a^{\dagger}]=1$. The guiding center operators are intra-Landau level operators defined as:

[TABLE]

They satisfy:

[TABLE]

The single particle Hilbert space can be decomposed into a tensor product $|n\rangle\otimes|m\rangle$, where $|n\rangle$ are the Landau level indices on which ${a,a^{\dagger}}$ act and $|m\rangle$ are intra-Landau level indices on which $R_{a}$ acts. Now the projected interaction into the Landau level of interest is obtained by imagining we have two flavors of particles with different mass tensors ($g_{1}=Q_{1}^{T}S_{1}^{2}Q_{1}$ and $g_{2}=Q_{2}^{T}S_{2}^{2}Q_{2}$) which interact via a potential that depends only on interparticle distance:

[TABLE]

where $r_{1}$ and $r_{2}$ above are understood to be operators. Using Eq. (35) we decompose the position of each particle as:

[TABLE]

where we are using matrix notation for the two component levi-civita symbol $\epsilon$. We have a similar expression for $r_{2}$. Using this we get:

[TABLE]

In the above expression the terms containing the operator $\pi_{1}$ produce inter-Landau level mixing and we proceed by projecting them into the zeroth landau level which is defined as $a|0\rangle=0$. By using Eq. (34) in combination with the BCH formula one can show that:

[TABLE]

where $\mathbf{q}=q_{x}+iq_{y}$. Then one obtains that the projected Hamiltonian is:

[TABLE]

Notice that $q_{1,2}$ are linear functions of $q$ described in Eq. (39).

For higher Landau levels we need to modify the density form factors. Using Eq. (40) and the algebra of raising and lowering operators one can show the following identities:

[TABLE]

where $n$ is the Landau level of interest and $L_{n}$ are the Laguerre polynomials. Therefore the interaction projected to the $n$ Landau level is:

[TABLE]

Appendix B Landau levels and density form factors for anisotropic Dirac cones

The derivation is very similar to that of Galilean electrons. We start from an anisotropic Dirac Hamiltonian:

[TABLE]

where $g$ is a $3\times 2$ tensor to which we apply a singular value decomposition of the form $g=R^{T}SQ$, where $R\in SO(3)$ describes a rotation of Pauli matrices $\sigma_{a}$ in pseudospin space, $Q\in SO(2)$ describes the transformation of the principal axes in real space, and $S$ is a $3\times 2$ matrix whose upper $2\times 2$ block is diagonal and characterizes the anisotropy of velocities $S={\rm diag}\{(v_{x}/v_{y})^{1/2},(v_{y}/v_{x})^{1/2}\}$, and whose lower row has zero entries, and $v=(v_{x}v_{y})^{1/2}$. Using the same definition of ladder operators as in Eq. (34), the Hamiltonian and the spectrum are:

[TABLE]

where $s=\pm$ labels negative and positive energy Landau levels and $k$ is the guiding center intra-Landau level index. From these expressions we can obtain the form factors of the interaction projected to the Dirac Landau levels. One can show that the zero Landau level is identical to the Galilean case, and for the excited Landau levels one gets the same expressions discussed in the previous section but with the modification:

[TABLE]

Appendix C Landau levels and density form factors in tilted Dirac cones

To the Hamiltonian, $H$, appearing in Eqs. (44), (45) we add a perturbation describing the tilt of the Dirac cone along the x-direction:

[TABLE]

where $a$ is defined in Eq. (34). The (un-normalized) perturbed zero energy Landau level to first order in $\delta v_{x}$ is found to be:

[TABLE]

where $\tau=\delta v_{x}/(\sqrt{2}v_{y})$. This verifies that the zero Landau level does not carry dipole as expected since particle-hole times inversion is a symmetry. Generically the $n$-th Landau level will break inversion symmetry however. The valence and conduction first Landau levels have a perturbed form:

[TABLE]

As described in the main text we imagine now two Dirac valleys at opposite momenta and with tilts $\delta v_{x}$ of opposite sign and equal magnitude. Let us assume we are filling only one of the two first Landau levels. The form factor obtained from the above perturbative expression is:

[TABLE]

where $f_{nm}(q)\equiv e^{-l^{2}|q|^{2}/4}\langle n|e^{il^{2}q\cdot\pi^{\prime}}|m\rangle$. To leading order in $\tau$ we can write the form factor as:

[TABLE]

One valley will have a form factor corresponding to $\tau$ and the other corresponding to $-\tau$ since the tilts have opposite signs. We label them by $\alpha=\{+,-\}$. If one performs the calculation of the exchange energy with these form factors then one finds that the exchange integrals in Eq. (4):

[TABLE]

The exchange energy can then be shown to be:

[TABLE]

with:

[TABLE]

Clearly for a repulsive interaction with $v_{q}>0$ we get that $X_{\tau}$ is strictly positive. Therefore the ground state is the Ising nematic ferroelectric state: $n_{z}=\pm 1$.

Appendix D DMRG implementation

Numerically we consider $N$ electrons moving along the surface of torus with a magnetic field perpendicular to its surface. $L_{x}$ ($L_{y}$) represents the circumference of the torus along $x$ ($y$) direction and they satisfy the relation $L_{x}L_{y}=2\pi N_{\phi}$, where $N_{\phi}$ represents the number of orbitals. Here, we set the magnetic length $l_{B}\equiv\sqrt{\hbar c/eB}$ as the unit of length. We choose a square torus with aspect ratio 1, i.e., $L_{x}=L_{y}$. Periodic boundary condition requires $k_{y}\equiv{2\pi j}/L_{y}$ with $j=0,1,...,N_{\phi}-1$. The Coulomb interaction on finite size system has the form

[TABLE]

with $V(q)=2\pi{e^{2}}/{L_{x}L_{y}\varepsilon q}$, and the wave-vectors are chosen to guarantee the periodicity of the interaction $V\left(x+L_{x},y\right)=V\left(x,y+L_{y}\right)=V(x,y)$. The projected Coulomb interaction between particles $i,j$ into the Lowest Landau level of a valley with mass anisotropic dispersion has the form:

[TABLE]

where $q_{1}$ is a linear functions of $q$ controlled by the anisotropic mass tensors described in Eq. (39). The matrix elements take the explicit form,

[TABLE]

Here, the Kronecker delta with the prime means that the equation is defined modulo $N_{\phi}$. We also consider a uniform and positive background charge so that the Coulomb interaction at $q=0$ is absent Yoshioka84 .

For the case of two anisotropic valleys, the Coulomb interaction includes both the intra-component and inter-component interactions. After projecting onto the lowest Landau level, the effective Hamiltonian reads as:

[TABLE]

Here, the first two terms are intra-component Coulomb interaction, and the last term is inter-component Coulomb interaction, and $q_{1,2}$ are the linear functions of $q$ controlled by the anisotropic mass tensors of each valley in described in Eq. (39). A similar expression to Eq. (57) holds for the explicit matrix elements.

Simulations are realized by mapping the single particle orbitals into one dimensional lattice. The advantage of the density matrix renormalization group (DMRG) method, compared to exact diagonalization, is that it can achieve reliable results for larger system sizes. In this manuscript, we impose momentum-space DMRG method on torus, where each site in $x$ direction corresponds to different states labeled by the momentum in $y$ direction. The standard procedure is similar to real-space DMRG, but since the non-zero Coulomb interaction exists between the orbitals satisfying the momentum conservation $j_{1}+j_{2}=j_{3}+j_{4}$ mod $N_{\phi}$, only small number of momentum sectors in each DMRG block contributes to the ground state in the initial process. To ensure the convergence, one has to save more momentum sectors with keeping a small number of states in these sectors during the initial process, which may become important in the later DMRG sweep, while the remaining states are selected following the standard DMRG algorithm to reach the required truncation accuracy. For our calculation, we keep around $2000\sim 8000$ states and perform measurement during each DMRG sweep until the results converge.

We found that in the ground state the charge distributes uniformly into the two orbitals for isotropic case 777The method selects a single sector even in this case in which there are SU(2) related degenerate copies., i.e., $m_{x}/m_{y}=1$, while the ground state is realized in a sector in which all electrons polarize into a single orbital with the same flavor for anisotropic case, as shown in the Fig. 5 for $N=N_{\phi}=24$ system with $(m_{x}/m_{y})^{1/4}=1.5,2$.

To study charged excitations we search for ground states with one particle added to those described above. The behavior of the charged excitations is markedly distinct for large and small anisotropic ratio $m_{x}/m_{y}$. When the anisotropic ratio is large enough, as shown in Fig. 6 for $N=N_{\phi}+1=21$ system with $(m_{x}/m_{y})^{1/4}=2$, $N_{\phi}$ electrons will fully occupy the orbitals with the same flavor, while the extra one will stay at the orbital with another flavor. The position of this extra electron is determined by the total momentum $K/(2\pi/N_{\phi})$ targeted in DMRG simulation, but different momentum sectors have the same energy due to translational invariance. This is consistent with the behavior of the conventional quasiparticle predicted by Hartree-Fock analysis. However, for smaller $m_{x}/m_{y}$ ratio, as shown in Fig. 7 with $(m_{x}/m_{y})^{1/4}=1.2$, the total number of electrons in the majority flavor will be smaller than $N_{\phi}$, while he number of particles in the minority flavor will be larger than one, indicating that the quasiparticle has a nontrivial pseudo-spin texture. Here, the number of particles in each flavor depends on both anisotropic ratio $m_{x}/m_{y}$ and total number of particles $N$, and therefore extrapolations to $N\rightarrow\infty$ are required. Fig. 8 shows the finite size scaling of the charge gap in the quasiparticle regime. We calculate the charge gap $\Delta_{gap}$ for the systems with $N_{\phi}=12$ to $N_{\phi}=28$ and the anisotropic ratio $(m_{x}/m_{y})^{1/4}=1.6$ and $(m_{x}/m_{y})^{1/4}=1.8$. We fit the data as a quadratic function of $1/N_{\phi}$, and estimate the thermodynamic limit of the gap by extrapolating $1/N_{\phi}\rightarrow 0$.

Appendix E Comparison between Hartree-Fock and DMRG quasiparticle gap

The charge gap in Hartree-Fock theory can be shown to coincide with minus twice the energy per particle of the ground state, which, for the Coulomb interaction, can be found to be:

[TABLE]

Here $E^{0}_{HF}$ is the Hartree-Fock energy of the ground state with no quasiparticles and $K$ is the elliptic integral of the first kind. This formula corrects a typo in Ref. Abanin, . Figure 9 illustrates the Hartree-Fock charge gap for the conventional quasi-electron and quasi-hole pair computed within Hartree-Fock theory and compared to the results from DMRG. After extrapolation to infinite sizes the results show good agreement in the large mass anisotropy regime where quasi-particles are no longer expected to be skyrmions.

The physical reason for the decreasing gap at large anisotropies is that as the anisotropy increases the overlap between orbitals becomes smaller and their exchange energy gain is reduced. Because the energy gain per particle is reduced it is also less energetically costly to add or remove particles, making the gap a decreasing function of mass anisotropy.