Nonlinear Hall effect with time-reversal symmetry: Theory and material realizations

The appearance of a Hall conductance necessarily requires breaking of time-reversal symmetry, either by an external magnetic field or by the internal magnetization of a material. However, as a second response, Hall dissipationless transverse currents can appear even in time-reversal symmetric conditions in non-centrosymmetric materials. Moreover, this non-linear effect has a quantum origin: it is related to the geometric properties of the electronic wavefunctions and encoded in the dipole moment of the Berry curvature. Here we review the general theory underpinning this effect and discuss various material platforms where non-linear Hall transverse responses have been theoretically proposed and experimentally verified. On the theoretical front, the link between the non-linear Hall effect and the Berry curvature dipole is discussed using Boltzmann transport theory. On the material front, different platforms, including topological crystalline insulators, transition metal dichalcogenides, graphene, and Weyl semimetals are reviewed.


I. INTRODUCTION
component that corresponds to the Hall conductivity can be written as with αβ the two-dimensional Levi-Civita symbol. The constraints on σ H due to crystalline symmetries can be identified by first recalling that a generic point group symmetry is represented by an orthogonal matrix O. Since the current and the electric fields transform as vectors under a generic coordinate change, the conductivity tensor must transform according to O T σ O. This transformation rule implies that the antisymmetric Hall conductivity transforms as a pseudoscalar As a consequence of the equation above, already the presence of a single point-group symmetry with det{O} = −1, such as a mirror line, forces the Hall conductance to vanish 22 .
Let us now consider the current response at second-order in electric fields. This defines a non-linear conductivity tensor In strict analogy with the linear response regime, we can separate the non-linear tensor χ αβγ into components that contribute to the electrical power dissipation j α E α , and dissipationless Hall components, which, contrary to σ H , are not forced to vanish only by time-reversal invariance. This Hall conductance can be identified by antisymmetrizing the first index with either the second or the third -both choices are equivalent because the tensor is by construction symmetric in the last two indices. In the particular case of two-dimensional systems, there are thus two independent components of the non-linear dissipationless Hall conductance reading and transform as a pseudovector under a point group symmetry, i.e.
The presence of a mirror line in this case only forces the non-linear Hall pseudovector to be orthogonal to the mirror plane. Therefore, χ H can be non-vanishing. On the contrary, the presence of two or more mirror lines require χ H to correspond to the null vector. Note also that the different transformation properties of the linear (Eq. 1) and nonlinear (Eq. 2) conductivity tensors imply that in certain crystalline structures non-linear dissipationless Hall currents can exist in the complete absence of linear Hall currents, even if time-reversal is explicitly broken. As an example, in crystals with C v point group symmetry the linear Hall conductance σ H vanishes independent of the presence of time-reversal symmetry. On the contrary χ H will be only orthogonal to the unique mirror line of the crystal. Following similar arguments, it is possible to derive the symmetry constraints on the (non)linear conductivity tensor in a three-dimensional material. The three independent components of the linear Hall conductance given by σ H γ = γαβ σ αβ /2 transform as a pseudovector. Therefore, precisely as the non-linear Hall conductance of twodimensional systems, the Hall vector of a three-dimensional material must be normal to any mirror plane. The presence of two independent mirror planes therefore forces all components of σ H γ to vanish. The three-dimensional non-linear Hall conductivity has instead nine independent components that transform as a rank-two pseudotensor. This tensor can be decomposed into a symmetric and an antisymmetric part 14 . The antisymmetric part transforms as a vector, and is potentially non-vanishing in the ten polar point groups C n and C nv where n = 1, 2, 3, 4, 6. In this case, the polar axis determines the direction of the non-linear Hall vector. The existence of a finite symmetric part of the three-dimensional non-linear Hall conductance is instead strongly dependent on the presence of mirror symmetries. All non-centrosymmetric crystal point groups without left-handed symmetries can potentially have a non-vanishing nonlinear Hall conductance. The only crystal point groups with mirror symmetries that allow for a non-vanishing symmetric part of the non-linear Hall tensor are instead C 1v , C 2v and S 4 . As we will see in the next section, the non-vanishing non-linear Hall conductance can be linked to the BCD using a Boltzmann equation approach in the constant relaxation time approximation.
B. Non-linear Hall effect due to Berry curvature dipole We now employ the semiclassical Boltzmann transport framework in a simple single band model to relate the nonlinear Hall conductance introduced in the preceding section to the BCD. We start by recalling that in the absence of externally applied magnetic fields, the semiclassical equations of motion accounting for the anomalous velocity term 9,23,24 can be written asṙ where is the energy dispersion of the metal in question, E is the driving electric field, and Ω k is the Berry curvature defined as Ω ka = abc ∂ k b A kc with A kc ≡ −i u k | ∂ kc |u k the Berry connection associated to the Bloch waves |u k with crystalline momentum k. To proceed further, we use that the electronic distribution function f k,r,t satisfies the semiclassical Boltzmann equation ∂ t f r,k,t +k · ∇ k f r,k,t +ṙ · ∇ r f r,k,t = I coll {f } where the collision integral in the relaxation time approximation reads In the equation above τ indicates an intraband scattering time and f 0 k the equilibrium distribution function. We are interested in a stationary and homogeneous solution to the Boltzmann equation. Using Eqs. 3, 4 we have To proceed further, we expand the distribution function in a series as f k = f 0 where the term f n k is understood to vanish to the E n order. One then finds the following structure for the linear and non-linear terms in the distribution function We can now compute the current using the relation j = e d d k/(2π) dṙ f k . At linear order in the electric field we have The first term on the right hand side of the equation above is the usual semiclassical contribution to the conductance that is expressed in terms of the electronic group velocities. This term can provide a finite contribution to the transverse resistance in low-symmetric crystals. However, such contributions are clearly symmetric in the α ↔ β exchange and therefore cannot contribute to the Hall conductance. The second term on the right hand side corresponds instead to the well-known intrinsic contribution to the anomalous Hall conductance controlled by the Berry phases accumulated by the particle motion on the Fermi surface 8 . Because the Berry curvature Ω k is an odd function of the momentum when time-reversal is preserved, this Hall conductance vanishes in time-reversal symmetric conditions, in perfect agreement with Onsager relations 21 . Let us now instead consider the non-linear response of the current ∝ E 2 . It is simple to show that the non-linear current takes the following form The first term, which is of entirely semiclassical origin, vanishes in time-reversal symmetric condition since it involves the three index tensor ∂ kα (k)∂ k β ∂ kγ f 0 k that is odd under time-reversal. After integration by parts, the non-vanishing non-linear conductivity tensor can be therefore written as We therefore have that in time-reversal symmetric conditions the dissipationless Hall non-linear current is regulated by the first moment of the Berry curvature, the BCD 14 , over the occupied states The BCD is subject to the same symmetric constraints of the non-linear Hall conductance introduced in the preceding section. In two-dimensional systems, for instance, the Berry curvature is a pseudoscalar and therefore the BCD is a pseudovector precisely as χ Hγ . In our foregoing discussion, we have considered the dc limit. When accounting for an ac driving electric field, nonlinear Hall currents yield a current response at twice the driving frequency, with the addition of a rectified current. Specifically for a driving electric field E α (t) = Re E α e iωt , with E ∈ C, the resulting current at twice the frequency j 2ω α = χ αβγ E β E γ while the rectified current j 0 α = χ αβγ E β E γ . The ac non-linear response function takes the following form It is interesting to note that at frequencies above the width of the Drude peak ωτ 1 but below the interband transition threshold, the non-linear response function becomes independent of the scattering time and therefore provides a direct measure of a geometric property of the electronic wavefunctions.

