Reconstruction of the polarization distribution of the Rice-Mele model

TL;DR

This paper reconstructs the polarization distribution in the Rice-Mele model using gauge invariant cumulants, revealing how the distribution shifts and distorts under various parameters and topological cycles.

Contribution

It introduces a method to reconstruct the polarization distribution from moments in the Rice-Mele model, linking topological properties with distribution shape changes.

Findings

Distribution shifts towards the midpoint as on-site interaction decreases.

Shape of the distribution becomes asymmetric away from the fully dimerized limit.

Distribution moves across unit cells during topological cycles with significant shape distortion.

Abstract

We calculate the gauge invariant cumulants (and moments) associated with the Zak phase in the Rice-Mele model. We reconstruct the underlying probability distribution by maximizing the information entropy and applying the moments as constraints. When the Wannier functions are localized within one unit cell, the probability distribution so obtained corresponds to that of the Wannier function. We show that in the fully dimerized limit the magnitude of the moments are all equal. In this limit, if the on-site interaction is decreased towards zero, the distribution shifts towards the midpoint of the unit cell, but the overall shape of the distribution remains the same. Away from this limit, if alternate hoppings are finite, and the on-site interaction is decreased, the distribution also shifts towards the midpoint of the unit cell, but it does this by changing shape, by becoming asymmetric around the maximum, as well as by shifting. We also follow the probability distribution of the polarization in cycles around the topologically non-trivial point of the model. The distribution moves across to the next unit cell, its shape distorting considerably in the process. If the radius of the cycle is large, the shift of the distribution is accompanied by large variations in the maximum.

Figures

Tables

Table 1. Table 1: − log 10 ⁡ χ 2 subscript 10 superscript 𝜒 2 -\log_{10}\chi^{2} rounded to the first digit shown for the reconstructed probabilities in Figs. 4 and 5 .

Fig. 4 panel 1	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$
$J^{\prime }/J=0$	$-2$	$8$	$-1$	$7$
$J^{\prime }/J=0$	$-0.5$	$8$	$-0.4$	$6$
$J^{\prime }/J=0$	$-0.3$	$7$	$-0.2$	$6$
$J^{\prime }/J=0$	$-0.1$	$7$	$0$	$6$
Fig. 4 panel 2	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$
$J^{\prime }/J=0.3$	$-2$	$6$	$-1$	$5$
$J^{\prime }/J=0.3$	$-0.5$	$6$	$-0.4$	$4$
$J^{\prime }/J=0.3$	$-0.3$	$5$	$-0.2$	$4$
$J^{\prime }/J=0.3$	$-0.1$	$5$	$0$	$4$
Fig. 4 panel 3	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$
$J^{\prime }/J=0.5$	$-2$	$5$	$-1$	$4$
$J^{\prime }/J=0.5$	$-0.5$	$5$	$-0.4$	$4$
$J^{\prime }/J=0.5$	$-0.3$	$5$	$-0.2$	$4$
$J^{\prime }/J=0.5$	$-0.1$	$5$	$0$	$4$
Fig. 4 panel 4	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$
$J^{\prime }/J=0.7$	$-2$	$5$	$-1$	$4$
$J^{\prime }/J=0.7$	$-0.5$	$5$	$-0.4$	$4$
$J^{\prime }/J=0.7$	$-0.3$	$4$	$-0.2$	$4$
$J^{\prime }/J=0.7$	$-0.1$	$4$	$0$	$4$
Fig. 5	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$	$\Delta /J$	$-{\log }_{10}\left({\chi }^2\right)$
$J^{\prime }/J=1.0$	$-1$	$8$	$-0.6$	$7$
$J^{\prime }/J=1.0$	$-0.9$	$8$	$-0.5$	$7$
$J^{\prime }/J=1.0$	$-0.8$	$7$	$-0.4$	$7$
$J^{\prime }/J=1.0$	$-0.7$	$7$

Table 2. Table 2: Values of the parameters according to the parametrization used in Figs. 2 - 5 are shown. Also shown are values of − log 10 ⁡ χ 2 subscript 10 superscript 𝜒 2 -\log_{10}\chi^{2} rounded to the first digit for probability distributions corresponding to the points in Figs. 6 and 7 .

	Fig. 6			Fig. 7
	$\Delta /J$	$J^{\prime }/J$	$-{\log }_{10}\left({\chi }^2\right)$	$\Delta /J$	$J^{\prime }/J$	$-{\log }_{10}\left({\chi }^2\right)$
