Diego García Ovalle

Armando Pezo

Aurélien Manchon
email: aurelien.manchon@univ-amu.fr

Influence of the surface states on the nonlinear Hall effect in Weyl semimetals

We investigate the influence of surface states on the nonlinear Hall response driven by the Berry curvature dipole in non-centrosymmetric time-reversal invariant Weyl semimetals. To do so, we perform a tomography of the Berry curvature dipole in a slab system using a minimal two-band model. We find that in the type-I phase, the nonlinear Hall response is not particularly sensitive to the presence of Fermi arcs or other trivial surface states. However, in the type-II phase, we find that these surface states, be they topologically trivial or not, contribute substantially to the Berry curvature dipole, leading to a strong thickness dependence of the nonlinear Hall response. This feature depends on the nature of the surface states and, henceforth, on the slab geometry adopted. In order to assess the validity of this scenario for realistic systems, we performed Berry curvature dipole calculations by first principles on the WTe2, confirming the dramatic impact of surface states for selected slab geometries. Our results suggest that surface states, being topological or not, can contribute much more efficiently to the nonlinear Hall response than bulk states. This prediction is not limited to topological semimetals and should apply to topologically trivial non-centrosymmetric materials and heterostructures, paving the way to interfacial engineering of the nonlinear Hall effect.

I. INTRODUCTION

Conventional wisdom inherited from Hall's foundational work [START_REF] Hall | On a New Action of the magnet on Electric[END_REF][START_REF] Hall | On the "Rotational Coefficient" in nickel and cobalt[END_REF] states that Hall currents flowing transverse to the injected charge direction are only permitted as long as time-reversal symmetry is globally broken. In other words, ordinary, anomalous or topological Hall effects only exist either in the presence of an external magnetic field or in magnetic materials displaying a net magnetization [START_REF] Nagaosa | Anomalous Hall effect[END_REF]. However, this long-lived statement has been recently challenged by two important observations. First, it has been realized that anomalous Hall effect does not necessitate the presence of a net magnetic moment to exist. For instance, in an antiferromagnet the Hall effect can be nonzero as long as no crystal symmetry reinstalls the time-reversal effectively. This is particularly true for certain classes of antiferromagnets with a non-collinear magnetic configuration [START_REF] Chen | Anomalous Hall Effect Arising from Noncollinear Antiferromagnetism[END_REF][START_REF] Kubler | Non-collinear antiferromagnets and the anomalous Hall effect[END_REF][START_REF] Nakatsuji | Large anomalous Hall effect in a non-collinear antiferromagnet at room temperature[END_REF][START_REF] Nayak | Large anomalous Hall effect driven by non-vanishing Berry curvature in non-collinear antiferromagnet Mn3Ge[END_REF] (see also Ref. [START_REF] Bonbien | Topological Aspects of Antiferromagnets[END_REF]). Second, the time-reversal symmetry breaking necessary to obtain Hall effect is not required anymore at the second order in the electric field. Under certain conditions, anomalous Hall effect can appear in non-magnetic materials to the second order of the electric field [START_REF] Sodemann | Quantum nonlinear hall effect induced by berry curvature dipole in time-reversal in-variant materials[END_REF][START_REF] Ma | Observation of the nonlinear hall effect under time-reversal-symmetric conditions[END_REF][START_REF] Kang | Nonlinear anomalous hall effect in few-layer WTe2[END_REF]. Recently, it has been proposed that such a second-order Hall effect exists in collinear antiferromagnets where the linear Hall effect is absent, e.g., CuMnSb [START_REF] Shao | Nonlinear anomalous hall effect for néel vector detection[END_REF] .Whereas the intrinsic contribution to the linear anomalous Hall effect in ferromagnetic and antiferromagnetic compounds is associated with the Berry curvature of the material's ground state, the second order, nonlinear anomalous Hall effect (NLHE) is rather associated with the Berry curvature dipole (BCD). From a symmetry standpoint, the minimal requirement is inversion symmetry breaking which ensures that the Berry curvature does not vanish, but this is not sufficient: mirror symmetry also needs to be broken to obtain a finite BCD. Sodemann and Fu [START_REF] Sodemann | Quantum nonlinear hall effect induced by berry curvature dipole in time-reversal in-variant materials[END_REF] identified the crystallographic point groups that possess the minimal requirements for the observation of NLHE, a study recently refined by Du et al. [START_REF] Du | Quantum theory of the nonlinear Hall effect[END_REF] to include both intrinsic (BCD related) and extrinsic mechanisms [START_REF] Nandy | Symmetry and quantum kinetics of the nonlinear hall effect[END_REF] allowed by symmetry.

Different material candidates have been explored experimentally and theoretically as suitable options to obtain NLHE. From an experimental point of view, quadratic responses in the electric field have been detected, among others, in bilayers and few layers of WTe 2 [START_REF] Ma | Observation of the nonlinear hall effect under time-reversal-symmetric conditions[END_REF][START_REF] Kang | Nonlinear anomalous hall effect in few-layer WTe2[END_REF] which have C 2v point group. From a theoretical perspective, other possibilities for large BCD values have been proposed, including transition metal dichalcogenides [START_REF] Du | Band signatures for strong nonlinear hall effect in bilayer WTe2[END_REF][START_REF] You | Berry curvature dipole current in the transition metal dichalcogenides family[END_REF][START_REF] He | Giant nonlinear hall effect in twisted bilayer WTe2[END_REF][START_REF] Singh | Engineering weyl phases and nonlinear hall effects in t d -MoTe2[END_REF][START_REF] Zhang | Electrically tuneable nonlinear anomalous hall effect in twodimensional transition-metal dichalcogenides WTe2 and MoTe2[END_REF][START_REF] Zhou | Highly tunable nonlinear hall effects induced by spin-orbit couplings in strained polar transition-metal dichalcogenides[END_REF], graphene [START_REF] Battilomo | Berry curvature dipole in strained graphene: A fermi surface warping effect[END_REF][START_REF] Pantaleón | Tunable large berry dipole in strained twisted bilayer graphene[END_REF][START_REF] Zhang | Giant nonlinear hall effect in strained twisted bilayer graphene[END_REF], and especially Dirac [START_REF] Samal | Nonlinear transport without spin-orbit coupling or warping in two-dimensional dirac semimetals[END_REF] and Weyl semimetals (WSM) [START_REF] Zeng | Nonlinear transport in weyl semimetals induced by berry curvature dipole[END_REF] (See also Ref. [START_REF] Facio | Strongly enhanced berry dipole at topological phase transitions in bitei[END_REF]). Ab initio simulations [START_REF] Zhang | Berry curvature dipole in weyl semimetal materials: An ab initio study[END_REF] have stimulated further theoretical studies of BCD in 3D WSMs, because of their rich geometrical features and their potential benefits to create highly efficient electronic transport devices. NLHE in WSMs has also been verified analytically by applying perturbation theory [START_REF] Gao | Second-order nonlinear hall effect in weyl semimetals[END_REF], suggesting that the transport is sensitive to intraband processes, the chemical potential and the tilting of the Weyl nodes. Besides, similar optical effects such as second harmonic generation can be confirmed by applying Floquet [START_REF] Morimoto | Semiclassical theory of nonlinear magneto-optical responses with applications to topological dirac/weyl semimetals[END_REF] and manybody quantum formalisms [START_REF] Rostami | Nonlinear anomalous photocurrents in weyl semimetals[END_REF]. Last but not least, a full Green function theory of the NLHE has been proposed recently [START_REF] Du | Quantum theory of the nonlinear Hall effect[END_REF], pointing out the differences between extrinsic and intrinsic contributions. The former is associated with higher moments of the impurity potential whereas the latter is solely associated to the band structure. In this article, we focus on the nonlinear response arising on the intrinsic mechanism driven by the BCD. Nonmagnetic WSMs such as TaAs or WTe 2 are particularly interesting platforms for the realization of NLHE because inversion symmetry is necessarily broken and Berry curvature diverges at the Weyl nodes. WSMs possess pairs of doubly degenerate linearly dispersive states, forming Weyl cones at Fermi level [START_REF] Wan | Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridates[END_REF]. According to the Nielsen-Ninomiya theorem [START_REF] Nielsen | The adler-bell-jackiw anomaly and weyl fermions in a crystal[END_REF], each pair of nodes carries Berry curvature monopoles of opposite chirality which are connected via Fermi arcs lying at opposite surfaces of the slab [START_REF] Armitage | Weyl and dirac semimetals in three-dimensional solids[END_REF]. Type I WSMs, such as elemental Tellurium [START_REF] Rodriguez | Two linear regimes in optical conductivity of a type-i weyl semimetal: The case of elemental tellurium[END_REF], Janus superlattices [START_REF] Meng | A type of robust superlattice type-i weyl semimetal with four weyl nodes[END_REF] and Ta or As compounds [START_REF] Lv | Experimental discovery of weyl semimetal TaAs[END_REF][START_REF] Xu | Discovery of a weyl fermion semimetal and topological fermi arcs[END_REF][START_REF] Sun | Topological surface states and fermi arcs of the noncentrosymmetric weyl semimetals TaAs, TaP, NbAs, and NbP[END_REF], are characterized by point-like Fermi surface in the bulk and vanishing density of states. Type II WSMs, such as MoTe 2 [START_REF] Deng | Experimental observation of topological fermi arcs in type-ii weyl semimetal MoTe2[END_REF][START_REF] Tamai | Fermi arcs and their topological character in the candidate type-ii weyl semimetal MoTe2[END_REF] and WTe 2 [START_REF] Bruno | Observation of large topologically trivial fermi arcs in the candidate type-ii weyl semimetal WTe2[END_REF][START_REF] Wu | Observation of fermi arcs in the type-ii weyl semimetal candidate WTe2[END_REF] but also the magnetic candidate Co 3 Sn 2 S 2 [START_REF] Guin | Zero-field nernst effect in a ferromagnetic kagome-lattice weyl-semimetal Co3Sn2S2[END_REF][START_REF] Morali | Fermi-arc diversity on surface terminations of the magnetic weyl semimetal Co3Sn2S2[END_REF][START_REF] Li | Surface states in bulk single crystal of topological semimetal Co3Sn2S2 toward water oxidation[END_REF][START_REF] Wang | Large intrinsic anomalous hall effect in half-metallic ferromagnet Co3Sn2S2 with magnetic weyl fermions[END_REF][START_REF] Okamura | Giant magneto-optical responses in magnetic weyl semimetal Co3Sn2S2[END_REF][START_REF] Ikeda | Critical thickness for the emergence of weyl features in Co3Sn2S2 thin films[END_REF], offer a slightly different paradigm as the Weyl cone spectrum is tilted in momentum space, breaking Lorentz invariance. As a result, the Weyl points arise at the boundary between electron and hole pockets. Notice that certain compounds can support type I as well as type II Weyl nodes [START_REF] Meng | Ternary compound HfCuP: An excellent weyl semimetal with the coexistence of type-i and type-ii weyl nodes[END_REF][START_REF] Zhang | Coexistence of type-i and type-ii weyl points in the weyl-semimetal OsC2[END_REF].