C. Disorder-induced contributions to the non-linear Hall effect
The BCD does not completely determine the NLHE in a quantum material. Using either a quantum Boltzmann transport approach 19 or a semiclassical Boltzmann theory beyond the constant relaxation time approximation [25][26][27] , it can be shown that other disorder-mediated contributions to the non-linear Hall conductance exist. These "extrinsic" contributions are the non-linear counterparts of the side-jump and skew-scattering contributions appearing in the linear AHE 9 . In this section, we review the Boltzmann semiclassical transport framework by following Ref. 27 , to which we refer the readers interested in further details.
Let us start out by rewriting the collision integral appearing in the kinetic equation Eq. 4 as where the label l is a composition a band and momenta indices (n, k) whereas l = n d d k/(2π) d . In addition, ω ll is the disorder averaged scattering rate between the Bloch waves with quantum numbers l and l . The scattering rate ω ll can be related to the T-matrix element The scattering T matrix is defined by T l l = l|V imp |ψ l withV imp indicating the impurity potential operator whereas |ψ l represents the eigenstate of the full Hamiltonian H = H 0 +V imp that satisfies the Lippman-Schwinger equation For weak disorder one can approximate the scattering state |ψ l by a truncated series in powers of V ll = l|V imp |l as Inserting this expression in the expression of the T-matrix element of the disorder potential, we can therefore expand the scattering rate in powers of the disorder strength as The scattering rate ω can be neglected since they only renormalize the second-order scattering rate ω (2) ll . On the contrary, the antisymmetric contributions to ω (3,4) ll yield the non-linear skew-scattering and side-jump contributions to the NLHE in time-reversal symmetric conditions. To show this, we first go back to the Boltzmann equation and rewrite explicitly the collision integral using the symmetric and antisymmetric contributions as follows: where, now, ω (a) l l contains the antisymmetric contributions of both ω (3) ll , and ω (4) ll . The equation above, however, does not account for the microscopic displacement experienced by a wavepacket when scattering from a generic state l to l , δr ll . Assuming smooth impurity potentials, the gauge-invariant expression for this coordinate shift [28][29][30] , usually referred to as side-jump, reads as where arg is the phase of the complex number and we introduced the operatorD k k = ∂ k + ∂ k . In the presence of the external driving electric field, the microscopic displacement δr ll yields an energy shift ∆U ll = −eE · δr ll , the scattering rate should account for. The symmetric second-order scattering rate, specifically, has to be therefore modified as where we have used the antisymmetric property [c.f. Eq. 11] of the side-jump δr ll = −δr l l . The Boltzmann equation Eq. 10 is then modified as where we introduced the quantity Having in our hands the anomalous and the skew-scattering linear and non-linear distribution functions, we can now evaluate the response current to a driving electric field j = e lṙ f l . In doing so, we note that the side-jump velocity introduced above explicitly enters into the semiclassical equation of motion 9,24,27-31 This implies that in the linear response regime there are three disorder-mediated contributions beyond the semiclassical and anomalous Hall conductance of Eq. 6. There is a first contribution due to the anomalous side-jump distribution that yields a conductivity In addition, the side-jump velocity appearing in the equation of motion contributes with an additional conductivity We point out that we have neglected terms ∝ v sj l g (adis,sk),1 since they are of higher-order power in the scattering rate ω. Finally, the skew-scattering contribution to the distribution function leads to the conductivity In systems with broken time-reversal symmetry, Eqs. 18,19,20 provide the side-jump and skew-scattering contributions to the anomalous Hall conductance. We note that Eq. 20 can be split into two different contributions corresponding to the antisymmetric scattering rates ω (3) and ω (4)9 . In time-reversal symmetric conditions these extrinsic contributions to the Hall conductance are forced to vanish, precisely as the intrinsic Berry-curvature mediated anomalous Hall conductance. This is because the coordinate shift δr l l , and consequently the side-jump velocity v sj l , are even under time-reversal. Since the conventional group velocity ∝ ∇ k is instead odd under time-reversal, one finds that sidejump conductivities σ sj, (1,2) are forced to vanish. The same holds true for the skew-scattering contribution since the antisymmetric scattering rate is odd under time-reversal.
We now evaluate the non-linear disorder-induced conductivity. Using the aforementioned symmetry constraints, it is easy to show that the terms quadratic in the driving electric field where the anomalous velocity Ω is directly coupled to the anomalous distribution g adis,1 l and the skew-scattering distribution g sk,1 l are forced to vanish by time-reversal. As a result, we have that the only Berry curvature mediated term corresponds to the dipole of Eq. 8. As for the linear conductivity, there are three disorder-induced non-linear terms. The first contribution comes from the second-order anomalous distribution function g adis,2 l . From Eq. 15 and after simple manipulations it can be recast as where we introduced the tensor The contribution due to the side-jump velocity can be instead put as Finally, the skew scattering contribution to the non-linear distribution function yields the additional conductivity Eqs. 21,22,23,24 provide the disorder-induced contributions to the non-linear Hall effect in the dc limit and assuming a constant relaxation time. As mentioned above, these results can be easily generalized assuming a driving a.c. electric field and a general relaxation time (see Ref. 27 for further details). In closing this section, we point out that disorder-induced contributions to the non-linear Hall conductance have a dependence in the impurity concentration n i given by This also implies that precisely as the BCD contribution Eq. 8, the disorder-mediated corrections to the non-linear Hall conductance grow linear with the relaxation time. Nevertheless, the different contributions to the non-linear Hall conductance can be distinguished in the experimental realm by using the scaling between the non-linear Hall signal and the conventional longitudinal (linear) resistivity, in strict analogy with the scaling used in the AHE 32,33 .