$A^{\ast }$	$-1.9875$	$0.9875$	$8$	$-.3995$	$.9975$	$6$
$B$	$-1.154$	$.333$	$7$	$-.3149$	$.8181$	$6$
$C$	$-.5359$	$.072$	$7$	$-.1705$	$.7047$	$6$
$D$	$0$	$0$	$6$	$0$	$.6666$	$6$
$E$	$.5359$	$.072$	$8$	$.1705$	$.7047$	$6$
$F$	$1.154$	$.333$	$7$	$.3149$	$.8181$	$6$
$G$	$2$	$1$	$7$	$.4$	$1$	$6$
$H$	$3.4641$	$3$	$6$	$.3849$	$1.2222$	$6$
$I$	$7.4641$	$13.928$	$8$	$.2419$	$1.4189$	$6$
$J$	$0$	$\infty$ $\left(1/0\right)$	$7$	$0$	$1.5$	$6$
$K$	$-7.4641$	$13.928$	$7$	$-.2419$	$1.4189$	$6$
$L$	$-3.4641$	$3$	$6$	$-.3849$	$1.2222$	$6$

Equations

$${\varphi }_{Zak}=\text{mbox}Im\ln \prod _{I=0}^{M-1}⟨{\Psi }_0(K_I)|{\Psi }_0(K_{I+1})⟩,$$

$${\left[\prod _{I=0}^{M-1}⟨{\Psi }_0(K_I)|{\Psi }_0(K_{I+1})⟩\right]}^{\Delta K}=\exp \left(\sum _{n=1}^{\infty }\frac{(i\Delta K{)}^n}{n!}{\tilde{C}}_n\right),$$

$${\tilde{C}}_1=i\sum _{I=0}^{M-1}\Delta K{\gamma }_1(K_I)$$

$${\tilde{C}}_2=-\sum _{I=0}^{M-1}\Delta K[{\gamma }_2(K_I)-{\gamma }_1(K_I{)}^2]$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$${\mu }_C^{(1)}$$

$${\mu }_C^{(2)}$$

$${\mu }_C^{(3)}$$

$${\mu }_C^{(4)}$$

$$C_{M,\alpha }=\frac{L}{2\pi }{\int }_{-\frac{\pi }{L}}^{\frac{\pi }{L}}dK\prod _{i=1}^d⟨u_{nK}|{\partial }_K^{m_i}|u_{nK}⟩,$$

$$u_{nK}(x)=\sum _{p=-\infty }^{\infty }\exp (iK(pL-x))a_n(x-pL),$$

$$\frac{L}{2\pi }{\int }_{-\pi /L}^{\pi /L}dK{\int }_0^{L}dx|u_{nK}(x){|}^2={\int }_{-\infty }^{\infty }dx|a_n(x)|=1.$$

$$⟨u_{nK}|{\partial }_K^m|u_{nK}⟩=\sum _{\Delta p=-\infty }^{\infty }\exp (-iK\Delta pL){\int }_{-\infty }^{\infty }dx$$

$$a_n^{\ast }(x-\Delta pL)(-ix{)}^m{a}_n(x).$$

$$C_{M.\alpha }=\sum _{\Delta p_1=-\infty }^{\infty }...\sum _{\Delta p_d=-\infty }^{\infty }\delta [\Delta P,0]$$

$$\prod _{j=1}^d\left\{{\int }_{-\infty }^{\infty }d{x}_j(-i{x}_j{)}^{m_j}{a}_n^{\ast }(x_j-\Delta p_{j}L)a_n(x_j)\right\},$$

$$C_2=\frac{\hslash }{\pi q_e^2{n}_0}\int \frac{d\omega }{\omega }\bar{\sigma }(\omega ),$$

$$\chi \le \frac{2q_e}{V{\Delta }_g}{C}_2.$$

$$S[P(x)]=-\int dxP(x)\ln P(x),$$

$${\mu }_P^{(k)}=\int dxP(x)x^k,$$

$$\frac{\delta }{\delta P(x)}\left[S[P(x)]-\sum _k{A}_k({\mu }_P^{(k)}-{\mu }_C^{(k)})\right]=0,$$

