Measurement of the Berry curvature of solids using high-harmonic spectroscopy

Berry phase and Berry curvature have become ubiquitous concepts in physics, relevant to a variety of phenomena, such as polarization, various Hall effects, etc. Studies of these phenomena call for characterization of Berry phase or curvature which is largely limited to theory, and a few measurements in optical lattices. In this work, we report polarimetry of high-harmonic emission from solids and exploit this novel capability to directly retrieve the Berry curvature of α-quartz. We show that the two manifestations of broken inversion symmetry in solids lead to perpendicular or parallel polarization of even harmonics with respect to the driving field. Using semiclassical transport theory, we retrieve the Berry curvature from spectra measured in perpendicular polarization, the results being supported by ab initio calculation. Our work demonstrates an approach for the direct measurement of Berry curvature in solids, which could serve as a benchmark for theoretical studies.


SUPPLEMENTARY NOTE 1: THE ROLE OF BROKEN SYMMETRY IN NONLINEAR RESPONSES Perturbative response -expansion to different broken inversion symmetries
Although these derivations have been given before, it is convenient to present them here in a consistent formalism to provide the reader with a coherent overview of the topic (and prepare the later extension of semiclassical transport theory in solids).

Non-inversion symmetric media
Assuming that the electric field has the general form E(t) = E 0 H(t) sin(ω L t) where E 0 , H(t), ω L are the peak, envelope, angular frequency of the laser electric field, the perturbative response derived from perturbation theory [1] is P medium (t) = 0 ∞ n=1 χ (n) E n (t) where n runs over all odd and even positive integers because of non-vanishing even-order susceptibilities in non-isotropic media [1]. In the long pulse regime, this formula can be expanded using De Moivre's formula and binomial theorem as follows where F [X] is the Fourier transform of X: with a = (2m + 1 − 2k)ω L , b = (2m − 2k)ω L .
Clearly, the Fourier transform of P medium (t) contains both odd and even harmonics of the fundamental carrier frequency ω L . If the trigonometric function on the right-hand side of Supplementary Eq. 1 is of cosine type (the same electric field, with only a change of carrier-envelope phase of π/2), the expansion can be performed similarly and the result stills includes both even and odd harmonics.
Inversion-symmetric media with the electric field as the source of symmetry breaking In this case, certainly the even-order susceptibilities vanish completely, only the odd orders χ (2n+1) remain. However, the addition of the second harmonic to the fundamental electric field leads to a non-inversion-symmetric electric field: E(t) = E 0 sin(ω L t) + E 1 sin(2ω L t) . Similarly, we can also expand the perturbative response of this electric field as follows: The A term in Supplementary Eq. 5 is very similar to the right-hand side of Supplementary Eq. 3 except for the running index which is always odd. Therefore, it is straightforward to see that the A term includes all odd harmonics in the spectral response. Similarly, the B term includes 2a and 2b thus all the even harmonics are included in this term. As a result, P Et (ω) consists of both even and odd harmonics.
The above consideration shows that regardless of the exact cause of broken symmetry (medium or electric field), as soon as the inversion symmetry in the process is broken, the polarization response will include both even and odd harmonics, although with different spectral intensities.
Origin of even and odd harmonics in the semiclassical model without anomalous term for solids We start our discussion with the simplest model describing HHG from solids: high-energy photons are emitted due to Bragg scattering of a time-dependent electron wavepacket in a conduction band (beyond nearest-neighbor approximation) of solids. Here we follow previous works [2][3][4][5][6][7] and slightly extend recent work [8] to show how even and odd harmonics could be emitted in the framework of this semiclassical model and its implications on the quantum-mechanical model. The incident electric field is a linearly polarized 30 fs pulse at the carrier wavelength of 800 nm.
From the definition of the group velocity and the "acceleration theorem" [9,10], which is very similar to the main text, without the anomalous terms: where v ν (k), E ν (k) are group velocity and band dispersion of a given band ν, E and B are electric, magnetic fields. In the non-relativistic regime, v × B E thus we arrive at The solution of this differential equation is with the vector potential Here we emphasize that we do not take for granted that k(0) = 0 and we keep k(0) in the formula since it plays an important role in the harmonics, as will be revealed later. Since we can always decompose the band dispersion into Fourier series as where n max is the maximum number of distant neighbors considered in our model and a is the lattice constant. We should emphasize that the trigonometric function on the right-hand side of Supplementary Eq. 10 can only be of cosine type, not sine type, due to time-reversal symmetry [11], which is very important for later expansions. Inserting Supplementary Eq. 10 into Supplementary Eq. 7 we have Taking Supplementary Eq. 9 into account, we arrive at which can be rewritten using the trigonometric identities If we assume , ω L , ϕ CE are the peak field strength, envelope, carrier angular frequency and carrier-envelope phase of the vector potential of the incident electric field, Supplementary Eq. 13 can be expanded as In the multi-cycle regime or monochromatic electric field, G(t) = 1, we can apply the Jacobi-Anger expansion to arrive at Here J l is the Bessel function of the first kind of order l. Now the associated spectrum of the intraband current can be calculated as proportional to Equations 15, 16, 17 comprise the full analytical solution of the semiclassical model under minimal approximation. Two conclusions can be drawn at this stage: The spectrum associated with the intraband current consists of modulated Dirac combs with both even and odd harmonics (represented in Supplementary Eq. 16 and 17) which are included in the first term and second term respectively, under general conditions. The spectral intensity of the harmonics depends on J l , n, ν,n , a, E0 ωL , but are independent of ϕ CE because |F [A]| 2 ω=lωL = |F [B]| 2 ω=lωL = π/2, which is expected for a long pulse. And the consequences of these equations are It is important to remember that time-reversal symmetry results in a symmetric (even function) band dispersion, thus the trigonometric function in Supplementary Eq. 12 is of sine type. If time-reversal symmetry was not in force, the trigonometric function in Supplementary Eq. 12 would be of cosine type hence the first term in Supplementary Eq. 15 would contain odd harmonics whereas the second term would contain even harmonics. The result is very similar to the case of even harmonics due to Berry curvature as discussed in [12].
In a typical semiclassical model, it is natural to start with k(0) = 0 thus sin nak(0) = 0, cos nak(0) = 1 hence first term in Supplementary Eq. 15 vanishes completely and the second term is maximized. This means that the even harmonics are canceled entirely, only the odd harmonics remain. Therefore, a typical semiclassical simulation in this case will result in a spectrum of purely odd harmonics.
If for any reason, the initial electron wavepacket in the semiclassical model does not start at 0 which means k(0) = 0, both the first term and the second term will exist, the emitted spectrum will inevitably contain both odd and even harmonics. This result is best illustrated in Supplementary Figure 1 using the semiclassical model [8].
Odd and even harmonics in the semiclassical model with anomalous term (Berry curvature) Here we consider an extended version of Supplementary Eq. 6 which is Eq. 1 in the main text.
where v ν (t) is the ordinary, parallel component of the current that has the full expansion described in Supplementary Eq. 15. Similarly we can expand the second term v ν (t) ⊥ using the Fourier expansion of Ω ν (k) (Eq. 2 in the main text) as: Now for simplicity, we assume that H(t) = G(t) = 1 and ϕ CE = 0, Supplementary Eq. 21 reduces to: We consider the limit where k(0) = 0 then Supplementary Eq. 23 reduces to: The associated emitted spectra can be calculated using F v ν (t) ⊥ as: where Ω B = aeE 0 /h is the Bloch frequency. Supplementary Eq. 26 is the general expression in the long pulse limit of anomalous current due to non-vanishing Berry curvature, which is very similar to [12]. Clearly, the associated spectrum contains only even harmonics of the fundamental laser frequency. The most important message from this derivation is that the broken symmetry of the medium results in an odd function of the Berry curvature, therefore its trigonometric function is of sine type. This together with the Jacobi-Anger expansion leads to Supplementary Eq. 25,26 which contain only even harmonics.

