Unraveling materials Berry curvature and Chern numbers from real-time evolution of Bloch states.

Materials can be classified by the topological character of their electronic structure and, in this perspective, global attributes immune to local deformations have been discussed in terms of Berry curvature and Chern numbers. Except for instructional simple models, linear response theories have been ubiquitously used in calculations of topological properties of real materials. Here we propose a completely different and versatile approach to obtain the topological characteristics of materials by calculating physical observables from the real-time evolving Bloch states: The cell-averaged current density reveals the anomalous velocities that lead to the conductivity quantum. Results for prototypical cases are shown, including a spin-frozen valley Hall and a quantum anomalous Hall insulator. The advantage of this method is best illustrated by the example of a quantum spin Hall insulator: The quantized spin Hall conductivity is straightforwardly obtained irrespective of the non-Abelian nature in its Berry curvature. Moreover, the method can be extended to the description of real observables in nonequilibrium states of topological materials.

When the Hamiltonian of a system is subject to adjustable periodic parameters, the eigenstates can acquire a nontrivial gauge-independent phase over the adiabatic evolution in the parameter space (1, 2). In various areas of physics the presence and importance of such geometrical phases have been recognized and characterized through the formulation attributed to Berry (1) and Chern and Simons (3). In particular, the geometrical phase resulting from a variation of the Bloch vector () of a periodic Hamiltonian () has recently attracted the most attention in the condensed-matter community (4–6). This Bloch geometrical property can also be discussed by a local field, known as Berry curvature (), defined on the parameter space (2), instead of the loop-integrated phase. Once a finite Berry curvature is present, a solid will attain an intrinsic anomalous Hall conductivity and, in this regard, the Berry curvature can be thought of as an intrinsic magnetic field of the material (2, 5, 7). Furthermore, as established by Chern and Simons (3), the integration of the Berry curvature over a closed surface is quantized, named a Chern integer, which has been recognized as an essential indicator of the band topological nature (8). For example, for 2D Bloch electrons, the Chern number can be determined by the integration of the Berry curvature over the whole Brillouin zone (BZ): .

In addition to the anomalous charge Hall conductivity, various other transport properties can reflect the effects of the Berry curvature, such as the valley Hall or spin Hall effect and others (7, 9–11). Except for instructionally designed simple models, the calculation of Berry curvature and the Chern number for real materials requisitely involves a perturbative approach within linear response theory, which in most cases requires the Wannierization technique to cope with the fine-grid integral over the BZ (12, 13). In the linear response formulation, a monochromatic field () is introduced to perturb the ground state, and the off-diagonal element () of the conductivity tensor is written in terms of the unperturbed eigenstates, by means of the Kubo formula (2, 9). In particular, the static Hall conductivity () has been assigned the most significant physical meaning owing to the work of Thouless et al. (8): The quantized Hall conductivity of the Landau-level insulator is exactly related to the integer indicator of the topology ( as ). The intrinsic Chern number of a material, as opposed to that of a quantum Hall state formed by an external magnetic field, was pursued in the search for quantum anomalous phases (14), and Haldane’s (15) suggestion for the spin-frozen Chern insulator ignited flourishing interest in topological states of matter: Kane and Mele (16, 17) suggested the topological quantum spin Hall (QSH) effect by considering two copies of Haldane’s model, constituting the time-reversal partner to each other. Bernevig, Hughes, and Zhang (18) presented the same concept of a QSH effect through the exploration of an intrinsically band-inverted semiconductor and its surface state (19). Extension to 3D and the relevant topological numbers has been pursued, and the topologically protected metallic surface bands at the boundary of distinct topological states have been observed and discussed in the perspectives of spin-resolved dissipationless carrier transport (4, 20–22). More comprehensive overviews of the topological phases of materials, and their historical developments, have been provided in a few recent review articles (2, 4, 22).