A remarkable aspect of WSMs is the nature of their surface states. As mentioned above, alike topological insulators WSMs possess topologically protected surface states in the form of spin-momentum locked Fermi arcs that connect bulk Weyl nodes of opposite chirality. In type I WSMs, the Fermi arcs coexist with the projection of electron (or hole) pockets when the chemical potential lies away from the neutrality point. In type II WSMs, the Fermi arcs coexist with projected electron and hole pockets irrespective of the value of the chemical potential, as well as with trivial closed loops called "track states" [START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF]. As a consequence, surface states of WSMs can be rich, resulting in enhanced Edelstein effect [START_REF] Johansson | Edelstein effect in Weyl semimetals[END_REF], and unconventional patterns in quantum oscillation experiments [START_REF] Potter | Quantum Oscillations from Surface Fermi-Arcs in Weyl and Dirac Semi-Metals[END_REF][START_REF] Bulmash | Quantum oscillations in Weyl and Dirac semimetal ultrathin films[END_REF][START_REF] Wang | Anomalous Phase Shift of Quantum Oscillations in 3D Topological Semimetals[END_REF] (see also Ref. [START_REF] Moll | Transport evidence for Fermi-arc-mediated chirality transfer in the Dirac semimetal Cd3As2[END_REF]). Previous works pointed out that topological materials defined in slab geometries can exhibit interesting transport properties due to finite size effects and the behavior of surface states inside the samples. In this context, it has been shown that the anomalous Hall conductivity is highly influenced by surface states such as Fermi arcs in Weyl systems without time-reversal symmetry [START_REF] Breitkreiz | Large contribution of fermi arcs to the conductivity of topological metals[END_REF], even in presence of disorder [START_REF] Lopez | Multiterminal conductance at the surface of a weyl semimetal[END_REF]. Additional studies have been performed in confined geometries to clarify, among others, the behavior of chiral magnetic effects [START_REF] Gorbar | Chiral separation and chiral magnetic effects in a slab: The role of boundaries[END_REF], the magnetoresistance [START_REF] Alekseev | Magnetoresistance of compensated semimetals in confined geometries[END_REF] and the quantum Hall effect in Dirac semimetals [START_REF] Schumann | Observation of the quantum hall effect in confined films of the three-dimensional dirac semimetal Cd3As2[END_REF].

In this work, we seek to understand how the surface states of nonmagnetic WSMs influence the NHLE response driven by the BCD. To do so, we consider a minimal 2-band model of a time-reversal invariant WSM with inversion symmetry breaking in a slab geometry, so that bulk and surface states are treated on equal footing [START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF]. This model exhibits four Weyl points: the minimum number of degeneracies due to time-reversal symmetry. These four points are associated with local divergencies of the Berry curvature, as depicted on Fig. 1. In this sense, after neglecting the vanishing components of the BCD tensor due to mirror symmetries in the 3D lattice, Zeng et al. [START_REF] Zeng | Nonlinear transport in weyl semimetals induced by berry curvature dipole[END_REF] recently reported that NLHE requires Weyl cone tilting and an asymmetric Fermi surface when the nodes lie at the same energy. The NLHE is also influenced by the distance between nodes. Accordingly, our study gives further insight about the implications of the Fermi arc configurations on the BCD. Remarkably, it also complements a recent study that comprises a surface BCD due to the projection of Fermi arcs in type I WSMs [START_REF] Wawrzik | Infinite berry curvature of weyl fermi arcs[END_REF].

Our paper is organized as follows. In Section II, we recall the general formalism for nonlinear Hall transport driven by BCD and basic physical considerations for our model of interest. In Section III, we provide and discuss our results regarding the connection between surface states and NLHE, its dependency on the Weyl cones tilting and its layer decomposition. In Section IV, we also collate these outcomes with realistic numerical simulations on a WTe 2 slab with different cuts. Finally, we summarize and state our main conclusions in Section V.

II. GENERAL THEORY AND MODEL

Let us begin by recalling basic elements of NLHE theory [START_REF] Sodemann | Quantum nonlinear hall effect induced by berry curvature dipole in time-reversal in-variant materials[END_REF]. We start with a nonmagnetic crystal and its corresponding Bloch Hamiltonian H k , whose eigenstates |u nk satisfy Schrödinger equation H k |u nk = nk |u nk with eigenenergies nk . If the crystal is subjected to a sufficiently small electric field E =Re[E exp(iωt)], E ∈ C 3 , such that the adiabatic limit is still valid, and assuming weak disorder, the first order Hall current is induced by the Berry curvature [START_REF] Nagaosa | Anomalous Hall effect[END_REF]

Ω nk = i m =n u nk |v k |u mk × u mk |v k |u nk ( nk -mk ) 2 , (1)
where vk = ∂ k H k is the velocity operator. Upon time reversal operation T we have T Ω nk T -1 = -Ω n-k , and thus the anomalous Hall effect vanishes. Nonetheless, when inversion symmetry P is further broken, then PΩ nk P -1 = Ω n-k , and one can show to the lowest order in scattering time τ , that the (rectified) second order Hall current reads [9]

