Giant c-axis nonlinear anomalous Hall effect in Td-MoTe2 and WTe2

While the anomalous Hall effect can manifest even without an external magnetic field, time reversal symmetry is nonetheless still broken by the internal magnetization of the sample. Recently, it has been shown that certain materials without an inversion center allow for a nonlinear type of anomalous Hall effect whilst retaining time reversal symmetry. The effect may arise from either Berry curvature or through various asymmetric scattering mechanisms. Here, we report the observation of an extremely large c-axis nonlinear anomalous Hall effect in the non-centrosymmetric Td phase of MoTe2 and WTe2 without intrinsic magnetic order. We find that the effect is dominated by skew-scattering at higher temperatures combined with another scattering process active at low temperatures. Application of higher bias yields an extremely large Hall ratio of E⊥/E|| = 2.47 and corresponding anomalous Hall conductivity of order 8 × 107 S/m.

I n ferromagnetic conductors, passage of charge current can generate an electric field in the transverse direction even without the application of an external magnetic field. This anomalous Hall effect (AHE) requires broken time reversal symmetry and originates both from topological aspects of the material's band structure and electron scattering coupled to the spin-orbit interaction 1 . While parsing the different contributions of the AHE has been a grand challenge in condensed matter physics for many decades, the overall strength of the AHE can be clearly quantified by the Hall ratio. Defined as the ratio between the transverse and longitudinal conductivities or fields, it ranges between~10 −2 in common magnets and~10 −1 in more exotic, giant AHE compounds 2, 3 .
In non-centrosymmetric materials retaining time reversal symmetry, it is possible to realize a nonlinear AHE (NLAHE), whereby the current induces an effective magnetization in the sample and establishes a transverse electric field that increases quadratically with applied electric field 4-12 . Experimentally, an AC voltage of frequency ω is applied across the sample and a voltage of either frequency 2ω or 0 (DC) is detected in the transverse direction. Such an effect was first demonstrated in twodimensional (2D) WTe 2 within the plane of the layers 13,14 , as well as in several subsequent material systems [15][16][17][18] . In contrast to the ordinary AHE, the strength of the NLAHE can be quantified by the transverse field relative to the square of the applied longitudinal field. This value was measured to be 1.4 × 10 −9 m/V for 2D WTe 2 at low temperature 14 .
Here, we report the observation of a new thickness-dependent, c-axis NLAHE in the same transition metal dichalcogenide family, where application of an in-plane current generates a Hall field perpendicular to the layers 6 . We focus mainly on T d -MoTe 2 , although we also verify that WTe 2 exhibits qualitatively similar behavior. The c-axis NLAHE strength can be as large as 0.43 m/V. Due to the quadratic increase of the transverse field with increasing current, a Hall ratio of E ⊥ /E || = 2.47 can be achieved at higher bias, the largest yet reported without edge states. The corresponding anomalous Hall (transverse) conductivity, 8 10 7 S=m, is furthermore the largest in any material. These effects arise primarily from skew-scattering at higher temperatures combined with another scattering process at lower temperatures that increases the NLAHE strength exponentially with linear longitudinal conductivity. By performing measurements across different temperatures and sample thicknesses down to the 2D limit, we furthermore uncover a general scaling behavior of the NLAHE strength with conductivity, which we parse and analyze in the two temperature regimes.