Semiclassical simulations
As a simple illustration for the consideration of the broken symmetry in the semiclassical model, we performed numerical simulations [8,13] for a linearly polarized 30 fs pulse at the carrier wavelength of 800 nm, peak electric field strength 0.5 V/Å and the results are shown in Supplementary Figure 1. There are several requirements for constructing the EUV polarizer which would be suitable for our experiments: The polarizer should be able to transmit a broad photon energy range, preferably much broader than the range of our experiments: 0 − 30 eV, well in the vacuum ultraviolet regime.
The polarizer should have a high contrast between the horizontal and vertical polarization.
The polarizer should not change the direction of the output beam as compared to the input beam so that the energy calibration of the whole spectrometer during polarizer's measurements would not be an issue.
The first requirement results in the use of metallic mirrors at grazing incidence. The grazing-incidence angle is optimized such that the total throughput and contrast are maximized (second requirement). In order to satisfy the third requirement, a symmetric assembly should be used. The final design consists of four bare gold mirrors, all at grazing angle of 20 degrees, assembled as illustrated in Supplementary Figure 2a and main text, Fig. 2a. The whole assembly is placed on a rotary stage that has the rotation axis identical to the propagation axis of the laser beam, which satisfies the third requirement.  [14]. The enhancement factor (or contrast) is calculated as the total reflectivity for S-pol divided by the total reflectivity for P-pol and plotted in Supplementary Figure 2d after four reflections. Here, the calculated contrast between the S-pol and P-pol beam at 17 eV (11 th harmonic) is ≈ 70 times.