j 0 a = e 3 τ 2 2 adc D bd E b E * c . (2)
Here, latin indices refer to the components of the usual cartesian basis and adc is the Levi-Civita tensor. D bd is the BCD defined as [9]

D bd = BZ d 3 k (2π) 3 n (v nk ) b (Ω nk ) d ∂f nk ∂ nk , (3)
with v nk = ∂ k nk the eigenvalues of the velocity operator and f nk the Fermi distribution function. In the zero-temperature limit, ∂ nk f nk → -δ( nk -µ), setting µ as the chemical potential. Equation ( 2) neglects the intrinsic contribution to NLHE that appears when timereversal is broken as well as higher order extrinsic contributions (i.e., side-jump and skew scattering) [START_REF] Du | Quantum theory of the nonlinear Hall effect[END_REF][START_REF] Nandy | Symmetry and quantum kinetics of the nonlinear hall effect[END_REF].

The nonmagnetic WSM slab is built on the spinless 2-band model in a cubic lattice introduced in Ref. [START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF]. The bulk Hamiltonian reads

H B = γ(cos 2k z -cos k 0 )(cos k x -cos k 0 )σ 0 -m(1 -cos 2 k x -cos k y ) + 2t x (cos k z -cos k 0 ) σ1 -2t sin k y σ2 -2t cos k x σ3 , (4)
with σi , i = 1...3 being the 2 × 2 Pauli matrices, σ0 = 1 2×2 and the cubic Brillouin zone is C = [-π, π] 3 . The parameters of Eq. ( 4) are set to

k 0 = π 2 , m = 2, t x = 1 2
and t = 1. This WSM possesses four Weyl nodes of zero energy located at k * = ± π 2 (x + ẑ) (see Fig. 1). Importantly, for the parameters adopted in this work, γ = 2 sets the transition point between type I (γ < 2) and type II (γ > 2) WSM phases. In the remaining of this paper, the transport properties will be investigated as a function of the chemical potential µ as well as of the tilting of the Weyl cones controlled by γ. Interestingly, Eq. ( 4) is constrained by the following mirror symmetries [START_REF] Zeng | Nonlinear transport in weyl semimetals induced by berry curvature dipole[END_REF]

M † x H B (k x , k y , k z )M x = H B (-k x , k y , k z ), (5)
M † z H B (k x , k y , k z )M z = H B (k x , k y , -k z ), (6)
which impose that the Fermi arcs connecting the Weyl nodes lie on the (z, x) surfaces. In addition, these symmetries constrain the BCD tensor, so that the only nonvanishing elements are D zx and D xz . In other words, a second order Hall current can only be obtained in the (y, z) (≡ D zx ) and (x, y) (≡ D xz ) planes. Consequently, the Fermi arcs cannot contribute to the NLHE response. Nevertheless, as discussed below, other trivial surface states can substantially impact the second order response.

We now design the slabs by discretizing Eq. ( 4) along a given direction, x, ŷ or ẑ, of the cubic Brillouin Zone (see, e.g., Ref. [START_REF] Manchon | Semirealistic tight-binding model for spin-orbit torques[END_REF]). The Hamiltonian loses periodicity along the chosen axis. The new Hamiltonian H n S of size 2L × 2L, with n = x, ŷ, ẑ depending on the growth direction, and L the number of layers, is given by

H n S =             H n 0 H n 1 H n 2 0 ... 0 H n † 1 H n 0 H n 1 H n 2 . . . . . . H n † 2 H n † 1 . . . . . . . . . 0 0 H n † 2 . . . . . . . . . H n 2 . . . . . . . . . . . . . . . H n 1 0 . . . 0 H n † 2 H n † 1 H n 0             . ( 7
)
For the block matrices in Eq. ( 7), H n 0 is the intralayer Hamiltonian that retains in-plane periodicity after the cut, while H n 1 and H n 2 are the nearest neighbor and second nearest neighbor interlayer Hamiltonian, respectively. For the cut along x, the block matrices are given by

H x 0 = -γ cos k 0 (cos 2k z -cos k 0 )σ 0 -m 1 2 -cos k y + 2t x (cos k z -cos k 0 ) σ1 -2t sin k y σ2 , (8)
H x 1 = γ 2 (cos 2k z -cos k 0 ) σ0 -tσ 3 , (9)
H x 2 = m 4 σ1 . (10)
A cut along ŷ leads to

H ŷ 0 = γ(cos 2k z -cos k 0 )(cos k x -cos k 0 )σ 0 -m sin 2 k x + 2t x (cos k z -cos k 0 ) σ1 -2t cos k x σ3 , (11)
H ŷ 1 = m 2 σ1 + itσ 2 . (12)
Finally, a cut along ẑ gives We illustrate the slab band structures for each cut in Fig. 2, for γ = 1 (type I WSM) and γ = 3 (type II WSM). Since the slab Hamiltonians retain a Berry curvature perpendicular to the plane of the slab, we define D x zx as the BCD associated with H x S when the electric field is applied along ẑ, and D ẑ xz is the BCD associated with H ẑ S when the electric field is along x. These are the only relevant BCD components for the slab geometries presented above due to the symmetry restrictions. We now move forward to the method for calculating the layer decomposition of transport coefficients in the slabs. Since the periodic part of the Bloch function can be written in terms of the complete layer basis as |u nk =

H ẑ 0 = -γ cos k 0 (cos k x -cos k 0 )σ 0 -m(1 -cos 2 k x -cos k y ) -2t x cos k 0 σ1 -2t sin k y σ2 -2t cos k x σ3 , (13)
H ẑ 1 = -t x σ1 , (14)
H ẑ 2 = γ 2 (cos k x -cos k 0 ) σ0 . (15)
D lk = - 1 π Im Tr S l G R k , (16)
with the retarded Green's function

G R k = lim η→0 + ( + iη)1 -H n S -1 . ( 17
)
This decomposition provides an insight about the influence of the surface states for each growth direction. Additionally, in order to obtain the contribution of layer l to the NLHE response along, say, ŷ, one simply needs to perform the substitution vy → S l vy into Eq. ( 1) and then into (3).

III. RESULTS AND DISCUSSION

A. Band structure and surface states

Let us first consider how the Weyl cone tilting impacts the slab band structure and its surface states. Since only D x zx and D ẑ xz are non-vanishing, and because the Weyl nodes are located in the (k x , k z ) plane, we focus on slabs composed of L = 25 layers and normal to the x and ẑ directions. Without loss of generality, we fix the chemical potential to µ = 0.2 and select the exemplary cases γ = 1, 3 to illustrate the differences between a type I and a type II WSM, respectively. To better understand how bulk and surface states evolve upon tilting the Weyl cones, we represent the density of states in the 2D Brillouin zone and projected on a set of representative layers, ranging from the bottommost surface to the topmost surface. Our results are depicted in Fig. 3(a,c) for the x-cut and Fig. 3(b,d) for the ẑ-cut.

Figures 3(a,c) and 3(b, d) show different behaviors of the bulk and surface states when tuning the tilting parameter. When n = x [Fig. 3(a,c)], notice two distinct situations depending on the value of γ: if γ = 1 [type I, Fig. 3(a)], the density of states is dominated by the Weyl nodes across the whole slab, with a surface state composed of the projected nodes connected by degenerate Fermi arcs. When γ = 3, [type II, Fig. 3(c)], the density of states of the central layer is composed of electron and hole pockets touching at the Weyl nodes, as expected. These pockets result in large trivial surface Fermi pockets enclosing the nodes, whose connections change direction from ẑ (type I) to x (type II) WSMs [START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF]. A map of the Fermi arcs' reconnection and surface states' evolution as a function of the tilting γ can be found in Appendix A.

