Generalized Wilson loop method for nonlinear light-matter interaction

Nonlinear light–matter interaction, as the core of ultrafast optics, bulk photovoltaics, nonlinear optical sensing and imaging, and efficient generation of entangled photons, has been traditionally studied by first-principles theoretical methods with the sum-over-states approach. However, this indirect method often suffers from the divergence at band degeneracy and optical zeros as well as convergence issues and high computation costs when summing over the states. Here, using shift vector and shift current conductivity tensor as an example, we present a gauge-invariant generalized approach for efficient and direct calculations of nonlinear optical responses by representing interband Berry curvature, quantum metric, and shift vector in a generalized Wilson loop. This generalized Wilson loop method avoids the above cumbersome challenges and allows for easy implementation and efficient calculations. More importantly, the Wilson loop representation provides a succinct geometric interpretation of nonlinear optical processes and responses based on quantum geometric tensors and quantum geometric potentials and can be readily applied to studying other excited-state responses.


INTRODUCTION
Nonlinear light-matter interaction plays a pivotal role in ultrafast optics 1 , bulk photovoltaics 2 , nonlinear optical sensing and imaging 3 , optical transistor 4 , efficient generation of entangled photon pairs for quantum computing 5 , etc. In particular, noncentrosymmetric materials are known to hold even-order nonlinear photocurrent responses under an external electromagnetic field. For example, the wave packet of charge carriers can be displaced in real space upon photon excitation via a second-order process, resulting in shift current 2,6 that accounts for the shift mechanism for the bulk photovoltaic effect.
Field-dependent nonlinear photocurrent can be obtained by solving the quantum kinetic equation of density matrix using perturbation theory. Subsequently, it can be calculated by firstprinciples methods such as density-functional theory 7-10 and Wannier interpolation 11,12 with sum rules. However, the sum-overstates approach involves an ad hoc cutoff that induces divergence at band degeneracy and optical zeros. Moreover, it suffers from the convergence issue with respect to the number of states. A direct approach is largely underexplored.
The Wilson loop method was originally proposed by Wilson 13 in 1974 for computing gauge field on a closed path. It is ubiquitous to gauge theories and has been widely used to calculate Berry curvature, Chern number, and other topological invariants in condensed matter physics, which are the hallmarks of a rich set of low-energy transport phenomena governed by the linear response of intraband process, including quantum Hall effect 14 , quantum anomalous Hall effect 15 , spin Hall effect 16 , and quantum spin Hall effect 17,18 . Unlike the above linear responses, shift current involves interband transitions and its conductivity tensor is proportional to the quadratic electric field E 2 . Young and Rappe 19 reformulated the shift vector using a gauge-invariant discrete expression similar to the King-Smith and Vanderbilt formalism of electric polarization 20 .
Recently, Shi et al. represented the photon-drag shift vector with the Wilson loop formalism, an important geometric quantity in shift current tensor and photon-drag shift current tensor 21 . These motivate us to develop a general approach for nonlinear optical (NLO) responses by representing interband Berry curvature, quantum metric, and shift vector in a generalized Wilson loop.
Here we present a physically intuitive gauge-invariant Wilson loop approach for direct and efficient calculations of NLO responses with Wilson loop representation, using the shift vector and shift current conductivity tensor as examples. In the Wilson loop picture, the geometrical nature of the shift current can be viewed as the difference in the spontaneous polarization determined by the interband Berry connection between the valence and conduction bands upon direct optical transition. Unlike the standard sum-ofrule method, our Wilson loop approach is free of the convergence issue with respect to the number of states, and avoids the cumbersome divergence at band degeneracy and optical zeros where the dipole matrix element is zero, that is, r a mn k ð Þ mr a j jn h i¼ 0. This quantum geometric approach can be easily implemented and allow for efficient calculations.
We demonstrate this powerful approach in two representative cases, including (1) a Rice-Mele model with an extra staggered onsite potential and (2) monolayer ferroelectric GeS with firstprinciples Wannier tight-binding Hamiltonian. The results calculated by the Wilson loop approach are in excellent agreement with the exact analytic solutions of the Rice-Mele model and the numerical results of monolayer GeS with the sum-over-states approach. Furthermore, we provide a Wilson representation of the geometrical shift vector where the integral of the Wilson loop results in polarization difference between two bands upon optical transition, illustrating the geometrical nature of the shift current. In general, gauge-invariant geometric quantities, e.g., quantum metric, Berry curvature, and shift vector, can be all represented by the Wilson loop naturally. The generalized Wilson loop approach developed here can be readily applied to other linear and nonlinear responses and allow for direct geometric interpretation of these quantities.