This work is motivated by the question of whether the Berry curvature can be extracted from physical observables, such as the total current or others, that can be obtained from the time-evolving states of a solid () through the time-dependent Schrödinger equation. For actual time propagation, in the present work, we perform first-principles time-dependent density functional theory (TDDFT) calculations, and we show that the Berry curvature and Chern numbers for bands can be obtained from the time profiles of the total current. This method does not require the commonly used Wannierization technique, and physical observables (charge or spin current) are directly computed as expectation values from the time-evolving solid states. Exemplary real-material systems are presented, including a trivial insulator, a valley Hall insulator, and a quantum anomalous Hall insulator (QAHI). Using an example of a quantum spin Hall insulator (QSHI) we explicitly demonstrate that this dynamical approach offers a natural standpoint for spin Hall conductivity when the Berry curvature is non-Abelian and thus gauge dependent.

Theoretical Framework and Computation Method

The computation method consists of the use of a spatially uniform electric field (E field) in the form of a time-dependent vector potential: . Time-evolving wavefunctions are computed through the time-dependent Kohn–Sham (KS) equation, as derived within TDDFT (23, 24),where indicates the scalar-type external potential including atomic pseudopotential, and and are Hartree exchange-correlation density-functional potentials and the spin-orbit coupling term, respectively. The Kohn–Sham wavefunctions are two-component spinors or single-component wavefunctions (the spin-frozen or spin-polarized ones). Detailed parameters related to standard DFT calculation (25) and the time-integration algorithms are described in .

Once the time-evolving solid states () are computed, the expectation values of an observable can be directly evaluated. For example, the real-time–dependent cell-averaged current density can be obtained by taking the expectation value of the gauge-invariant mechanical momentum,where and are the occupation factor and the volume of the D-dimensional unit cell (e.g., the area of the 2D unit cell is denoted by hereafter), respectively. Here the gauge-invariant mechanical momentum is defined as , to which the velocity operator is related by . In this computation, the angular momentum-dependent atomic potential is described by a separable form of a nonlocal pseudopotential (), and thus the nonzero commutator with the position operator () should be taken into account, as explained elsewhere (26). Later for the QSHI, our main interest will be focused on the spin current operator () (see below) instead of the velocity operator () used in Eq. . If the parameters for the Hamiltonian, such as the E-field strength, are suitably chosen to make the adiabatic evolution, up to a first-order variation of the wavefunctions, the velocity expectation can be written as (2)where is the instantaneous band energy, and is the Berry curvature with being the periodic part of the time-dependent KS wavefunctions . In the adiabatic regime, the wavefunction (), the band energy (), and the Berry curvature () converge to the corresponding stationary values of the Bloch state with the Bloch vector of . To ensure that the time-propagating electron states adiabatically follow the ground-state band structure, the E field in the present calculations is smoothly turned on over a period of 10–25 fs, and the strength of the E field () is well below the Zener tunneling criterion (27). Since the integral of the group velocity () over the BZ is zero, the integrated Berry curvature of the band can naturally be deduced from the integrated velocity: .

Up to now, most first-principles topological studies of real materials have used the Wannierization method using maximally localized Wannier functions (MLWFs) to compute the Berry curvature and to derive the topological numbers (12). From the standpoint of equality of Thouless et al., the topological nontrivialness of the insulators can be discussed in terms of physical quantity, such as charge or spin Hall conductivity (28, 29). However, such trials have mostly been hampered by the nonintuitiveness of spin Hall current, particularly when the Hamiltonian is furnished with a measurable spin-orbit coupling (SOC) of real materials (16, 17, 29). Here, as the dynamical approach, noted in Eq. , provides a better computational versatility in dealing with experimentally accessible observables, we focus on the direct evaluation of the physical gauge-invariant velocity operator, as given in the left-hand side of Eq. . This method requires only the occupied orbitals and thus avoids the large sum over states of the response formula (13). The required fine sampling of the BZ is here realized by the adiabatic sweep with a fine time step. Furthermore, even when a strong field perturbs the system, expectation values can be evaluated, which enables us to capture the Berry curvature effects in the nonadiabatic or nonequilibrium regime.