Calibration
In order to properly examine the polarization of the HHG spectra from solids, the polarizer contrast should be calibrated. Its calibration is done by performing measurements of HHG from gases under identical experimental conditions. It is known that HHG from interaction of noble gases with linearly polarized multi-cycle pulses is linearly polarized, parallel to the polarization of the input pulse. We performed polarization measurement of HHG from xenon under high contrast linearly polarized input pulses and the result is shown in Supplementary Figure 3a.
The cross-cut of the spectrogram at 11 th harmonics (17 eV) shows the contrast, spectral intensity at parallel polarization divided by spectral intensity measured in perpendicular polarization, of ≈ 60. Because the polarization purity of the incident laser electric field is very high, more than three orders of magnitude, we can conclude that the contrast of our polarizer at 17 eV is 60 : 1.
Since the experiments in the main text show the contrast of 30 : 1 at 17 eV, it means that the remaining perpendicular polarization accounts for not more than 2/60 which is ≈ 3.3% of the total spectral intensity. Therefore all polarizations measured in the main text are linearly polarized, with a depolarization ratio of less than 4%.

Methodologies
Although the ideas related to Berry phase and curvature have been developed over the last three decades, actual calculations of Berry phase and curvature of condensed-matter systems only recently gained interest due to development in studies of electronic charge and spin transport of electrons. Here we outline three methodologies for calculating the Berry curvature that are used in our work. Interested readers could consult [15] for a more comprehensive review on this topic.

Direct implementation
In solid-state physics, the parameter space R on which the quantum system evolves is k-space. Thus the generic Berry curvature Ω n (R) is represented by Ω n (k) where n is a quantum number specifying the band index and k is the crystal momentum. The definitions of the Berry connection A n (k) and Berry curvature Ω n (k) [15][16][17] are given by , where u n,k (r) is the cell-periodic part of the Bloch wavefunction Ψ n,k (r). The integral in Supplementary Eq. 28 is performed over the whole unit cell, spanned by the (â 1 ,â 2 ,â 3 ) primitive vectors. Since the derivative operation is done in k-space, the direct implementation suffers from the fact that a large number of ab-initio calculations are needed in order to converge the calculation of Supplementary Eq. 28, which in extreme cases can reach up to one million k-points [15].