On the other hand, for n = ẑ [Fig. 3 the slab, for γ = 3 [type II, Fig. 3(d)], a surface contribution appears due to the emergence of track states [START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF]. Nonetheless, the surface Fermi pockets associated with the projection of the bulk electron and hole pockets remain very small and the surface states are dominated by the Weyl nodes, in sharp contrast with the x-cut discussed in Fig. 3(a,c). We therefore expect the surface states to have completely different impact on the NLHE signal in these two different slab geometries.

B. Nonlinear Hall response

Let us now turn our attention towards the NLHE response for the two slabs. We compute the BCD components given by Eq. ( 4) in both slab geometries, (a) D x zx /L (x-cut) and (b) D ẑ xz /L (ẑ-cut), as well as in the 3D bulk structure, D 3D zx(xz) (black symbols), in order to better identify the impact of the surface states. Notice that the BCD calculated in the slab geometry is normalized by the number of layers L to allow for a quantitative comparison with the BCD calculated in the bulk. By definition, the BCD is therefore unitless. In the following, we set the chemical potential to µ = 0.2 and represent the BCDs as a function of the tilting parameter γ, as displayed in Figs. 4(a-b). In addition, the corresponding ratio between the BCD and the density of states is reported on Figs. 4(c-d) and discussed further below. For the x-cut [Fig. 4(a)], the bulk BCD (D 3D zx , black symbols) displays a nonlinear dependence as a function of the tilting parameter, as observed by Zeng et al. [START_REF] Zeng | Nonlinear transport in weyl semimetals induced by berry curvature dipole[END_REF]: it first increases smoothly with γ, reaches a maximum and decreases for large γ. We attribute this tendency to the contribution of the Fermi surface between the Weyl nodes, as also mentioned by Zeng et al. [START_REF] Zeng | Nonlinear transport in weyl semimetals induced by berry curvature dipole[END_REF], and the upper limit that should reach the tilting parameter regarding the inclination of the Weyl nodes. In contrast, the slab BCD (D x zx , colored symbols) increases sharply and reaches a plateau at large γ. Furthermore, the value of this plateau is substantially larger than the maximum value obtained in the bulk although it slightly decreases when increasing the number of layers, suggesting that surface states play a major role in the NLHE response. Since the transition between type I and type II regimes produces a change in the configuration of the surface states, with the apparition of large Fermi pockets [Fig. 3(c)], it is clear that the mismatch between the slab and bulk responses arise from the influence of the Fermi pockets' projections on the BCD across the slab.

Conversely, for the ẑ-cut [Fig. 4(b)], the bulk and slab BCDs show a similar behavior, increasing continuously and displaying local plateaus at the type I-type II transition, as well as for large tilting parameter γ. Based on Fig. 3(d), we infer that the surface track states and other trivial states do not affect qualitatively the behavior of BCD in finite slab samples. In fact, in the ẑ-cut, the slab BCD converges faster towards the bulk BCD when increasing the number of layers than in the x-cut configuration.

Since the density of states of carriers also changes when tuning the band structure, we report the corresponding ratio between the BCD and the density of states in Figs. 4(c-d). In Fig. 4(c), for a slab normal to x, the renormalized BCD exhibits a qualitatively similar trend across the type I-to-type II transition, although the slab BCD is markedly larger than the bulk BCD when approaching the type II regime. Therefore, the Fermi pockets' projections on this slab enhance the value of the BCD due to the reconfiguration of the surface states and the presence of more local states within the unit cell. In contrast, for a slab normal to ẑ in Fig. 4(d) shows a smaller discrepancy between the slab and bulk calculations suggesting that the BCD is weakly impacted by the track states and other trivial states.

The calculations displayed in Fig. 4 have been performed at µ = 0.2. Looking a the band structures in Fig. 2, one expects the relative contribution of bulk and surface states to the NLHE to vary when increasing the chemical potential. In particular, when µ leaves the regions of the Weyl cones, the bulk states should dominate the transport. Nonetheless, as further discussed in Appendix B, we find that as long as the chemical potential lies close to the Weyl cones, the distinct behavior of the BCD in the two different slabs, x-cut and z-cut, is qualitatively the same in the type II regime: the surface states substantially contribute to the nonlinear transport for the x-cut, whereas they are negligible in the z-cut.

To complement this analysis, we compute the layer decomposition of the NLHE in slab geometries, for γ = 1, 3 and L = 25. The result is reported in Fig. 5 for (a) γ = 1 and (b) γ = 3 for x-(black) and ẑ-cut (red). As we can see from Fig. 5(a), in a type I WSM the NLHE response strongly differs at the edges, i.e, at the top (l=1) and bottom (l=25) layers. From Figs. 3(a) and (b) we attribute this behavior to the surface states driven by the degenerate Fermi arcs that only appear in the x-cut. In the case of the type II WSM, Fig. 5(b), we note that the magnitude of the central layer (l = 13) in the x-cut is larger than that in the ẑ-cut. Actually, the main contributions stem from layers located underneath the surfaces (l = 5, l = 21) rather than on the surface layers. This fact corroborates our claim that in type II WSMs a decisive factor to enhance NLHE is the presence of Fermi pockets' projections at the surface rather than track states. It is worth mentioning that as long as the surface states are topologically protected (i.e., Fermi arcs), a weak influence of surface disorder is expected. Nonetheless, topologically trivial states, such as the ones that contribute the most to the NLHE in our study, are more sensitive to surface disorder, which should impact their overall response. We leave this observation to further studies.

Finally, let us comment on the experimental signature of the surface contribution to the NLHE. Based on the analysis provided above, it is clear that the surface states contribute substantially to the overall signal in the x-cut, whereas it contributes marginally to the signal in the ẑcut. As a consequence, one expects a strong difference between these two slabs when measuring the NLHE as a function of the thickness. Figure 6 displays the thickness dependence of the NLHE in a type II WSM calculated at µ = 0.2 and γ = 3, in both ẑ-cut (blue) and x-cut (red). Whereas the NLHE of both slabs converge towards the bulk values at large thicknesses (dashed lines), they exhibit a distinct behavior at small thicknesses. In the x-cut where the surface states contribute massively to the NLHE (see Fig. 4), the signal increases substantially at small thicknesses, whereas in the ẑ-cut no such behavior is reported. This distinct feature can be used as an indication of surface state-driven NLHE in experiments.

IV. NONLINEAR HALL EFFECT IN WTE2 SLABS

Let us now consider a realistic system, WTe 2 in its orthorhombic phase, and compute the NLHE from first FIG. 6. (Color online) Thickness dependence of the BCD in a type II WSM for the ẑ-cut (blue) and x-cut (red), as compared to the bulk values (dashed lines). The x-cut shows a very strong deviation from the bulk value at small thicknesses, illustrating the importance of the surface states in this case.

principles. WTe 2 is a well-known type II WSM [START_REF] Bruno | Observation of large topologically trivial fermi arcs in the candidate type-ii weyl semimetal WTe2[END_REF][START_REF] Wu | Observation of fermi arcs in the type-ii weyl semimetal candidate WTe2[END_REF], in which NLHE has been originally reported [START_REF] Ma | Observation of the nonlinear hall effect under time-reversal-symmetric conditions[END_REF][START_REF] Kang | Nonlinear anomalous hall effect in few-layer WTe2[END_REF]. For the density functional theory simulations [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF], we used the Perdew-Burke-Ernzerhof [START_REF] Perdew | Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation[END_REF][START_REF] Perdew | Generalized gradient approximation made simple[END_REF] exchangecorrelation functional. The geometry optimizations were performed using a plane-wave basis as implemented in the Vienna Ab-initio Simulation Package (VASP) [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]. We have employed 400 eV for the plane-wave expansion cutoff with a force criterion of 5 µev/ Åand a reciprocal space sampling containing 16 × 16 × 14 k-points within the Brillouin zone. The ionic potentials were described using the projector augmented-wave (PAW) method [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF], post-processing calculations were performed using Wan-nierTools [START_REF] Wu | Wanniertools : An open-source software package for novel topological materials[END_REF]. The band structure is displayed in Fig. 7 with the inset showing the unit cell. The band crossings are located within the X -Γ path in momentum space, such that the Fermi level was set to zero near this region. With this consideration, we can project the Weyl points on selected surfaces, i.e., ẑ (corresponding to the (001) direction) and x (corresponding to the (100) direction).

Let us first look at the density of states in the bulk and at the surface. Fig. 8 displays the projected density of states in the bulk (a-c) and at the surface (b-d) for a slab cut along the ẑ direction (a-b), and for a slab cut along the x direction (c-d). As can be observed from Figs. 8(ab), in the ẑ-cut the resulting density of states of the slab is larger at the surface than in the bulk, especially near the origin ( Γ point). On the other hand, the opposite situation happens when the slab is normal to x, as can be noticed from Figs. 8(c-d). From these results and comparing to the lattice model in Fig. 3, we deduce that We now move on to the computation of the surface contribution of the NLHE for the two cuts considered above. For the calculation of the BCD, we have used Eq. ( 3) performing the sum over a sample Brillouin zone containing 500 × 500 k-points. Figures 9(a) and (b) display the band structure and the BCD, respectively, computed for a slab cut along the ẑ direction. The color bar in (a) represents the projections on the top and bottom surfaces of a slab containing 25 layers. One can distinguish both surfaces by the dark and yellow lines near the Fermi level. In panel (b), the BCD is computed as a function of the energy for slabs containing an increasing number of layers, from 15 to 25. It is clear that the energy profile of the BCD strikingly depends on the slab thickness. This thickness dependence reflects the influence of the surface states. For the sake of comparison, we also reported the value of the BCD computed in the bulk (dashed line). For an infinitely thick slab, the peaks present below the Fermi level and associated with the surface states disappear. Similarly, Figs. 9(c,d) display the band structure and BCD, respectively, calculated for a slab cut along the x direction. In contrast to the ẑcut discussed above, the surface states cannot be clearly identified in the band structure that is instead dominated by bulk states. In panel (d), the BCD is computed for three different slab thicknesses. Interestingly, the qualitative behavior of the BCD is similar, displaying a peak whose position is weakly influenced by the slab thickness. For an infinitely thick slab, the bulk BCD (dashed line) conserves the overall structure, with a slight reduction in magnitude. These calculations show that surface states do substantially impact the NLHE in realistic materials. An analogous behavior is obtained for the minimal model in the type II regime (γ = 3), see Appendix B.

V. CONCLUSIONS

We have investigated the influence of surface states on the NLHE response of non-centrosymmetric timereversal invariant WSMs. Using both a model Hamiltonian and realistic first principles calculations, we have demonstrated that depending on the direction of the cut, surface states emerge that can substantially contribute to the NLHE of the slab. Notice that the topological nature of the surface state (topologically protected arcs, or track states) has no impact on the overall BCD, only the rela- [START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF], and we assign this behavior to the absence of mirror symmetry along ŷ in the original Hamiltonian given by Eq. ( 4).

tive number of states occupying the surface and the bulk matters. We emphasize that the relative contribution of the surface states with respect to the bulk states is in fact very large, leading to dramatic thickness-dependence of the NLHE response, in particular in type-II WSMs. This observation, confirmed by first principles calculations on WTe 2 slab geometries, suggests that surface states can contribute much more efficiently to NLHE than their bulk counterpart. We expect that the experimental signature of the surface-driven NLHE can be identified upon varying WSMs film thickness, depending on the growth orientation of the slab. These conclusions are not limited to WSMs and should apply to other topological systems and topologically trivial classes of non-centrosymmetric materials.

reported on Fig. 10. We observe a clear transition from Fermi arcs to Fermi pockets on the surface states (top and bottom panels). Indeed, a strong tilting can severely enlarge the Fermi pockets surrounding the Weyl nodes, leading to the merging between them and the generation of a trivial pocket. This is in fact what is happening at γ = 3 in Fig. 3(d). In other words, if one increases the tilting of the Weyl cones, the connectivity of the Fermi arcs is modified and can even induce the emergence of a trivial Fermi pocket from the merging of two Fermi pockets surrounding two opposite nodes.

FIG. 1 .

 1 FIG. 1. (Color online) Berry curvature of the Hamiltonian (4) in the (kx, kz) plane with ky = 0 and intrinsic parameters k0 = π 2 , m = 2, tx = 1 2 , γ = 1 and t = 1. Note that the Berry curvature of the two-band model is not sensitive to the value of γ [51]. The solid lines represent the Berry curvature vector field (Ωx, Ωz), whereas the blue (red) dots represent the local positive (negative) divergencies taking place at the Weyl nodes.

FIG. 2 .

 2 FIG. 2. (Color online) Band structures of the WSM slab of 25 layers, with γ = 1 (a-c) and γ = 3 (d-f) along the high symmetry path of the cubic lattice for different cuts: (left panels) x for ky = 0 and kz ∈ [-π, π], (center) ŷ along the path A1 = (-π, 0, -π) → Γ = (0, 0, 0) → A2 = (π, 0, π) and (right) ẑ for ky = 0 and kx ∈ [-π, π]. The central panels clearly displays type I (b) and type II (e) Weyl cones. The surface states at the top and bottom surfaces are represented with a solid red line. The flat red bands in (a), (c), (d) and (f) are surface states at zero energy that are highly localized at the corresponding surfaces.

L

  l=1 |ω lk ω lk |u nk , one can extract information about layer l by applying the projection operator S l = |ω lk ω lk |, l ∈ [1, L] to an observable O. In this way, it is easy to see that O = L l=1 S l O and L l=1 S l = 1 2L×2L . We use this description to separate the contributions of the density of states for a slab with a fixed number of layers,

  (b, d)], we obtain a different behavior. Whereas for γ = 1 [type I, Fig. 3(b)], the Weyl nodes remain disconnected throughout