RESULTS
Geometrical shift current response Shift current originates from the difference between the realspace charge center of the valence and conduction bands upon optical transition. It is a bulk effect as the separation of photoexcited electrons and holes does not rely on a p-n junction with a built-in electric field. Unlike conventional photovoltaics, the generated open-circuit voltage can go beyond the bandgap, hence the power conversion efficiency is not limited by Shockley-Queisser limit for a single p-n junction. Under homogeneous linearly polarized light illumination, shift current can be written as 2,6,7,19 (2) where r nm ¼ i n ∂ k j jm h ifor n ≠ m and A n ¼ i n ∂ k j jn h iare interband and intraband Berry connection for states m j i and n j i, respectively. ℏω nm corresponds to the band energy difference. f is the mn k ð Þ is the well-known shift vector described by the derivative of the phase and the difference of intraband Berry connection, also known as the quantum geometric potential 22 . Although the difference of Berry connections is gauge dependent, the shift vector is gauged invariant. We will discuss the gauge transformation property later.
The geometric aspect of the shift current is related to quantum metric and Berry curvature through the Christoffel symbols of the second kind 23,24 . The local quantum geometric tensor Q ab mn ¼ ∂ a mQ ∂ b n was originally proposed by Provost and Vallee 25 , whereQ ¼ 1 ÀP andP is the ground-state projection operator 26 . It indicates that the geodesic quantum distance between two quantum states in the Hilbert space can be viewed as absorption strength in the interband optical process, e.g., Q ab mn ¼ r a mn r b nm . Q ab mn consists of a symmetric quantum metric g ab mn and an antisymmetric Berry curvature Ω ab mn , i.e. Q ab mn ¼ g ab mn À i 2 Ω ab mn . Quantum metric g ab mn and Berry curvature Ω ab mn play quite different roles. For example, the off-diagonal g ab mn and diagonal Ω ab mm contribute to the linear response coefficients of interband and intraband processes, respectively. In contrast, both quantum metric g ab mn and Berry curvature Ω ab mn play a crucial role in second-order responses. As we will show below, gauge-invariant geometric quantities such as quantum metric g ab mn , Berry curvature Ω ab mn , and shift vector R a;b mn can be all represented by the Wilson loop naturally.