Kubo formula
In order to overcome the expensive calculation cost induced by the direct implementation, there have been multiple methodologies proposed [15,17]. One widely used methodology is the use of a Kubo-like formula [18,19]: Here,v is the usual velocity operator, ε i,k is the eigen-energy of the band index i and crystal momentum k, and for compact representation, the Dirac notation is used. Equation 29 replaces the computationally demanding derivative operation in Supplementary Eq. 28 by evaluation of the velocity matrix element and the loop through all the band indices {i}. Theoretically, all the band indices (occupided and unoccupied states) have to be included in the evaluation of Supplementary Eq. 29. However, practically only a small number of bands are needed as we will show in the next section.

Fukui formula
Another good method that we employed in our calculation is a technique introduced by Fukui, Hatsugai, and Suzuki [20,21] and being used later [22,23] by defining a so-called link variable: whereμ is the vector in the direction µ = (k x ,k y ,k z ) with the fractional amplitude 1/(N µ − 1) with N µ being the number of k-points used in discretizing the Brillouin zone in µ direction. Next, the lattice field strength is defined by and for large N µ , the corresponding Berry curvature will be calculated as with the area of the plaquette A plaquette =k x ·ky (N1−1)(N2−1) where (k x ,k y ,k z ) are primitive vectors in the reciprocal lattice. This method has been proven [20,21] to be extremely stable even for a very small number of k-points.

Numerical considerations
Since SiO 2 (α−quartz) has a hexagonal lattice, not the simple cubic lattice with orthogonal primitive vectors, thus (k x ,k y ) = 60 • and (k x ,k z ) = (k y ,k z ) = 90 • as shown in the main text, Fig. 2b, and vectorial k is used throughout the calculations. Here, the direct implementation and Kubo formula result in a full vectorial Berry curvature Ω n (k) = (Ω x,n (k), Ω y,n (k), Ω z,n (k)) whereas the Fukui formula gives the result in only one predefined direction Ω i,n (k) where i = (x, y, z). Within the scope of our work, the quantity of interest is Ω z,n (k) for all k points making the directions Γ − M and Γ − K with the corresponding fractional k-points (0, 0, 0) − (0, 0.5, 0) and (0, 0, 0) − (1/3, 1/3, 0). Because there are 3 silicon atoms and 6 oxygen atoms in a unit cell, there are 48 valence electrons accordingly contributing to the valence bands. Due to spin, each band is doubly degenerate thus there are in total 24 valence bands. Our focus is on the bands close to the Fermi level which has the quantum number n close to 24 (top valence band) and 25 (bottom conduction band). The focus of our calculation is the first conduction band in Γ − M direction where there is no degeneracy with other bands, leaving a smooth band dispersion and accordingly no sharp jump in the Berry curvature.
For the direct implementation, we used 6 th order central finite difference for the evaluation of the derivation. Combining this with three different k directions, it results into 18 different sets of u n,k (r) and corresponding Bloch wavefunctions Ψ n,k (r) for each k point and a given band index. The step size for the derivative evaluation is different (much smaller) than the step size used for the evaluation of the Berry curvature on the whole k points spanning the direction of interest.
For the Kubo formula, there is no derivation in k-space, but instead an evaluation of the velocity matrix elements which can be easily carried out using Fast Fourier Transform. In addition, all of the band indices need to be included thus the total number of Bloch wavefunctions required is not small.
For the Fukui formula, only one band index is needed but due to Supplementary Eq. 30, 31, 8 different sets of u n,k (r) are required for evaluation of Ω i,n at one k point.
All calculations are performed at room temperature (300K).