$$P(x)=C\exp \left(-\sum _k{A}_k{x}^k\right),$$

$${\chi }^2=\sum _k({\mu }_P^{(k)}-{\mu }_C^{(k)}{)}^2,$$

$${\hat{H}}_{SSH}=-J\sum _{i=1}^{N/2}{c}_{i,A}^{†}{c}_{i,B}-J^{\prime }\sum _{i=1}^{N/2}{c}_{i,B}^{†}{c}_{i+1,A}+\text{mbox}H.c.,$$

$${\hat{H}}_{\Delta }=-\Delta \sum _{i=1}^{N/2}{c}_{i,A}^{†}{c}_{i,A}+\Delta \sum _{i=1}^{N/2}{c}_{i,B}^{†}{c}_{i,B}.$$

$$J=\frac{t}{2}+\frac{\delta }{2},J^{\prime }=\frac{t}{2}-\frac{\delta }{2}.$$

$$\hat{H}={\hat{H}}_{SSH}+{\hat{H}}_{\Delta }.$$

$$\hat{H}=\sum_{k}\left(\begin{array}[]{cc}\Delta&-\rho_{k}\\ -\rho_{k}^{*}&-\Delta,\end{array}\right)$$

$${\rho }_k=J{e}^{ik}+J^{\prime }{e}^{-ik}.$$

$$\left(\begin{array}[]{c}\alpha_{k}\\ \beta_{k}\end{array}\right)=\left(\begin{array}[]{c}\sin\left(\frac{\theta_{k}}{2}\right)\\ e^{-i\phi_{k}}\cos\left(\frac{\theta_{k}}{2}\right)\end{array}\right),$$

$${\theta }_k$$

Full text

I Introduction

One way to derive the Berry phase Berry84 ; Pancharatnam56 ; Xiao10 is to form a product of scalar products between quantum states at different points of the space of external parameters (Bargmann invariant Bargmann64 ) and to take the continuous limit along a cyclic curve. An extension Souza00 ; Hetenyi14 of this derivation, keeping higher order terms, leads to gauge invariant cumulants (GIC) associated with the Berry phase. One is lead to ask two questions. The GICs give information of the distribution of what physical quantity? Can one reconstruct the probability distribution from the GICs?

The answer to the first question depends on the physical context in which the Berry phase is defined. In a crystalline solid the Berry phase (or Zak phase Zak89 , in this context) corresponds to the macroscopic polarization. Zak showed Zak89 that the phase itself corresponds to the expectation value of the position over a Wannier function. For the higher order GICs Souza, Wilkens and Martin Souza00 showed that they only correspond to the cumulants of the distribution of the position associated with Wannier functions, if the Wannier functions themselves are localized within the unit cell (non-overlapping among different unit cells). Indeed, in the construction of tight-binding based lattice models, one starts with a continuum description, and assumes a localized basis of non-overlapping Wannier functions (see for example, Ref. Essler05 ). In practice, however, constructing such a localized basis is not trivial Marzari97 .

The distribution of the polarization gauges the extent to which the system is localized in the full configuration space, a criterion Kohn64 which distinguishes an insulator from a conductor. The second GIC was shown Kudinov91 ; Souza00 to be proportional to the integrated frequency dependent conductivity (sum rule). A gauge dependent definition of the spread (similar to the second GIC) was used to define the maximally localized Wannier function Marzari97 . Also, the second cumulant was proposed Resta99 to distinguish conductors from insulators. In Ref. Hetenyi14 the simplest system with a Berry phase, an isolated spin-$\frac{1}{2}$ particle in a magnetic field, was considered, and it was shown that (based on calculating the first four cumulants) the moments of this underlying distribution are all equal.

The Zak phase was measured in Ref. Atala13 in an optical lattice setup which corresponds to the experimental realization of the Su-Schrieffer-Heeger (SSH) model Su79 and its extension the Rice-Mele (RM) model Rice82 . The RM model is a lattice model with an alternating on-site potential, and hoppings with alternating strengths, depending on whether a given bond is odd or even. An interesting characteristic Vanderbilt93 ; Xiao10 of the RM model is its topological behavior which manifests when an adiabatic cycle in the parameter space of the Hamiltonian encircles the point ($\Delta=0$, $J=J^{\prime}=1$). Due to the fact that the polarization as a function of the parameters of the Hamiltonian is not single valued, the polarization in such a process changes by a “polarization quantum.” A recent related study Nakajima15 realized quantized adiabatic charge pumping Thouless83 , also in the RM model.