Wilson loop approach of shift vector and shift current
Gauge-invariant single band Berry curvature in a discretized Brillouin zone can be calculated by Fukui-Hatsugai-Suzuki method 27 , Ω c n k ð Þ ¼ Im ln W n k ð Þ, with (4) where q a is an infinitesimal displacement vector along the corresponding a direction near k point and ϵ abc is the Levi-Civita symbol. Now we derive a Wilson loop formula for the shift current response by reformulating shift vector with a strategy similar to that of King-Smith and Vanderbilt 20 . The intraband Berry connection contrast, the interband Berry connection between states m j iand n j i is given by r mn k ð Þ ¼ i n ∂ k j jm h i¼ r mn ðkÞ j je iϕ mn , where ϕ mn is the phase of interband Berry connection: ϕ b mn ðkÞ ¼ Im ln r b mn ðkÞ ð Þ . Thus, the shift vector can be reformulated as The first term on the right-hand side can be evaluated using finite difference around k along direction a, ∂ ka Im ln r b mn ðkÞ Thus, Since hn; kjr b jm; ki does not depend on q, we can rewrite We then arrive at the Wilson loop representation of shift vector (9) where W mn k; q a ; r b ; r b À Á denotes general absorption on a Wilson loop, i m; k þ q a r b n; k þ q a n; k þ q a jn; k h i n; k r b m; k ; Herein, the interband Berry curvature contributing to nonlinear injection current 10 can also be represented by the Wilson loop In fact, W mn k; q ¼ 0; r a ; r b À Á defined on the Wilson loop is identical to a quantum geometric tensor, and its real and imaginary part gives the symmetric quantum metric g ab mn and antisymmetric Berry curvature Ω ab mn at finite crystal momentum k, respectively.
We can further rewrite the shift current conductivity tensor using the Wilson loop representation as H. Wang et al.
We then obtain The above equation is a Wilson loop formula for shift current, which is only dependent on the imaginary part of the Wilson loop. This formula avoids the ambiguity in the definition of the argument around optical zeros where r a mn k ð Þ hm; kjr a jn; ki ¼ 0 and at band degeneracies where E m = E n or ω nm = 0. The Wilson loop representation of Berry curvature and geometrical shift vector are illustrated in Fig. 1a, b. In fact, the Wilson loop approach provides an equivalent expression as the Young-Rappe formula 19 , which involves a more complicated Wilson loop, including indirect optical transition matrix elements as shown in Supplementary Fig. 1.
Gauge invariance is guaranteed on each local Wilson loop in the discretized Brillouin zone. Under an arbitrary local gauge transformation u n k ð Þ ! e iϕ n k ð Þ u n k ð Þ; or; n; k j i ! e iϕ n k ð Þ n; k j i, Berry connections transform as Shift vector is clearly gauge-invariant Hence, the quantum geometric tensor W mn k; q; r; r ð Þ is also gauge-invariant. Figure 1c shows the geometrical Berry curvature and shift vector using a two-band model. The geometrical meaning of the shift vector in Wilson loop representation is illustrated in Fig. 1d, which clearly shows it is related to the difference of the real-space charge center or spontaneous polarization for the valence and conduction bands upon direct optical transition. It is also known that the geometrical shift vector contributes to nonreciprocal Landau-Zener tunneling 28 .