Convergence
For all calculations reported in this manuscript, the self-consistent calculation was performed first using an existing software package [24,25], then the calculation of the Bloch wavefunctions is done using k-meshes of different sizes to test for convergence. The results are presented in Supplementary Figure 4.
As it is expected from the nature of the direct implementation and the Kubo formula, there is no dependence of the Berry curvature calculation on the size of the grid that spans Γ − M, Supplementary Figure 4a,c. However, if the fractional dk is 10 −4 or larger, the calculation does not converge, as shown in Supplementary Figure 4b. Therefore, convergence of the direct implementation is reached when the fractional dk is 10 −6 or smaller which is very compatible with the requirement mentioned in [15]. For the Kubo formula, as mentioned in the previous section, evaluation of Supplementary Eq. 29 is dependent on the number m of the bands taken into account. Figure S 4d shows that as long as more than 7 bands below and 7 bands above the desired band are taken into account, the calculation is converged. Nevertheless, for the Fukui formula, the standard deviation changes as a function of the grid size. Reliable results can only be achieved if the Brillouin zone is discretized by 100 or more k points. For speed and convenience reasons, the direct approach is the method of choice for any further calculations.
Dependence on exchange-correlation energy functional approximations Supplementary Figure 5a,b,c show how different exchange-correlation approximations affect the bandstructure and consequently the calculated Berry curvature (Supplementary Figure 5d). While it was well known that the introduction of the modified Becke Johnson exchange potential with correlation effects [26] should result in a more accurate bandgap. However, the bandgaps retrieved in all approximations here are different from the experimentally measured one of about 9 eV [27,28]. This reveals the shortcommings of the LCAO approach as compared to the full-potential (linearized) augmented plane-wave and local orbitals (FP-(L)APW + lo) method as used in [26,29]. For this reason, our experiment and follow-up developments provide a new approach to benchmark ab-initio calculations.  Figure  6b. In addition, if we consider the fact that VB3 is more likely to be filled than CB1 to be populated, this leaves CB1 as the most effective band, influencing Berry phase effects in this direction which would correspond to the Ω effective (k) in the main text. Furthermore, Supplementary Figure 6c shows that the Berry curvature of the same band is very much (100 times) weaker in the Γ − K direction. For symmetry reason, it should be zero. In the semiclassical model, this translates into a spectral intensity that is about four orders of magnitude less intense. This is certainly too weak to be seen within the precision of our measurements.
A direct comparison between the experimentally retrieved Berry curvature and the theoretically calculated value is shown in the main text, Fig. 3b. Near-quantitative agreement is obtained with the MGGA functional (shown) and the LDA functional (not shown), whereas the GGA functional appears to overestimate the Berry curvature by a factor of ≈ 2. Discrepancy of the values can be attributed to: (i) imperfections of the experimental measurements, (ii) first order derivation of the adiabatic transport theory utilized in determining the retrieved value from experimental results, and (iii) limitations inherent to the ab-initio calculations: even in the DFT:LCAO technique, different functionals used already lead to amplitude variations up to a factor of ≈ 2. (Supplementary Figure 5d Figure 7 shows the calculated density of states in SiO 2 crystal using different approximations. Because the density of states expands more than 30 eV with a relatively constant amplitude for both crystal directions, this is not in favor of the picture of multiple plateaus due to higher-lying conduction bands [30,31]. As a result, more extensive experiments and simulations should be carried out in order to fully explain the second plateau behavior in HHG from quartz. In many cases where the absolute spectral intensity is not needed, only the relative intensity is of importance, for instance in Fig. 2b,c,d of the main text, we do not need to calibrate the spectral intensity perfectly. However, for the application of retrieving the Berry curvature, accurate spectral intensity is needed. In order to determine the spectral intensity as precisely as possible, we took care of the following issues:

Wavelength calibration
Since the spectrometer is grating-based, we perform the calibration in wavelength domain first, then the spectrum is converted to frequency domain for plotting. By measuring the spectrum of the incident light, the carrier central wavelength is determined to be λ 0 ≈ 798.9 nm. We recorded HHG spectra from gases and solid samples. The formed high-harmonics are calibrated using the grating equation. Note that the value in Γ − K direction (should be zero due to symmetry reason) has been multiplied by 100 to be visible in the same plot.

Intensity calibration
Multi-channel plate and phosphor's screen are assumed to have a flat spectral response in our measurement range.
Polarizer reflectivity: for S-pol and P-pol are calculated using [14] as discussed before.
Grating reflectivity: the grating efficiency for the two different polarizations (S-pol and P-pol) are obtained from the manufacturer and applied to the corresponding cases.
Effect of slit on the total spectral intensity: since the different wavelengths diffract differently, the spatio-spectral distribution of the beam after the slit could be illustrated in Supplementary Figure 8. Therefore, the beams carrying high-energy photons have a smaller divergence, thus the percentage of transmission is higher than for low-energy photons. As a first order calibration to this effect, we use the following formula: S(ω) averaged after slit = S(ω) averaged before slit · ω 2 First order propagation approximation and beam profile averaging effect: because we measured only the spectra, the phase was not measured, thus we could not perform the back propagation to get the exact spatially resolved EUV spectra at the focus. Therefore, we assume the averaged spectral intensity measured before the slit is the  signal at the center of the beam (the degree of nonlinearity is high). First order propagation inside the crystal is approximated by: S(ω) propagated = S(ω) microscopic · ω 2 . As a consequence, we have: Electric field strength: is calculated using the measurement of the beam profile at the focus for various intensity settings and the temporal characterization of the electric field using Transient Grating -Frequency-Resolved Optical Gating.

Fitting of the Berry curvature
Basis of the fitting: Associated spectra of different Fourier coefficients of the Berry curvature In order to examine the reliability of the fitting, we make use of the analytical derivation of the associated spectra of the Berry curvature: Supplementary Eq. 25, 26. Although the derivations apply for the long pulse regime, because we are only interested in the relative intensity of the harmonic peaks, not their shapes, these equations are of great applicability. Indeed, we could assign Supplementary Eq. 25 as v ν (t) ⊥ = eE0 h ∞ n=1 γ n A n with A n = − ∞ m=1 (−1) m J 2m−1 (na ē h E0 ωL ) sin [2mω L t] + sin [(2m − 2)ω L t] . Now we could plot associated spectra of different {A n } as illustrated in Supplementary Figure 9 for the electric field strength of E 0 = 1.0 V/Å. Due to the properties of the Bessel functions, the summation over m converges at m = 30 for our energy range. From Supplementary Figure 9, it is clear that: (i) all A n exhibit only even harmonics, with different number of harmonics as well as their relative intensity; (ii) the higher the number n is, the higher the associated photon energy is. The first feature is obvious from the formula while the second feature is a consequence of the Bessel functions. Therefore, {A n } forms a basis set of HHG spectra from Berry curvature. This basis set is possibly complete but redundant because the harmonics of A n at high n are very similar. In addition, within our photon energy range, all A n with n > 20 contribute very little and thus play almost no role in the fitting result.
In conclusion, the above analysis shows that the basis set formed by the series {A n } could be well utilized for retrieving the Berry curvature from our measured HHG spectra, especially for the range of n = 1 → 10.

The fitting
We first simulated the k(t) using the measured electric field and a given set of dipole matrix elements. The fitting is done using Eq. 2 (main text) and the trust-region reflective algorithm in a commercial software package. The fitting error is defined as: Since the basis set used above is not orthonormal, there are different sets of Fourier coefficients that give similar fitting results. Therefore, in order to quantify more accurately the retrieved Berry curvature, we performed a large number (10 4 ) of fittings using random-number generation for the initial guesses. The fitted results with errors lower than a given threshold are then used to illustrate the probability density of the retrieved Berry curvature given in Fig.  3b main text. More accurate measurements and calibration of this experiment could be carried out in the future, for instance with atomically-thin quartz samples, which could help deriving a more accurate effective Berry curvature in quartz. Unfortunately such samples are not presently available.