Design principles for shift current photovoltaics

While the basic principles and limitations of conventional solar cells are well understood, relatively little attention has gone toward maximizing the potential efficiency of photovoltaic devices based on shift currents. In this work, we outline simple design principles for the optimization of shift currents for frequencies near the band gap, derived from the analysis of a general effective model. The use of a novel sum rule allows us to express the band edge shift current in terms of a few model parameters and to show it depends explicitly on wavefunctions via Berry connections in addition to standard band structure. We use our approach to identify two new classes of shift current photovoltaics, ferroelectric polymer films and single-layer orthorhombic monochalcogenides such as GeS. We introduce tight-binding models for these systems, and show that they exhibit the largest shift current responsivities at the band edge reported so far. Moreover, exploring the parameter space of these models we find photoresponsivities that can exceed $100$ mA/W. Our results show how the study of the shift current via effective models allows one to improve the possible efficiency of devices based on this mechanism and better grasp their potential to compete with conventional solar cells.

C ost-effective, high-performing solar cell technology is an essential piece of a sustainable energy strategy. Exploring approaches to photo-current generation beyond conventional solar cells based on pn junctions is worthwhile given that their performance is in practice constrained by the Shockley-Queisser limit 1 . One of the most promising alternative sources of photocurrent is the bulk photovoltaic effect (BPVE) or 'shift current' effect, a nonlinear optical response that yields net photocurrent in materials with net polarization [2][3][4][5][6][7][8][9][10][11] . Contrary to conventional pn junctions, the BPVE is able to generate an above band-gap photovoltage 12 , potentially allowing the performance of BPVE-based photovoltaics to surpass conventional ones. However, closed-circuit currents generated via the BPVE reported in the literature have typically been small compared with those generated in pn junction photovoltaics [13][14][15] . Recent interest in the BPVE also stems from the proposal that it may be at work in a promising class of materials for photovoltaics known as hybrid perovskites 13 , an extremely active field of research [16][17][18][19][20][21][22][23][24][25][26][27][28][29] .
The fundamental requirement for a material to produce a current via the BPVE is that it breaks inversion symmetry, allowing an asymmetric photoexcitation of carriers. But despite considerable case-by-case study of the BPVE, the necessary ingredients to optimize a BPVE-based solar cell are not sufficiently well understood. As with conventional solar cells, band gaps in the visible (1.1-3.1 eV) 15,30 and large electronic densities of states 14,31 are always beneficial. In addition, to produce a solar cell that responds to unpolarized sunlight, a highly anisotropic material must be used, since otherwise there is no preferred direction for the current to flow. But beyond these natural requirements, our only guiding knowledge is that the shift current depends explicitly on the nature of the electronic wavefunctions 31,32 and that it is not correlated with the material polarization in any obvious way 15 despite the fact that both shift currents and polarization originate from inversion symmetry breaking.
In the current situation, a more generic understanding of what makes the BPVE strong is highly desirable. When tackling complex material science problems, stripping off all complications and optimizing the simplest model that captures the relevant physics often proves the best strategy, as shown for example in thermoelectricity studies [33][34][35] . In this work, we present simple design principles for BPVE optimization based on the study of an effective model for the band edges. With this model, the band edge shift current is given by the product of the joint density of states (JDOS) and a matrix element, both given by simple expressions in terms of a few model parameters. The simplicity of the model allows us to derive the main principle that band edges with semi-Dirac type of Hamiltonians are the best starting point to obtain large band edge prefactors. In addition, by relating the effective model parameters to realistic tight-binding models, we can predict that several materials with the required band structure have larger shift currents than any reported so far.