FIG. 3 .

 3 FIG. 3. (Color online) Projected density of states (DOS) for (a, c) n = x and (b, d) n = ẑ as a function of the momentum coordinates for selected layers l for the type I (γ = 1) (a, b) and type II (γ = 3) phase (b, d). For the x-cut in the type I phase (a), the Weyl nodes dominate in the bulk and get connected by degenerate Fermi arcs at the surface, whereas in the type II phase (c), large trivial Fermi pockets appear across the slab. For the ẑ-cut in the type I phase (b), the Weyl nodes remain disconnected across the slab, whereas in the type II phase (d), trivial Fermi pockets dominate in the bulk with track states and projected Weyl nodes at the surface. The track states' projections are centered around kx = 0, ±π, which are indicated with blue arrows on each corresponding layer.

FIG. 4 .

 4 FIG. 4. (Color online) BCD components as a function of the tilting parameter γ for slab and bulk systems. (a) Slab (D x zx /L) and bulk (D 3D zx ) components for the x-cut, for different slab thicknesses. (b) Slab (D ẑ xz /L) and bulk (D 3D xz ) components for the ẑ-cut, for different slab thicknesses. (c, d) Corresponding ratio between the BCD coefficient and the corresponding density of states for bulk and slab systems for the x-cut (c) and ẑ-cut (d).

FIG. 5 .

 5 FIG. 5. (Color online) Layer dependent contribution for the BCD normalized by the number of layers (in this case, L = 25). We show the tomography for a type I WSM with γ = 1 in (a) and for a type II WSM with γ = 3 in (b).

FIG. 7 .

 7 FIG. 7. (Color online) Bulk band structure of WTe2 type II WSM obtained from density functional theory simulations. The inset displays the unit cell.

FIG. 8 .

 8 FIG. 8. (Color online) WTe2 density of states at Fermi level, projected on the bulk (a,c) and top surfaces (b,d) for a slab geometry containing 25 layers. Panels (a,b) correspond to a cut along ẑ whereas panels (c,d) correspond to a cut along x.

FIG. 9 .

 9 FIG. 9. (Color online) Band structure of a WTe2 slab containing 25 layers and cut along ẑ (a) and x (c). For this case the color bar represents the projections on the bottom (-1) and top (+1) layers. (b,d) Corresponding BCD for bulk and slabs containing 15, 20 and 25 layers.

FIG. 10 .

 10 FIG. 10. (Color online) Layer-resolved density of states of the topmost (top panels), central (middle panels) and bottommost (bottom panels) surface of a system with 25 layers and µ = 0.2. The red arrows indicate the Fermi arcs, and the blue arrows point to the center of the track states. Finally, indexes (a-d) show the cases γ = 1, 2, 2.5, 3, respectively, i.e., tracking the transition from type I to type II WSM. The mismatch between the top and bottom surfaces is in agreement to what is reported in[START_REF] Mccormick | Minimal models for topological weyl semimetals[END_REF], and we assign this behavior to the absence of mirror symmetry along ŷ in the original Hamiltonian given by Eq. (4).