A. Tilted massive Dirac cones
Beside the general symmetry constraints discussed in Sec. II A, the Berry curvature-mediated contribution to the non-linear Hall conductance is subject to other point-group symmetry restrictions, which are of primary importance in identifying materials that possess substantial BCD. Let us first consider specifically the role played by rotational symmetries. As long as time-reversal symmetry is preserved, all systems with an evenfold rotation symmetry C n (with n = 2, 4, 6) cannot have a finite BCD. This is because the composed symmetry C 2 Θ, Θ being the time-reversal operator, is an antiunitary symmetry that squares to one. Moreover this symmetry acts locally in momentum space since it brings the two-dimensional momentum k back to itself. These two properties imply that the Berry curvature is forced to vanish for all momenta in the two-dimensional BZ 34 . Consequently, a non-vanishing BCD can only appear in two-dimensional crystals where there is either a threefold rotation symmetry or all rotation symmetries are broken. For the latter and as discussed in Sec. II A, a single mirror line can still exists, in which case the BCD is forced to be orthogonal to the mirror line. In the former case, instead, the BCD will not be pinned to any specific direction. The time-reversal symmetric NLHE in a C 3 -symmetric setting has been not proposed or experimentally realized so far. Therefore, in the remainder we will limit ourselves to discuss systems with a single mirror line symmetry, i.e. with point group symmetry C s . In order to have a sizable Berry curvature, the low-energy electronic properties of such systems must be described by massive Dirac cones. Furthermore, the absence of rotational symmetries implies that there is no fermion multiplication theorem 35 , and the minimum number of massive Dirac cones allowed by symmetry is two 36 . To preserve the mirror symmetry these two massive Dirac cones will be located specularly with respect to the mirror-symmetric line of the two-dimensional BZ. To make things concrete, consider for instance the reflection symmetry to map a point with coordinates (x, y) to (x, −y). The massive Dirac cones will then be centered around two valleys Λ 1,2 = k x , ±k y . The effective Hamiltonian close to these valleys can be derived by accounting for all symmetry-allowed terms in a k·p expansion. Time-reversal symmetry gives the constraint on the effective Hamiltonian Here Θ is the time-reversal operator that is represented as Θ = iσ y τ x K where the Pauli matrix vectors σ and τ act in spin and valley space respectively, whereas K is the complex conjugation. Similarly, the mirror symmetry gives the constraint M −1 The mirror symmetry operator takes the form M y = −iσ y τ x since it exchanges the valleys, as time-reversal does, and acts in spin space as exp{−iπσ y /2}. Retaining terms up to linear order in momentum k and neglecting intervalley mixing terms -this condition amounts to keep only symmetry allowed terms ∝ τ 0 , τ z -the low-energy continuum Hamiltonian gets the following form: In the equation above we have neglected a term of the form ασ x τ z since this term shifts in a time-reversal symmetric manner the two valleys Λ 1,2 , and can be thus reabsorbed in the effective Hamiltonian by a proper redefinition of the momentum k y . Eq. 25 corresponds to the low-energy theory of a massive Dirac cone except for the term proportional to a, which produces a tilt as shown in Fig. 1(a). The presence of this tilt term is key for the occurrence of a finite BCD 14,19,27 , as we now discuss. It is straightforward to show that the Berry curvature is independent of the tilt strength and can be written as where the ± sign refers to the conduction and valence bands respectively, and τ = ±1 for the two valleys Λ 1,2 . In each valley the Berry curvature is an even function of the momentum k measured relatively to the Dirac points. Therefore, the BCD density ∂ ky Ω τ kz is an odd function of the momentum. Consequently, the contribution to the BCD coming from each valley is identically zero in the absence of a tilt term. The situation is different for a = 0. In this case the Fermi lines become elliptical, and therefore the contribution to the BCD in each valley is different from zero [c.f. Fig. 1(b-e)]. In addition, the contribution from the two different valleys is same since the elliptical distortion of the Fermi lines is opposite in the two valleys. As a final result, we therefore have that the BCD density is directly proportional to the tilt parameter a. Another way to understand the tilt-driven occurrent of a finite BCD is to transform D y as a line integral over the Fermi line of the system as In the absence of a tilt term, this integral vanishes in each valley since the Berry curvature is constant on the Fermi lines whereas the group velocity in theŷ direction is equal but opposite on the opposite sides of the Fermi surface. Tilting the Dirac cones allows for a mismatch between left and right movers in each valley [see Fig. 2(a),(b)] leading to a finite contribution to the BCD. Furthermore, the contribution to the dipole coming from the two valleys is equal since the opposite mismatch between left and right movers is compensated by the opposite values of the Berry curvature. Using Eq. 27 it is possible to provide an analytical expression for the BCD in systems with tilted massive Dirac cones. Consider for simplicity the case in which v x = v y = v, and assume the tilt parameter a is small as compared to the Fermi velocity v. Expanding the Dirac delta as δ ( − F ) = δ 0 − F + aτ k y ∂ 0 δ 0 − F , with 0 indicating the energy of the untilted massive Dirac cone, the BCD can be found to be where we considered the Fermi energy to be in the conduction band. It is important to note that as a function of the chemical potential, the BCD displays a characteristic non-monotonous behavior 14 [see Fig. 2(c)]. A similar feature is also found when explicitly computing the disordered-averaged contributions to the NLHE discussed in the previous section 19,27 . The overall contribution does not cancel out and is still finite.