Results
Density of states in one-and two-dimensions. In our search for materials we should look for large JDOS in systems where the band edge is closely aligned with the peak of the solar spectrum, around 1.5 eV. Since the band edge always induces a Van Hove singularity in the density of states, the requirement of a large peak in the photoresponse can be naturally better satisfied by lowdimensional materials, which generically present stronger singularities 36 . Materials of one and two dimensions are therefore the focus of this work. Among one-dimensional materials, ferroelectric polymers are suitable candidates for shift-current photovoltaics: they strongly break inversion symmetry, some have suitable band gaps for photovoltaics applications [37][38][39][40] , and they can be produced in macroscopically oriented samples. For these reasons, we consider solar cells consisting of such polymer films, shown in Fig. 1a. Two-dimensional materials 41  potential for photovoltaics, as shown by demonstration of a pnjunction photovoltaic effect in dichalcogenide heterostructures [42][43][44] , and in few-layer black phosphorus 45 . However, these well known two-dimensional (2D) semiconductors have vanishing shift currents because of either inversion or rotation symmetry. Group IV monochalcogenides have emerged in the past years as a new familiy of inversion-breaking, anisotropic 2D materials with fascinating properties [46][47][48][49][50] , and interest in growing as thin films of all four members of the family, GeS 51-54 , GeSe 53,54 , SnS 55,56 and SnSe 57-59 has now been isolated experimentally. In this work, we show that GeS is ideally suited to realize high values of the BPVE. Their GeS structure is shown in Fig. 1c.
To understand how to optimize the photoresponse, we first discuss how the shift current can be computed for a tight-binding model, and then we proceed to apply this formalism to describe a generic band edge and the response of particular materials.
Shift current. In this work we consider the shift current contribution to the BPVE and we shall use both terms interchangeably (note the BPVE can have other contributions as well 6 ). With electric field E b (o) at frequency o and linearly polarized in the b direction, the shift current is a DC response of the form 6 Defining an intensity for each polarization, , we define the photoresponsivity k abb as the current density generated per incident intensity J a ¼ k abb I 0,b , which gives k abb ¼ 2s abb =cE 0 . Note that in conventional solar cells the current is also linear with intensity. For a D-dimensional system, k abb takes the form 7,9 where C ¼ 4g s pe 3 =' 2 E 0 c, with c being the speed of light, E 0 the vacuum permittivity and g s ¼ 2 accounts for the spin degeneracy.
In what follows we set : ¼ 1. Summation of indices is explicitly indicated using the summation symbol. The sum is over all Bloch bands, with o nm ¼ E n À E m the energy difference between bands n and m and f nm ¼ f n À f m the difference of Fermi occupations, which we take at zero temperature. The integrand is where r a nm are the inter-band matrix elements of the position operator (or inter-band Berry connections), defined as r a nm ¼ i nj@ ka m h ifor na m and zero otherwise, where n j i is the eigenstate of band n. A semicolon denotes a generalized derivative Generic two-band model. With the aim of describing the shift current response of the band edge of a semiconductor, next we consider the shift current of a generic two-band model. The Fourier transform of the real space Hamiltonian is performed with the choice of phases c m;k x ð Þ ¼ 1 is a localized orbital and x i is the position of site i in the unit cell. This choice is made in order to naturally incorporate the action of the position operator, see refs 60-62. The Hamiltonian matrix takes the form where s 0 is the identity matrix, s i ¼ s x , s y , s z are the Pauli matrices and E 0 and f i ¼ f x , f y , f z are generic functions of momenta k (the momentum label is omitted to simplify notation). The conduction and valence bands are given by Note that this basis choice implies that the Hamiltonian matrix elements are not periodic in the Brillouin Zone, H ij (k þ G)aH ij (k) with G a reciprocal lattice vector.
To compute the shift current, the direct use of equation (3) requires the evaluation of derivatives of Bloch functions, which can be difficult to compute numerically. Previous works 4,7,9 have addressed this problem with the use of identities that replace wavefunction derivatives with sums over all states of matrix elements of Hamiltonian derivatives. These identities are known as sum rules and rely on the fact that momentum and velocity operators are proportional in the plane wave basis p ¼ mv, which is not true in the tight-binding formalism. In this work we derived a generalized sum rule appropriate for tight-binding models (see Methods section), from which the integrand equation (3) can be evaluated for any two-band model in terms of the Hamiltonian derivatives only. The result is where the compact derivative notation f i;a @ ka f i and E ;b @ k b E is used. Equation (5) is one of the main results of this work. Several general principles to maximize the band edge shift current can be derived from this expression. A straightforward one is that, since this expression does not depend on E 0 , particle-hole asymetry does not influence the shift current at all. Therefore E 0 is set to zero from now on. The additional term that appears only for tight-binding models in this more general sum rule is f m f i,b f j,ab , which is absent in previous formulations. For a direct band gap, this term dominates the response exactly at the band edge, since to lowest order in k the first term always has constant contribution, while the second one is at least linear in k for any model due to the energy derivative E ;b . For this term to be finite, the three Pauli matrices in the Hamiltonian must have constant, linear and quadratic coefficients, in any order. Satisfying this low-energy constraint can be taken as another general principle in the search for materials with large shift current. More explicit guidelines can be obtained by considering an explicit low-energy model with a direct band gap at a time reversal invariant momentum. Expanding the Hamiltonian around it we get Time reversal symmetry H*( À k) ¼ H(k) prevents quadratic terms in s y , and we have taken the linear term to be in the x direction without loss of generality. Note this type of linear term requires the breaking of any C n rotation symmetry with n42.
The band gap of this model is E g ¼ 2E k¼0 . Evaluating equation (5) we get Also note that in order to have a non-zero shift current quadratic terms in s x or s z are required. In 2D, the fact that I xyy is in general non-zero means that the current need not be in the direction of the electric field polarization.
The shift current close to the band edge can now be obtained by substituting equations (7) and (8) into equation (2), which gives (9) provides an analytical formula for o close to the band edge for a very general class of models. This simple expression allows one to disentangle the contributions of the shift current integrand and the JDOS and hence to optimize them independently.
To maximize the response we therefore require band structures where the JDOS has a strong singularity. It is well known that in the one-dimensional (1D) case, the generic JDOS diverges as a square root, N(o)p(o À E g ) À 1/2 . 1D systems such as polymers or nanowires or systems in the quasi-1D limit will in general have a large response. In 2D, the band edge JDOS has a finite jump of N(o) ¼ (m x m y ) 1/2 /2p, where m i are the average effective masses for valence and conduction bands. A singular N(o) thus occurs in 2D when the inverse effective mass vanishes. In the effective model in Equation (6), this happens when d ¼ 0, which realizes what we may call a gapped semi-Dirac dispersion 63 , since the coefficients of s y and s x are linear and quadratic in momentum, respectively. In such a case we have N(o)p(o À E g ) À 1/4 (full expressions for N(o) may be found in the Methods section).
For materials with large JDOS, the current can be further enhanced by appropriately tuning the parameters in Equations (7) and (8). This is most easily discussed if these parameters can be related to microscopic lattice models. In the next section, we discuss tight-binding models for simple materials that realize the described types of band structures.
Material realizations and lattice models. As a realization of the 1D case, we consider ferroelectric polymers that break inversion symmetry such as polyvinylidene fluoride or disubstituted polyacetilene 39,40,64 . This system is described by the tight-binding model schematically shown in Fig. 1b, defined in terms of two types of hoppings, t 1 and t 2 , alternating on-site potentials ± D, and orbital centres at x ¼ 0 and x ¼ x 0 . With our choice of basis functions, the Hamiltonian is specified by f where a ¼ 10 Å is the lattice constant and the distance between closest neighbours is ref. 64 x 0 ¼ 0.48a. For estimates of the tight-binding parameters, we consider the example of disubstituted polyacetilene that was experimentally realized in ref. 39, with a band gap of 2.5 eV. For regular polyacetilene, where D ¼ 0, the hopping parameters and band gap have been estimated as ref. 64 Assuming the same hopping for the disubstituted version, we use D ¼ 1.0 eV to match the observed band gap. Note that the dispersion does not depend on x 0 .
Using Equations (2) and (5) we can now compute the shift current for this 1D model. Expanding about the low energy momentum k x ¼ p/a and performing a constant rotation of the Pauli matrices, we obtain an effective model as equation (6) with parameters k y ¼ 0 and To be able to compare the responsivity of these materials to that of a three-dimensional system, we consider a stack of polymers as depicted in Fig. 1a, separated by a distance d which we take to be equal to the lattice constant of the polymer d ¼ a.
The photoresponsivity is then k abb The typical photoresponsivity spectrum of this model with this convention is shown in Fig. 2a.
For the 2D case, we require a layered material that breaks both inversion and rotational symmetries. The most popular of the recently isolated 2D semiconductors break either inversion (BN, MoS 2 ) or rotational symmetries (black phosphorus 65 , ReS 2 (ref. 66)), but not both. An inversion symmetry breaking version of the strongly anisotropic black phosphorus, a group V element, can be obtained combining elements of the IV and VI groups.
These group IV monochalcogenides, such as GeS, are predicted to be stable in the monolayer form with the orthorhombic structure of black phosphorus 46,47 .
These materials can be described with a tight-binding model similar to the one used for black phosphorus [67][68][69] . While the GeS unit cell contains two Ge-S pairs at different heights, a unit cell with a single Ge-S pair can be used when the physics to be probed is insensitive to the heights of the atoms (see Methods for a detailed explanation). The two band Hamiltonian is specified by Ák Þ, and f z ¼ D. a 1 and a 2 are the lattice vectors. See Fig. 1d for the definition of the hopping integrals. Again note the dispersion is independent of x 0 . The specific values of the tight-binding parameters for GeS have been obtained by fitting an ab-initio calculation as described in the Methods section, where the coefficients of the low-energy model near the band edge are also shown. Note in this lattice structure there is a mirror symmetry y-À y, which is represented as the identity, and restricts a xy ¼ b xy ¼ 0. (This is so because both conduction and valence bands are even under the symmetry, as it also happens in black phosphorus. This is also the result of our ab-initio calculation.) This symmetry still allows a linear term of the form k x s y , crucial for the semi-Dirac type of band structure. In this model, the semi-Dirac limit is realized when We consider a stack of monolayers separated by d ¼ a, as shown in Fig. 1c. In this case, we consider an inert spacer layer between the GeS layers to avoid the restoration of inversion symmetry that would occur if we were to stack GeS into its natural bulk form. The three-dimensional photoresponsivity of this model, given by k abb 3D ¼ k abb 2D =d, is computed using equations (2) and (5). To make contact with the 1D case we consider a stacking distance d ¼ a(|a 1 | þ |a 2 |) 1/2 and x 0 ¼ 0.18a. The results are shown in Fig. 2b. We see that both k xxx and k xyy are in general finite, and the polarization average is also finite due to the strong anisotropy.
The response of the monochalcogenides is large because they are close in parameter space to the gapped semi-Dirac Hamiltonian. This is best illustrated by considering the evolution of a fictitious system where the hoppings are tuned (with t 3 ¼ 0 for simplicity) to the semi-Dirac case |t 1 |/t 2 ¼ 2, where the divergence of the response is clearly appreciated. This evolution is shown in Fig. 2c.
Further optimization. After describing the representative tightbinding models with large JDOS, we may now address a more systematic analysis of the photoresponsivity. First, we consider exploring the phase diagram of the monochalcogenides by sweeping |t 1 |, t 2 in parameter space while the band gap is fixed at 1.89 eV by choosing D appropriately and t 3 ¼ 0 for simplicity. Figure 3a shows the polarization averaged photoresponsivity, , for the parameters x 0 ¼ 0.18a and y ¼ 0.69. This phase diagram summarizes nicely the most physically relevant regimes where the shift current is large due to a divergent JDOS, namely the 1D limit where t 1 j j ( t 2 , and the semi-Dirac regime where t 1 j j$ 2t 2 . In this phase diagram, the point corresponding to t 1 and t 2 of GeS is shown as a white circle with blue outline. Next we illustrate a very important feature of the behaviour of the shift current integrand. Equations (7) and (8) depend generically on the hoppings and lattice parameters. The energy does not depend on the parameter x 0 , but the wavefunctions do. In Fig. 3b, we show the peak photoresponsivity as a function of |t 1 |/t 2 and x 0 . A large response is observed in the semi-Dirac limit t 1 j j=t 2 $ 2. However, a very strong dependence on x 0 and even a sign change is also observed. The dependence on x 0 dramatically illustrates the fact that the shift current depends not only on the band structure but also on the wavefunctions. This can be seen explicitly in the fact that the effective mass m À 1 x ¼ 4a 2 x t 1 t 2 =E g is independent of x 0 , but the combination v F a x appearing in the shift current integrand is not. In particular a x vanishes for , which means that regardless of the JDOS, the band edge response can actually be zero. This behaviour is characteristic of Berry connections, which depend explicitly on the positions of the sites in the unit cell.

Discussion
In this work, we have shown how an effective model for the band edge enables a clean separation of the two factors that contribute to a large shift current: the standard JDOS and the shift current matrix element. This model also allows us to readily identify materials with semi-Dirac-like Hamiltonians as those where both factors can be made large. Several other general conclusions can be drawn from the form of the effective shift current integrand in equations (7) and (8). First, since the 1/o 3 factor becomes 1/E 3 g at the band edge, materials with smaller gaps are expected to have larger shift currents. A second conclusion is that while looking for materials with large JDOS is a good guiding principle, the shift current integrand depends on other microscopic details that can change the response dramatically. Within our simple model, the shift current can be maximized by bringing the two sites of the unit cell closer together, which is a requirement that the monochalcogenides satisfy well. Materials that may perform even better than GeS may be searched for exploring different chemical compositions, alloying or by strain engineering.
Our results were made possible by the derivation of a new sum rule appropriate for tight-binding models. With this sum rule, our work can be easily extended to tight-binding models with more than two bands, or systems where the minimum direct gap is not at a time-reversal invariant momentum. We expect that the formalism developed here will provide the necessary link to combine ab-initio methods with effective models, allowing for more in-depth, systematic study of shift current photovoltaics. i i i i i ii i i i i i) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) Our results should be compared with known ferroelectric materials that have been recently studied. In the visible range of frequencies, op3 eV, we find peak values of 0.1 mAW À 1 in BiFeO 3 (ref. 30), 1 mAW À 1 in hybrid perovskites 13 and a maximum 10 mAW À 1 in BaTiO 3 (ref. 14) or NaAsSe 2 (ref. 15). The realistic materials that we propose present larger responsivities, with the additional advantage that the peak is by construction at the band edge. Moreover, as Figs 2c and 3b show, peak responses on the order of several hundreds of mAW À 1 could be achieved with materials closer to the semi-Dirac regime. To compare with conventional photovoltaic mechanisms, the total current per intensity of a crystalline Si solar cell exposed to sunlight is about 400 mAW À 1 (ref. 71).
Given these numbers, our work is a sign that shift current photovoltaics capable of surpassing conventional solar cells may be close at hand, and a push to investigate their full potential using methods discussed in this work-along with established techniques-is warranted. We believe that the simple principles derived in our work will serve as a guide for both theory and experiment in the development and optimization of the next generation of shift current photovoltaics.