In this paper we calculate the leading GICs associated with the Zak phase for the RM model. Based on the GICs (or associated moments, GIMs) we approximately reconstruct the distribution associated with the polarization. The RM model is a lattice model, which implies that the underlying Wannier functions are non-overlapping among different unit cells, and that the GICs correspond to the distribution associated with the Wannier function. Hence, our reconstructed probabilities correspond to the squared modulus of the Wannier function. We show that in the fully dimerized limit the GIMs should all have the same magnitude, and that the sign of odd GIMs switch sign with respect to the direction of the polarization. We also focus on the line of the parameter plane where the polarization shows a line of discontinuity (see Fig. 3, lower panel, left inset). We also present two model calculations in which the evolution of the probability distribution is followed around the topologically nontrivial point of the RM model. As expected, the distribution migrates to the next unit cell, although its shape varies considerably during the cycle.

Reconstructing a probability distribution from knowledge of a finite set of moments is an ill-posed mathematical problem which already has a long history Smoluchowski17 , although there has been a renewed interest in the last decades John07 ; DeSouza10 . The scientific applications are also quite broad; image processing Sluzek05 , calculating magnetic moments Berkov00 , or molecular electronic structure Bandyopadhyay05 . In our study, we opt for a reconstruction based on maximizing the entropy Bandyopadhyay05 ; Collins77 ; Mead84 of the underlying probability distribution.

This paper is organized as follows. In the next section we introduce the GICs associated with the Zak phase. We then discuss their connection to the distribution associated with the Wannier functions. In section IV we discuss the connection of the cumulants to response functions, after which the reconstruction procedure is presented. In section VI the Su-Schrieffer-Heeger and Rice-Mele models are introduced. Subsequently, the behavior of the moments for the fully dimerized limit is studied. Section VIII contains our results and analysis before concluding our work.

II Gauge invariant cumulants associated with the Zak phase

Consider a one-dimensional system whose Hamiltonian which is periodic in $L$. We take Bloch functions parametrized by the crystal momentum, $\Psi_{0}(K)$ on a grid of $M$ points $K_{I}=2\pi I/(ML)-\pi/L$, with $I=0,...,M-1$. The Zak phase can be derived from a product of the form

[TABLE]

by taking the continuous limit ($M\rightarrow\infty$). The product in Eq. (1) is known as the Bargmann invariant Bargmann64 . We will derive the Zak phase, as well as the associated gauge invariant cumulants (GIC). We start by equating the product in Eq. (1) to a cumulant expansion,

[TABLE]

with $\Delta K=2\pi/M$. We now expand both sides and equate like powers of $\Delta K$ term-by-term, mindful of the fact that the left-hand side includes a product over $I$. For example, the first-order term will be

[TABLE]

the second will be

[TABLE]

with $\gamma_{i}(K)=\langle\Psi_{0}(K)|\partial^{i}_{K}|\Psi_{0}(K)\rangle$. Straightforward algebra and taking the continuous limit ($\Delta K\rightarrow 0$, $M\rightarrow\infty$) gives up to the fourth order term,

[TABLE]

The quantities $C_{n}$ in Eq. (5) are the GICs associated with the Zak phase (the Zak phase itself being equal to $C_{1}$). The difference between $\tilde{C}_{i}$ and $C_{i}$ is the multiplicative factor $L/2\pi$, which is also how the phase is defined by Zak Zak89 . This assures that the first moment corresponds to the average position associated with square modulus of the Wannier function (Eq. (10) in Ref.  Zak89 ). When the underlying probability distribution is well defined the associated moments can be defined based on the cumulants. Following this standard procedure we also define a set of moments. For the first four moments the expressions are

[TABLE]

As discussed below, when the Wannier functions of a particular model are localized within the unit cell, these moments correspond to the moments of the polarization, alternatively, to the distribution of the Wannier functions themselves.

We remark that in general, the Berry phase is a physically well-defined observable, which is thought not to correspond to an operator acting on the Hilbert space. The Zak phase, however, is known to correspond to the total position, and is the basic quantity in expressing the polarization in the modern theory King-Smith93 ; Resta94 ; Resta98 .

III Connection to the distribution of Wannier centers