Rice-Mele model of a one-dimensional ferroelectric system
To demonstrate the generalized Wilson loop approach, we first use the one-dimensional Rice-Mele (RM) model of ferroelectric systems with broken inversion symmetry. The tight-binding Hamiltonian is illustrated in Fig. 2a, which is given by where t is the hoping parameter and δ t denotes the dimerization of the chain related to the distortion with respect to the centrosymmetric structure with t i ¼ t 2 þ À1 ð Þ i δt 2 . Δ t is the staggered onsite potential between two sites. c y i and c i are the fermion creation and annihilation operators, respectively. The inversion symmetry is broken when δ t ≠ 0 and Δ t ≠ 0. It leads to the following Bloch Hamiltonian where a is the lattice parameter. We use the following parameters for GeS, which yields a bandgap of 1.9 eV 9 , t ¼ À1:0eV, δ t ¼ À0:83eV, and Δ t ¼ À0:45eV. The shift vector of the twoband RM model 29 has an analytical solution, The conduction and valence band energies E v;c , as shown in Fig. 2b, are given by It is clear that the shift vector is reversed when δ t or Δ t changes sign, enabling ferroelectric-driven shift photocurrent switching 10 . This is also verified by the Wilson loop approach as shown in Fig. 2d. Numerically, the shift vector and shift current conductivity are usually calculated by the sum rule with mass or diamagnetic term for a tight-binding model. The generalized derivative of interband Berry connection for shift vector and shift current can be expressed by using the sum rule as 7,30 where Δ b nm ¼ v b n À v b m is the group velocity difference and hω nm ¼ hω n À hω m is the band energy difference. The mass term w ab nm ¼ n ∂ ka ∂ k b H j j m h icannot be neglected for the tight-binding model because the interband Berry connection is not gauge-covariant and its generalized derivative in the Hamiltonian gauge involves the second-order derivative of Hamiltonian 31 . In a two-band RM model with m = 1, n = 2, the shift vector using the sum rule at optical nonzero k-points (i.e.
It shows that the shift vector and shift current are vanishing in the two-band RM model without considering the mass term. The calculated shift vector and shift current conductivity with an effective area of 9.37 Å 2 by different methods are shown in Fig. 2c, d. The results demonstrate that the shift vector and the shift current calculated by the Wilson loop approach are in excellent agreement with the analytical solution and the sum-over-states approach.
Wannier tight-binding model of monolayer GeS Next, we demonstrate the Wilson loop method for real materials with a symmetrized Wannier tight-binding Hamiltonian. The details of the first-principles calculations and the constructions of symmetrized Wannier tight-binding Hamiltonian are described in Methods. Taking the ferroelectric monolayer GeS as an example, it has a C 2v point group with a mirror plane M x perpendicular to the x-axis. The crystal structure and band structure of monolayer GeS are shown in Fig. 3a, b, respectively. From group theory analysis, the components σ xxx (ω) and σ xyy (ω) vanish under linearly polarized light with in-plane polarization, which was verified in our calculation. Here, we focus on σ yyy (ω) and the corresponding shift vector R y;y cv . Figure 3c shows k-resolved shift vector R y;y cv between the top valence band and the bottom conduction band across the bandgap. The shift vector away from optical zeros can be~10 Å, much larger than its lattice constant. This is very different from a spontaneous electric polarization vector constrained within the lattice constant. Berry curvature of the top valence band is also calculated by the Wilson loop approach as shown in Fig. 3d. Given the mirror symmetry M x , we have verified the symmetry properties R y;y cv k x ; k y A Berry curvature dipole along x-direction is generated, which can induce a similar ferroelectric nonlinear Hall effect [32][33][34] in monolayer GeS. The intraband Berry curvature of the bottom conduction band and the interband Berry curvature across the gap are presented in Supplementary Fig. 2. Figure 4a shows the calculated frequency-dependent shift current conductivity σ yyy (ω). To convert the sheet conductivity σ 2D to bulk conductivity σ 3D , we set the effective thickness l to be x y z Ge S ky k x   .56 Å by σ 3D = σ 2D /l. We have verified the identity W mn k; q ¼ 0; r b ; r b À Á ¼ r b mn r b nm . The k-resolved absorption strength between the top valence band and the bottom conduction band is shown in Fig. 4b. The white region indicates optical zeros and has no contribution to shift current conductivity. To investigate the origin of large responses, we calculate the k-resolved shift current strength I a;b mn k; ω ð Þ at ω = 2.0 and ω = 2.8 eV using the Wilson loop approach, defined as The results are shown in Fig. 4c, d. The convergence was carefully checked with respect to the number of bands and k-point sampling in the first Brillouin zone, as shown in Supplementary Fig. 3, for shift current conductivity. The calculated shift current conductivity is well converged with a k-point mesh of 200 × 200 × 1 and two bands below 3 eV. In addition, frequencydependent σ yxx (ω) for monolayer GeS is shown in Supplementary  Fig. 4. Furthermore, we performed similar calculations with the Wilson loop approach for a different 2D materials monolayer WS 2 , and the results are shown in Supplementary Fig. 5. Our results clearly demonstrate that the generalized Wilson loop approach is not only efficient and generally applicable to both effective models and realistic materials but also avoids the summation over a large number of intermediate valence and conduction bands, making it valuable for computing NLO responses.
It should be noted that the geometrical shift vector at optical zeros cannot be calculated using the sum-over-states approach. Furthermore, while the large shift vector at optical zeros has no contribution to shift current response with vertical transitions, it can contribute to the shift current response when taking into account the photon-drag effect 21 involving indirect transitions or strong electron-phonon coupling effect. Our demonstration of geometrical shift vector in real materials will allow for theoretical investigations of a wide range of geometric effects induced by quantum geometric potential.

DISSCUSSION
A common challenge with a perturbation theory for NLO responses in the length gauge is the treatment of the position operator r for the extended Bloch states. The intraband part of the position operator is represented by m; k r i j jn; k 0 h i¼δ mn δ k À k 0 ð ÞAk ð Þþ ð i∂ k δ k À k 0 ð ÞÞ. The real-space coordinate r of the wave packet made from the Bloch wavefunctions is represented by r i ¼ i∂ k À A k ð Þ. NLO responses involve the matrix element of the commutator Þ plays a central role in many other nonlinear responses. The derivation of the commutator relation can be found in Supplementary Information. For example, the generalized derivative of dipole matrix element is written as which is a key physical quantity for second harmonic generation 7,35 . Hence, the Wilson loop approach developed here can be readily applied to other nonlinear optical effects such as second and third harmonic generation and linear and quadratic electro-optic effects. The spin-orbit interaction is weak in monolayer GeS, thus it is not considered in the present calculations. Nevertheless, the expression can be easily extended to include spin-orbit coupling, and generalized to the degenerate and near degenerate cases by considering connected and disconnected band sets and using k Á p perturbation theory to obtain a smooth variation of the Bloch states between k and k þ δk 36 . Furthermore, although all the calculations in this work are performed within the independent particle approximation, the Wilson loop approach can also be developed to include many-body effect.
In summary, we presented a gauge-invariant generalized approach for efficient and direct calculations of NLO responses with pure Wilson loop representation. This generalized Wilson loop method avoids the cumbersome issues of the commonly used sumover-states approach and allows for easy implementation and efficient calculation. The Wilson loop representation provides an elegant geometric interpretation of nonlinear optical processes and responses based on quantum geometric tensors and quantum geometric potentials responsible for shift current and Landau-Zener tunneling. The generalized Wilson loop method developed here can be readily applied to study other linear and nonlinear responses such as second and third harmonic generation, linear and quadratic electro-optic effect, as well as magnetic injection current and magnetic shift current 37 .

First-principles calculations of atomistic and electronic structure
First-principles calculations for structural relaxation and quasiatomic Wannier functions were performed using density-functional theory 38,39 as implemented in the Vienna Ab initio Simulation Package (VASP) 40 with the projector-augmented wave method for treating core electrons 41 . We employed the generalized-gradient approximation of exchange-correlation functional in the Perdew-Burke-Ernzerhof form ref. 42 , a plane-wave basis with an energy cutoff of 400 eV, and a Monkhorst-Pack k-point sampling of 10 × 10 × 1 for the Brillouin zone integration.

Generalized Wilson loop approach of shift current using firstprinciples tight-binding Hamiltonian
To compute the Wilson loop related quantities, we first construct quasiatomic Wannier functions and symmetrized first-principles tightbinding Hamiltonian from Kohn-Sham wavefunctions and eigenvalues under the maximal similarity measure with respect to pseudoatomic orbitals 43,44 . A total of 16 quasiatomic Wannier functions were obtained for monolayer GeS. Using the developed tight-binding Hamiltonian we then compute Berry curvature, shift vector and shift current using a dense kpoint sampling of 200 × 200 × 1. Sokhotski-Plemelj theorem is employed for the Dirac delta function integration with a small imaginary smearing factor η of 0.02 eV. We checked the convergence of shift current conductivity tensor with respect to the number of bands and the k-point sampling as well as different q values of the reciprocal lattice used in the Wilson loop method. (see Supplementary Figs. 3, 4).

Symmetrization of the tight-binding Hamiltonian
The construction of Wannier functions for crystals may not preserve space group symmetries. To avoid artificial symmetry breaking, we performed symmetrization of the tight-binding Hamiltonian. The Hamiltonian is invariant under symmetry operation g in the group G where D k g ð Þ is k-dependent representation of the symmetry and t g is the translation vector of the symmetry. We define the symmetrized Hamiltonian using the group averagẽ To apply the group average to the above tight-binding Hamiltonian with all crystalline symmetry constraints in real space, we rewrite the hopping matrices 45 where S g is the real-space rotation matrix and T ml ij ¼ S g r m À r l ð ÞÀr j À r i , r i is the position of the localized orbitals in the unit cell. The band structure of monolayer GeS with and without symmetrization of the Hamiltonian is shown in Supplementary Fig. 6.