Methods
Shift current. To make contact with previous work, we note the shift current integrand in equation (3) is sometimes expressed in terms of the phase of the inter-band matrix element r is known as the shift vector. The response to a natural light source such as sunlight, which is unpolarized, is obtained by averaging k abb over polarization. Taking E y ð Þ ¼ E j j cos y; sin y ð Þwe have Sum rule. The expression for the shift current presented in the main text can be obtained by the use of a sum rule for the quantity r a nm;b , which is obtained from the identity Evaluating both sides explicitly for nam, the identity can be expressed as where In the evaluation, we used v a nm ¼ ir a nm o nm : n 6 ¼ m ð16Þ The first equality follows from @ k n m j h i¼0 if man, while the last two follow from @ ka n H j jm h i¼ d nm @ ka E n . Note this sum rule contains the extra term w ab nm compared with ref. 9, where H ¼ p 2 /2m þ V(x) and w ab nm ¼ d nm d ab /m, which has no off diagonal component. Quite importantly, the term w ab nm in tight-binding models is the one responsible for all band edge contributions. Also note that it has been argued before that I xxx ¼ 0 for a two-band model 4 , which is actually only true if w ab nm ¼ 0.
Two-band model. For the case of two bands, m ¼ 1, n ¼ 2 the use of the sum rule for the shift current integrand in equation (3) leads to the simplified expression To evaluate this expression we compute the wave functions of H with n ¼ 1, 2, Z ¼ ( À 1) n , and f k ¼ arctan(f y /f x ). The required matrix elements are where the off diagonal matrix element s i ¼ c 1 s i j jc 2 h iis and the diagonal velocity matrix elements are computed from equation (15). The imaginary part in equation (17) can be taken using Im s Ã i s j Â Ã ¼ À P m E ijm f m =E and this leads to equation (5) in the main text.
Joint density of states. To compute the JDOS, we first start with the 1D case. Close to the band edge, we expand the energies of conduction and valence bands as x /2m x where the total effective mass m À 1 x ¼ jm 1;x j À 1 þ jm 2;x j À 1 is given by and solve for k o ð Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . Rescaling 2m x we get where we get the expected 1D singularity. For the generic 2D case, again we expand o 12 EE g þ k 2 x /2m x þ k y 2 /2m y , where m x is still given by equation (22) and We consider the case when m x 40, m y 40, so that the minimum does lie atk ¼ 0. By rescaling 2m x and 2m y we get in polar coordinates which is the expected constant result. Finally, the semi-Dirac case occurs in 2D when m À 1 y ¼ 0, which in the absence of second neighbour hopping occurs exactly at d ¼ 0. In this case, we keep the complete expression for o 12 In polar coordinates we have We now rescale a x , a y instead, solve for k and get Ab-initio calculation and tight-binding fit for GeS. Owing to the lack of tight-binding models for monochalcogenide materials 46,47 , we have derived the tight-binding parameters by fitting the electronic structure of GeS ab initio. We used the PBE 72 approximation to the exchange correlation functional, ultrasoft pseudopotentials 73 , Quantum-ESPRESSO 74 and Wannier90 (ref. 75) computer packages. The cutoff for electron wavefunction is set to 40 Ry and cutoff for electron density to 200 Ry. Internal coordinates and in-plane lattice constants were fully relaxed. Vacuum region between repeating images of GeS monolayers is 17 Å. Wannier functions were constructed from a 12 Â 12 regular k-mesh grid. The maximally localized Wannier functions were constructed in a standard way by projecting into hydrogenic s-like and p-like orbitals on both Ge and S atoms along with two s-like orbitals in the vacuum region that are needed to represent the vacuum states. The frozen window for the disentanglement procedure spans up to 6.2 eV above the Fermi level. The crystal structure of GeS is orthorombic with space group Pnma (No. 62) and lattice vectorsl 1 ¼ (l 1 , 0) andl 2 ¼ (0, l 2 ), with l 1 ¼ 4.53 Å and l 2 ¼ 3.63 Å and contains two Ge and two S atoms. The structure can be seen as two GeS zigzag chains separated by a height of h ¼ 2.32 Å. The ab-initio results for the conduction and valence bands near the G point are shown in Fig. 4 and have mostly p z character. This system can be effectively described with a two site tight-binding model. This can be done because the lattice structure has glide symmetries with mirror reflection z-À z and translationsã 1 ¼ (a x , a y ) andã 2 ¼ (a x , À a y ), with a x ¼ l 1 /2 and a y ¼ l 2 /2. When the out of plane positions of the atoms are not relevant for the problem of interest, one can define a smaller two site unit cell where the glides play the role of lattice vectors (as it is done in black phosphorus 69 ). The Ge and S sites in this effective tight-binding model are located at (0, 0) and (x 0 , 0), with x 0 ¼ 0.62 Å. This is the tight-binding model employed in the main text. The parameters of this model are obtained from the ab-initio calculation as follows.
Since our aim is to model faithfully only the low-energy bands around the Gamma point, it will suffice to consider a single p z orbital per site in the tightbinding model. The minimal model parameters are the on-site potential difference D between Ge and S p z orbitals and the three nearest neighbours hoppings t i , with i ¼ 1, 2, 3, which are all between Ge and S atoms. In addition, to reproduce the small particle-hole asymmetry of the gap, we also consider two further neighbour hoppings t 0 1 and t 0 2 , which connect Ge-Ge or S-S pairs (we assume the same values for both species to simplify).
The tight-binding Hamiltonian takes the form H ¼ E 0 þ S i s i f i (k) with coefficients where, as defined in the text, F k ð Þ ¼ e ia1Ák þ e ia2Ák À Á . Our tight-binding fit is intended to reproduce faithfully the bands and wavefunctions close to the band edge, where the effective low-energy model applies. This model is given by where a constant term is omitted as it can be absorbed in the chemical potential. The effective model parameters are related to the tight-binding parameters as The key to obtain a reliable tight-binding parametrization is that, since the shift current depends sensitively on the actual wavefunctions, the tight-binding model should be fitted to wavefunction-dependent quantities in addition to the band energies. The simplest gauge invariant quantity that depends on wavefunction phases is the bracket of two covariant derivatives with D m ¼ q m À iA m , with A m ¼ i u k j@ m u k the Berry connection. The real and imaginary parts of this tensor are known as the Berry curvature and the quantum metric. A fit that reproduces this tensor correctly in addition to band energies ensures that the wavefunction structure around the G point is correctly accounted for, so that any other gauge invariant quantity computed in the effective model should be the same as that computed ab initio.
The Berry curvature O(k) is defined as O k ð Þ ¼ E mn Im @ m u k @ n u k j Â Ã ¼ rÂA: The Berry curvature around G for the tight-binding model is given by Since O vanishes at the origin, we take q y O as one extra input for the fit. The quantum metric is defined as The only non-vanishing component of the quantum metric at k ¼ 0 is given by so we take g xx as another extra input for the fit.
In summary, we take as ab-initio input parameters the gap, the four effective masses and the lowest order Berry curvature and quantum metric, q y O and g xx . The difference in effective masses for electron and hole bands, accounted for the term E 0 , can be fitted independently with the hoppings t 0 1 and t 0 2 . Since E 0 has no impact in the shift current response, the hoppings t 0 1 and t 0 2 are not considered in the main text. The rest of the input is fitted with t 1 , t 2 and t 3 , the on-site potential D and x 0 , and the results of the fit are shown in Table 1. While x 0 is in fact known from the lattice structure of GeS to be 0.62 Å, obtaining it independently from the tight-binding fit, which gives a close value of 0.52 Å provides an additional check of the validity of the model.
Data availability. The data that support the findings of this study are available from the corresponding author upon request.