Cumulants of the type described in the previous section appear in the theory of polarization Souza00 . In this section we connect the cumulants to the distribution of Wannier centers. We consider a typical term contributing to cumulant $C_{M}$, which can be written in the form

[TABLE]

where $\sum_{i=1}^{d}m_{i}=M$ and where we have used the periodic Bloch functions $u_{nK}(x)$ as a basis. The periodic Bloch functions can be written in terms of Wannier functions,

[TABLE]

where $a_{n}(x)$ denote the Wannier functions. With this definition it holds that

[TABLE]

We can rewrite a scalar product appearing in Eq. (7) as

[TABLE]

Substituting Eq. (10) $C_{M,\alpha}$ and integrating in $K$ results in

[TABLE]

where $\Delta P=\sum_{j=1}^{d}\Delta p_{j}$ and $\delta[\Delta P,0]$ is a Kronecker delta.

We note that if the Wannier functions are localized in one unit cell, then the summation in the scalar product of Eq. (10) will be restricted to the term $\Delta p=0$. In this case, the cumulants $C_{M}$ will correspond to those of the Wannier centers.

IV Relation to response functions

The second GIC associated with the polarization gives a sum rule for the frequency-dependent conductivity. This was shown for a finite system by Kudinov Kudinov91 , and the derivation was extended to periodic systems by Souza, Wilkens, and Martin Souza00 , by replacing the ordinary matrix elements of the total position operator by their counterparts valid in the crystalline case. Their result is

[TABLE]

where $q_{e}$ denotes the charge, $n_{0}$ the density, and $\bar{\sigma}(\omega)=(V/8\pi^{3})\int d{\bf k}\sigma^{\bf k}(\omega)$.

For an insulating (gapped) system one can show that the second cumulant provides an upper bound for the dielectric susceptibility, $\chi$. This was shown by Baeriswyl Baeriswyl00 for an open system. This derivation is also easily extended to periodic systems by the appropriate replacement of the total position matrix elements, resulting in,

[TABLE]

In this equation $\Delta_{g}$ denotes the gap, $V$ denotes the volume of the system.

For higher order cumulants, the derivation of relations such as Eq. (13) are not possible. However, in the classical limit, the cumulants correspond exactly to the response functions of their respective order ($C_{2}$ gives $\chi$, $C_{3}$ gives the first non-linear response function, etc.).

V Reconstruction of the probability distribution

If the Wannier functions can be assumed to be localized within a unit cell, the moments calculated based on the GICs correspond to the actual moments associated with the Wannier orbitals. If all the moments are known, the full probability distribution can be reconstructed. However, in practice, usually only a finite number of cumulants are available. In this case the cumulants can be used as constraints to improve the form of the probability distribution. The first and second cumulants give the average and the variance, and if only these two are available, the best guess for the probability distribution is a Gaussian. Higher order cumulants refine this guess. The third cumulant (skewness) provides information about the asymmetry of the distribution around the mean, while the fourth order one, (kurtosis) represents how sharp the maximum of the distribution is approached from either side.

Below we calculate the GICs of the Rice-Mele model, which is a lattice model (in other words, the Wannier functions are completely localized on particular sites), and approximately reconstruct the probability distribution of the polarization. Our reconstruction is based Bandyopadhyay05 ; Collins77 ; Mead84 on maximizing the information entropy under the constraints provided by the moments calculated. The expression for the entropy we use is

[TABLE]

minimized as a functional of $P(x)$ under the constraints

[TABLE]

as well as the constraint that $P(x)$ is normalized. The functional minimization of Eq. (14) under the constraints results in the functional differential equation

[TABLE]

where $\mu_{C}^{(k)}$ are the moments obtained from the cumulants of the Berry phase (see Eq. (6)), and $A_{k}$ are Lagrange multipliers. The solution of Eq. (16) is

[TABLE]

where $C$ is the normalization constant. We determine the constants $A_{k}$ by numerically minimizing the quantity

[TABLE]

as a function of $A_{k}$. As our initial guess in all cases studied below, we take the Gaussian distribution defined by the first two cumulants obtained for the particular case. The minimization procedure we applied is the simulated annealing technique Kirkpatrick83 . Below our reconstructions are based on calculating the first six GIMs in all cases.

VI Su-Schrieffer-Heeger and Rice-Mele models