Results and Discussion

Trivial Atomic Insulator.

To demonstrate the efficiency in obtaining topological characteristics of materials from real-time evolving band states, we first consider the time evolution of a trivial insulator. As a representative example of an atomic insulator, here we devised a Gedanken system by placing a He atom in the simple cubic (SC) lattice with a given fixed lattice constant of = 3 Å, as schematically depicted in Fig. 1, Inset. The valence band is derived from the atomic orbital, and the band dispersion is presented along the x direction in Fig. 1. In this system, simultaneous presence of time reversal and inversion symmetry enforces that at every point for the spin-frozen bands. Onto this self-consistently converged ground state, we applied a weak static electric field along the x direction, Å, and then the KS states were allowed to evolve over time. The selected KS states are depicted with symbols in Fig. 1 in units of : , , , and . The time-evolving band energy of these KS states is presented in Fig. 1 up to 200 fs. This model of an atomic insulator by construction possesses a large band gap of eV, and the real-time evolutions of the KS states strictly follow the ground-state energy surface. Fig. 1 demonstrates that these evolving Bloch states () move across the BZ along the ground-state band with a time-evolving Bloch vector: given by .

Fig. 1.

Time propagation of a KS band of a 3D atomic insulator in response to an applied E field. (A) The static band structure of the solid with the He atom in a simple cubic cell. (B) Time variation of the band energy of four selected KS states. (C) Time variation of longitudinal () and transverse ( and ) velocities of the KS state starting from the point. A, Inset depicts the unit cell with the indication of the E field. In A and B, the four selected Bloch states are denoted by a square, a circle, a triangle, and a diamond.

The time profile of the velocity of the state starting from the point (denoted by the circles in Fig. 1 ) is shown in Fig. 1. It is noteworthy that the calculated transverse velocities always remain at zero (), which is consistent with the absence of Berry curvature in this system. The other three sampled Bloch states also carry completely vanishing transversal velocities, as summarized in . The calculated longitudinal velocity, using the formula given in the left-hand side of Eq. , is indeed consistent with the instantaneous curvature of the band at : . This periodic oscillation of the longitudinal velocity, with the period of fs as shown in Fig. 1, originates from the periodic nature of the energy band in the reciprocal cell, which can be compared with the Bloch oscillation (30).

Valley Hall System: Inversion Symmetry-Broken Graphene.

We now investigate a system which possesses a locally nonvanishing Berry curvature, but the band Chern number is absent owing to the presence of time-reversal symmetry. The simplest example in this perspective can be achieved from the spin-frozen bands of the graphene by intentionally breaking its inversion symmetry. To that end, we performed standard DFT+U calculation by adding an asymmetric U potential of eV and 0 eV for A and B sites, respectively (Fig. 2, Upper). As a result, the Dirac cones in the K and K′ valleys develop a band gap of eV, as shown in Fig. 2, Upper. The Berry curvature of this artificially inversion-broken graphene band is calculated and presented in Fig. 2, Lower (7). Onto the self-consistently converged ground state, we applied a constant and uniform static E field along the x direction ( V/Å) and performed the time propagation. The E field was gradually turned on over the initial 20 fs (). The evolutions of three selected k points near the K valley (K and ) are displayed in Fig. 2: The initial location of these three points is denoted in the BZ in Fig. 2, and the time profiles of the band energies [the valence band maximum (VBM) and the conduction band minimum (CBM)] are summarized in Fig. 2. Note that the state starting from the point (red line in Fig. 2) arrives at the top of the K valley at fs. This time evolution of band energies is schematically illustrated on the cone surface in Fig. 2.

Fig. 2.