B. Surface states of topological crystalline insulators
The surface states of the topological crystalline insulator SnTe 37-39 represent a prime example of tilted massive Dirac cones. This makes SnTe a paradigmatic material platform for the realization of the quantum NLHE in timereversal symmetric conditions. Generally speaking, the group-IV tellurides have a high-temperature rocksalt crystal structure with a face-centered cubic Brillouin zone (BZ) [c.f. Fig. 3]. The BZ is bounded by six square faces and eight hexagonal faces 37 . The centers of the latter, commonly denoted by L, represent the equivalent high-symmetry points of the BZ where the fundamental gap of the SnTe or PbTe is located. The band structure in the immediate vicinity of the four equivalent L points can be captured using a four-band k · p model 40 that span spin space and the p-orbital space of the anion (Te) and cation (Pb,Sn). In PbTe, the valence band is derived from the p-orbital of the anion while the conduction band from the cation Pb. On the contrary, SnTe has an inverted band ordering with the valence band derived from the cation Sn and the conduction band from Te. The inverted ordering of the bands of SnTe as compared to PbTe establishes the former as a topological crystalline insulator. The bulk crystalline topology of the material comes about the non-trivial mirror Chern number associated to the ΓL 1 L 2 plane of the BZ [see Fig. 3], which is left invariant by the (110) mirror symmetry of the rocksalt crystal structure. The non-trivial value of the mirror Chern number then implies the presence of counterpropagating midgap modes on the surface projections of the mirror invariant plane 37,41 . This, for instance, applies to theΓ −X 1 line of the (001) surface BZ. Since the mirror Chern number associated to the ΓL 1 L 2 plane n M = 2, there will be two pairs of counterpropagating edge modes along theΓ −X 1 line. In addition, since the left and right edge modes belong to the different ±i mirror sectors, their crossings cannot be gapped. On the contrary, away from theΓ −X 1 line the in-gap edge modes can be gapped since they are not protected by mirror symmetry. The crossings on theΓ −X 1 line of the (001) surface BZ are the Dirac points of two topologically protected surface Dirac cones. By rotational symmetry, such surface Dirac cones also appear along theΓ −X 2 line. One can derive the electronic characteristics of these four surface Dirac cones by using a two-dimensional k · p theory close to theX 1,2 points of the surface BZ 42,43 . It can be found by simply noticing that theX 1,2 points of the surface BZ are invariant under a twofold rotation symmetry C 2 , and the two mirror line symmetries M x,y -from here onwards we denote withx the (110) direction perpendicular to the mirror plane. As in our former analysis, these symmetry operations, together with time-reversal symmetry, constrain the allowed terms in the k · p expansion for the surface Dirac cones. Furthermore, the representation of the four symmetries can be derived by noticing that at theX 1,2 points of the (001) surface BZ two out of the four equivalent L points of the BZ are projected [see Fig. 3]. Imagine now that instead of considering an atomically sharp interface between the vacuum and the material bulk, one would consider a smooth interface separating SnTe from PbTe. Then, one would expect Dirac-like domain wall states, each of which coming about from the change in the mass gap in one specific L point of the BZ. At theX 1,2 one would then expect two different flavours of Dirac-like states. Consequently, the effective surface k · p theory should account for such a flavor degree of freedom beside the internal spin degree of freedom. This holds true even when the smooth interface is substituted by an atomically sharp interface. With this in mind, we can proceed to find the representations of the symmetry operations that constraint the surface low-energy theory. Let us explicitly consider the low-energy theory for the surface states projected onto theX 1 point of the surface BZ. The theory close to theX 2 point can be obtained following the same strategy.
The internal time-reversal symmetry does not act on the flavor degree of freedom. Therefore it can be represented, as usual, by Θ = iσ y K where K is the complex conjugation and the Pauli vector σ acts in spin space. The twofold rotation operator C 2 interchanges the L 1,2 valleys while it can be represented in spin space by exp{−iπσ z /2}. As a result, the twofold rotation operator is represented by C 2 = −i τ x σ z . Similarly to the twofold rotation operator C 2 , also the mirror symmetry M y interchanges the L 1 point with L 2 . Hence, this mirror symmetry can be represented as M y = −iσ y τ x . On the contrary, the mirror symmetry M x = −iσ x τ 0 as it does not interchanges the two projected L points. With this, we can write the surface k · p theory away from theX 1 point of the surface BZ. Specifically, by retaining flavor mixing terms up to zeroth order in momentum the low-energy theory can be found to be As expected, the surface states are gapped specifically on the k y = 0 line, i.e. along theX 1 −M direction of the surface BZ. On the contrary the M x symmetry guarantees the presence of two zero-energy Dirac points for k y = ±v y / √ m 2 + δ 2 on the k x = 0 line. In perfect agreement with the foregoing general symmetry analysis, the in-gap surface states realize two Dirac cones located at Λ 1,2 with an effective Hamiltonian where we introduced the momentum δk = k − Λ 1,2 and the renormalized velocity v x = v x δ/ √ m 2 + δ 2 . As mentioned above, this form of the surface states for SnTe has been found by neglecting flavour mixing term ∝ τ x,y,z linear in momentum. If the latter were to be considered, the following perturbation should be added to the original surface k · p Hamiltonian: Projecting this terms in the low-energy Hamiltonian reveals that the surface Dirac dispersion of Eq. 30 is modified in two different ways. First, it leads to an additional renormalization of the Fermi velocity v x . Second, it yields a tilt of the surface Dirac cones. The corresponding low-energy theory indeed reads where the tilt parameter a is directly proportional to the k · p parameters α y , β whereas χ = ±1 distinguishes the Λ 1,2 valleys. The absence of a surface energy gap does not allow for a finite Berry curvature. The latter, however, naturally arises 44,45 when considering a ferroelectric distortion whereby the Sn and Te atoms are displaced along different directions. SnTe is known to undergo a structural transition at low temperatures 46-48 that involves precisely such a ferroelectric distortion. Besides yielding a finite ferroelectric polarization, this structural distorsion breaks completely the rotational symmetry 49 . Moreover if the displacement vector is along thex direction, the structural distortion breaks the M x mirror symmetry, whereas a displacement in the orthogonal direction leads to a loss of the M y reflection symmetry. Using the preceding symmetry analysis, it can be shown that breaking the M x symmetry changes the crossings of the states along the Γ −X 1 line into avoided level crossings, and ultimately the Dirac cones centered at Λ 1,2 acquire a mass gap. Breaking instead the M x symmetry preserves massless Dirac fermions close to theX 1 point of the surface BZ 43 . The situation is of course reversed when considering the surface states located close to thē X 2 point of the surface BZ. All in all, we therefore have that a ferroelectric distortion with the displacement vector oriented along a principal crystallographic direction of the surface BZ will result in the gapping of two out of the four total (001) surface Dirac cones in SnTe.
To concretely show how a mass gap for the surface Dirac cones at Λ 1,2 appears when the M x reflection symmetry is broken, it is sufficient to notice that in this situation an additional symmetry-allowed term has to be included in the k · p Hamiltonian of Eq. 29. It reads Projecting this term onto the low-energy surface state Hamiltonian, we find that the ferroelectric distortion yields two opposite masses for the two surface Dirac fermions centered at Λ 1,2 . Therefore, the full Hamiltonian can be recast as where χ = ±1 distinguishes the Λ 1,2 valleys. We have therefore reached the massive tilted Dirac fermion model introduced in the preceding section and characterized by a non vanishing BCD.