The SSH model was first introduced Su79 to understand the properties of one-dimensional polyacetylene. The RM model is an extension of the SSH model, it includes an additional term, consisting of an alternating on-site potential, added in order to extend the SSH model to diatomic polymers. In recent decades it has been studied extensively due to the wealth of interesting physical phenomena it displays: topological soliton excitation, fractional charge, and non-trivial edge statesTakayama80 ; Su80 ; Jackiw76 ; Heeger88 ; Ruostekoski02 ; Li14 . It was also realized as a system of cold atoms trapped in an optical lattice in one dimension recently Atala13 . The Berry phase in the RM model was studied by Vanderbilt and King-Smith Vanderbilt93 . In that study the point of the model in parameter space of the model which is metallic (and which is responsible for the topologically nontrivial behavior) was encircled in parameter space. This leads to the increase of $C_{1}$ (the Berry phase, or the polarization) by one polarization quantum, consistent with the quantization of charge transport Thouless83 ; King-Smith93 .

The hopping part of the SSH Hamiltonian reads:

[TABLE]

where $N$ denotes the number of sites, the on-site potential has the form

[TABLE]

The model is shown schematically in Fig. 1. This figure shows the one-dimensional lattice, including sublattices, the alternating hoppings, and the on-site potential. The unit cell is indicated in shaded yellow. Also shown is the continuous variable $x$, which runs from $-\infty$ to $\infty$, and will serve as the axis for the reconstructed probability distributions of the polarization calculated below.

The hoppings can also be expressed in terms of the average hopping $t$ and the deviation $\delta$ as

[TABLE]

The total Hamiltonian we consider is

[TABLE]

The parameters $J$ and $J^{\prime}$ are hopping parameters corresponding to hopping along alternating bonds. We take the lattice constant to be unity (the unit cell is two lattice constants). The parameter $\Delta$ denotes the on-site potential, whose sign alternates from site to site. This model is metallic for $J=J^{\prime}$ and $\Delta=0$ but is insulating for all other values of the parameters. In reciprocal space this Hamiltonian becomes

[TABLE]

where

[TABLE]

At a particular value of $k$ we can write the eigenstate for the lower band as

[TABLE]

where

[TABLE]

The cumulants can now be written in terms of the eigenstates. For example,

[TABLE]

and the other cumulants can be constructed accordingly (note that the unit cell is $L=2$).

VII Fully dimerized limit

Here we show that in the fully dimerized limit the GIMs should all have the same magnitude. In Ref. Hetenyi14 we pointed out that the Berry phase can be related to an observable $\hat{O}$ fixed by requiring that

[TABLE]

This definition does not uniquely fix the operator $\hat{O}$. For example, for the magnetic field example the matrix $\sigma_{z}/2$ or $(\sigma_{z}+I)/2$ both satisfy Eq. (28). This arbitrariness causes a shift in the first cumulant. However, only the operator $(\sigma_{z}+I)/2$ will give a distribution in which all moments are equal, since this matrix has the form

[TABLE]

and is equal to itself when raised to any power.

In the case of the RM model we first write the Hamiltonian with the parameter $K$ explicitly as

[TABLE]

The operator $\partial_{K}\hat{H}(K)$ is the current,

[TABLE]

We now write a form for the operator $\hat{O}$ as

[TABLE]

Evaluating the commutator gives

[TABLE]

For the case $J^{\prime}=0$ we can chose $x_{j}=0$ and $y_{j}=1$, so that $i[\hat{H}(K),\hat{O}]$ corresponds to the current. This is not the only choice, but with this choice the operator $\hat{O}$ when written in $k$-space corresponds to

[TABLE]

which gives equal moments. Clearly, the choice for the spatial coefficients $x_{j}$ and $y_{j}$ are due to the fact that in this case the system consists of a set of independent dipoles. When $J$ is taken to zero, and $J^{\prime}$ kept finite, then the appropriate choice to fix $\hat{O}$ is $x_{j}=0$ and $y_{j}=-1$. If instead the sign of $\Delta$ is changed $\hat{O}$ is again defined by the $x_{j}=0$ and $y_{j}=-1$. These results are clearly due to the reversal of the direction of the dipole moment within the unit cell. The results presented in Fig. 2 corroborate our derivation.

VIII Results and Analysis