Time propagation of KS states in an inversion symmetry-broken graphene system. (A) The static band structure (Upper) and the Berry curvature (Lower). (B) Schematics of the BZ and the Dirac cone with three selected k points in the K valley with , where is the real-space lattice constant. (C) Time variation of the band energies of the VBM and CBM states at the selected three k points in the K valley. (D) Time profile of the longitudinal velocities of the VBM states of the three k points in the K valley. D, Inset shows the same longitudinal velocity of the point in the K′ valley. (E) Time profile of the transverse velocity of the VBM states starting from the K + δ and K′ + δ points. (F) The same as E for the state starting from the K + δ point with various strengths of the E field: . Here, indicates the strength of the E field used for C–E. The critical E field in the Zener tunneling model () is defined in the text. In A, Inset, atomic symbols schematically depict the inversion-broken sublattices of two C atoms.

The time evolution of the longitudinal velocities of these three selected VBM states in the K valley is presented in Fig. 2. The calculated velocities coincide with the instantaneous slope of the band energy dispersion: For example, the states starting from the exact K point have zero initial velocity. It is noteworthy that, after around 70 fs, the velocities of all three states reveal a similar linear trend. After this time, all three states move down the linear surface of the Dirac cone. We also performed the same calculation for the point in the K′ valley, which showed the same longitudinal behavior as that in the K valley: Fig. 2, Inset shows the velocity of the state starting from K′ + δ, which is almost identical to that starting from K + δ. On the other hand, the transverse velocities of the states in the K valley are contrastingly different from those of the states in the K′ valley. Fig. 2 shows that the transverse velocities of the states starting from and K′ + δ points have opposite sign. This demonstrates that the spin-frozen bands of the inversion-broken graphene can indeed reveal the valley Hall transport: The carriers in the K valley deflect with an opposite anomalous velocity from that in the K′ valley (7). In Fig. 2 , we demonstrate that the two states arrive at the top of the valley at around 46 fs, leading to zero-longitudinal velocity. At this point the instantaneous velocities consist only of the anomalous velocity given by , which directly monitors the Berry curvature on the peak point (7). The Berry curvature obtained this way is , which is comparable to the value of obtained from the Kubo formula with Wannierized bands (Fig. 2, Lower) (13).

This example thus suggests that the real-time propagation of the KS equations can indeed be used to evaluate the Berry curvature, Hall conductivity, and consequently the band Chern number, by choosing a reasonable parameter, such as a modestly weak E field. To make the time-evolving states follow well the adiabatic energy surfaces, the strength of the E field needs to be sufficiently weak compared with the band gap (). This feature can be discussed quantitatively in terms of a Zener tunneling model (27): The critical strength of an E field to induce tunneling across the band gap is estimated to be , where is the reduced effective mass between the CBM and the VBM. The E field used for Fig. 2 ( V/Å) corresponds to . In Fig. 2, we compare the simulation results with four different strengths of the E field: with . It is noteworthy that, for , an abrupt oscillatory feature emerges as the Bloch state passes the Dirac cone region, which indicates the hybridization of CBM and VBM states and consequently a breakdown of the adiabaticity.

Quantum Anomalous Hall Insulator.

Our method of time propagation can also be directly applied to the QAHI case. Since the first experimental realization of the QAH phase on a magnetized topological insulator (14), various studies ensued to find a new material that can preserve the intrinsic quantum Hall conductivity in an elevated temperature. Several works, in this perspective, have focused on the possible transitions of a 2D hexagonal lattice from the QSH to the QAH (21, 31). One of the most intriguing examples in this direction is the single layer of Sn (32), named stanene, with various hydrogen or halogen coverages: The full coverage of halogens reveals a QSH phase, whose band structures can be described with the Bernevig–Hughes–Zhang (BHZ) model Hamiltonian, while the ideal half coverage on one side results in a QAH phase (31). This half-hydrogenated Sn (HHS) is depicted in Fig. 3. Before the evaluation of the time-evolving state, we calculated the static ground-state band structure and also the Berry curvature using MLWFs, as summarized in Fig. 3. In agreement with the previous study, HHS exhibits a highly localized peak in the Berry curvature () around the point, which is integrated to a single quantum of the Chern number () (31).