ACKNOWLEDGMENTS D.G.O thanks to Andrés Saúl and Carlos Quezada for their help and fruitful discussions related to the computational implementation of the work. Besides, D.G.O and A.M acknowledge support from the Excellence Initiative of Aix-Marseille Université -A*Midex, a French "Investissements d'Avenir" program.

Appendix A: Fermi arcs reconnection upon tilting

When transiting from type I to type II WSM, increasing the tilting of the Weyl cones induces a reconnection of the Fermi arcs. This is clearly observed in the y-cut, as Appendix B: Dependence on the Fermi level All our calculations were performed away from the neutrality point, for µ = 0.2. One might wonder how the Fermi level impact the relative contribution of the surface states and the bulk states to the NLHE. In Fig. 11, we show the dependence of the BCD as a function of the Fermi level for the (a, c) x-and (b, d) z-cut, for type I (γ = 1) and type II (γ = 3) WSM. We find that the scenario proposed in the main text qualitatively holds: the x-cut, which possesses a strong surface contribution, exhibits a very strong difference between bulk BCD and slab BCD, regardless of the Fermi level, contrary to the z-cut whose surface states are much less prominent. This effect is similar to that observed in the realistic calcula-tions of WTe 2 , see Fig. 9.


Naif Hadadi

Adel Belayadi

Ahmed Alrabiah

Ousmane Ly

Collins Ashu Akosa

Michael Vogl

Hocine Bahlouli

Aurelien Manchon

Adel Abbout

Pseudo electric field and pumping valley current in graphene nanobubbles

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

An interesting and powerful property of graphene is its ability to be stretched elastically up to 25%, and to control the induced strains in different ways [START_REF] Masir | Pseudo magnetic field in strained graphene: Revisited[END_REF][START_REF] Guinea | Generating quantizing pseudomagnetic fields by bending graphene ribbons[END_REF]. In addition, the unique coupling between mechanical deformation and electronic structure along with the possibility of deformation makes graphene attracts considerable attention [START_REF] Settnes | Graphene nanobubbles as valley filters and beam splitters[END_REF][START_REF] Juan | Aharonov-bohm interferences from local deformations in graphene[END_REF]. This interplay was examined by observing the Landau levels that can form in graphene due to an induced strain [START_REF] Juan | Aharonov-bohm interferences from local deformations in graphene[END_REF]. Indeed, deformed graphene (Fig. 1) can generate an effective gauge field ⃗ A [START_REF] Juan | Aharonov-bohm interferences from local deformations in graphene[END_REF], to which one can associate a pseudo magnetic field (PMF), ⃗ B = ⃗ ∇ × ⃗ A. This stimulated pseudo magnetic field would allow electrons to behave as if they were subjected to a strong real magnetic field with strength a few orders of magnitude higher than the one generated by superconducting magnets. [START_REF] Guinea | Energy gaps and a zero-field quantum hall effect in graphene by strain engineering[END_REF].

Currently, there is a great interest in generating and controlling the valley degree of freedom of electrons in semiconductors. Indeed, the ability to manipulate valley electrons can potentially enable advanced valleyresolved electronic devices. In particular, this valley controllability opens up the possibility of using the momentum state of electrons, holes, or excitons as a completely new paradigm in information processing [START_REF] Vitale | Val-leytronics: opportunities, challenges, and paths forward[END_REF]. Moreover, pure valley currents are interestingly non-dissipative currents with no accompanying net charge flow, akin to pure spin currents [START_REF] Shimazaki | Generation and detection of pure valley current by electrically induced berry curvature in bilayer graphene[END_REF]. This property is very useful in seeking ultra-low power devices. Different routes have been proposed for generating valley currents such as using quantum pumps with a Dirac gap [START_REF] Wang | Pure valley current generation in graphene with a dirac gap by quantum pumping[END_REF], and using electrically induced Berry curvature in bilayer graphene [START_REF] Shimazaki | Generation and detection of pure valley current by electrically induced berry curvature in bilayer graphene[END_REF].

In the quantum pumping technique, a cyclic change of the AC voltage of the gates leads to the variation in the scattering matrix of the device and the generation of a DC current [START_REF] Wang | Pure valley current generation in graphene with a dirac gap by quantum pumping[END_REF]. This method of generating a charge DC current without bias voltage between two electrodes [START_REF] Wang | Pure valley current generation in graphene with a dirac gap by quantum pumping[END_REF][START_REF] Brouwer | Scattering approach to parametric pumping[END_REF] can be generalized to account for the spin or even the valley degree of freedom. In this work, a time-dependent nano-bubble in graphene, that can be created by a time-dependent voltage of an Atomic Force Microscopy (AFM) tip, capacitively coupled to the device [START_REF] Klimov | Electromechanical properties of graphene drumheads[END_REF], will be utilized as a quantum pumping device to examine the possibilities of generating a nonzero valley current with a zero net charge current. The deformation itself can be initially created by different techniques such as corrugated substrate engineering [START_REF]Strain superlattices and macroscale suspension of graphene induced by corrugated substrates[END_REF][START_REF] Viola Kusminskiy | Pinning of a two-dimensional membrane on top of a patterned substrate: The case of graphene[END_REF] or gas inflation [START_REF] Bunch | Impermeable atomic membranes from graphene sheets[END_REF]. The pseudo magnetic field associated with the deformation in graphene is valley dependent and changes its sign between the two valley points K and K'.This feature allows the PMF to preserve time-reversal symmetry unlike a real external magnetic field. Consequently, a time manipulation of the PMF will induce different behaviour for the K and K' electrons and thus a valley-resolved study of the current is necessary [START_REF] Nedell | Deep learning of deformation-dependent conductance in thin films: Nanobubbles in graphene[END_REF]. It is worth mentioning that different studies have looked at the effect of other dynamical strains on the electronic transport. Ref. [START_REF] Vaezi | Topological electric current from time-dependent elastic deformations in graphene[END_REF] for example, considered in-plane strain in a gapped graphene to induce a topological current, Ref. [START_REF] Jiang | Generation of pure bulk valley current in graphene[END_REF] looked at the symmetry requirements for pumping valley current, while Ref. [START_REF] Sela | Quantum hall response to time-dependent strain gradients in graphene[END_REF] investigated the pseudo quantum Hall regime due to strain. Other shapes or mechanisms to get the valley filtering were also investigated [START_REF] Amasay | Transport through dynamic pseudogauge fields and snake states in a corbino geometry[END_REF][START_REF] Belayadi | Valley-dependent tunneling through electrostatically created quantum dots in heterostructures of graphene with hexagonal boron nitride[END_REF] 2. MODEL Fig. 2. Tight-binding model for a quasi-1D graphene system with zigzag edges. a) shows the local deformation as a local bump in the center of the graphene sheet, altering the translational invariance of the quasi 1D waveguide. b) shows a face side of the Gaussian strain bump which is centered in the middle of the system. The green links in the zoom of c) represent the interface through which the pumped current is calculated. The red parts in c) are the leads (reservoirs) that indicate that the system is infinite in that direction. An animation sketching the system can be found in [START_REF]Supplemental material on website[END_REF] Graphene consists of carbon atoms, that are arranged in an hexagonal honeycomb structure. This arrangement is a triangular lattice with a basis of two atoms per unit cell where the lattice vector for this structure are ⃗ a 1 = a 2 (3,

√

3) and ⃗ a 2 = a 2 (3, -√ 3) (Fig. 2) and the constant a is carbon-carbon distance with value d 0 ∼ 1.42 Å [START_REF] Carrillo-Bastos | Gaussian deformations in graphene ribbons: Flowers and confinement[END_REF]. The electronic structure of strained graphene is studied by utilizing a first nearest neighbor tight-binding model approximation, H = -∑ ⟨i j⟩ t i j c † i c j where the sum runs over the nearest neighboring atoms. Figure 2 shows the semi-infinite graphene system, where the plane waves are coming along the zigzag direction from left and right leads. Experimentally, the deformation can be created in a graphene sheet by using an atomic force microscope's tip [START_REF] Klimov | Electromechanical properties of graphene drumheads[END_REF][START_REF] Lee | Measurement of the elastic properties and intrinsic strength of monolayer graphene[END_REF][START_REF] Nemes-Incze | Preparing local strain patterns in graphene by atomic force microscope based indentation[END_REF]. In the presence of a nano-bubble the nearest neighbors hopping parameters become all different. Thus, the strain can be included in the system by altering the hopping parameter to the strained bond length. The adjustment of the hopping in the presence of the strain can be described by