C. Transition metal dichalcogenide monolayers
Having established the occurrence of a finite BCD from the topologically protected surface states of SnTe, we next show the occurrence of this phenomenon also in two-dimensional crystalline materials characterized by a strong spinorbit coupling. Due to the presence of massive Dirac cones, transition metal dichalcogenides (TMDs) in their so-called 1H form [see Fig. 4(a)] have been proposed as promising platforms for the occurrence of sizable BCD 50 . Because of the threefold rotation symmetry of the crystalline structure, 1H-TMDs do not satisfy the symmetry constraints discussed above. Consequently an external perturbation, namely uniaxial strain, needs to be applied.
Generally speaking, 1H monolayers of group-VI dichalcogenides MX 2 (with M=Mo,W and X=S, Se, Te) have a direct band-gap in the visible frequency range 51,52 with conduction and valence bands edges that are located at the corners of the two-dimensional hexagonal BZ. As a result, the low-energy properties of these systems are very similar to graphene but with two important differences. First, inversion symmetry is explicitly broken thus implying a non-vanishing Berry curvature 53,54 . Second, TMDs have a strong spin-orbit coupling resulting from the d-orbitals of the heavy metal atom. The appearance of massive Dirac cones can be immediately understood by first neglecting spin-orbit coupling all together. Density functional theory calculations [55][56][57][58] show that the electronic band structure consists of partially filled d-bands of the heavy metal ion lying between M-X s − p bonding and antibonding bands. The trigonal prismatic coordination of the transition metal ions [see Fig. 4(a)] splits its d-bands into a single one-dimensional representation -the d z 2 orbital -and two different two-dimensional representations that differ from each other in the ±1 eigenvalue of the reflection symmetry with respect to the (horizontal) mirror plane. Specifically, the d x 2 −y 2 and d xy orbitals form the so-called irreducible E representation, whereas the d xz and the d yz orbitals give rise to the so-called E two-dimensional irreducible representation 59 . At the Γ point of the BZ the valence and conduction states close to the Fermi level are given by the single d z 2 orbital and twofold degenerate states composed of the d xy and d x 2 −y 2 orbitals. These two degenerate states split along the Γ − K line thus allowing to describe the low-energy properties in terms of an effective two-band model. The states forming the effective low-energy doublet at the K point are given by the combined orbitals d z 2 and d xy ± id x 2 −y 2 as can be simply rationalized by noticing that the little group at the K point is C 3h . We emphasize that it is the absence of a vertical mirror symmetry σ v in the little group that guarantees the splitting of the E doublet at the K point. On the other hand, time-reversal symmetry guarantees that the states at the K and K points, which are related to each other by time-reversal symmetry, are given by d xy + iτ d x 2 −y 2 with τ the valley index. The effective Hamiltonian away from the K, K of the BZ can now be derived using two-dimensional k · p theory. To do so, we notice that the threefold rotation symmetry operator can be represented as C 3 = exp{−iπσ z τ z /3} whereas the time-reversal operator Θ = τ x K. Finally, the mirror symmetry interchanging the valleys can be simply represented as M = τ x . The presence of these symmetries allows the presence of a Dirac mass ∆ σ z /2 where ∆ corresponds to the TMD direct gap, while the triad of Pauli matrices σ x,y,z now acts in the orbital space. Using the transformation of momenta under the threefold rotation symmetry one finds to linear order in momentum a isotropic Dirac theory of the form The strong spin-orbit coupling of the metal d orbitals require the inclusion of the atomic spin-orbit coupling L · S. The ensuing effective Hamiltonian 53 can be then recast in the form with s z indicating the spin degrees of freedom, and λ the halved spin splitting of the valence band edge caused by spin-orbit coupling. Note that the eigenstates of the effective k · p Hamiltonian are common eigenstates of the spin in the direction perpendicular to the two-dimensional TMD sheet due to the presence of the horizontal mirror symmetry M z . The low-energy Hamiltonian in Eq. 36 corresponds to two copies (but with different gaps) of the effective Hamiltonian for the surface states of topological crystalline insulators without the tilt term. The absence of the latter is enforced by the threefold rotation symmetry, which, in turn, implies the existence of three vertical mirror symmetries, and forces the BCD to vanish. However, the application of uniaxial strain can lower the point group symmetry making the tilt term symmetry allowed and the BCD finite. The presence of a strain-induced tilt term cannot be captured in an effective k · p theory considering strain terms at lowest order. However, a microscopic model based on ab initio derived Wannier functions with strain effects incorporated by varying the interatomic bond length 60 shows the appearance of a finite BCD in complete agreement with the symmetry analysis. The size of the BCD is, however, much smaller than the one theoretically predicted to occur on the surface states of SnTe. The latter has been estimated to be in the nm range 14 , whereas BCDs of the order of 10 −2Å are expected in the 1H structure of both WSe 2 50 [see Fig. 4(b)] and MoS 2 61 . Although the hexagonal structure is the most energetically favorable one, TMD monolayers are present also in a different structural form: the so-called 1T d phase. It corresponds to the bulk T d phase in which WTe 2 and MoTe 2 realize Weyl semimetals [63][64][65][66] . The monolayer 1T d phase intrinsically breaks bulk inversion symmetry. However, at sufficiently low-temperature the monolayers undergo a slight structural distortion thanks to which inversion symmetry is recovered. This structure is the so-called 1T structure, which, as we now review, possess gapless tilted Dirac cones if spin-orbit coupling is removed all together. To show this, let us consider the relevant crystalline symmetries 62 of the 1T structure [see Fig. 5(a)]: a modified mirror symmetryM x corresponding to a conventional mirror line symmetry M x followed by a translation by half of a lattice vector t (e x /2), and a screw rotationC 2x that is the product of a twofold rotation C 2x around thex axis followed by the same fractional translation t (e x /2). Note that the combined presence of these two symmetries implies the presence of the inversion symmetry I =C 2x ×M x . In the two-dimensional BZ of the system, there are two high-symmetry lines, i.e. k y = 0 and k y = π where the bands can be labelled by the eigenvalues ofC 2x . The latter fall into two momentum-dependent branches ±e −ikx/2 as follows from the fact that C . This topological semimetal protected by screw rotation becomes a quantum spin-Hall insulating state when the intrinsic spin-orbit coupling is taken into account [67][68][69][70][71] . In this case, the Berry curvature is identically zero due to the concomitant presence of inversion and time-reversal symmetry. However, breaking inversion symmetry via either an externally applied electric field or considering a structural distortion to the non-centrosymmetric 1T d phase, all the conditions for the appearance of a BCD are met. Both these situations have been theoretically studied using density functional theory calculations 50,72 . In particular, it has been shown that WTe 2 in the 1T d possess BCDs of the order of 10 −1Å . Larger values have been instead found using the tunability provided by the external electric field in the 1T phase both in WTe 2 and MoTe 2 monolayers. Even more importantly, a BCD of the order of 10 0Å has been estimated using circular photogalvanic measurements in WTe 2 73 . Similarly to the density functional theory calculation studies, these experiments have highlighted a large tunability of the BCD.