Fig. 3.

Time propagation of the KS states of the ideal half-hydrogenated Sn in a quantum anomalous Hall phase. (A) Schematic geometry. (B) The band structure and Berry curvature. (C) The time-averaged Hall conductivity calculated from the time-evolving states. Red and blue colors in lines in B represent the texture of the spin polarization in the positive and the negative z direction, respectively. C, Inset shows the real-time profile (gray) and the time-averaged profile (red) of the Hall conductivity after the turning-on period.

To evaluate the charge Hall conductivity from the expectation value of the velocity operator, as noted in Eqs. and , we applied a static E field in the x direction and performed the time propagation of the uniformly sampled Bloch states. Since the band gap in this case is quite small, to achieve a reasonable adiabatic quality, we chose smaller strengths of the E field which was turned on gradually over a 25-fs period. The transverse current is calculated from the time-evolving Bloch states, and the Hall conductivity is given byThe calculation results are presented in Fig. 3, which shows that, after the initial turning-on period ( fs), the time-averaged profile converges well to the single quantum of the conductance. In this QAH phase of HHS, the Berry curvature is not dominated by a single band, but rather distributed over multiple bands around the point (31), where the inverted VBM and CBM states are intricately hybridized. However, the Hall current calculated as a sum of occupied valence bands reveals a good convergence. Detailed electronic band structures for HHS are given in .

Quantum Spin Hall Insulator.

We now extend our example of the time-evolution study to a QSHI. As an example, we calculated the single-layer bismuthane in which the hydrogenated Bi atom locates in the hexagonal lattice (21, 33), as shown in Fig. 4. When the SOC is intentionally turned off, the band structure reveals the typical Dirac cone structure of a honeycomb bipartite lattice, as depicted by the dashed line in Fig. 4. When we include the SOC of Bi atoms, the bands near the Dirac cones open a gap of 0.8 eV at K and K′ points, presented as solid lines (Fig. 4). The optimized geometry of the layer has small buckling but preserves the inversion symmetry, and thus the spinor states in the K and the K′ valley constitute the time-reversal and inversion partners, enforcing the spin-up and spin-down states to be degenerated over the whole BZ (16, 21). Results for a similar system without inversion symmetry are presented in . In Fig. 4, the valence bands are labeled from to , where the subscript up and down arrows (e.g., ) indicate the doubly degenerate bands. It should be noted that the spin in this case is not well polarized in the z direction, but has varying textures depending on the k points as a result of the strong SOC.

Fig. 4.

Time propagation of the KS states of bismuthane in a quantum spin Hall phase. (A) Schematics of the geometry and the BZ. (B) Band structures with (solid line) or without (dashed line) SOC. Doubly degenerate valence bands are labeled from to from the highest band. (C) The spin Hall conductivity and the charge Hall conductivity calculated from the time-evolving states. (D) The band-resolved contributions to the spin Hall conductivity from the doubly degenerate second valence band () and from the others (). C, Inset depicts the bias static E field and the spin-resolved current.

We computed the electron dynamics of this QSHI system by gradually ramping the E field over fs toward a static value of V/Å. As we above enumerated the charge current as an expectation value from the time-evolving wavefunctions, the spin value and the spin current can also be computed as an expectation value from these spinor wavefunctions. There have been many theoretical works to devise a better form of the spin current operator (9, 34–36). The definition and the practicality of the operator used in the present work are explained in detail in . Using the obtained spin current, the spin Hall conductivity is evaluated as

Fig. 4 depicts the time profile of the spin Hall conductivity and the charge Hall conductivity.