t i j = t 0 e -β ( d i j d 0 -1) , ( 1
)
where t 0 = 2.7 eV, β = 3.37 for graphene structure [START_REF] Stegmann | Current splitting and valley polarization in elastically deformed graphene[END_REF]. d i j represents the length of the strained lattice bonds whereas d 0 is the lattice constant in the absence of strain.

3. RESULTS

Pseudo Magnetic Field

In our model, the graphene out-of-plane deformation is local with a characteristic width σ = 5 nm, a maximum height h 0 = 3.5 nm [START_REF] Settnes | Graphene nanobubbles as valley filters and beam splitters[END_REF][START_REF] Juan | Aharonov-bohm interferences from local deformations in graphene[END_REF][START_REF] Pengfei | Programmable graphene nanobubbles with three-fold symmetric pseudo-magnetic fields[END_REF] and an overall Gaussian shape modeled by z = h 0 e -x 2 +y 2 2σ 2 . In the limit of low-energy carriers, the deformation is equivalent to an induced effective vector potential given by [START_REF] Masir | Pseudo magnetic field in strained graphene: Revisited[END_REF] ⃗ A = -hβ 2ea 0 ((ε xxε yy ) x -2ε xy ŷ), where the stress tensors ε i j can be calculated by

ε i j = 1 2 [(∂ i z)(∂ j z)].
A straightforward calculation leads to the following form of the vector potential,

⃗ A 0 = -hβ ea 0 σ 4 z 2 ( x 2 -y 2 x -2xy ŷ). ( 2
)
The corresponding pseudo magnetic field can be obtained by

⃗ B 0 = ⃗ ∇ × ⃗ A 0 : ⃗ B 0 = 4hβ ea 0 σ 6 z 2 (y 3 -3x 2 y)ẑ. (3)
The pseudo magnetic field shows alternating positive and negative regions in the deformation position as shown in Fig. 3. This magnetic field is valley dependent and changes its sign between the high symmetry points K and K ′ and is responsible for the valley filtering in graphene nano-bubbles [START_REF] Settnes | Graphene nanobubbles as valley filters and beam splitters[END_REF]. Using eq.1, we can demonstrate that the maximum relative change in the hopping parameter is proportional to β h 2 0 /σ 2 . This implies that for a small deformation width σ , it becomes necessary to reduce the maximum height to remain within the constraints of stretchable graphene [START_REF] Masir | Pseudo magnetic field in strained graphene: Revisited[END_REF][START_REF] Guinea | Generating quantizing pseudomagnetic fields by bending graphene ribbons[END_REF]. In our particular case, a numerical examination of max(∆t i j /t i j ) confirms that we are, indeed, within this specified limit.

Induced Pseudo Electric Field

The dynamics of the charge carriers in graphene under a dynamical external mechanical strain is a challenging problem. The generated frequencies of the perturbation projects into the sub-Terahertz domain with reported frequencies up to 200 GHz [START_REF]200 ghz maximum oscillation frequency in cvd graphene radio frequency transistors[END_REF]. We note in passing that such frequencies are much smaller than the bandwidth of graphene. Therefore, standard techniques like Floquet theory [START_REF] Eckardt | High-frequency approximation for periodically driven quantum systems from a floquet-space perspective[END_REF][START_REF] Bukov | Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering[END_REF][START_REF] Vogl | Light driven magnetic transitions in transition metal dichalcogenide heterobilayers[END_REF][START_REF] Abanin | Effective hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems[END_REF][START_REF] Rodriguez-Vega | Lowfrequency and moiré-floquet engineering: A review[END_REF][START_REF] Vogl | Analog of hamilton-jacobi theory for the time-evolution operator[END_REF][START_REF] Vogl | Effective floquet hamiltonian in the low-frequency regime[END_REF] are hard to apply, which makes a full numerical treatment necessary and which will be the subject of our work. Things become more interesting when the strain is a nano-bubble exhibiting a few hundred Teslas pseudo magnetic field. Indeed, an oscillating bump with a height h = h 0 + δ h(t) will change the vector potential to a new form

⃗ A = (1 + δ h(t) h 0 ) 2 ⃗ A 0 .
In the Weyl gauge, the Maxwell equations show an induced pseudo electric field in the limit δ h(t) h 0 ≪ 1, to the lowest order:

⃗ E(x, y) = -2 δ . h(t) h 0 ⃗ A 0 , (4)
where the dot represents the time derivative. This result accounts for the small perturbations due to a small amplitude in the out-of-plane vibration of the nano-bubble which allows us to keep our interpretation within the linear response theory. It is worth noting that although the pseudo magnetic field is strong (hundreds of Teslas), the induced electric field is only a few kV/cm. This is mainly explained by the fact that these large values of B 0 are rather due to the large gradient of ⃗ A 0 (variation on very short scales) than to its magnitude. The pseudo-electric field map is depicted in Fig. 3 where a threefold symmetry emerges and the system alternates between an inward and outward pseudo electric field for the K and K ′ valleys.

Fig. 4. The generated currents through the interface shown in Fig. 2. The two types of currents at K and K ′ are in phase and average to zero. The charge current is obtained for h 0 = 3.5 nm and δ h = 0.35 nm. E F = 0.28t 0 . T is the period of oscillation of the bubble.

Charge and Valley Currents

As, in the case of charge-pumping by the dynamics of a skyrmion deformation [START_REF] Brouwer | Scattering approach to parametric pumping[END_REF][START_REF] Abbout | Cooperative charge pumping and enhanced skyrmion mobility[END_REF] and spin-pumping by the precession of magnetic textures [START_REF] Tserkovnyak | Spin pumping and magnetization dynamics in metallic multilayers[END_REF][START_REF] Tatara | Consistent microscopic analysis of spin pumping effects[END_REF], an oscillating nano-bubble is expected to alter the charge density and induce a charge/valley current through the interface between the system and the lead [37]. To investigate this property, we adopt a nano-bubble oscillating at hundreds GHz with δ h(t) = δ h 0 sin(ω 0 t + φ ) with a height modulation strength h 0 = 3.5 nm, a characteristic width σ = 5 nm and a perturbation magnitude δ h 0 = 0.35 nm. The way to generate these oscillations is unimportant at this stage but can be initiated by a capacitively coupled atomic force microscope's tip or a laser on top of the strained part. The first step towards the calculation of the charge current is to solve the time-dependent Schrodinger equation to obtain the wave functions at different times and then calculate the different currents in the system in one-body formalism [START_REF] Kloss | Tkwant: a software package for timedependent quantum transport[END_REF][START_REF] Groth | Kwant: a software package for quantum transport[END_REF]. It is important to mention that a more accurate procedure involves an integration over the different contributions from all energies in the Fermi sea. This was safely avoided due to the small driving frequencies (compared to the Fermi level E F = 270 meV) and the small excitations considered in this study (hω/E F ≪ 1). This reduces the needed Fig. 5. The charge current (Left) and the valley current (Right) pumped through an interface with the lead (see figure 2). The average of each of them is zero. The charge oscillates at the same frequency as the nano-bubble, whereas the pumped valley current shows higher order harmonics.The blue dots represent the times at which the current maps were plotted in Fig. 7. T is the period of oscillation of the nanobubble numerical resources and keeps the problem at the Fermi surface (similar to the case of the spin-pumping theory). Figure 4 shows the result of the pumped charge current I charge = I K + I K ′ , which expectedly averages to zero as explained by the scattering approach to parametric pumping. Indeed, varying one parameter is not enough to generate a sustainable charge current [START_REF] Brouwer | Scattering approach to parametric pumping[END_REF].