We first look at the system with $J^{\prime}=0$. In this case, the band structure of the system is simply two flat lines in the Brillouin zone. The system can be thought of as a simple two-state system. We calculated the first four GICs, from which we obtained the corresponding GIMs. The results are shown in the uppermost panel of Fig. 2. The moments as a function of $\Delta/J$ all fall on the same curve in this case. If the hopping parameters $J$ and $J^{\prime}$ are switched (not shown), the sign of the odd moments changes, the even moments remain the same. These results are in accordance with section VII.

Fig. 2 also shows the cumulants for different ratios; $J^{\prime}/J=0.3,0.5,0.7$. The deviation of the cumulants from one another is more pronounced, and increases with an increase of $J^{\prime}/J$. However, the moments become equal for any $J^{\prime}/J$ when $\Delta\rightarrow\pm\infty$. In this case also, the system becomes an independent array of two state systems. The band energies in all these cases vary continuously with $k$ across the Brillouin zone.

The results for the case $J^{\prime}/J=1$ are also shown separately in Fig. 3, as well as the limits $J^{\prime}\rightarrow_{\pm}J$. For finite $\Delta/J$ the odd cumulants are zero, indicating an even probability distribution. The ratio of the second and fourth cumulants rule out a Gaussian. As $\Delta/J\rightarrow 0$ a discontinuity in the slope of the band develops. In this case, the cumulants $C_{2}$ and $C_{4}$ diverge. The lower panel in this figure shows what happens when $J^{\prime}$ is close to $J$ (bigger or smaller) but the two are not quite equal ($J^{\prime}=J+\epsilon$, $\epsilon$ a small number). We see that in this case the first moment is one or minus one, depending on the sign of $\epsilon$, and zero is not approached as $\epsilon\rightarrow 0$ from either side. The left inset in the lower panel shows the behavior of the first moment on the $\Delta-\delta$ plane, indicating a discontinuity along the line $\Delta<0,\delta=0$ (the well-known result of Vanderbilt and King-Smith Xiao10 ; Vanderbilt93 ). The moments and cumulants we find are consistent with the behavior shown in the left inset of the lower panel of the figure.

In Fig. 4 we show examples of reconstructed probability distributions for $J^{\prime}/J=0.0,0.3,0.5,0.7$, in each case for several values of $\Delta/J$. $\chi^{2}$ (defined in Eq. (18)) is tabulated in the appendix (Table 1). The most localized example ($J^{\prime}/J=0$ and $\Delta/J=-2$) shows a sharp peak around $x=1$; as $\Delta$ decreases the curves shift to the left and spread out, but their shape is always very similar (for smaller values of $\Delta/J$ this is emphasized in the inset). The maximum of the probability distribution is always between zero and one. These curves are all cases for which all the moments are equal. As the alternate hoppings ($J^{\prime}$) are turned on, the shifting occurs in a qualitatively different manner. Initially ($\Delta/J=-2$ in all cases) the curves are centered very near $x=1$. $\Delta/J=-2$ is for most cases well in the region where the moments are equal. As $\Delta$ decreases, the distributions shift, but they do this by becoming asymmetric about their mean, with the density increasing on the side left of the maxima of the distributions. The shape of the distributions changes considerably. This is clearly due to the fact that in these latter cases the moments vary as $\Delta$ is varied, and they are not all equal. The maxima for the cases for which $J^{\prime}/J\neq 0$ shift much less as $\Delta/J$ is varied. When $\Delta/J$ changes sign (results not shown), the polarization becomes centered around $x=0$ end of the unit cell and the probability distributions are reflections of the ones shown in Fig. 4 across $x=1/2$.

The probability distributions for the case $J^{\prime}=J=1$ are also shown separately in Fig. 5 (with $\chi^{2}$ tabulated in Table 1), as well as the case $J^{\prime}$ close to $J$. All of the $J^{\prime}=J$ distributions are symmetric around the origin. As $\Delta/J\rightarrow 0$ the distribution broadens, and it is clear that a conducting phase is approached Resta99 . If $\Delta<0$ then the polarizations are localized near $x=\pm 1$ depending on whether $J^{\prime}$ is smaller or larger than $J$. This is consistent with Fig. 3.