Remarkably, even though is not a good quantum number and varies over different k points (35), the spin Hall conductivity converges well to the conductivity quantum just after the initial turning-on period of the E field (). On the other hand, the charge Hall conductivity consistently remained at zero, which is a natural outcome of the time-reversal symmetry of the system. As depicted in Fig. 4, Inset, the longitudinal current of the insulator is summed to zero, but the Hall current of each spin sector was directed oppositely, as a result producing a finite spin Hall current.

The band-resolved contribution of each doubly degenerate valence band (denoted by , , , and in Fig. 4) to the spin Hall conductivity is presented in Fig. 4. We observe that the QSH phase of bismuthane is dominated solely by the second valence band (), while the effects of the others (, , and ) are marginal. This can be explained by the structure of the band inversions, which can be inferred from the comparison between the bands with and without SOC. Bands and are inverted at K and K′, and the VBM band () has two points of band inversion: one with the CBM band at K and K′ and the other with the band at . This double inversion renders the band topologically trivial, and the full occupation of the inverted pair of and bands makes their topology (the spin Chern number) cancel each other. As a result, the band , which is doubly degenerate and is solely inverted from the band at , remains as a single source of the quantized spin Hall state.

Here, we discuss more in depth the structure of the Berry curvature of the QSHI state. The degeneracy in the valence bands requires the Berry curvature to be defined as a matrix constructed in the degenerate subspace, putting in evidence its non-Abelian character (2, 36). The components of the Berry curvature matrix () and of the spin matrix () depend on the choice of the basis (gauge dependent) and thus cannot be observed in an experiment. Physical observables are gauge independent, and it is shown that the charge Hall and the spin Hall current can be written in terms of Tr[] and Tr[], respectively (36). In this regard, our approach based on the expectation value of the physical observables (charge or spin current operator, as given in Eqs. and ) is advantageous over other methods: The spin Hall current, as presented in Fig. 4, does not require any additional cost or separate treatment due to the degeneracy. Moreover, since the expectation value of an operator is independent of the unitary rotation within the degenerate subspace, the spin current can be written in terms of the velocity given by the diagonalized basis at each k point:The form of the last expression of Eq. is valid regardless of whether the spin is well defined as a good quantum number throughout the whole BZ or varies over k points (component-wise detailed derivation is provided in ). This suggests that the spin Hall conductivity, as written in Eq. , can be used to identify the spin Chern number, exactly in the same way as introduced in the simplified Kane–Mele model (16, 17), even when a strong SOC demotes the spin from the status of a good quantum number of the system (18, 29). In , we present the same results for an inversion symmetry-broken bismuthane, in which the degeneracy is modestly lifted, but the calculated spin Hall conductivity converges well into the quantized value, identifying its QSH phase. This strategy to deal with the 2D quantum spin Hall conductivity can be applied to the time-reversal planes in the BZ of a 3D topological insulator, such as () and () (37), allowing us to discern between strong and weak 3D topological insulators ().

Conclusions and Outlook

In summary, we demonstrated here that the calculated total time-dependent current from the real-time propagation of the Kohn–Sham Bloch states provides a completely alternative method to explore the topological character of solids. Results for exemplary cases were presented, including a trivial atomic insulator, a valley Hall system, and a quantum anomalous Hall system. On a prototypical example of the quantum spin Hall insulator, we discussed that this direct evaluation of the physical observables can serve as a natural platform for an adequate description of the non-Abelian Berry curvature. We highlight that the concept suggested here is not necessarily limited to the DFT-based single-particle scheme. One-body physical observables (such as charge or spin current), derivable from the time-evolving many-body states, can be used to gauge the anomalous behaviors rooted in the geometrical phase structures of the quantum mechanical wavefunctions. Moreover, the flexibility of the proposed computation scheme is clearly advantageous to deal with general time-dependent perturbations (38–40). To illustrate this feature, in , we present transiently emerging Hall current in a graphene nanoribbon when a circularly polarized driving force is applied as an externally driven time-reversal breaking mechanism (41, 42).