We investigate theoretically valley-resolved lateral shift of electrons traversing an n–p–n junction bulit on a typical tilted Dirac system (8-Pmmn borophene). A gauge-invariant formula on Goos–Hänchen (GH) shift of transmitted beams is derived, which holds for any anisotropic isoenergy surface. The tilt term brings valley dependence of relative position between the isoenergy surface in n region and that in the p region. Consequently, valley double refraction can occur at the n–p interface. The exiting positions of two valley-polarized beams depend on the incident angle and energy of incident beam and barrier parameters. Their spatial distance D can be enhanced to be ten to a hundred times larger than the barrier width. Due to tilting-induced high anisotropy of the isoenergy surface, D depends strongly on the barrier orientation. It is always zero when the junction is along the tilt direction of Dirac cones. Thus GH effect of transmitted beams in tilted Dirac systems can be utilized to design anisotropic and valley-resolved beam-splitter.
Xixuan Zhou, Jianlong Zheng(郑建龙), Feng Zhai. Anisotropic and valley-resolved beam-splitter based on a tilted Dirac system[J]. Communications in Theoretical Physics, 2022, 74(7): 075701. DOI: 10.1088/1572-9494/ac6fc2
1. Introduction
The common wave nature of light and electrons makes it possible to manipulate electrons by means of components inspired by geometrical optics [1–19]. For ballistic electrons traversing a boundary separating two uniform regions with different potentials, reflection and refraction happen in analogy to a light beam incident on an interface between materials with different optical indices. Both positive and negative refraction have been observed [14] experimentally in graphene junctions. It is well known that a totally reflected beam of light undergoes a lateral displacement relative to the incident point. This phenomenon is referred to as Goos–Hänchen (GH) shift [20] and has been generalized to partial reflections and also transmitted beams [3]. The GH effect in graphene p–n–p junctions has been predicted to result in an 8e2/h conductance plateau [8]. The GH shifts of transmitted electron beams have been utilized to design a spin beam splitter [6] based on a two-dimensional (2D) electron gas or valley beam splitter [11] based on strained graphene. There are extensive studies on the GH shift of electrons in materials with an isotropic Fermi surface, especially in Dirac materials with linear dispersion. Very recently, Ghadiri and Saffarzadeh [19] have calculated the GH shift of transmitted electrons on the surface of a topological insulator with a potential barrier and a warped Fermi surface. Their derivation on the GH shift is made only for a special barrier orientation which is straightforward but cumbersome. A general formula and numerical approach are necessary for the GH shift of electrons on an anisotropic Fermi surface.
Materials with a gapless dispersion and anisotropic Fermi surface are similar to anisotropic optical medium. Among these materials tilted Dirac systems have attracted great research interest. Tilted Dirac cones can appear in 8-Pmmn borophene [21, 22], organic conductor α -(BEDT-TTF)2I3 [23, 24], quinoid-type or partially hydrogenated graphene [25, 26], and 1T${}^{{\prime} }$ monolayer transition metal dichalcogenides [27, 28]. Among them, the organic 2D conductor α-(BEDT-TTF)2I3 has been experimentally verified to host a pair of tilted and anisotropic Dirac cones. Some unusual electronic and transport properties caused by the tilt of the Dirac cone have been demonstrated theoretically, such as strongly anisotropic plasmon [29] and optical conductivities [30], unique intervalley damping effect in magnetoplasmons [31], valley-dependent Weiss oscillation [32], nearly perfect valley polarization induced by a single ferromagnetic gate [33], and light-driven metal–insulator transition [34, 35]. Electron optics based on tilted Dirac systems has been discussed in several works [33, 36, 37] which concern valley-dependent electron retroreflection, oblique Klein tunneling, generalized Klein tunneling, and Veselago focusing. The GH shift in tilted Dirac systems has not been studied so far. In this work, we take 8-Pmmn borophene as an example to investigate the GH shift of electrons on an anisotropic Fermi surface.
2. Model and formalism
The 2D electron system in 8-Pmmn borophene under consideration is depicted in figure 1, which is modulated by an electric potential U(r) (yellow-shadowed region) due to a metallic gate on top. The electric potential in the borophene plane (XOY) varies only along the x axis and vanishes outside the region 0 < x < W. The X axis is perpendicular to the Σv mirror in pristine 8-Pmmn borophene. The angle between these two directions is denoted as α (0 ≤ α ≤ π/2). The motion of an electron in a given valley τ can be described by the low-energy effective Hamiltonian [22, 32, 37]
Here τ = ±1 stands for the kD and − kD valleys, ΣX, ΣY together with ΣZ are three Pauli matrices, Σ0 is a unit matrix, and $\hat{{\boldsymbol{p}}}=-{\rm{i}}$ℏ∇ is the momentum operator. Hereafter the velocity, length and energy are, respectively, in units of vF = 106 m s−1, L0 = 50 nm and E0 = ℏvF/L0 = 13.2 meV. Without specification, the three velocity parameters vX, vY and vt in units of vF are taken, respectively, as [22] 0.86, 0.69, and 0.32. For some 2D organic conductors [23, 24] or graphene-type materials [25, 26] with two degenerate and tilted Dirac cones, the low-lying excitations near Dirac points can be also described by the Hamiltonian equation (1).
Figure 1. Schematic of the considered device. The electric potential varies only along the x axis and vanishes outside the shadowed region. The x axis has angle α to a symmetry axis X of pristine 8-Pmmn borophene.
2.1. Scattering states
The momentum operators ${\hat{p}}_{X}$ and ${\hat{p}}_{Y}$ in equation (1) can be written as ${\hat{p}}_{X}={\hat{p}}_{x}\cos \alpha -{\hat{p}}_{y}\sin \alpha $ and ${\hat{p}}_{Y}={\hat{p}}_{x}\sin \alpha \,+{\hat{p}}_{y}\cos \alpha $. The velocity operator ${\hat{v}}^{x}=\tfrac{\partial {\hat{H}}^{\tau }}{\partial {\hat{p}}_{x}}$ and ${\hat{v}}^{y}=\tfrac{\partial {\hat{H}}^{\tau }}{\partial {\hat{p}}_{y}}$ read
The Hamiltonian (1) is translationally invariant along the y direction. Its eigenstate with energy E admits the form ${{\rm{\Psi }}}^{\tau ,q}({\boldsymbol{r}})=\exp ({\rm{i}}{qy}){\psi }^{\tau ,q}(x)$. Here ℏq is the conserved transverse momentum. In a region with label j and a constant U(r) ≡ Uj, the spinor wave function ψτ,q(x) can be written as
For a real root k of equation (7), one can calculate the group velocity ${v}_{j,k}^{x/y}$ and propagate angle φj,k for the plane wave ${\psi }_{j,k}^{\tau ,q}{{\rm{e}}}^{{\rm{i}}{kx}}$
In the case that equation (7) has real roots, kj1 is selected as the root k with ${v}_{j,k}^{x}\gt 0$. The group velocity can be also calculated from the dispersion E = E(k, q) determined by equation (7), ${v}_{j,k}^{x}={\left(\tfrac{\partial E}{\partial k}\right)}_{q}$ and ${v}_{j,k}^{y}={\left(\tfrac{\partial E}{\partial q}\right)}_{k}$. By means of the Maxwell relation on variables k, q, and E in equation (7), one gets another expression on the propagate angle
For the incidence of electrons from the left lead (x < 0), the wave function ${\psi }_{j=L}^{\tau ,q}(x)$ in the left lead has the form equation (3) with UL = 0 and ${C}_{L1}=1/\sqrt{{v}_{L,{k}_{L1}}^{x}}$. The wave function ${\psi }_{j=R}^{\tau ,q}(x)$ in the right lead (x > W) is also written as equation (3) but with UR = CR2 = 0. The reflection and transmission coefficients are
For brevity, we have suppressed the valley τ, energy E, and wave vector q dependence of kj1, kj2, Cj1, Cj2, ${v}_{j,k}^{x/y}$, φj,k, r, and t. As a function of these variables, the transmission coefficient can be written as t = tτ(E, q).
2.2. GH shift of outgoing beams
Consider a wavepacket consisting of incident electrons in valley τ with the same energy E but different transverse wave vector q
Here ${\psi }_{I}(q)={\psi }_{L,{k}_{L1}}^{\tau ,q}$. The normalizable real function f(q) determines the profile of the wavepacket. The outgoing beam can be written as
Here ${\psi }_{O}(x;q)={C}_{R1}{\psi }_{R,{k}_{R1}}^{\tau ,q}{{\rm{e}}}^{{\rm{i}}{k}_{R1}x}$. We calculate its average transverse position for a given x,
The profile of the incident beam is usually taken as a Gaussian with center q and width Δq, i.e. $f(p)\propto \exp \left[-{\left(p-q\right)}^{2}/(2{{\rm{\Delta }}}_{q}^{2})\right]$. For a narrow width Δq, we can approximate ∫F(p)f2(p)dp in equations (17) and (18) as F(q)∫f2(p)dp and then yield
The GH shift is defined as the difference in the transverse position between the outgoing beam (at the interface with abscissa x > W) and the incident beam (at the interface x = 0)
This is the general expression of the GH shift for a narrow Gaussian beam. It applies to any multicomponent electron scattering state. When ψO(x; q) = t(x; q)ψI(q) with t(x; q) the scalar transmission coefficient, equation (20) becomes ${s}_{\mathrm{GH}}(q;x)=-\tfrac{\partial \arg t(x;q)}{\partial q}$. In this special case, the GH shift is determined only by q-dependence of the transmission phase $\arg t(x;q)$, as in [3] and [11]. If several transmitted beams coexist [19], one can calculate the GH shift for each of them from equation (20) where ψO(x; q) should be the spinor of the corresponding right-propagating mode.
Consequently, the center of the outgoing beam moves along the ray defined by the outgoing angle ${\varphi }_{R,{k}_{R1}}$. Note that in [19] the propagate angle for the propagating mode ${\psi }_{j,k}^{\tau ,q}{{\rm{e}}}^{{\rm{i}}{kx}}$ is defined as $\arctan \tfrac{q}{k}$ (the orientation of phase velocity). From the view of GH shift, the propagate angle should be defined as equation (8) rather than $\arctan \tfrac{q}{k}$.
It can be also checked that sGH is invariant under any smooth scaling transformation ${\psi }_{I}(q)\longrightarrow {\psi }_{I}(q)\exp [{\rm{i}}\zeta (q)]c(q)$ with c(q) > 0. We can choose ζ(q) to satisfy
Under such a gauge, i.e. replacing ψI(q) with ${\psi }_{I}^{\mathrm{new}}(q)\,={\psi }_{I}(q)\exp [{\rm{i}}\zeta (q)]$ so that ${\left[{\psi }_{I}^{\mathrm{new}}(q)\right]}^{+}d{\psi }_{I}^{\mathrm{new}}(q)=0$, the second term in equation (20) vanishes. For the numerical calculation of sGH, this gauge choice can circumvent the continuity requirement of ψI(q). For a smoothing profile of U(x), the scattering state for the incidence ψI(q) can be calculated numerically by means of the scattering matrix method [38]. One then applies the central difference formula to calculate the derivative of ψI(q) and ψO(x; q) with respect to q. The GH shift in equation (20) is obtained directly from the scattering state and its derivative.
3. Results and discussions
For simplicity, we consider the electric potential with a rectangular shape [see figure 2(a)], i.e. U(x) = UΘ(x)Θ(W − x) with Θ(x) the Heaviside step function and U > 0 the barrier height. For an energy E > 0 and barrier orientation α = 15°, the isoenergy line in the lead region for electrons in valley τ is plotted in figure 2(b), which is an ellipse centered at the point with pX = 0 and ${p}_{Y}=-\tau {{Ev}}_{t}/({v}_{Y}^{2}-{v}_{t}^{2})$. Here pX and pY are (dimensionless) momentum components along the X and Y axis depicted in figure 1. The ellipse equation is given by equation (7) and denoted as fτ(k, q; E) = 0. Note that f−1(k, q; E) = f1( − k, − q; E). The major and minor axis of the ellipse is along the pY and pX axis and have length $2E/\sqrt{{v}_{Y}^{2}-{v}_{t}^{2}}$ and $2E/[{v}_{X}\sqrt{1-{\left({v}_{t}/{v}_{Y}\right)}^{2}}]$. One can see from equation (9) that the incident direction is perpendicular to the tangent of the isoenergy line in the (k, q) plane. As an example, in figure 2(b) we show a positive incident angle φc2 for electrons in valley τ = −1. Note that electrons in two valleys with the same incident angle φ have different transverse wave vectors q±. The isoenergy surfaces in the barrier region are given by fτ(k, q; E − U) = 0, which are similar to those in the lead region with a scaling factor ∣E − U∣/E. In the case E < U, the ellipse f−1(k, q; E − U) = 0 is below the ellipse f1(k, q; E − U) = 0, as shown in figure 2(b). In the (k, q) plane a horizontal (dashed) line q = q0 represents the conserved transverse wave vector. When it intersects with both the ellipse fτ(k, q; E) = 0 and fτ(k, q; E − U) = 0, refraction occurs at interface x = 0 for an incident electron in valley τ. The sector of allowed q0 depends on the topmost and bottommost points of the two ellipses. For 0 < 2E < U the critical values qc1 and qc2 together with the critical angle φc2 are depicted in figure 2(b) for valley τ = −1. In this case, a valley-unpolarized beam with incident angle φ ∈ ( − φc2, φc2) will split in the barrier region. The two refracted beams with valley index τ = ±1 propagate at different angles φ± , as depicted in figure 2(c). The refracted angles can be determined from the ellipse fτ(k, q, E − U) = 0. The existing position of refracted beam τ at the interface x = W is given by the GH shift sτ = sGH(qτ; x = W), which usually differs from ${s}_{\tau }^{\mathrm{cl}}=W\tan {\phi }_{\tau }$ predicted by the geometric ray.
Figure 2. (a) Rectangular potential barrier with height U and width W. (b) Isoenergy lines (for an energy E > 0) of electrons in valley τ = 1 (blue lines) and τ = −1 (red lines) in the lead and in the barrier. (c) Representation of valley-dependent refraction angle φ± and GH shift s±.
In the barrier region with index j = B and Uj = U, one can find two plane (or evanescent) waves ${\left(1,{m}_{B,{k}_{B1}}^{\tau ,q}\right)}^{{\rm{T}}}{{\rm{e}}}^{{\rm{i}}{k}_{B1}x}$ and ${\left(1,{m}_{B,{k}_{B2}}^{\tau ,q}\right)}^{{\rm{T}}}{{\rm{e}}}^{{\rm{i}}{k}_{B2}x}$ from equations (3)–(7). For an incident mode ${\left(1,{m}_{L,{k}_{L1}}^{\tau ,q}\right)}^{{\rm{T}}}$, one can work out the transmission coefficient
In figure 3 we plot the valley-resolved transmission probability Tτ = ∣tτ∣2 and GH shift sτ at the interface x = W as functions of the barrier width W under a fixed incident angle φ and energy E. The barrier height and orientation are fixed at U = 4 and α = 0. Four typical cases are considered and discussed as follows.
Figure 3. Transmission T±1 (in black) and GH shift s±1 (in red) at the existing position as functions of barrier width W for electrons in valley τ = 1 (solid lines) and τ = −1 (dashed lines). The blue lines represent ${s}_{\tau }^{{cl}}$ if corresponding refracted beam exists. The barrier height and orientation are U = 4 and α = 0. The incident energy E and incident angle are (a) E = 1, φ = 0; (b) E = 1, φ = 60°; (c) E = 1.2, φ = 60°; and (d) E = −1, φ = 60°.
The case of normal incidence (φ = 0) is presented for E = 1 (an n–p–n junction) in figure 3(a), where positive/negative refraction (see blue lines) happens in the barrier region for electrons in valley τ = +1/−1. Obviously, this observation violates Snell's law. The reason is that the anisotropy of the isoenergy line results in an orientation deviation between the momentum and group velocity. For α = 0, the isoenergy lines for the two valleys are symmetric about the horizontal k axis. Therefore, the transverse wave vectors q± and refraction angles φ± satisfy q− = −q+ and φ− = −φ+. One can also observe that the transmission probability for two valleys coincides. Actually, due to the fact that the inversion operator transforms the Hamiltonian ${\hat{H}}^{\tau =1}$ into ${\hat{H}}^{\tau =-1}$, the transmission coefficient satisfies a general relation
Since q− = −q+, this relation leads to T−1 = T+1 and s−1 = − s+1, as seen in figure 3(a). The oscillation of transmission probability with the width W arises from the Fabry–Perot interference. The GH shift sτ oscillates around the classical value ${s}_{\tau }^{\mathrm{cl}}$ with a small amplitude.
In figure 3(b), the incident angle and energy are chosen as φ = 60° and E = 1 so that φ is close to the critical angle φc2 depicted in figure 2(b) (but for α = 0). The two refraction angles are φ+ = 13° and φ+ = − 60°. In the barrier region, the longitudinal wave vector for valley τ = +1 is much larger than that for valley τ = −1. Consequently, the Fabry–Perot oscillation for valley τ = −1 has a longer period and larger peak-to-trough ratio than that for valley τ = +1. The wave effect on GH shift is weak for valley τ = +1 but drastic for valley τ = −1. At the first (second) resonant peaks of T−1, the amplitude of GH shift s−1 exceeds 48 (96), which is 20 times larger than the corresponding barrier width. As demonstrated in [11], the GH shift at a resonant peak increases quickly with the peak-to-trough ratio. When the energy E turns from 1 to 1.1, the critical angle φc2 is much close to φ = 60° so that ∣s−1∣ can exceed 800 > 100W. Near these resonant peaks, the valley spatial separation D = ∣s+1 − s−1∣ is much larger than the counterpart ${D}^{\mathrm{cl}}=| {s}_{+1}^{\mathrm{cl}}-{s}_{-1}^{\mathrm{cl}}| $.
Once the critical angle φc2 is smaller than the incident angle φ, the refracted beam for valley τ = −1 vanishes. Such a case is presented in figure 3(c) with φ = 60° and E = 1.2. One can see that the transmission T−1 and GH shift s−1 decays quickly with the barrier width. The corresponding transmitted beam is difficult to detect. The refracted beam for valley τ = −1 remains. In comparison with the case of normal incidence, the oscillation of ${s}_{+1}-{s}_{+1}^{\mathrm{cl}}$ with the width W is more remarkable due to the larger peak-to-trough ratio in the transmission T+1. Accordingly, the polarity of s+1 alternates with W. One can see offsets between peaks of T+1 and s+1, which are absent in strained graphene [11].
In the case E < 0 < U (a p–${p}^{{\prime} }$–p junction), the isoenergy line fτ(k, q; E) = 0 is enclosed in the ellipse fτ(k, q; E − U) = 0 so that the valley double refraction happens for all incident angles. For E = −1 , the refracted beams for φ = 60° are shown as blue lines in figure 3(d), whose orientation angles 25.4° and −7.66° are close to their limits (26.3° and −7.03°) under φ $\longrightarrow $ 90°. Now both ${s}_{+1}-{s}_{+1}^{\mathrm{cl}}$ and ${s}_{-1}-{s}_{-1}^{\mathrm{cl}}$ oscillate rapidly with W. The two oscillations have close periods and amplitudes. The valley separation D can vary from 0 to a value higher than 2Dcl.
Hereafter we consider only the n–p–n junction with barrier width W = 3 and barrier height U = 4. We plot in figures 4 and 5 the variation of transmission T±1 and GH shift s±1 with the incident angle φ. For the barrier orientation α = 0, the transmission (GH shift) for electrons in valley τ = +1 and that for electrons in valley τ = −1 are symmetric (antisymmetric) about the incident angle φ = 0, as shown in figure 4. This feature results from equation (29). It is absent for other barrier orientations (see figure 5 for α = 45°). For a given valley, in all cases the angular distribution of transmission (GH shift) is not symmetric (antisymmetric) about φ = 0. There are several angles where the transmission Tτ is perfect. For a resonant peak of Tτ, when its left/right (nearest) transmission minimum is lower than the right/left one, the corresponding maximum of ∣sτ∣ locates on its left/right side. The highest ∣sτ∣ depends on the incident energy E and barrier orientation α. For E ≤ U/2 (such as E = 1), the transmission T−τ is remarkable near the angle with the largest ∣sτ∣, enabling beam splitting with a large valley spatial separation. In contrast, For E > U/2 (such as E = 3), the angular sector with large transmission T−τ does not overlap with that with large Tτ (or ∣sτ∣), allowing valley filtering of incident beams. For E = U/2 and α = 0, both T−1 and T+1 reach 1 at φ = 0. The valley spatial separation in this case is almost twice the largest ∣s±1∣. For E = 1 (E = 2 and E = 3), the largest ∣sτ∣ for all incident angles under α = 45° is near half (twice) of that under α = 0. One can understand these features from the isoenergy surfaces as the way in figure 2(b).
Figure 4. Transmission T±1 (in black) and GH shift s±1 (in red) as functions of incident angle φ for electrons in valley τ = 1 (solid lines) and τ = −1 (dashed lines). The barrier height and width are U = 4 and W = 3. We take α = 0 and set E = 1 (a), 2(b), and 3(c).
Figure 5. Same as figure 4, but for the barrier orientation α = 45°.
Since isoenergy lines depcited in figure 2(b) are high anisotropic, one can expect that the valley spatial separation D depends strongly on the barrier orientation α. In figure 6 we plot the transmission and GH shift as functions of α for a fixed energy E = 1. For the transport along the tilted direction of the Dirac cone (α = 90°), the product of Σz and conjugation K transforms the Hamiltonian H+1 into H−1, which leads to t−1(E, q) = t+1(E, q). As a result, the transmission and GH shift are valley-independent at α = 90°. In the case of normal incidence (figure 6(a)), as α increases from 0° to 90°, the transmission T±1 > 0.6 is noticeable. The valley separation D moves up firstly to its global maximum ≈5 (at α ≈ 29.8°) and then decreases oscillatingly to zero. At the fixed incident angle φ = 60° (figure 6(b)), there are several values of α where T−1 approaches 1 and T+1 is remarkable. The valley separation D has a maximum ≈39 at α ≈ 32.3° near the first peak of T−1.
Figure 6. Transmission T±1 (in black) and GH shift s±1 (in red) as functions of barrier orientation α for electrons in valley τ = 1 (solid lines) and τ = −1 (dashed lines). The barrier height and width are U = 4 and W = 3. We take E = 1 and set φ at 0 in (a) and 60° in (b).
The tilt term in equation (1) is essential for the valley spatial separation of transmitted beams. In figure 7 we plot the transmission and GH shift as functions of tilt velocity vt for the barrier orientation α = 0 and a fixed energy E = 1. In the case vt = 0, one has T+1 = T−1 and s+1 = s−1 for all incidences. At normal incidence (figure 7(a)), when vt increases from 0 to 0.35 the transmission T±1 varies slightly while the valley separation D increases gradually from 0 to 4.9. The valley separation D exceeds 77 (16) at large vt ≈ 0.51 (0.46). As for the incident angle φ = 60° (figure 7(b)), the transmission T−1 has two sharp peaks located at vt ≈ 0.27 and 0.34 where the valley separation D is larger than 27 and 150.
Figure 7. Transmission T±1 (in black) and GH shift s±1 (in red) as functions of tilt velocity vt for electrons in valley τ = 1 (solid lines) and τ = −1 (dashed lines). The barrier height, width, and orientation are U = 4, W = 3, and α = 0. We take E = 1 and set φ at 0 in (a) and 60° in (b).
4. Conclusions
In summary, we have derived a general expression for the GH shift of transmitted beams. The formula is independent of the phase choice of incident states and applies for any anisotropic isoenergy surface. We take 8-Pmmn borophene as an example to demonstrate valley beam-splitter based on tilted Dirac systems. The tilt term leads to valley-dependent isoenergy surfaces. Valley double refraction can happen at the n–p interface because the relative position between isoenergy surfaces in the lead and in the barrier region depends on the valley index. A valley-unpolarized beam traversing a potential barrier may split into two valley-polarized beams. When the refracted angle for one valley is close to its critical value, it is shown that the valley spatial separation D of transmitted beams can be ten to a hundred times larger than the barrier width. The high anisotropy of the isoenergy surface leads to a drastic variation of D with the barrier orientation α. D is always zero when the barrier is along the tilted direction (α = 90°) and usually reaches a maximum at α ≠ 0 for a given incident angle and energy. It is also shown that D can be enhanced by the tilt velocity. The remarkable GH shift of valley-polarized beams can be measured by the transverse magnetic focusing technique as in [14].
This work was supported by the National Natural Science Foundation of China (Grant No. 11774314).
GhadiriHSaffarzadehA2022 Electron beam splitting at topological insulator surface states and a proposal for electronic Goos–Hänchen shift measurement Phys. Rev. B105 085415
OsadaTKiswandhiA2020 Possible current-induced phenombena and domain control in an organic Dirac fermion system with weak charge ordering J. Phys. Soc. Jpn.89 103701