In Figs. 6 and 7 we show the evolution of the reconstructed probability distributions along two cyclic paths which encircle the topologically non-trivial point of the RM model, one with radius unity, the other with radius $0.2$ in the $\Delta/t$, $\delta/t$ plane. In these calculations the parametrization was different from the previous ones, here $t$ was set to unity, rather than $J$ (see Eq. (21)). For the points $A^{*}$, $B$, … in Figs. 6 and 7 the values of $\Delta/J$ and $J^{\prime}/J$ are shown in Table 2. The upper panel in both figures show the evolution of the different GIMs(GICs). The even moments are single-valued, the odd ones are not. This follows from gauge invariance properties of the cumulants (Eq. (5)). The first cumulant is only gauge invariant modulo $2\pi$ times an integer Xiao10 ; Vanderbilt93 , the others do not change at all due to a gauge transformation. The odd GIMs depend on combinations of the GICs which involve odd combinations of the cumulants, therefore they are not multivalued in general. In both sets of figures (6 and 7) the points $A^{*}$ are not exactly on the $\phi=-\pi/2$ axis, but instead we numerically realize the limit $\phi=\lim_{\delta\Phi\rightarrow 0^{+}}(-\pi/2+\delta\Phi)$. In the actual calculation we took $\delta\Phi=2\pi/1000$. Also, the point $\phi=-\pi/2$ or $\phi=3\pi/2$ is excluded from the curves shown in the upper panels of the two figures.

The example with radius unity (Fig. 6) remains mostly in the fully dimerized limit, as can be seen in the upper panel of the figure. The odd moments and even moments are always equal. Except for a small region near $\phi/\pi=0.5$ the absolute values of the moments are equal. The lower panel shows the evolution of the probability distribution along the path. Starting from a relatively sharp distribution localized near $x=1$, the maximum moves to the left. Before reaching half the unit cell, the distribution spreads. After passing through the midpoint the system, where the maximum is the smallest, the distribution begins to localize again until $x=0$. From there this tendency is repeated. Indeed, the distribution ends up at $x=-1$ at the end of the process: the Wannier function “walked” to an equivalent position in the next unit cell. For the case of the smaller radius (0.2, Fig. 6) the initial distribution is broader, and as the cycle is traversed, the maximum of the distribution oscillates with a smaller amplitude, but the “walking” to a new equivalent position still occurs.

In both Figs. 6 and 7 it is clear that the odd moments do not correspond to single-valued functions. The values of the odd moments depend on whether we approach the original point from which the cycle begins ($\delta=0,\Delta<0$) from the left or the right. At the same time, the probability distributions for some cases with $\delta=0,\Delta<0$ are shown in Fig. 5; they are centered around zero and they spread as $\Delta/J$ is decreased. This suggests the limiting cases from either direction give different results from the result for fixing the Hamiltonian parameters such that $\delta=0,\Delta<0$.

IX Conclusion

We studied the gauge invariant cumulants associated with the Zak phase. We have shown that for localized Wannier functions they correspond to the cumulants of the Wannier centers. They are also related to the dielectric response functions of a given system. We calculated the cumulants for the Rice-Mele model. In the limit of isolated dimers, all the moments (extracted from the gauge invariant cumulants) are equal. This can be justified for this case by constructing the operator which corresponds to the Berry phase explicitly. Deviations from this behavior come about when the hopping parameters are both finite. For a system with equal hopping parameters the odd cumulants vanish. We have also reconstructed the full probability distribution of the polarization based on the gauge invariant cumulants and have studied how they evolve as functions of different parameters of the Hamiltonian. In particular we calculated the evolution of the distribution around the topologically non-trivial point of the model. We anticipate that detailed experimental measurements can also provide a probability distribution of the polarization for comparison with our predictions.

Acknowledgments

The authors acknowledge financial support from the Turkish agency for basic research (TÜBITAK, grant no. 113F334). We also thank L. G. M. de Souza for helpful discussions on the topic of probability reconstruction from moments.

Appendix

In Table 1 values of the negative base ten logarithm of $\chi^{2}$ rounded down to the first digit (defined in Eq. (18)) is shown for the reconstructed probabilities in Figs. 4 and 5. In all cases $\chi^{2}$ decreased at least eight orders of magnitude from its initial value during the simulated annealing calculation.

In Table 2 the values of the parameters according to the parametrization used in Figs. 2-5 are shown. Also shown are values of the negative base ten logarithm of $\chi^{2}$ rounded down to the first digit for probability distributions corresponding to the points in Figs. 6 and 7.