The frequency of the signal is the same as that of the oscillating nano-bubble. So far, these results look like the same as for the spin pumping of I x , I y and I charge in Ferromagnets (precessing around the z-axis): zero average and the same frequency as the driving perturbation. Since the pseudo electric field is valley dependent, the contribution coming from each K-point to the charge current is investigated and plotted in Fig. 4. The two currents, I K , and I K ′ have amplitudes of the same order of magnitude and both oscillate at the same angular frequency ω 0 . The two currents are not exactly the same, because for a Fermi energy slightly higher than zero, the modes in the two different valleys are non-symmetric and have slightly different wavenumbers. The valley current, which is defined as I valley = I K -I K ′ is of huge interest in low energy consumption devices. The reason behind it is that the net charge transferred can be zero whereas at the same time the valley current is not vanishing. This eliminates the dissipation of energy as heat while transporting information through the valleys' degrees of freedom. The valley current is well defined in the interface region (away from the bubble), because the lead is uniform and the conducting modes (incoming wave functions from the lead) do not mix. One has to compute the time-dependent wave function for each mode, using T-kwant package [START_REF] Kloss | Tkwant: a software package for timedependent quantum transport[END_REF][START_REF] Groth | Kwant: a software package for quantum transport[END_REF], and identify the different modes in each valley. The contribution of each mode to the current between sites a and b reads:

I ab = i ψ † b (H ab ) † ψ a -ψ † a H ab ψ b (5)
where, ψ a(b) is the wave function of a given mode at site a(b). H ab is the matrix element between sites a and b (in our case, it represents the hopping parameter between a and b in graphene). The contribution of modes of a given valley are summed together to finally express I valley . Inside the scattering region, this procedure can still be used, but the current is not exactly a pure valley current. Indeed, the scattering between the valleys is possible, yet it remains small (the points K and K' are far apart in the Brillouin zone) [START_REF] Nedell | Deep learning of deformation-dependent conductance in thin films: Nanobubbles in graphene[END_REF]. For this reason, the approximation of the local current in the scattering region is still valid [START_REF] Nedell | Deep learning of deformation-dependent conductance in thin films: Nanobubbles in graphene[END_REF]. Another approach is to project the wave functions on the valley-resolved ones of pristine graphene [START_REF] Stegmann | Current splitting and valley polarization in elastically deformed graphene[END_REF]. Fig. 5 shows the result for the valley current pumped through the same interface. The average is again zero and the amplitude is much smaller than that of the charge current. The most important thing to notice is that the frequency of the I valley is not exactly the same as that of the nano-bubble. As we can see in Fig. 7, the valley current exhibits a main frequency of 3ω 0 with small contributions from ω 0 , 2ω 0 and 5ω 0 . The appearance of higher-order harmonics was recently reported for spinpumping in the presence of spin-orbit interaction [START_REF] Ly | Spin-orbit coupling induced ultrahigh-harmonic generation from magnetic dynamics[END_REF] or in the presence of non-collinear magnetic structure [START_REF] Ly | Noncollinear antiferromagnetic textures driven high-harmonic generation from magnetic dynamics in the absence of spin-orbit coupling[END_REF]. In our case, none of the discussed causes in [START_REF] Ly | Spin-orbit coupling induced ultrahigh-harmonic generation from magnetic dynamics[END_REF][START_REF] Ly | Noncollinear antiferromagnetic textures driven high-harmonic generation from magnetic dynamics in the absence of spin-orbit coupling[END_REF] are present in our system. The charge and valley pumping depend on the parameters of the system. As we can see it from the expression of the pseudo electric field, the latter is proportional to the perturbation δ h o h 0 . Naturally, this is expected to be the same for the charge and valley currents. Fig. 7 b) shows the amplitude of the I charge and I valley which are clearly linearly varying with the change δ h o h 0 . For stronger perturbations, higher non-linear effects are expected to manifest. Another way to increase the different pumped currents is to play with the frequency and use a Terahertz spectrum. This is a challenging task, but very promising results suggest that this projection is reachable [START_REF]200 ghz maximum oscillation frequency in cvd graphene radio frequency transistors[END_REF]. The zero-average of the valley current makes it very difficult to exploit in devices since no accumulation can be obtained and only oscillating current redistribution of the valleys is achieved. With proper engineering, we still can think of it as the basis for a potentially new Field-effect transistor [START_REF] Choi | Ferromagnet-free all-electric spin hall transistors[END_REF].

Local currents

Since the pseudo electric field is local and mainly appears in the central system, investigating the current at different regions than the interface with the lead might be of a great interest. The plot of the current map will show how the regions of the maximum pseudo magnetic field influence the current direction and generation. The wave functions are calculated at all sites of the systems which allows us to express the local currents. Fig. 6 illustrates the map obtained at different times during one period. The current lines for the charge and the valley are symmetric which is just the reflection of the symmetry of the Gaussian bump considered in this work. This will be different if a triaxial strain was adopted. In fact, it is known that a Gaussian deformation filters the valleys whereas the triaxial one splits them [START_REF] Settnes | Graphene nanobubbles as valley filters and beam splitters[END_REF].The Fermi energy is chosen low to reduce the number of modes due to the limited computing resources in this kind of 2D maps (current is calculated between each two links, at each time) and this doesn't change the general conclusions. The maps for the valley current and the charge current are intrinsically different: for the charge current, it is hard to see the positions of the valley-dependent maxima of the pseudo-magnetic field. The situation is different for the valley current. We clearly notice loops of different orientations showing the dependence of the pseudomagnetic field on the valley of the carriers. The size of the current's loops can be changed via the parameter σ ie, creating very peaked bumps or smoothly varying ones over a larger region. The number of modes for the chosen Fermi energy is very small. In order to increase the current, one needs to increase it either by changing the width of the system or by increasing the Fermi energy. It's worth mentioning that we did investigate deformations with elliptical symmetry, deliberately tilting them to disrupt left-right symmetry. Unfortunately, this approach did not lead to a sustained average valley current.

CONCLUSION

The pseudo magnetic field generated by a strained graphene sheet was analyzed. In the presence of a timedependent oscillating strain, a corresponding pseudo electric field is stimulated giving rise to a redistribution of the quantum states in momentum space, leading to charge and valley pumping. We demonstrated that the valley current exhibits extra-harmonics and that on average, the pseudo electric field did not sustain a direct current, and rather a zero average was obtained. This work opens the door to further interesting studies like the pseudo-electric field generated by the propagation of mechanical waves in graphene or the distribution of local electric currents due to phonons and temperature effects on graphene.

Fig. 1 .

 1 Fig. 1. A mechanical deformation of a graphene sheet. The bump which can have different shapes (spherical, Gaussian, ...) induces an out-of-plane stretching of the hoppings between neighboring atoms.

Fig. 3 .

 3 Fig. 3. The color map indicates the threefold symmetric pseudo-magnetic field caused by the circularly symmetric Gaussian deformation. The vector field illustrates the stimulated electric field by the time-dependent pseudo-magnetic field. 'a' is the lattice constant.

Fig. 6 .Fig. 7 .

 67 Fig.6. The current map of the pumped charge current (top) and the valley pumped current at different times for 6 modes with E F = 0.07t 0 . The nano-bubble size is σ = 5 nm and the height δ h 0 = 0.1h 0 . The position is expressed in nm.The times at which the maps are plotted are inside one period of oscillation and increase from left to right. (the 4 different phases of an oscillation). The blue dots in Fig.5show the time where the maps are plotted. ( for these figures, the time increases from left to right.)

ACKNOWLEDGMENT

The authors acknowledge computing time on the supercomputer SHAHEEN at KAUST Supercomputing Centre and the team assistance. A.A. and M.V. gratefully acknowledge the support provided by the Deanship of Research Oversight and Coordination (DROC) at King Fahd University of Petroleum and Minerals (KFUPM) for funding this work through exploratory research grant No. ER221002.