D. Bilayer transition metal dichalcogenides
Sizable BCDs have been theoretically predicted 74 and experimentally verified in bilayer 75 and few-layer 76 WTe 2 . TMD bilayers, in particular, display dipoles of the order of a few nanometers, one order of magnitude larger than the values predicted and experimentally observed in the electrically activated 1T phase of TMD monolayers, and four order magnitudes than the one predicted in 1H monolayers. Moreover the size of the dipole can be additionally tuned with the application of an out-of-plane electric field. The origin of such a large BCD in TMD bilayers can be understood by considering a simple mechanics of hybridization between Dirac fermions appearing in two isolated monolayers 62 . Assume for simplicity that on the screw lines of the TMD monolayers the Dirac cones have a vanishingly small tilt term. In the unrealistic case in which the two layers are completely decoupled the bilayer will then feature two massless Dirac cones one on top of each other [see Fig. 6(a)]. Consequently, the electronic band structure of the bilayer system will be twofold degenerate at all wavevectors. Next, we can assume that the two layers are coupled but stacked without any relative displacement in the x − y plane. The interlayer coupling will energetically split the Dirac cones, and thus remove the twofold degeneracy of the bands. However, the Dirac points on the screw lines will still be present [see Fig. 6(b)] since the screw rotation is preserved. Finally, one can consider the stacking experimentally realized in, e.g. WTe 2 . The relative displacement between the two layers breaks the screw rotation symmetry. Therefore, the Dirac cones acquire a mass and Berry curvature [see Fig. 6(c)]. Moreover, since the original Dirac cones have opposite chirality the Berry curvature of the resulting massive Dirac cones are opposite. After they become gapped, these untilted Dirac cones yield a strong BCD since the right moving states and the left moving states are characterized by an opposite Berry curvature. Adding now the effect of the spin-orbit coupling and the original tilt of the Dirac cones does not change qualitatively any of the conclusion.

IV. BERRY CURVATURE DIPOLE IN THE ABSENCE OF SPIN-ORBIT COUPLING
The mechanism at work in bilayer TMDs does not strictly require the presence of titled massive Dirac cones, neither of a sizable spin-orbit coupling. This feature highlights the possibility to observe the non-linear Hall effect even in two-dimensional materials made of light elements. Graphene has been consequently put forward as a paradigmatic example of spin-orbit-free materials a non-vanishing BCD 77 . Importantly, signatures of a non-linear Hall effect with time-reversal symmetry have been recently found in strained bilayer graphene 78 . As we will discuss below, the origin of the non-linear Hall effect in graphene does not stem from the Dirac cone shifting mechanism of bilayer WTe 2 . It is instead the warping of the Fermi surface the physical characteristic, which is responsible for the non-vanishing of the BCD.

A. Uniaxially strained monolayer graphene
The wallpaper group of monolayer graphene is p6mm. It is generated by the point group C 6v and in-plane translations. The point group is generated by a threefold rotation symmetry, a mirror symmetry and a twofold rotation. Since spin-orbit coupling can be neglected all together in graphene, the twofold rotation coincides with inversion symmetry. As a result, the Berry curvature identically vanishes. A possible way to open up a gap and break inversion symmetry is to induce a sublattice imbalance 79 , i.e. a charge density wave instability. This can be achieved for instance by placing graphene on a lattice-matched substrate. Density functional theory calculations 80 first identified hexagonal boron nitride as a possible substrate inducing a band gap of the order of 25 meV. The small nominal lattice mismatch between graphene and hexagonal boron nitride can also lead to a large moire' superlattice characterized by the generation of mini Dirac cones at the expense of a full band gap opening [81][82][83][84] . Signatures of a commensurateincommensurate transition for graphene on top of hexagonal boron nitride 85 have revealed that both these situations can occur in practice: moire' gapless regions are separated by the gapped commensurate ones of interest for our discussion below. As for the case of TMD monolayers in the 1H phase, the point group C 3v of gapped graphene does not allow for a non-vanishing BCD. However, the application of uniaxial strain lowers the symmetry of the system to C v ; the presence of a single mirror line then allows for the onset of a non-linear Hall effect. Let us now examine the electronic characteristic of a graphene layer in the presence of uniaxial strain. To do so, we start out from the graphene tight-binding Hamiltonian 86 in its simplest form, i.e. considering only hopping processes between nearest neighbor atomic sites. It reads: where a † and b † (a, b) are creation (annihilation) operators on the A and B sublattices respectively. In the equation above, ∆ indicates the substrate-induced band gap, the subscript i runs over all unit cell positions, and we introduced the three nearest neighbor vectors Furthermore, the presence of strain implies that the hopping amplitudes t n explicitly depend on the nearest neighbor vectors. Specifically, t n = t 0 (1 − β δu n ) where t 0 is the hopping amplitude for the pristine threefold rotation symmetric honeycomb lattice, the lattice parameter β can be determined by Raman spectroscopy whereas the relative distance changes δu n can be expressed in terms of the strain tensor components i j as To proceed further, we go to momentum space and write the Bloch Hamiltonian as where we have rewritten the momenta as k = K ( ) + q since we are interested in the electronic properties close to the K or K valleys of the BZ given by K = 4π 3 √ 3a , 0 and K = − 4π 3 √ 3a , 0 . The Bloch Hamiltonian can be then FIG. 7. Berry curvature and dipole density for unstrained (left panels) and strained (right panels) monolayer graphene. When strain is present the deformation of the Fermi surface leads to a finite value of the dipole. Note that we have considered for simplicity the trigonal warping term that respects the threefold rotation symmetry, i.e. λ1 = λ2 = λ3. The distortion of the Fermi surface is therefore entirely due to the anisotropic Fermi velocity.
expanded to linear order in the small momenta q. Using simple vector identities and assuming an anisotropic biaxial strain with xx = yy = 0 and xy ≡ 0 the continuum low-energy Hamiltonian near the two valleys of the BZ can be then recast as where ξ = ±1 is the valley index. In the equation above, A x = √ 3 β( xx − yy )/(2a) is the well-known straininduced "pseudo"-gauge field [87][88][89] whereas v F = √ 3t 0 a/2 is the Fermi velocity of the Dirac carriers in unstrained samples. In addition, v x = v F [1 − β(3 xx + yy )/4] and v y = v F [1 − β( xx + 3 yy )/4] are renormalized Fermi velocities that become anisotropic when explicitly considering the momentum-strain coupling 90 . The Dirac cones explicitly acquire an anisotropic character, consistent with the reduction of the point group symmetry. The warping of the Fermi surface can be instead captured by keeping terms up to quadratic order in the momentum q 2 . They read ∆H ef f (q) = (λ 1 q 2 y − λ 2 q 2 x )σ x + 2ξq x q y λ 3 σ y , where the warping coefficients λ 1,2,3 are explicitly renormalized by strain. In the absence of it, λ 1 ≡ λ 2 ≡ λ 3 and we reach the trigonal warping of pristine graphene 86 . By limiting to homogeneous strain, the presence of the pseudo-gauge field can be simply reabsorbed by a proper redefinition of the small momentum: The pseudo-gauge field in fact only shifts the Dirac cones away from the K ( ) valley. The anisotropy in the Fermi velocity combined with the warping of the Fermi surface instead produces a distortion of the Fermi lines [see Fig. 7 ] that endows the system with a non-zero BCD even if an explicit tilt term is absent 77 . Precisely as for the case of TMDs both the gapped Dirac cones contribute in an equal manner to the total BCD 14 . Note also that a finite value of the BCD occurs even if the strain-induced renormalization of the warping coefficients is disregarded [c.f. Fig. 7] . This finite BCD is pinned to the direction orthogonal to the surviving mirror line of the system, and thus parallel to the zigzag direction of the honeycomb lattice. For electronic densities n el ∼ 10 10 cm −2 , and a substrate-induced gap ∆ 20 meV, the size of the BCD has been calculated to lie in the 10 −3 nm range assuming at 5 % strain. This is the same order of magnitude of the BCD in the 1H phase of monolayer TMDs discussed above.