Results
Nonlinear susceptibility tensor and measurement geometry. The left side of Fig. 1a shows the crystal structure of T d -MoTe 2 . The a and b axes lie within the plane of the layers, while the c-axis points out-of-plane. In addition to a mirror plane that is normal to the a-axis, there is a glide plane that is normal to the b-axis and a two-fold screw axis along the c-axis, putting the bulk crystal in the non-centrosymmetric space group Pmn2 1 19 . The corresponding point group mm2 has the nonlinear susceptibility tensor: along the c-axis, while nonlinear currents in-plane strictly vanish: In contrast, since the glide plane and screw axis symmetries are broken for the surface layers, ultrathin samples understood to be in space group Pm 21 , with allowed nonlinear currents both in-plane and out-of-plane.
To measure this c-axis NLAHE, we have fabricated samples with vertical as well as in-plane contacts. In particular, our devices consist of underlying gold electrodes with MoTe 2 crystals of various thicknesses transferred on top. To maintain consistency, all MoTe 2 flakes were exfoliated from a single piece of bulk crystal grown by the flux method (see Methods section). Fewlayer graphene was used as a vertical top electrode with insulating hexagonal boron nitride (h-BN) blocking all but the tip of the vertical contacts to prevent mixing with in-plane currents and fields. In order to protect the sample from degradation, the entire transfer process was performed within a nitrogen-filled glovebox, while a final layer of h-BN was used to cover the whole structure. A device schematic is shown (without the top h-BN layer for clarity) on the right of Fig. 1a and an optical image of a representative device is shown in the inset.
Since current applied along the a-and b-axis of the crystal generally give rise to different nonlinear Hall currents or fields as per the different elements of the susceptibility tensor d 31 and d 32 , we have controlled for the sample orientation by selecting MoTe 2 flakes that were rectangular in shape and aligning them with the circular electrode pattern. This allows for current injection along both the long and short directions of the flake in the same device, which in almost all cases correspond to the a-and b-axis, respectively. For either current I||a or I||b, our device geometry yields measurement of the following linear and nonlinear voltages: in-plane longitudinal (V xx ), in-plane Hall (V xy ), and out-of-plane Hall (V xz ).
Determination of crystalline axes. In order to confirm the crystal orientation of the flakes, we first measured the longitudinal magnetoresistance, ð Þ for constant current, along the two current directions. Fig. 1b shows representative data from a 47-nm-thick device at 1.4 K: I||a (I||b) exhibits quadratic (linear) magnetoresistance up to 12.5 T, consistent with previous measurements on bulk crystals 22 . Overall, MR is large and non-saturating due to close electron-hole compensation combined with relatively high carrier mobilities [23][24][25] . In addition, we have performed optical second harmonic generation (SHG) measurements, which can also distinguish between the a-and b-axis of MoTe 2 . Fig. 1c shows SHG intensity measured on the same sample at 80 K with incident and scattered light polarization in parallel as a function of the polarization angle relative to the a-axis, at which a node can be seen, consistent with previous results 26 . Unlike the electrical second harmonic measurements we are to present, optical SHG is mostly sensitive to the symmetry properties of the surface layers, which allow for inplane responses (see Supplementary Note 1).
Observation of c-axis nonlinear anomalous hall effect. Fig. 2a shows representative NLAHE measurements taken on a 127-nmthick sample at 2 K. The upper two panels show V xz measured at DC and second harmonic (2ω ¼ 354Hz) frequencies as a function of the first harmonic (ω ¼ 177Hz) V 2 xx for I||a and I||b. All four traces show linear behavior at low bias as expected for the NLAHE 14 , although the slope is larger for current along the aaxis for reasons which we shall discuss. At higher bias the curves become slightly sublinear, possibly due to sample heating. Furthermore, for each current direction, the DC and 2ω amplitudes are comparable for a given applied bias. In contrast, as shown in the lower two panels, relatively small second harmonic voltage is measured for V xx or V xy for comparable bias, consistent with the allowed symmetries discussed earlier for bulk-like crystals. Although, V 2ω xy for I||a is in principle allowed for thin or surface  For the remainder of this work, we shall focus only on the second harmonic component of the c-axis NLAHE. The strength of this effect can be defined as the slope of the linearized plots presented at low bias. In order to understand and demonstrate the reproducibility of these results, we have performed similar measurements across five samples of different thicknesses down to the 2D limit: 127, 70, 47, 32, and 9 nm. The NLAHE strength is shown as a function of sample thickness for both I||a and I||b in the upper panel of Fig. 2b. We have used electric field values, E 2ω z = E ω x À Á 2 , instead of voltage in order to account for the different dimensions of the samples. The local E ω x at the vertical contact area was obtained from finite element method simulations of each individual device (see Methods and Supplementary Note 4). For every sample, the a-axis shows larger strength than the b-axis, while for both axes, the strength decreases substantially with decreasing thickness. In particular, for I||a in the thickest (127 Scaling with sample conductivity. The differences observed between the samples and two crystalline axes can be attributed to their differences in conductivity. In the lower panel of Fig. 2b, we have plotted the residual longitudinal conductivity σ xx0 for each sample and axis, which shows a similar trend. This reflects the three-dimensional electronic character of MoTe 2 -decreasing thickness increases surface scattering, which lowers the residual conductivity by over an order of magnitude in our thinnest sample, consistent with previous results 25 . In Fig. 3a, we further show the temperature dependence of E 2ω z = E ω x À Á 2 and conductivity measured for the 127-nm-thick sample-both decrease with increasing temperature. In the lower inset, we show the temperature-dependent longitudinal resistivity for I||a in the same sample, which conversely decreases with decreasing temperature. Specifically, the resistivity scales linearly with temperature at higher temperatures, but begins to saturate in the crossover region roughly centered at T L $ 18K . Taken together, these results indicate that conductivity (or resistivity) may be the fundamental parameter upon which the NLAHE strength depends. To show this dependence explicitly, in the main panel of Fig. 3b, we have plotted, in log-log scale, E 2ω z = E ω x À Á 2 for both axes of each sample versus temperaturedependent conductivity normalized to the residual conductivity, σ xx =σ xx0 , and a general trend appears to emerge. Specifically, at lower conductivity (higher temperature) for a given sample and axis, E 2ω z = E ω x À Á 2 scales closely with σ 2 xx as compared with the guide-to-eye in gray. Above (below) a certain conductivity (temperature), it increases at an even faster rate. To see this crossover point more clearly, we have extracted the local slope υ from the data in the main panel as a function of temperature T (normalized to the T L of each sample/axis, which ranges between 13-18 K), yielding the relationship E 2ω The inset of Fig. 3b shows υðTÞ for all the different traces in the same color scheme as that used in the main panel. Above T L , υ is near 2 for every trace, while it rises continuously below this temperature. We have also performed similar measurements on bulk-like WTe 2 (which shares the same crystal structure) as well as MoTe 2 Hall bar devices for I||a (see Supplementary Note 5 and 6, respectively), both of which show qualitatively similar behavior. This indicates that the strength and scaling of the c-axis NLAHE is not unique to MoTe 2 nor a result of the circular electrode geometry.

Discussion
In order to understand the mechanism behind the NLAHE and the observed scaling with conductivity, we turn to the theory behind the ordinary AHE, as it is understood that the NLAHE strength, which can also be written as (σ zz being the linear c-axis conductivity and σ ðNLÞ AH is the nonlinear anomalous Hall conductivity element zx), scales in the same way as anomalous Hall conductivity σ AH in the ordinary (linear) AHE 9,11,14 . For the latter, σ AH consists of intrinsic Berry curvature, skewscattering, and side-jump(-like) scattering contributions 1 . The intrinsic part does not depend on longitudinal conductivity (scales with σ 0 xx ), while upon changing temperature, skewscattering scales with σ 2 xx and side-jump contains terms that scale with σ 0 xx , σ 1 xx , and σ 2 xx 9,11,27,28 . It is therefore not possible to distinguish between the various mechanisms by scaling analysis a priori. Nonetheless, it is generally understood that skewscattering dominates for highly conducting samples with large residual conductivity 1 . Our MoTe 2 samples show a conductivity comparable to the most highly conducting AHE systems, and so we shall assume that the υ ¼ 2 scaling at higher temperatures is attributed primarily to skew-scattering.
To further our understanding, we have fit all our data for E 2ω z = E ω x À Á 2 above T L to the functional form: Aσ 2 xx þ B, and the extracted A and B values are plotted in the lower panels of Fig. 4a as a function of the residual conductivity σ xx0 of the sample/axis. The upper panel shows a representative fitting for the 127-nmthick sample with I||a. The constant term B is expected to contain contributions both from the intrinsic Berry curvature and extrinsic side-jump(-like) events and be independent of residual conductivity 11,27,28 . Although this value varies over two orders of magnitude between samples, we indeed do not observe any clear trend with σ xx0 . We have further calculated the upper limit for the predicted intrinsic contribution to B based on previous theory (see Supplementary Note 7) 6 , and the result is marked by the gray line. This theoretical value is one order of magnitude less than our smallest experimental value, suggesting that intrinsic Berry curvature only plays a minor role in the large NLAHE observed here. Although, the intrinsic contribution may in principle be enhanced by external perturbations, such as strain, or symmetrybreaking by the surface layers [29][30][31] .
On the other hand, A initially decreases with increasing σ xx0 , but saturates at higher σ xx0 . For the AHE, the skew-scattering term is given by A ¼ ασ À1 xx0 9,11,27,28 , where α is termed the skewscattering coefficient and the residual conductivity is taken as a measure of the disorder from impurities. This general scaling cannot be directly applied for our samples, since we are not tuning the impurity concentration by changing thickness, while measurements along the different crystalline axes for the same sample show different residual conductivity due to differences in band structure, not disorder. By reducing thickness, we have instead increased the contribution of surface scattering to the total resistivity, which by Matthiessen's rule can be written as: xx0sðurfaceÞ . Taking into account the different symmetries for the bulk and surface layers, we further assign different skew-scattering coefficients between the two and write: In the thick limit, we expect σ xx0 $ σ xx0b , and so A $ α b σ À1 xx0b , consistent with the nonzero constant value we observe at high conductivity. If we further approximate that σ xx0b is independent of thickness and equal to the measured σ xx0 of our thickest sample, we can fit our data across the entire thickness (residual conductivity) range for each axis to the formula above (see colored lines in the middle panel of Fig. 4a). The extracted skew-scattering coefficients are shown in the inset. Overall, α is larger for bulk skew-scattering and similar between the two axes.
For temperatures below T L , E 2ω z = E ω x À Á 2 increases further with conductivity beyond the σ 2 xx dependence seen at higher temperatures (see Fig. 3b). Such behavior has not been previously observed experimentally for the AHE, and so it remains an open theoretical question as to the microscopic scattering mechanism responsible for this effect. A possible candidate is side-jump, as it has been recently predicted that, for high-purity samples below the Debye or Bloch-Gruneisen temperature, scaling parameters for side-jump(-like) contributions may actually be temperaturedependent 11 . We proceed to fit our data in this regime in order to determine an empirical scaling relationship. The continuous increase of υ below T L suggests an exponential dependence on σ xx , and so we have fit the additional contribution to: I||a. While C sharply with increasing σ xx0 , D decreases. The net effect, however, is such that the NLAHE strength at low temperature is increased by a roughly constant factor above that expected from skew-scattering alone. In the lower panel of Fig. 4b, we have plotted as a function of σ xx0 , the ratio between the measured E 2ω z = E ω x À Á 2 value at 2 K and Aσ 2 xx0 þ B, the extrapolated skew-contribution, and this factor falls between 3 and 5 for all but one sample/axis.
The combination of the different scattering processes, together with the nonlinear dependence of the Hall field on applied bias, allows us to obtain the largest Hall ratio in electric field observed to date outside the quantum (anomalous) Hall regime. In Fig. 5a, we have plotted E 2ω z =E ω x , the ratio between the Hall field generated and the applied longitudinal field, as a function of E ω x measured for the 127-nm-thick sample at 2 K, which shows the largest NLAHE strength amongst all our samples. The ratio is always larger for I||a, for which it rises initially with increasing E ω x , reaching a maximum of 2.47, an order of magnitude larger than the analogous quantity measured in giant linear AHE compounds 2,3 . At larger bias, E 2ω z =E ω x decreases as E ω x itself becomes sublinear with applied current (see Supplementary Note 8), likely due to heating.
We would further like to quantify the nonlinear anomalous Hall conductivity at this peak Hall ratio, which can be expressed as: σ x σ zz . We have measured both the a-and c-axis resistivity of a separate bulk MoTe 2 crystal (grown under the same conditions as that used for our devices) and further calculated the conductivity anisotropy as a function of the position of the Fermi level using density functional theory (see Supplementary Note 9). σ zz =σ xx $ 0:25ðexp:Þ À 0:6ðth:Þ for I||a near the charge neutrality point where electrons and holes compensate, yielding σ ðNLÞ AH $ 8 10 7 S=m. Assuming a similar anisotropy ratio for WTe 2 , we estimate σ ðNLÞ AH $ 5 10 7 S=m at the peak Hall ratio for our bulk-like WTe 2 device (see Supplementary Note 5). For comparison, we have plotted these values in Fig. 5b, together with those measured in-plane for ultrathin MoTe 2 and WTe 2 , as well as the linear anomalous Hall conductivities measured in various magnetic systems as a function of the residual longitudinal conductivity 2,3,14,32-43 . Overall, σ AH is larger in more highly conducting systems. The values for the c-axis response in MoTe 2 and WTe 2 , in particular, are the largest yet reported.
In conclusion, we have observed an extremely large c-axis NLAHE in T d -MoTe 2 and WTe 2 that is dominated by asymmetric electron scattering. Our work reaffirms the importance of extrinsic contributions to the (NL)AHE, especially in highly conducting metals, and opens a new direction for obtaining giant Hall ratios and conductivities in non-centrosymmetric systems without breaking time reversal symmetry.

Methods
Crystal synthesis. 1T′-MoTe 2 (T d -WTe 2 ) single crystals were grown by the flux method using Te as a solvent. Mo (W) [Alfa Aesar, 99.9%], Te [Alfa Aesar, 99.99%] powders were ground and placed into alumina crucibles in a ratio of 1:25 (1:99) and sealed in a quartz ampoule. After the quartz ampoule was heated to 1050°C (1000°C) and held for 2 days (10 h), the ampoule was slowly cooled to 900°C (460°C ) over 120 h (100 h) and centrifuged. Shiny and plate-like crystals with lateral dimensions up to several millimeters were obtained.
Device fabrication. Au (45 nm)/Ti (5 nm) electrodes in a circular geometry were pre-patterned on Si wafers with 285-nm-thick SiO 2 using conventional photolithography and electron-beam deposition. Graphene (Coorstek), h-BN (HQ graphene), and (Mo/W)Te 2 flakes were exfoliated onto blank SiO 2 /Si wafers. After the desired flakes were identified, a polymer stamp coated with polycarbonate was used to sequentially pick up the entire h-BN/graphene/h-BN/(Mo/W)Te 2 /h-BN heterostructure to avoid contamination between the layers. The heterostructure was then aligned and transferred onto the pre-patterned electrodes. The entire exfoliation and transfer process were performed within a nitrogen-filled glovebox to avoid degradation of (Mo/W)Te 2 in air.
Transport measurements. Both the magnetotransport and NLAHE measurements were primarily carried out in a pumped Helium-4 cryostat with a base temperature of 1.4 K. The latter was further cross-checked in an optical cryostat with base temperature of 5 K. For the second harmonic NLAHE measurements, an AC current with frequency between 17 to 277 Hz was passed along either the a-or b-axis of the (Mo/W)Te 2 crystal, and V xx , V xy , and V xz voltages were measured at both the first harmonic (X channel) and second harmonic (Y channel) frequencies using an SR830 lock-in amplifier. A Keithley 2450 source measure unit was further used to measure the DC voltage response of the NLAHE.
Optical second harmonic generation measurements. Rotational anisotropy SHG measurements were taken in a normal incidence geometry using a pulsed laser with a pulse duration of~325 fs, a repetition rate of 200 kHz, and an incoming fundamental wavelength of 800 nm. The MoTe 2 samples were held at 80 K inside an optical cryostat. Using a single-photon sensitive detector, the intensity of the reflected SHG was measured as a function of the angle between the incident polarization and the x-axis in the lab coordinate frame. The incident fundamental and the reflected SHG polarizations can be selected to be either parallel or crossed, forming two polarization channels for the SHG measurements.
Finite element method simulations. Simulations of the potential and current distributions in our devices were carried out using the static current conduction solver in Elmer, an open source multiphysical simulation software utilizing the finite element method. For each circular device, the geometry of the (Mo/W)Te 2 flake and electrodes were constructed and meshed using FreeCAD and imported into Elmer. For a given direction and magnitude of current, the MoTe 2 conductivity tensor was iteratively adjusted to match the potential difference measured across the longitudinal voltage leads, V xx . The local longitudinal electric field at the vertical contacts, E x ,was then extracted from the simulation.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.