Anisotropic Excitons and their Contributions to Shift Current Transients in Bulk GaAs
Shift current transient are obtained for near band gap excitation of bulk GaAs by numerical solutions of the semiconductor Bloch equations in a basis obtained from a 14 band kp model of the band structure. This approach provides a transparent description of the optically induced excitations in terms of interband, intersubband, and intraband excitations which enables a clear distinction between different contributions to the shift current transients and fully includes resonant as well as off-resonant processes. Using a geodesic grid in reciprocal space in our numerical solutions, we are able to include the electron-hole Coulomb attraction in combination with our anisotropic three-dimensional band structure. We obtain an excitonic absorption peak and an enhancement of the continuum absorption and demonstrate that the excitonic wave function contains a significant amount of anisotropy. Optical excitation at the excitonic resonance generates shift current transients of significant strength, however, due to the electron-hole attraction the shift distance is smaller than for above band gap excitation. We thus demonstrate that our approach is able to provide important information on the ultrafast electron dynamics on the atomic scale.
pacs:72.40.+w, 78.47.J-, 78.55.Cr
The optical excitation of non-centrosymmetric crystals can be used to generate photocurrents on ultrafast time scales without the need of an external bias. As shown by Sipe and co-workers, the lack of inversion symmetry in zincblende III-V semiconductors results in a non-vanishing zero-frequency second-order optical susceptibility which corresponds to photocurrents that can be generated by optical excitation with a single frequency.Sipe and Shkrebtii (2000); Nastos and Sipe (2006) One can distinguish between three types of photocurrents: (i) injection currents originating from non-symmetric electronic distributions in k-space after resonant above band gap excitation, (ii) shift currents which are due to the spatial motion of optically excited carriers in real space after above band gap excitation, and (iii) rectification currents that result from the non-resonant polarization generated for below band gap excitation.
Here, we investigate the ultrafast dynamics of bulk GaAs following the near-band gap excitation by femtosecond laser pulses and focus on analyzing shift currents which are responsible for the bulk photo-voltaic effect. Experimentally, shift currents have been investigated in bulk semiconductorsZhang et al. (1992); Dalba et al. (1995); Chen et al. (1997); Côte et al. (2002); Somma et al. (2014) and semiconductor quantum wellsLaman et al. (2005); Bieler et al. (2006, 2007); Priyadarshi et al. (2009); Duc et al. (2016). Previous theoretical research on shift currents was mainly performed in the frequency domain using a perturbative analysis of the light-mater interaction to derive analytical expression for the considered nonlinear response in terms of matrix elements and resonance denominators.Nastos and Sipe (2006); Young and Rappe (2012); Young et al. (2013); Král (2000); Brehm et al. (2014) Using the semiconductor Bloch equations (SBE) in the basis kp wave functions it is possible to obtain photocurrents directly in the time-domain as was shown for the case of injection currents Priyadarshi et al. (2010); Duc et al. (2010, 2011) as well as for shift and rectification currents Podzimski et al. (2015, 2016); Duc et al. (2016). This approach provides a transparent description of the optical excitations in terms of interband, intersubband, and intraband excitations, allows to treat the light-matter interaction non-perturbatively, and provides good agreement with experimental results on injection and shift currents of GaAs based quantum well systems, e.g., Priyadarshi et al. (2010); Duc et al. (2016).
Due to the tremendous numerical demands, excitonic effects have been neglected in most previous theoretical investigations of photocurrents in semiconductors. Whereas this is justified for high above band gap excitations where excitonic effects have negligible contributions, for near band gap excitations, the many-body Coulomb interaction, in particular the electron-hole attraction, strongly modifies the optical response and therefore needs to be incorporated into the theoretical approach. Often when many-body effects are considered simplified models for the band structure and electronic states, e.g., isotropic and/or parabolic band structures and a small number of bands, are used which significantly reduces the numerical requirements when solving the SBE.Zimmermann and Trallero-Giner (1997); Rudin and Reinecke (2002); Bhat and Sipe (2005); Turkowski and Ullrich (2008); Turkowski et al. (2009); Haug and Koch (2004) However, since the photocurents in general originate from anisotropies of the band structure and/or of the optical matrix elements, they cannot be described properly by models using isotropic band structure models. Furthermore, shift currents, in particular, involve non-resonant excitations and can therefore not be adequately described in models that consider only bands that are present at or near the band gap.Podzimski et al. (2015); Duc et al. (2016)
For the theoretical analysis of excitonic contributions to shift currents in bulk GaAs we therefore apply our approach of the SBE formulated in a basis of electronic eigenstates obtained from a 14 band kp model. The anisotropic electronic band structure and the electronic states are well suited to describe shift currents transients.Podzimski et al. (2015, 2016); Duc et al. (2016) For the incorporation of excitonic resonances we extend our previous analysis and include the Coulomb interaction in time-dependent Hartree-Fock approximation. Due to the anisotropic band structure and matrix elements and the need to incorporate off-resonant excitations, which means that the rotating-wave approximation cannot be applied, the solutions of the resulting SBE are numerically very demanding. To obtain converged results we employ a geodesic grid in k-space which provides results of reasonable accuracy with a significantly smaller numerical effort than a cartesian grid.
In Sec. II we describe the fundamentals of our theoretical approach and the idea behind the use of the geodesic grid in our simulations. The excitonic absorption spectrum and the anisotropic exciton wave function are presented in Sec. III.1. Numerical results on shift current transients including excitonic effects are shown and discussed in Sec. III.2.
Ii Theoretical Approach & Numerical Challenges
The kp-based extended Kane model is used for the calculation of the semiconductor band structure. The extended Kane model is represented by a 14-band Hamiltonian,
which describes the band structure of zincblende crystals near the -point, GaAs being one prominent example.Mayer and Rössler (1991); Winkler (2003); Terbin et al. (1979) The Hamiltonian includes the split-off band , the highest valence band , the lowest conduction band , and higher conduction bands and . The band structure is obtained by solving the eigenvalue equation
which is accomplished by a matrix diagonalization. The coupling between the lowest conduction band and the higher conduction bands and ,
is responsible for the shift current and consequently the respective bands have to be included in the simulations.Podzimski et al. (2015) At the band gap is while the distance between the valence band and the higher conduction band is .Aspnes and Studna (1973); Cardona et al. (1988) Thus when numerically solving the SBE significantly small time steps have to be used to resolve the rapid oscillations arising from these energetic differences of the involved bands.
The time evolution of the photoexcited system is described by the SBEDuc et al. (2010); Sternemann et al. (2013); Duc et al. (2016), i.e., the Heisenberg equations of motion for representing the coherences () and the occupations () of the system in k-space, which read
with and being the energies of band and at point in reciprocal space. The generalized Rabi-frequency
contains the light-matter interaction, here written in the velocity gauge , and the Coulomb interaction in time-dependent Hartree-Fock approximation.Duc et al. (2010) In the velocity gauge the light-matter interaction is described by the electrodynamic vector potential and the momentum matrix elementsChang and James (1989); Enders et al. (1995)
The Coulomb matrix element is given byDuc et al. (2010)
with , the dielectric constant of GaAsRössler and Strauch (2001), and is the volume of the reciprocal unit cell.
So, the Coulomb matrix elements depend on four band indices () and two three-dimensional k-space vectors (). Evaluating the Coulomb matrix for all 14 bands results in possible combinations of band indices. This number can be drastically reduced by considering only the resonantly excited bands near the band gap, i.e., and , which leaves band index combinations but still is sufficient to properly describe the excitonic absorption. Additionally, the Coulomb matrix couples all k-space vectors and which leads to a numerical effort when calculating the matrix elements and solving the SBE which grows quadratically as function of the total number of -points . Using a standard Cartesian grid with , with being the number of k-points in one direction, results in a numerical effort that grows with . To achieve converged results an unreasonably high memory and computer time would be required. A favorable alternative to a Cartesian grid is the so-called geodesic grid which is a spherical grid that is often used in climate simulations.Heikes and Randall (1995); Mahdavi-Amiri et al. (2015)
The kp band structure Hamiltonian can written as a sum which distinguishes the involved symmetriesLipari and Baldereschi (1970)
For small -vectors the spherical -terms have the largest contribution to the Hamiltonian while for large -vectors higher order -terms representing cubic and tetrahedral symmetry become more relevant. Thus at the -point the energy contributions to the total Hamiltonian can be ordered asLuttinger (1956)
which is the reason why the band structure at the -point can be reasonable well be approximated by a parabolic dispersion and the analytic solutions of Wannier excitons, which have a spherical -symmetry known from the hydrogen atom, are in good agreement with experimental results.Wannier (1937); Luttinger (1956) The application of a spherical grid takes advantage of the predominant spherical symmetry at the -point reducing the amount of wasted simulation space and consequently saving computational effort. In a spherical grid based on a spherical coordinate system, , the point density on a sphere surface is inhomogeneous and increases towards the pols. This numerically undesirable property of spherical coordinates is circumvented in a geodesic grid in which the sphere’s surface is constructed almost completely from equilateral hexagons and thus has a very homogeneous point density. The total number of points in a geodesic grid is with being the number of points per sphere and the number of spheres. Consequently, using the geodesic grid the numerical evaluation of the SBE with Coulomb interaction can be reduced from a scaling to a more moderate scaling which is a huge improvement regarding the required numerical resources. Further details can be found in Podzimski (2016).
Iii Numerical Results
Here, we present and discuss results on the excitonic absorption and the exciton wave function in Sec. III.1 and on the shift currents in Sec. III.2. In our calculations we use a geodesic grid with spheres and points per sphere and we thus assign more importance to the radial resolution than to the angular resolution. This can be justified by the fact that spherical symmetry is the dominant symmetry the of kp-Hamiltonian at the -point and that this resolution results in a converged exciton resonance. For numerical stability we introduce a weak screening and evaluate the Coulomb matrix element as
As screening constant we use nm which is much smaller than the inverse of the exciton Bohr radius nm nm and thus only weakly influences the spectra. In our simulations we use a dephasing time of fs which is chosen sufficiently long in order to spectrally separate the exciton peak properly from the band edge.
iii.1 Linear absorption and exciton wave function
To calculate linear absorption spectra, we numerically solve the SBE for excitation with a weak linearly polarized ultrashort pulse and obtain the time-dependent total optical polarization . Here, is the transition frequency between the valence band and the conduction band at wave vector and the term corresponds to the transition dipole matrix element. The linear absorption is then taken to be proportional to the imaginary part of the Fourier transformed polarization in the direction of the linearly excitation (), i.e., . Figure 1(a) shows the calculated absorption spectrum of bulk GaAs with and excitonic effects. Without Coulomb interaction the absorption starts at the band gap and follows a -dependence which originates from the density of states in three dimensions for a parabolic dispersion. When excitonic effects are included the absorption spectrum contains a well defined exciton peak at meV, thus the obtained exciton binding energy is in good agreement with literature values. Furthermore, Fig. 1(a) shows the Coulomb enhanced continuum absorption which is basically constant for energies slightly above the band gap. The increased absorption direction below the band gap originates from excited exciton states which are not resolved individually and merge with the continuum absorption.
Unlike a parabolic band structure, the kp band structure contains anisotropies, e.g., small splittings caused by the spin-orbit interaction that appear in the (111) and (100) directions. To investigate the influence of this anisotropy on the exciton wave function, we solve the SBE with a vector potential corresponding to a weak and slowly rising electric field, which is given by
where the excitation frequency chosen to be resonant with the exciton. We integrate the SBE up to ps and due to the narrow spectral width of the excitation only the exciton is significantly excited. Thus for the chosen exciting field, the microscopic polarization in reciprocal space is proportional to the exciton wave function, i.e.,
The microscopic polarization is determined as the -component of , which is the component in the direction of the linearly polarized excitation pulse. Since each point on a sphere of the geodesic grid contributes proportional the square of the radius when integrating over k-space, we multiply by in Fig. 1(b). We see, that the exciton wave function has a peak at approximately which corresponds to the inverse Bohr radius , with nm being the effective Bohr radius of the exciton in GaAs.Chow and Koch (1999) For an isotropic exciton the wave function does not depend on the direction in k-space and only a single line should be in Fig. 1(b). In our kp approach, however, the anisotropy of the band structure leads to a dependence of the exciton wave function on the direction in k-space, as is evidenced by the three different lines visible in Fig. 1(b). (A) corresponds to the in reciprocal space, where is the so called golden ratio which is used in the construction of the geodesic grid. Very similar results are also obtained in , , and , so we can write as a shorthand notation for these directions. (B) and (C) correspond to the 4 directions and , respectively. The difference between the wave function amplitudes in the different directions is significant and amounts to approximately 20% and thus clearly demonstrates the anisotropy of the exciton wave function.
iii.2 Excitonic and near-band-gap shift currents
The photocurrents induced by optical interband excitation are determined from solutions of the SBE byDuc et al. (2016)
The shift current is a response and as such represents the components of Eq. (13). Therefore is Fourier transformed and to the resulting a frequency filter is applied around . Considering a filter function which cuts off frequency components with meV, we obtain . By Fourier transforming back to the time domain we finally obtain . As is given by a double summation over the 14 bands, i.e., , one can distinguish different contributions arising from transitions between different valence and conduction bands. Thus one can write , as a sum
where () originates from intersubband transitions between different valence (conduction) bands and originates from interband transitions between a valence and a conduction band.
In Fig. 2(a) we show shift currents with and without Coulomb interaction which are generated by excitation with a pulse with a Gaussian envelope of 500 fs duration (FWHM of the pulse intensity) at the center of the pulse as function of the excitation frequency. The pulse is linearly polarized in (110)-direction which leads to a shift current in (001)-direction. Without Coulomb interaction the shift current vanishes for excitation frequencies below the band gap and follows the square root like behavior of the three-dimensional density of states for above band gap excitation. With Coulomb interaction the shift current is strongly enhanced. At a peak is visible in the total shift current which corresponds to excitation of the exciton resonance. For excitations energies above the band gap the total current initially grows with excitation energy and then approaches a basically constant value.
Analyzing the different sub currents, introduced above, offers an explanation for this behavior. and , the shift currents created by the conduction band and the interband polarization, have both a peak at the exciton and remain mostly constant for above band gap. , the valence band shift current, displays a very different behavior. At the exciton energy displays a negative current and remains negative up until , where the current changes its direction and then changes to a positive current. The reason for this behavior is two-fold. For excitation at the exciton and at the band gap the holes are correlated with the electrons and in average flow in the same directions. Due to the opposite sign of their charges, the resulting signs of corresponding currents are opposite. For excitation energies higher in the band the holes are less correlated with the electrons and therefore, as without Coulomb interaction, both currents have a positive sign, which corresponds to electrons and holes moving in opposite directions. In addition to the shift current contribution, also includes the intersubband coherence between the heavy hole and the light hole band. Since both bands are degenerate at the -point, their energy difference starts with and thus cannot be removed from the valence band shift current via Fourier filtering. This coherence is also present without Coulomb interaction and leads to a very small negative current when exciting at the band gap in .
The calculated time evolution of shift currents after excitation with a pulse with a Gaussian envelope of 500 fs duration that has its maximum at is shown in Fig. 3. For excitation at the band gap, see Fig. 2(b), and have a Gaussian shape, centered at . This is to be expected since the shift current involves is an off-resonant excitation process and in the limit of a purely off-resonant excitation the signal should strictly follow the envelope of the optical excitation Nastos and Sipe (2006). This behavior is also in agreement with previous results obtained with the SBEPodzimski et al. (2015, 2016); Duc et al. (2016) where excitonic effects have been neglected. , however, shows a weak temporal delay to positive times which originates from optically induced coherences between the heavy- and light-hole bands. When exciting above the band gap, see Fig. 2(c), and remain basically unchanged. , however, clearly displays a slow oscillatory behavior which is caused by the now dominant influence of the coherences between the heavy- and light-hole bands whose average energetic separation increases with increasing distance from the -point and thus with excitation energy.
When exciting at the exciton resonance, see Fig. 2(a), the total current and the sub currents have basically a Gaussian shape. In this case, however, a small delay of about to positive times is visible. This temporal shift does not significantly change when we use different dephasing times in our simulations. It can be interpreted by the following effect: When the excitation is tuned at or above the band gap, already the linear optical polarization corresponds to a superposition of a continuum of transitions with different frequencies which interfere destructively and thus lead to a decay of the total linear polarization on the time scale of the exciting pulse. This situation is different when we tune the excitation to the exciton resonance which is a single discrete optical transition whose linear polarization increases with the integral over the envelope of the exciting pulse, at least as long as dephasing is neglected. The above explained difference in the dynamics of the linear optical polarization depending on the spectral position of the excitation could be responsible for the delayed maximum of the excitonic shift current seen in Fig. 2(a).
Nastos and Sipe introduced the concept of the shift distanceNastos and Sipe (2006)
where is the photoexcited density, which describes the average distance that the electrons shift in space when optically excited from the valence to the conduction band. Eq. (15) was derived in a frequency-domain approach considering continuous wave excitation. Our time-domain simulations model transient situation with pulsed excitation and dephasing and relaxation of the coherences and occupations which makes Eq. (15) rather unsuitable. Based on simple electrodynamic consideration we define a time-dependent shift distance by
Due to the pulsed excitation and the dephasing and relaxation processes considered in our approach Eq. (16) can only approximate the shift distance and in particular for long times unphysically increases since decreases. In our numerical evaluations we find a plateau of the time-dependent shift distance during early times of the excitationsPodzimski (2016) and therefore determine in the (001) direction using Eq. (16) by averaging over the rising part of the incident optical pulse using the time interval fs to fs. For exciting at the exciton we obtain a shift distance of . When shifting the photon energy to above the band gap, the shift distance increases, e.g., and . This agrees with the interpretation that the attraction between holes and electrons causes a reduction of the shift distance for excitations near the exciton. Without Coulomb interaction we obtain a shift distance which is quite close to the average shift distance of calculated in Nastos and Sipe (2006) for GaAs considering continuous wave excitation high above the band gap. Since the 14 band kp band structure is not well suited to describe k vectors far away from the -point, we, however, cannot analyze excitations very high above the band gap and thus not directly compare with the results of Nastos and Sipe (2006).
Using the semiconductor Bloch equations in a basis obtained from a 14 band kp model of the band structure we calculate linear absorption spectra including excitonic effects and shift current transients of bulk GaAs for near band gap excitation. The electron-hole Coulomb attraction results in an excitonic resonance with a binding energy which is in good agreement with other approaches and experiments. We show that the anisotropy of the band structure leads to a significant anisotropy of the exciton wave function. The calculated frequency-dependence of the shift current has a peak at the exciton resonance. However, the enhancement at the exciton is weaker than the excitonic enhancement in the linear absorption and correspondingly the shift distance when exciting at the exciton is weaker compared to excitations above the band gap. The reduced excitonic shift distance as well as the sign change of the valence band current for below band excitation originates from the electron-hole attraction. In addition, for above band gap excitations coherences between the heavy- and light-hole bands lead to oscillatory signature in the current transients. Our findings demonstrate that our approach is able to provide important information on the ultrafast electron dynamics on the atomic scale. It is also applicable to study recently predicted Berry-curvature-induced currents arising from bound excitons Morimoto and Nagaosa (2016) which we plan to investigate for pulsed excitation in future work.
We are grateful to John Sipe for stimulating and fruitful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) through the projects ME 1916/2 and SFB/TRR 142 (Project A02). H. T. Duc acknowledges the financial support of the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under Grant No. 103.01-2017.42. We thank the PC (Paderborn Center for Parallel Computing) for providing computing time.
- J. E. Sipe and A. I. Shkrebtii, Phys. Rev. B 61, 5337 (2000).
- F. Nastos and J. E. Sipe, Phys. Rev. B 74, 035201 (2006).
- X.-C. Zhang, Y. Jin, K. Yang, and L. J. Schowalter, Phys. Rev. Lett. 69, 2303 (1992).
- G. Dalba, Y. Soldo, F. Rocca, V. M. Fridkin, and P. Sainctavit, Phys. Rev. Lett. 74, 988 (1995).
- Z. Chen, M. Segev, D. W. Wilson, R. E. Muller, and P. D. Maker, Phys. Rev. Lett. 78, 2948 (1997).
- D. Côte, N. Laman, and H. M. van Driel, J. Appl. Phys. 80, 905 (2002).
- C. Somma, K. Reimann, C. Flytzanis, T. Elsaesser, and M. Woerner, Phys. Rev. Lett. 112, 146602 (2014).
- N. Laman, M. Bieler, and H. M. van Driel, J. Appl. Phys. 98, 103507 (2005).
- M. Bieler, K. Pierz, and U. Siegner, J. Appl. Phys. 100, 83710 (2006).
- M. Bieler, K. Pierz, and U. Siegner, Phys. Rev. B 76, 161304 (2007).
- S. Priyadarshi, M. Leidinger, K. Pierz, A. M. Racu, U. Siegner, M. Bieler, and P. Dawson, J. Appl. Phys. 95, 151110 (2009).
- H. T. Duc, R. Podzimski, S. Priyadarshi, M. Bieler, and T. Meier, Phys. Rev. B 94, 085305 (2016).
- S. M. Young and A. M. Rappe, Phys. Rev. Lett. 109, 116601 (2012).
- S. M. Young, F. Zheng, and A. M. Rappe, Phys. Rev. Lett. 110, 057201 (2013).
- P. Král, J. Phys.: Condens. Matter 12, 4851 (2000).
- J. A. Brehm, S. M. Young, and A. M. Rappe, J. Chem. Phys. 141, 204704 (2014).
- S. Priyadarshi, A. M. Racu, K. Pierz, U. Siegner, M. Bieler, H. T. Duc, J. Förstner, and T. Meier, Phys. Rev. Lett. 104, 217401 (2010).
- H. T. Duc, J. Förstner, and T. Meier, Phys. Rev. B 82, 115316 (2010).
- H. T. Duc, J. Förstner, T. Meier, S. Priyadarshi, A. M. Racu, K. Pierz, U. Siegner, and M. Bieler, Phys. Status Solidi (c) 8, 1137 (2011).
- R. Podzimski, H. T. Duc, and T. Meier, Proc. SPIE 9361, 93611V (2015).
- R. Podzimski, H. T. Duc, and T. Meier, Proc. SPIE 9746, 97460W (2016).
- R. Zimmermann and C. Trallero-Giner, Phys. Rev. B 56, 9488 (1997).
- S. Rudin and T. L. Reinecke, Phys. Rev. B 65, 121311 (2002).
- R. D. R. Bhat and J. E. Sipe, Phys. Rev. B 72, 075205 (2005).
- V. Turkowski and C. A. Ullrich, Phys. Rev. B 77, 075204 (2008).
- V. Turkowski, A. Leonardo, and C. A. Ullrich, Phys. Rev. B 79, 233201 (2009).
- H. Haug and S. W. Koch, Quantum Theory of the Optical and Electronic Properties of Semiconductors, 4th ed., (World Scientific, 2004).
- H. Mayer and U. Rössler, Phys. Rev. B 44, 9048 (1991).
- R. Winkler, Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems (Springer, 2003).
- H.-R. Terbin, U. Rössler, and R. Ranvaud, Phys. Rev. B 20, 686 (1979).
- D. E. Aspnes and A. A. Studna, Phys. Rev. B 7, 4606 (1973).
- M. Cardona, N. E. Christensen, and G. Fasol, Phys. Rev. B 38, 1806 (1988).
- E. Sternemann, T. Jostmeier, C. Ruppert, H. T. Duc, T. Meier, and M. Betz, Phys. Rev. B 88, 165204 (2013).
- Y.-C. Chang and R. B. James, Phys. Rev. B 39, 12672 (1989).
- P. Enders, A. Bärwolff, and M. Woerner, Phys. Rev. B 51, 16695 (1995).
- U. Rössler and D. Strauch, eds., Landolt-Börnstein - Group III Condensed Matter, Lattice Properties, Vol. 41A1a (Springer-Verlag Berlin Heidelberg, 2001).
- R. Heikes and D. A. Randall, Monthly Weather Review 123, 1862 (1995).
- A. Mahdavi-Amiri, E. Harrison, and S. Faramarz, International Journal of Digital Earth 8, 750 (2015).
- N. O. Lipari and A. Baldereschi, Phys. Rev. Lett. 25 (1970).
- J. M. Luttinger, Physical Review 102, 1030 (1956).
- G. H. Wannier, Physical Review 52, 191 (1937).
- R. Podzimski, PhD thesis: Shift Currents in Bulk GaAs and GaAs Quantum Wells Analyzed by a Combined Approach of kp Perturbation Theory and the Semiconductor Bloch Equations (Universität Paderborn, 2016).
- W. W. Chow and S. W. Koch, Semiconductor-Laser Fundamentals (Springer-Verlag Berlin Heidelberg, 1999).
- T. Morimoto and N. Nagaosa, Phys. Rev. B 94, 035117 (2016).