B. Strained bilayer graphene
The appearance of a non-vanishing BCD in the complete absence of spin-orbit coupling and tilted Dirac cones is not specific of monolayer graphene. As we review below, a finite non-linear Hall voltage can appear in gated bilayer graphene subject to mechanical deformations. Even more importantly, the BCD has been shown to lie in the nm range and therefore comparable with the one of bilayer WTe 2 . To understand how such a large BCD appears in gated bilayer graphene, we first discuss strain effects on the electronic properties in Bernal stacked bilayer graphene. To do so, we start out with the effective continuum low-energy Hamiltonian 91-93 for the electrons on the A 1,2 and B 1,2 sublattices of the two 1, 2 carbon sheets. It reads: where π = k x + ik y and v 0/3 are two strain independent parameters with the dimension of a velocity. Specifically, v 0 is the Fermi velocity of the Dirac electrons in each monolayer and is thus related to the intralayer hopping amplitude γ 0 . By contrast, the velocity v 3 is related to the so-called "skew" interlayer hopping amplitude γ 3 connecting the A 1 and B 2 sites. The A 2 and B 1 lie on top of each other and are coupled by the intercoupling γ 1 . The effects of strain appear as two pseudo-gauge fields A 0/3 94 . Assuming the stresses are applied along the principal crystallographic directions, the pseudogauge fields A 0/3 = 3 4 γ 0/3 ( xx − yy )β 0/3 with β 0/3 two different dimensionless parameters related to the elastic properties of the material. Finally, the different on-site potential in the two layers ±∆/2 is due to an externally applied electric field 95 . For energies near the Fermi level the states coupled by γ 1 can be eliminated through a Schrieffer-Wolff transformation 96 . After a gauge transformation k x → k x − A 0 , one obtains the effective 2 × 2 Hamiltonian, where we introduced the effective mass m = γ 1 /2v 2 0 and the strain coupling w = 3 4 γ 3 ( xx − yy )(β 3 −β 0 ). The inclusion of the skew interlayer hopping has a twofold effect. First, it gives an explicit trigonal warping to the Fermi surface. Second, and most important, it changes the topology of the Fermi surface at low energies. In the absence of strain for instance, the low-energy quadratic band crossing point occurring for v 3 = 0 is changed in favour of the appearance of four Dirac cones, one of which is located at the K and K points of the BZ. The remaining three "leg" Dirac cones, with a Berry phase opposite to the central cone one, are arranged in a three-fold rotation symmetric fashion [see Fig. 8(a)] and have a distance in momentum given by κ L = mv 3 . The Lifshitz transition where the Fermi surface changes topology 93 occurs at L = mv 2 3 /2 in the ∆ = 0 inversion-symmetric case. The appearance of a non-vanishing BCD in the system can be understood considering the Berry curvature properties already in this strain-free (w = 0) case 77 . As shown in Fig. 8(a) the central Dirac cone has a C 3 symmetric Berry curvature and thus it gives a vanishing contribution to the BCD. When taken by themselves, the leg Dirac cones instead have only a mirror symmetric Berry curvature profile and therefore their contribution to the BCD is nonzero. This is because each of the leg gapped Dirac cones is described by the low-energy Hamiltonian introduced for monolayer graphene with a warping term such that λ 1 = λ 2 = λ 3 . However, the contributions of the leg Dirac cones cancel each other due to the threefold rotation symmetry. Once strain is introduced, the perfect cancellation of the leg Dirac cone contributions is lost and a finite overall BCD appears. In addition, the central massive Dirac cone also yields a non-zero contribution to the BCD for w = 0. As explicitly shown in Fig. 8  for a strain term w = −5 L , which corresponds roughly to a 1% strain, and a for gate-induced gap ∆ = 10 meV, the BCD D x 1 nm assuming a carrier density n 10 11 cm −295 . This value, which, as mentioned above, is comparable to the one experimentally found in bilayer WTe 2 can be boosted by an order of magnitude for smaller values of the gate-induced gap. And indeed dipoles even in the tens of nanometer scale have been recently found in the strained bilayer graphene nanoarchitectures of Ref. 78 .

V. THREE-DIMENSIONAL MATERIALS
Although the non-linear Hall effect in time-reversal symmetric conditions finds its natural realization in twodimensional materials, it can also arise in three-dimensional bulk crystals. TaIrTe 4 , a type-II Weyl semimetal ternary compound [97][98][99] , has been recently shown to possess a non-linear Hall effect surviving even at room temperature 100 . Generally speaking, the fact that a finite BCD can appear in Weyl semimetals with tilted cones [101][102][103] , be them either of type-I or of strongly overtilted type-II, can be understood considering k · p theory. Let us analyze, for simplicity, a Weyl semimetal with an isotropic Fermi velocity. Its effective low-energy Hamiltonian reads where χ = ±1 is the chirality or the topological charge of the Weyl cone 104 , and the tilt direction has been chosen to be theẑ direction. This continuum theory describes a type-II Weyl semimetal for |u x /v 0 | > 1, whereas for |u x /v 0 | < 1 the cone is of the type-I form. Since the total topological charge in the full BZ of a generic system has to vanish, each Weyl cone has to come with a partner of opposite chirality. In addition, time-reversal symmetry implies that each Weyl cone has a time-reversed partner of the same chirality. The total number of Weyl cones in a time-reversal symmetric topological semimetal is therefore 4n, with n integer. It is possible to show that each Weyl cone yields a finite BCD. This can be seen already in the untilted u x ≡ 0 limit. Using that the Berry curvature of a single Weyl cone is Ω = χk/(2 |k| 3 ) and that the carrier velocity is v = v 0k , the components of the BCD pseudotensor read as whereas the off-diagonal components identically vanish. We emphasize that there are inconsistencies in the literature as it concerns the BCD of a single Weyl cone 105 . We refer the reader to Ref. 101 for a discussion on this point. With the BCD that depends exclusively on the topological charge of the Weyl nodes, it immediately follows that in the absence of a tilt term the D ii 's summed over all nodes of a Weyl semimetal are identically zero. Assuming instead a sizable tilt term, the situation changes qualitatively. The diagonal non-vanishing components of the BCD indeed acquire a specific dependence on the dimensionless parameter δ = u z /v 0 . This has been evaluated in Ref. and reads The equations above reduce to Eq. 45 in the δ → 0 limit. In addition, it is interesting to note that, independent of the tilt term, the trace of the BCD pseudotensor is a universal quantity Tr{D} = χ/(4π 2 ). As in the untitled case, the off-diagonal components of the pseudotensor vanish as explicitly shown in Fig. 9(a). The δ-dependence of the dipole components D ii can make the total BCD finite in a generic Weyl semimetal. This is because the dipole tensors of the Weyl nodes related by time-reversal symmetry are same and are not entirely cancelled by the tensors corresponding to the Weyl nodes with opposite topological charge. This occurs in materials, such as SrSi 2 106 , where Weyl nodes of opposite topological charge are not mapped into each other by point group symmetries. According to this analysis, however, materials with additional mirror symmetries, TaAs being a paradigmatic example [107][108][109][110] , are forced to have a vanishing total BCD for the simple reason that Weyl nodes of opposite charge are mapped by the mirror symmetry, and therefore constrained to have opposite Berry curvatures.
Although this feature is in agreement with the fact that a three-dimensional material with mirror symmetries cannot have a symmetric BCD pseudotensor, it does not justify why the antisymmetric part of the BCD should be identically zero. Density functional theory calculations have shown that TaAs has a large dipole 112 with D xy = −D yx = 0 that is symmetry-allowed. The resolution of this paradox comes from the fact that higher-order terms in momentum k can change the properties of the Berry curvature. This can be demonstrated considering a k · p model for a pair of Weyl nodes in the presence of a mirror symmetry 111 . It reads This models predicts the presence of two Weyl nodes at k y = ± √ λ with a corresponding Fermi velocity v y = λ/m. Fig. 9(b) shows the corresponding D xy dipole surface density assuming a perfectly isotropic Fermi velocity for the two Weyl cones. As compared to the result obtained using the linear Weyl effective Hamiltonian [c.f. Fig. 9(a)], it is clear that the model Eq. 47 yields a finite antisymmetric BCD [c.f. Fig. 9] in agreement with the density functional theory calculations results for the type-I TaAs and the type-II MoTe 2 Weyl semimetals.
It is interesting to note that in the λ < 0 parameter range the model in Eq. 47 predicts an insulating state with a gap ∆ = (|λ| /m). Even in this regime there is a finite dipole (surface) density antisymmetric component [see Fig. 9(c)], provided the system can be doped with charge carriers and thus possesses a finite Fermi surface. The Rashba semiconductor BiTeI 113 has been predicted to become a strong three-dimensional topological insulator 114 at moderate pressures > 3 GPa. These two phases are separated by an intermediate Weyl semimetal phase 115 . These "topological" phase transition can be therefore captured by the model in Eq. 47 assuming that changing the pressure the parameter λ, or equivalently the bulk gap ∆, is varied. Fig. 9(d) shows the behavior of the antisymmetric BCD d z = (D xy − D yx )/2 computed in a recent density functional theory calculation study 111 for various carrier density. It highlights how the BCD gets strongly enhanced while approaching the Weyl semimetal phase. This is in agreement with the fact that the model of Eq. 47 predicts a divergent antisymmetric BCD in the insulating regime when the ∆ → 0 limit is taken [c.f. the inset in Fig. 9(d)]. This demonstrates how non-linear Hall measurements can provide a unique signature of the pressure-induced topological phase transition 116 in BiTeI, even if local order parameters are absent.

VI. CONCLUSIONS
In this concluding section, we summarize what has been achieved with the recent studies on the NLHE in timereversal symmetric conditions, and also point out possible directions for future research on this recently discovered phenomenon. Being regulated by the Berry curvature dipole, the time-reversal symmetric NLHE can be used as a direct probe of the geometry of the Bloch states in non-magnetic systems. This intrinsic quantum property can be only unveiled in crystals with unusually low crystalline symmetry. Bilayer WTe 2 has been the first investigated material where such symmetry requirements are fulfilled. It has led to the direct observation and measurement 75 of the BCD in a pristine crystal. The crystalline symmetry requirements for a non-vanishing BCD can be however also achieved using the capability of externally applied fields to lower the point group symmetry of a crystal. Breaking inversion symmetry, externally applied electric fields, for instance, allow for a time-reversal symmetric NLHE even in crystal with a native Berry curvature that identically vanishes. This scheme has been successfully applied in the quantum spin Hall 1T phase of transition metal dichalcogenides 73 . Mechanical deformations are yet another knob that can be used to decrease the crystalline symmetry of materials and eventually allow for a NLHE. Contrary to external electric fields, strain does not break the centrosymmetry of a crystal. However, it can reduce the rotational symmetry to make this effect visible. In principle, also externally applied magnetic field could be used to meet to lower the symmetry and meet the NLHE requirements. A planar magnetic field breaks the time-reversal invariance. However, the absence of Lorentz force combined with the presence of a residual mirror symmetry would guarantee the complete absence of a linear dissipationless Hall conductance thus making the non-linear Hall effect the only present transverse resistance. This effect has been recently dubbed anomalous non-linear planar Hall effect 34 .
Independent of the presence and nature of the external perturbation, the appearance of the time-reversal symmetric NLHE was originally believed to necessitate a substantial spin-orbit coupling. However, the experimental demonstration of non-linear transverse currents in strained (bilayer) grahene 78 has shown that this effect could be also pursued in systems made of light elements, where the presence of spin-orbit coupling could be neglected all together. This has opened new material avenues for studies of the NLHE with time-reversal symmetry. Twisted bilayer graphene -a material structure where correlated insulating behavior 117 and superconductivity 118 have been shown to arisehas been recently suggested to host a substantial BCD 119,120 on the order of tens of nanometers. In twisted bilayer graphene strain with an opposite sign in the two graphene layers has been imaged in STM experiments 121 . Moreover, when encapsulated with hexagonal boron nitride the inversion symmetry of twisted bilayer graphene is broken. As a result, all symmetry requirements for a finite BCD are fulfilled.
We conclude by mentioning that non-linear transverse currents in the presence of time-reversal symmetry have been suggested to occur also in the strong three-dimensional topological insulators of the Bi 2 Te 3 material class 122 . The anomalous single surface Dirac cones in these materials are characterized by an intrisic chirality that enables a novel skew scattering 123 similar to the one predicted in unstrained bilayer graphene 17 . However, it is important to note that the non-linear response found in these materials does not represent a dissipationless non-linear Hall conductance, and is therefore compatible with the trigonal symmetry of the Bi 2 Te 3 surface states 124 . Finally, non-linear responses using the topological surface states of Bi 2 Te 3 have been demonstrated to occur also in the presence of external planar magnetic fields 125 .