text:
# Berry-phase treatment of the homogeneous electric field perturbation in insulators

## I Introduction

Until the early nineties, the formulation of a quantum-mechanical theoretical framework for the study of the physics of electric polarization in solids had remained a challenging problem.  Even the definition of the polarization itself as a bulk quantity, independent of surface termination, was the subject of heated debate.

This picture changed when King-Smith and Vanderbilt (KS-V)  proposed a formulation (the modern theory of polarization - MTP), which resolved the conceptual difficulties associated with the definition of this quantity for continuous, periodic, charge distributions. In their work, the electric polarization of an insulating crystal is related to a Berry phase  computed from the valence wavefunctions. The existence of a band-structure Berry phase had already been discussed by Zak and coworkers,  before its connection with the electronic polarization was established by KS-V. Besides settling the important conceptual question related to the definition of polarization as a bulk quantity, the KS-V theory provided an entirely new framework for the computation of the polarization of a crystal maintained at vanishing homogeneous electric field. Since its formulation, the theory has been examined in greater detail by KS-V  and by Resta , and extended to many-body systems by Ortiz and Martin . The relation between polarization and the phases of the wavefunctions has also led to a reexamination of the role this quantity plays in the Density Functional Theory (DFT)  formulation of the ground-state properties of extended systems.

Of no less importance is the conceptual relationship between the spontaneous polarization and the centers of charge of the Wannier functions (WF) of the occupied bands, which was also discussed by KS-V  and previously by Zak.  This connection was later generalized by Nunes and Vanderbilt  (NV) to deal with insulators placed in non-zero external homogeneous electric fields: they introduced field-dependent “polarized” WF’s and a method for their computation. NV argued that, in the static-response regime, the state of an insulator under an external homogeneous electric field is one in which the periodicity of the charge density is retained, despite the fact that the perturbation lacks the lattice-translational symmetry of the unperturbed crystal. Such state is actually a long-lived resonance of the system, as rigorously demonstrated by Nenciu .

It is well known that within DFT , ground-state properties of condensed-matter systems, such as equilibrium lattice parameters, bond lengths and bond angles, among others, can be obtained with an accuracy of a few percent in comparison with experimental results. Within DFT , the NV method has recently been applied to the computation of the polarized WF’s and the dielectric constant of silicon and gallium arsenide . The latter quantity is related to the change of polarization due to a change of homogeneous electric field, in the linear regime, or equivalently, to the second-derivative of the total energy with respect to the homogeneous electric field.

Specific treatments have been developed (noticeably within DFT) for the study of the response of crystals to “external” perturbations, like phonons, stresses or homogeneous electric fields. The latter, on which we focus exclusively, can be taken as a homogeneous field or as the limit of long-wavelength perturbations.

In the so-called direct approach , supercell calculations are employed to study both the unperturbed and perturbed systems, with the response functions being obtained by numerical finite-difference analysis of the changes induced by a long-wavelength perturbation applied to elongated supercells. The non-linear response regime is directly accessible, although it must be disentangled from the linear response of the system. However, because of the use of supercells, the computational cost of this approach is rather high.

Alternatively, the specific response to a homogeneous electric field was considered within perturbation theory, already in the sixties. In the Random-Phase Approximation  (no local-field effects  included), the response of the wavefunctions is obtained through a sum over states, involving matrix elements of the position operator between valence and conduction states . This technique was generalized to the computation of second- and third-order susceptibilities . The need to compute many unoccupied bands is the bottleneck of this method.

Local-field effects can be reintroduced on top of such a sum-over-states approach either in a matrix-inversion framework , or in an iterative approach . In their calculation of linear susceptibilities, Levine and Allan included a so-called “scissor-operator” correction, that was understood later to compensate some deficiencies of local-density approximation computation of long-wavelength response functions . Also, Levine and coworkers proposed rather involved expressions for the second-order and third-order DFT susceptibilities .

To a large extent, Density Functional Perturbation Theory (DFPT)  overcome the limitations of the previously mentioned approaches, at the price of non-negligible additional coding. At the lowest order in the homogeneous-electric-field perturbation, this method was introduced by Baroni and collaborators . It is based on an iterative solution for the first-order change in the wavefunctions, which allows for the self-consistent inclusion of local-field effects, besides eliminating the cumbersome sum over conduction bands. It does not employ supercells, and can be applied to perturbations of arbitrary wavelengths. In the DFPT, the computational workload involved in the computation of linear-response functions is of the same order of that involved in one ground-state calculation.

DFPT is part of a class of formalisms in which perturbation theory is applied to a variational principle . This interesting combination leads to a generic “2n+1” theorem , as well as variational properties of even-order derivatives of the energy . For example, one can compute the third-order derivative of the energy from the first-order derivative of wavefunctions, and the fourth-order derivative of the energy is variational with respect to the second-order derivative of wavefunctions. The expressions derived in this framework are surprisingly simple and can be formulated at all order of perturbations.

However, the treatment of homogeneous electric fields in this variation-perturbation framework is plagued by difficulties similar to those encountered in the theory of polarization. Shortly after the appearance of the MTP, Dal Corso and Mauri , building upon the NV work, proposed a very concise expression for the second-order susceptibility which was later applied successfully to compute this quantity for a variety of systems .

In the present work, we formulate a perturbation theory of the static response of insulating crystals to homogeneous electric fields that combines the conceptual ideas of the MTP with the variation-perturbation theory. A major achievement of our work is the presentation of formulas valid for unrestricted order of perturbation theory. We also examine the low-order expressions in some detail, and recover the expression proposed by Dal Corso and Mauri .

The theory is worked out directly in reciprocal space, in terms of the Berry phase associated with the occupied bands of the perturbed crystal, in the manner of the MTP. The conceptual issues involved in the definition of a perturbative approach for the problem are addressed. The Berry phase is argued to remain a valid concept in the presence of the periodicity-breaking electric field. The periodicity of the charge density is assumed to survive the application of the field, and the Berry phase is obtained from the associated polarized periodic wavefunctions. By working out the perturbative approach in terms of these polarized states, we obtain very compact expressions for the high-order dielectric-response functions of the crystal. These can be numerically obtained on the basis of iterative equations for the second- and higher-order terms in the perturbation expansion of the wavefunctions of the system, as in the DFPT approach for other types of perturbation. We will not deal explicitly with the exchange and correlation parts of the DFT functional: the main difficulties that we want to address in this paper are not related to them. The formalism can be extended to include exchange and correlation terms in a self-consistent fashion, in the manner presented in Refs .

Any application of the MTP involves a discretization of the Berry-phase expression, in terms of a series of wavevectors for the electronic wavefunctions. We have discovered that such a discretization, that appears naturally also in the present framework, can be performed at two different conceptual levels when merged with perturbation theory: either after the derivation of formal expressions at different orders of perturbation theory, starting from a continuous Hamiltonian, or at the level of the field-dependent Hamiltonian itself, before any perturbation expansion is performed. We will refer to the first approach as the discretization after perturbation expansion (DAPE) formulation, and to the second as the perturbation expansion after discretization (PEAD) formulation.

In order to judge the relative merits of these two approaches, and also the correctness of the global framework, we analyze the behavior of a model tight-binding Hamiltonian, for which analytical responses to an electric-field perturbation have been obtained up to the fourth order.

The paper is organized as follows. In the next section, we address the conceptual issues associated with the definition of a Hamiltonian and its perturbative expansion for the electric-field problem. Sec. III summarizes the main results of the variational perturbation theory which are used in this work. In Sec. IV, we work out the continuous formulation of the problem and its perturbation expansion, from which we obtain the DAPE version of our theory. In Sec. V, we work from the start using a discretized expression for the polarization, which leads to the alternative PEAD formulation. The theory is applied to a model one-dimensional system in Sec. VI.

## II Insulators in an electric field: conceptual considerations

### A The modern theory of polarization

In the MTP , the change in electric polarization per unit volume induced by an adiabatic change in the crystalline potential (the self-consistent Kohn-Sham potential in the context of DFT) is written

$`𝚫𝐏={\displaystyle _{\lambda _1}^{\lambda _2}}{\displaystyle \frac{𝐏}{\lambda }}𝑑\lambda =𝐏(\lambda _2)𝐏(\lambda _1),`$ (1)
with $`𝐏(\lambda )`$ given in terms of a Berry phase associated with the occupied bands of the system

$`𝐏(\lambda )={\displaystyle \frac{2ie}{(2\pi )^3}}{\displaystyle \underset{n}{}}{\displaystyle _{BZ}}𝑑𝐤u_{n𝐤}^{(\lambda )}\left|_𝐤\right|u_{n𝐤}^{(\lambda )},`$ (2)
where $`e`$ is the electron charge, $`\lambda `$ is a parameter representing the adiabatic change in the potential, and the factor of 2 in the numerator accounts for spin (in this work we will consider only spin-unpolarized systems).

The gauge relation between periodic functions $`u_{n𝐤+𝐆}(𝐫)=e^{i𝐆𝐫}u_{n𝐤}(𝐫)`$ is established by requiring that the Bloch eigenstates be periodic in reciprocal space, i.e. $`\psi _{n𝐤}=\psi _{n𝐤+𝐆}`$, where $`𝐆`$ is a reciprocal-lattice vector. With this choice of gauge, the polarization changes given by Eq. 1 are defined to within a factor of $`(2e/\mathrm{\Omega })𝐑`$, where $`𝐑`$ is a lattice vector.  Eq. 2 was derived under the restriction that the macroscopic electric field inside the crystal vanishes. Moreover, it also requires that the set of wavefunctions be differentiable with respect to k

The actual evaluation of the polarization in Eq. 2 is carried out on a discrete mesh of points in reciprocal space. Because this expression depends on the phase relationships of wavefunctions at different k-points, the following discretized version was proposed by KS-V:

$`P_{}(\lambda )={\displaystyle \frac{e}{4\pi ^3}}{\displaystyle 𝑑𝐤_{}\underset{j=1}{\overset{N_k}{}}\mathrm{Im}\left\{\mathrm{ln}det\left[u_{n𝐤_j}^{(\lambda )}|u_{m𝐤_{j+1}}^{(\lambda )}\right]\right\}};`$ (3)
where $`P_{}`$ is the component of the polarization along the direction of a short reciprocal-lattice vector, $`𝐆_{}`$, and $`N_k`$ is the number of k-points sampling the Brillouin zone along that direction for each value of $`𝐤_{}`$, with $`𝐤_j=𝐤_{}+j𝐆_{}/N_k`$.

From a calculational point of view, this discretized expression ensures that the final result is unaffected by random numerical phases which may be introduced in the wave functions at different k-points, when these are independently determined by the diagonalization routine. However, Resta has taken the view that the discretized expression is to be regarded as more fundamental than the continuous form.  For the formulation of the electric-field response that we develop in the present work, discretization is crucial in order to define a numerically-stable minimization procedure. We will come back to this point in Sec. II D.

The Berry-phase expression can be transformed into a real-space integral involving the Wannier functions of the occupied bands, leading to a physically-transparent expression for the polarization in terms of the centers-of-charge of the Wannier functions:

$`𝐏(\lambda )={\displaystyle \frac{2e}{\mathrm{\Omega }}}{\displaystyle \underset{n}{}}{\displaystyle 𝐫\left|w_n^{(\lambda )}\right|^2𝑑𝐫},`$ (4)
where $`\mathrm{\Omega }`$ is the unit-cell volume.

In principle, the above expressions are valid only at vanishing electric field. However, it was soon realized  that Eq. 4 could be extended to the non-zero field problem, by introducing the so-called polarized Wannier functions. Polarization effects were then related to the field-induced shifts in the Wannier-function centers of charge. In the present work, Eq. 2 is argued to apply to the non-zero-field problem, thus defining a field-dependent Berry phase containing the polarization effects. This allows us to work out a perturbative approach for the finite-field problem.

### B Definition of the energy

The study of the problem of insulators under an external electric field has traditionally met with conceptual difficulties, related to the non-analyticity of the perturbation, as discussed in detail by Nenciu . Upon the application of the external field, the spectrum of electronic states changes non-analytically, with the band structure of the insulator at zero field being replaced by a continuum of eigenvalues spanning the entire energy axis (from $`\mathrm{}`$ to $`+\mathrm{}`$), even for a field of infinitesimal strength. From a mathematical point of view, the unbound nature of the perturbation term, $`e𝐫`$ (hereafter, we use $``$ for the magnitude of the electric field), hinders the straightforward application of perturbation theory, since the diagonal elements of the position operator in the basis of the unperturbed Bloch states are ill-defined . Strictly speaking, an infinite crystal in the presence of an external electric field does not have a ground state .

From the physical point of view, in the limit of weak to moderate fields, the tunneling currents (which destroy the insulating state at sufficiently high fields) can be neglected, and only the polarization of the electronic states by the external field is considered. One is left with a picture of the problem in which the insulating state of the unperturbed crystal is preserved, hence the band structure and the periodicity of the charge density are retained. The problem is thus physically well defined. The theory we develop below will concern this periodic polarized insulating state, resulting from the application of the electric field. Avron and Zak  have discussed the stability of the band structure in the present context, while, as mentioned in the introduction, Nenciu  has rigorously shown such state to be a long-lived resonance of the problem, in the regime of weak fields.

Building on the conceptual framework set by the MTP, NV introduced a practical real-space method to handle the problem on the basis of the WF formulation of the polarization. They noted that this polarized insulating state can be represented in real space by a set of field-dependent polarized WF’s. In this way, the relationship between the WF’s centers of charge and polarization can be extended to the non-zero field situation, and an energy functional is defined as follows

$$E[\{w^{}\},]=E^{(0)}\left[\{w^{}\}\right]\mathrm{\Omega }𝐏\left[\{w^{}\}\right],$$
(5)
where $`\{w^{}\}`$ is the set of field-dependent Wannier functions.

Because the state underlying the above expression is not a ground-state, rather a resonance, the energy functional is only well defined for WF’s of finite range.  The truncation of the WF’s provides a mathematical procedure for the regularization of the problem. We will come back to this pathology in section II D

As a corollary of the existence of polarized WF’s, we consider now the representation of the system in terms of polarized Bloch orbitals. In the following sections, we develop two alternative formulations in which the $`k`$-space MTP expressions for the electronic polarization, Eqs. 2 and 3, are extended to the non-zero-field problem. We will also have to regularize the ensuing expressions, which is possible by means of discretization of the $`k`$-space integrals in such a way that polarized valence bands are stable against mixing with the conduction bands.

For simplicity, we concentrate on a one-dimensional non-interacting spin-unpolarized system. The Hamiltonian for the unperturbed periodic insulator is given by

$`H^{(0)}=K+V_0,`$ (6)
where $`K`$ is the kinetic-energy operator, and $`V_0`$ is a periodic potential, i.e., $`[V_0,T_{\mathrm{}}]=0`$, where $`T_{\mathrm{}}`$ denotes a translation by a lattice vector. The zero-field function $`u_{nk}^{(0)}`$ is the periodic part of the unperturbed Bloch-orbital, obeying the eigenvalue equation $`H_k^{(0)}|u_{nk}^{(0)}=\epsilon _{nk}^{(0)}|u_{nk}^{(0)}`$, where

$`H_k^{(0)}=e^{ik\widehat{x}}H^{(0)}e^{ik\widehat{x}}`$ (7)
is the unperturbed cell-periodic Hamiltonian.

Under the action of an external electric field, the Hamiltonian becomes

$`H=H^{(0)}+e\widehat{x}.`$ (8)
As discussed above, this Hamiltonian is not amenable to a conventional perturbation treatment. We can arrive at an expression that is applicable to extended periodic systems, from the following considerations.

First, let us assume that we can define a set of field-dependent cell-periodic functions, representing the polarized state of the system. These are the Fourier transform of the field-dependent Wannier functions introduced by NV. Eq. 2 is thus extended to the non-zero field problem, with

$`P\left(\right)={\displaystyle \frac{ie}{\pi }}{\displaystyle \underset{n}{}}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑ku_{nk}^{}\left|{\displaystyle \frac{}{k}}\right|u_{nk}^{}.`$ (9)
defining a field-dependent polarization, including the spontaneous and induced parts of this quantity.

Next, combining this definition with Eq. 5 we write

$`E=E^{(0)}aP\left(\right)={\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}\left[u_{nk}^{}\left|H_k^{(0)}\right|u_{nk}^{}+u_{nk}^{}\left|\left(ie{\displaystyle \frac{}{k}}\right)\right|u_{nk}^{}\right],`$ (10)
for the total energy in terms of field-dependent cell-periodic functions, with the unit-cell volume $`\mathrm{\Omega }=a`$ for our 1D system.

Consider now the expansion of the set $`\{u_{nk}^{}\}`$ in terms of the complete set of zero-field periodic functions:

$`|u_{nk}^{}={\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}C_{nm}^{}(k)|u_{mk}^{(0)}.`$ (11)
In terms of this expansion, Eq. 10 is written

$`E={\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \underset{m,m^{}=1}{\overset{\mathrm{}}{}}}C_{nm}^{}(k)C_{nm^{}}^{}(k)\left[u_{mk}^{(0)}\left|H_k^{(0)}\right|u_{m^{}k}^{(0)}+u_{mk}^{(0)}\left|\left(ie{\displaystyle \frac{}{k}}\right)\right|u_{m^{}k}^{(0)}\right].`$ (12)
This expression is in the form of the matrix representation, in the basis $`\{u_{mk}^{(0)}\}`$, of the expectation value of the following “operator”:

$`H_k=H_k^{(0)}+\left(ie{\displaystyle \frac{}{k}}\right).`$ (13)
This suggests Eq. 13 as an ansatz for the cell-periodic Hamiltonian, now including the perturbation term $`U_k^{pert}=ie\frac{}{k}`$.

In subsection II C, we show that the perturbation Hamiltonian operator given in Eq. 12 by its matrix representation,

$`𝒰_{mn}(k)=u_{mk}^{(0)}\left|\left(i{\displaystyle \frac{}{k}}\right)\right|u_{nk}^{(0)},`$ (14)
is the periodic part of the $`\widehat{x}`$ operator.

Another way of arriving at this ansatz is by the following argument. Consider the first-order change in the total energy which can be obtained, without postulating the existence of the field-dependent functions, by combining Eqs. 2 and 5, as follows:

$`E^{(1)}=aP^{(0)}={\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(0)}\left|\left(ie{\displaystyle \frac{}{k}}\right)\right|u_{nk}^{(0)}.`$ (15)
This is simply the coupling of the spontaneous polarization of the system to the external field. From textbook perturbation theory, the first order change in the energy is given by the diagonal matrix elements of the perturbation term, leading again to the form of the Hamiltonian in Eq. 13.

Some remarks are needed about the application of perturbation theory to this Hamiltonian. Strictly speaking, $`𝒰_{mn}(k)`$ is not an operator (its transformation properties under unitary transformations will be discussed in subsection II C). However, it can be shown  that this is reflected only in the first-order change of the single-particle eigenvalues, when we analyze single-particle quantities obtained in the perturbation expansion of Eq. 13. All other single-particle quantities, such as wave-function derivatives and higher-order eigenvalue derivatives are actually gauge-invariant, and thus well defined. More importantly, in the following developments we will be interested only in quantities that are integrated over the Brillouin zone (the derivatives of the total-energy with respect to the applied field), which will be shown to be gauge invariant. As a final observation, it can also be shown  that under unitary transformations, the first-order eigenvalue acquires a change which, when integrated over the Brillouin zone, leads to a first-order energy derivative which is defined modulo the quantity $`e\mathrm{}`$, where $`\mathrm{}=Na`$ is a lattice vector ($`N`$ is an integer). This is consistent with the fact that, in the MTP, the zero-field polarization itself is defined modulo $`e\mathrm{}`$.

The continuous formulation of our theory is based on the application of a variational perturbation treatment to Eqs. 10 and 13. Alternatively, the PEAD formulation is derived by applying the variational principle to the total energy written in terms of the discretized form of the polarization. In this case, we combine Eqs. 3 and 5 to write

$`E[\left\{u_{nk_j}\right\};]={\displaystyle \frac{2}{N_k}}\left\{{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}u_{nk_j}\left|H_{k_j}^{(0)}\right|u_{nk_j}{\displaystyle \frac{e}{\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\left\{\mathrm{ln}det\left[S_{mn}(k_j,k_{j+1})\right]\right\}\right\},`$ (16)
where $`j`$ runs over the $`N_k`$ $`k`$-vectors in the discretized Brillouin zone, $`\mathrm{\Delta }k=2\pi /aN_k`$, and

$`S_{nm}(k_j,k_{j+1})=u_{nk_j}|u_{mk_{j+1}}`$ (17)
is the overlap matrix between states at adjacent points in the reciprocal-space mesh.

### C Position operator for periodic systems

In view of the above discussion, we examine now the action of the position operator in a space of periodic functions. Keeping in mind that we wish to retain the periodicity of the charge density, we seek to arrive at a consistent definition for the action of $`\widehat{x}`$ in that space. This problem has been recently tackled by Resta , who suggested an intrinsically many-body redefinition of $`\widehat{x}`$, in the context of periodic systems. In the spirit of retaining a single-particle picture, here we only offer a heuristic justification for the form of the perturbation term given in Eq. 13. For this, we use the crystal momentum representation (CMR), following the discussion in the paper by Blount.

Let $`f(x)`$ denote a square-integrable function. The full set of zero-field Bloch eigenstates of a periodic Hamiltonian forms a complete basis to expand $`f(x)`$:

$$f(x)=\frac{1}{2\pi }𝑑k\underset{n}{}\psi _{nk}^{(0)}(x)f_n(k)=\frac{1}{2\pi }𝑑ke^{ikx}\underset{n}{}u_{nk}^{(0)}(x)f_n(k).$$
(18)
The action of $`\widehat{x}`$ on $`f(x)`$ is given in the CMR by the expression

$$\widehat{x}f(x)=\frac{1}{2\pi }𝑑ke^{ikx}\underset{n}{}u_{nk}^{(0)}(x)\left[i\frac{f_n(k)}{k}+\underset{n^{}}{}𝒰_{nn^{}}(k)f_n^{}(k)\right].$$
(19)
where $`𝒰_{nn^{}}`$ is defined in Eq. 14.

Blount examined the transformation properties of the two terms appearing in Eq. 19, with respect to the choice of the phases of the Bloch orbitals. Consider the following CMR decomposition of $`\widehat{x}`$:

$$\widehat{x}=i\frac{}{k}+𝒰_{nn^{}}(k)=x_d+𝒰_{nn^{}}(k),$$
(20)
where $`x_d=i\frac{}{k}`$ is diagonal in the band index. He showed that when the Bloch orbitals are multiplied by a phase factor $`e^{i\varphi _n(k)}`$, the term $`x_d`$ transforms as $`x_d^{}=x_d\delta _{nn^{}}\varphi _n(k)/k`$, while a compensatory change occurs in the diagonal term $`𝒰_{nn}`$. So, the two terms in Eq. 20 do not transform separately like operators, while their sum does. In our formulation, we use the second term on the right, $`𝒰_{nn^{}}`$, to define a periodic Hamiltonian for the electric-field problem which, from this discussion, is not by itself an operator in the strict sense.

We show now that $`𝒰_{nn^{}}(k)`$ is translationally invariant, while $`x_d`$ (like $`\widehat{x}`$) is not. Consider a translation $`T_{\mathrm{}}`$ by a lattice vector $`\mathrm{}`$. The commutation relation for $`\widehat{x}`$ and $`T_{\mathrm{}}`$ is written.

$$[\widehat{x},T_{\mathrm{}}]=\mathrm{}T_{\mathrm{}}.$$
(21)
To obtain the commutation relation of the perturbation term in Eq. 13, we expand Eq. (21) in the CMR representation. Let $`g(x)=T_{\mathrm{}}f(x)=f(x\mathrm{})`$. From Eq. (18) and $`u_{nk}^{(0)}(x\mathrm{})=u_{nk}^{(0)}(x)`$ we get

$$g(x)=\frac{1}{2\pi }𝑑ke^{ik(x\mathrm{})}\underset{n}{}u_{nk}^{(0)}(x)f_n(k).$$
(22)
From this expression, it follows that

$$g_n(k)=e^{ik\mathrm{}}f_n(k).$$
(23)
Further, from Eq. (19) we obtain

$$T_{\mathrm{}}\widehat{x}f(x)=\frac{1}{2\pi }𝑑ke^{ik(x\mathrm{})}\underset{n}{}u_{nk}^{(0)}(x)\left[i\frac{f_n(k)}{k}+\underset{n^{}}{}𝒰_{nn^{}}(k)f_n^{}(k)\right];$$
(24)
while from Eqs. (19) and (23) we get

$`\widehat{x}T_{\mathrm{}}f(x)=\widehat{x}g(x)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑ke^{ikx}\underset{n}{}u_{nk}^{(0)}(x)\left[i\frac{g_n(k)}{k}+\underset{n^{}}{}𝒰_{nn^{}}(k)g_n^{}(k)\right]}`$ (25)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑ke^{ik(x\mathrm{})}\underset{n}{}u_{nk}^{(0)}(x)\left[\mathrm{}f_n(k)+i\frac{f_n(k)}{k}+\underset{n^{}}{}𝒰_{nn^{}}(k)f_n^{}(k)\right]};`$ (26)
where, in the last step, we use the result $`i\frac{}{k}g_n(k)=e^{ik\mathrm{}}\left[\mathrm{}f_n(k)+i\frac{}{k}f_n(k)\right]`$. Combining Eqs. (24) and (25) we arrive at the CMR expansion of the commutation relation in Eq. (21). Moreover, the above development immediately shows that

$`[x_d,T_{\mathrm{}}]=\mathrm{}T_{\mathrm{}}`$ (27)
$`[𝒰_{nn^{}}(k),T_{\mathrm{}}]=0.`$ (28)
From the above, we observe that the perturbation term in Eq. 13 is invariant under lattice translations.

### D Wannier-function cutoff in real space

The definition of an energy functional at finite fields requires careful analysis. Because the problem does not have a ground state, a regularization procedure is required for the definition of a numerically-stable functional capturing the physics of the state of the system after the electric field is turned on.

In the NV treatment of the problem in real-space, regularization is achieved with the introduction of truncated Wannier functions, which are constrained to zero beyond a real-space cutoff $`R_c`$ As discussed by NV, in the limit $`R_c\mathrm{}`$ the functional becomes pathological, with the property that a state having an arbitrary value for the polarization can be constructed without changing the value of the energy, when working at fixed polarization, or conversely with the development of a growing (infinitely many in the $`R_c\mathrm{}`$ limit) false local minima when working at fixed electric field.

In order to develop this analysis on a sound mathematical basis, one performs a Legendre transformation , from the $``$-dependent total energy $`E[]`$ to the $`P`$-dependent electric enthalpy $`\stackrel{~}{E}[P]`$:

$$\stackrel{~}{E}[P]=\underset{}{inf}\left\{E[]+aP\right\}.$$
(29)
The total energy was obtained previously \[see Eq.(5)\] thanks to the trial Wannier functions,

$$E[]=\underset{\{w\}}{inf}\left\{E[\{w\};]\right\}=\underset{\{w\}}{inf}\left\{E^{(0)}[\{w\}]aP[\{w\}]\right\},$$
(30)
while a constrained search alternatively gives its Legendre transform,

$$\stackrel{~}{E}[P]=\underset{\{w\}\mathrm{such}\mathrm{that}P[\{w\}]=P}{inf}\left\{E^{(0)}[\{w\}]\right\}.$$
(31)
The zero-electric field total energy functional of the Wannier functions is $`E^{(0)}[\{w\}]=_iw_i\left|H^{(0)}\right|w_i`$.

We aim at understanding the pathologies of $`E[]`$ by examining its expression as the inverse Legendre transform of $`\stackrel{~}{E}[P]`$,

$$E[]=\underset{P}{inf}\left\{\stackrel{~}{E}[P]aP\right\},$$
(32)
for which we need to characterize the minima of $`\stackrel{~}{E}[P]`$, as well as their local behavior.

We consider, for simplicity, the case of a single occupied band. For a given finite value of $`R_c`$, the electric enthalpy is a periodic function of $`P`$, and $`\stackrel{~}{E}[P_0]=E^{(0)}[\{w_0\}]=E_0`$ is the zero-field ground-state energy ($`w_0`$ is the zero-field valence-band Wannier function). For large values of $`R_c`$, it becomes possible to build a set of $`\mathrm{}`$-dependent functions (to be normalized),

$$|w=|w_0+P^{1/2}\mathrm{}^{1/2}|w_0^{cb}(\mathrm{}),$$
(33)
with arbitrary value of $`P`$, where $`w_0`$ is a zero-field valence-band Wannier function centered at the origin, and $`w_0^{cb}(\mathrm{})`$ is an empty conduction-band function centered at the site $`\mathrm{}`$ within the range of $`R_c`$, whose coefficient is on the order of $`\mathrm{}^{1/2}`$. We consider the lattice constant $`a`$ to be the unit of length. The energy for these states is

$$\frac{w\left|H^{(0)}\right|w}{w|w}=\frac{w_0\left|H^{(0)}\right|w_0+P\mathrm{}^1w_0^{cb}(\mathrm{})\left|H^{(0)}\right|w_0^{cb}(\mathrm{})}{1+P\mathrm{}^1}=E_0+P\mathrm{}^1(E_{cb}E_0),$$
(34)
where $`E_{cb}`$ is the expectation value of the energy for the conduction band Wannier function. Due to the exponential decay of Wannier functions for insulators, the value of the polarization is

$$\frac{w|x|w}{w|w}\frac{w_0|x|w_0+P\mathrm{}^1w_0^{cb}(\mathrm{})\left|x\right|w_0^{cb}(\mathrm{})}{1+P\mathrm{}^1}=P_0+P,$$
(35)
since $`w_0^{cb}(\mathrm{})\left|x\right|w_0^{cb}(\mathrm{})=\mathrm{}+P_0^{cb}`$. In the limit $`\mathrm{}\mathrm{}`$, these wavefunctions have an arbitrary value of $`P`$, and an energy infinitesimally close to $`E_0`$. The $`\stackrel{~}{E}[P]`$ curve becomes flat in this limit, and only derivatives in an infinitesimal region around the ground-state solution remain well defined. The development of multiple minima at finite fields corresponds to the same situation, as a growing number of minima with energies that become degenerate in the $`R_c\mathrm{}`$ are associated to states with different values of polarization. No global minimum as a function of polarization can be found. In the next subsection, we analyze the behavior of the energy functional for a model system in reciprocal space. We will show that the same pathology manifests itself in the limit $`\mathrm{\Delta }k0`$, where $`\mathrm{\Delta }k`$ is the discretization of the mesh of $`k`$-points in the Brillouin zone.

### E Reciprocal-space analysis of a model system

For the present analysis, as well as for the application of the perturbation expansions to be developed in the following sections, we chose a one-dimensional (1D) two-site periodic model defined by two parameters, the hoping integral $`t`$, and the on-site term which we choose as $`\mathrm{\Delta }/2`$ and $`\mathrm{\Delta }/2`$, for sites $`1`$ and $`2`$, respectively. The Hamiltonian can be rescaled by $`\mathrm{\Delta }`$ to become a one-parameter ($`\frac{t}{\mathrm{\Delta }}t`$) model, defined as

$`H={\displaystyle \underset{l}{}}\{{\displaystyle \frac{1}{2}}c_{2,l}^{}c_{2,l}{\displaystyle \frac{1}{2}}c_{1,l}^{}c_{1,l}+t[c_{1,l}^{}c_{2,l}+c_{2,l}^{}c_{1,l+1}+h.c.]\},`$ (36)
where $`l`$ runs over unit cells. Whenever we are concerned with the 1D model, we will consider all distances to be rescaled by the unit-cell period (i.e. we set $`a=1`$ in the present section, in Sec. VI, and in Appendix C), such that on each cell, denoted by the integer $`l`$, we have the basis functions $`\varphi _1(l)`$ and $`\varphi _2(l+1/2)`$.

We apply Bloch’s theorem to write the Schrödinger equation for the cell-periodic functions:

$`H_k^{(0)}|u_{nk}^{(0)}=\epsilon _{nk}^{(0)}|u_{nk}^{(0)},`$ (37)
where $`H_k^{(0)}`$ is the zero-field cell-periodic Hamiltonian. In the basis of periodic functions $`\chi _1=_l\varphi _1(l)`$ and $`\chi _2=_l\varphi _2(l+1/2)`$ we have

$`H_k^{(0)}=\left(\begin{array}{cc}{\displaystyle \frac{1}{2}}& 2t\mathrm{cos}\frac{k}{2}\\ 2t\mathrm{cos}\frac{k}{2}& \frac{1}{2}\end{array}\right).`$ (40)
The corresponding secular equation, $`det\left[H_k^{(0)}\epsilon _k^{(0)}1\right]=0`$, is easily solved for the eigenvalues

$`\epsilon _k^{(0)}={}_{}{}^{+}[{\displaystyle \frac{1}{4}}+4t^2\mathrm{cos}^2{\displaystyle \frac{k}{2}}]_{}^{\frac{1}{2}}={}_{}{}^{+}{\displaystyle \frac{1}{2}}[1+A^2\mathrm{cos}^2{\displaystyle \frac{k}{2}}]^{\frac{1}{2}},`$ (41)
where $`A=4t`$. Negative and positive eigenvalues correspond to valence and conduction bands, respectively. Because the Hamiltonian is real, we can use the following parameterization for the corresponding eigenstates

$`|u_{vk}^{(0)}`$ $`=`$ $`\left(\begin{array}{c}\mathrm{cos}\mathrm{\Theta }_k\\ \mathrm{sin}\mathrm{\Theta }_k\end{array}\right)\mathrm{e}^{i\alpha _{vk}}`$ (44)
$`|u_{ck}^{(0)}`$ $`=`$ $`\left(\begin{array}{c}\mathrm{sin}\mathrm{\Theta }_k\\ \mathrm{cos}\mathrm{\Theta }_k\end{array}\right)\mathrm{e}^{i\alpha _{ck}},`$ (47)
where $`\alpha _{vk}`$ and $`\alpha _{ck}`$ are real numbers, with no lack of generality. Coming back to the eigenvalue equation $`H_k^{(0)}|u_{vk}^{(0)}=\epsilon _{vk}^{(0)}|u_{vk}^{(0)}`$, we obtain

$`\mathrm{tan}\mathrm{\Theta }_k={\displaystyle \frac{\epsilon _{vk}^{(0)}+\frac{1}{2}}{2t\mathrm{cos}\frac{k}{2}}}.`$ (48)
Integrating Eq. 41 over the Brillouin zone gives the energy per unit cell:

$`E_0={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}𝑑k\epsilon _{vk}^{(0)}={\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\frac{\pi }{2}}}𝑑y\left[1+A^2\mathrm{cos}^2y\right]^{\frac{1}{2}}.`$ (49)
In order to discuss the pathology of the finite electric-field functional in $`k`$-space, we consider a set of trial cell-periodic functions

$`|u_k=\left(\begin{array}{c}\mathrm{cos}\mathrm{\Theta }_ke^{i\alpha _k}\\ \mathrm{sin}\mathrm{\Theta }_ke^{i\beta _k}\end{array}\right)`$ (52)
for $`k[\pi ,\pi ]`$, where $`\mathrm{\Theta }_k`$, $`\alpha _k`$, and $`\beta _k`$ are real numbers. Imposing the condition $`|u_{k+G}=e^{iGr}|u_k`$ we obtain

$`\mathrm{cos}\mathrm{\Theta }_{k+2\pi }e^{i\alpha _{k+2\pi }}`$ $`=`$ $`\mathrm{cos}\mathrm{\Theta }_ke^{i\alpha _k}`$ (53)
$`\mathrm{sin}\mathrm{\Theta }_{k+2\pi }e^{i\beta _{k+2\pi }}`$ $`=`$ $`\mathrm{sin}\mathrm{\Theta }_ke^{i\beta _k},`$ (54)
which leads to $`\alpha _{2\pi }\alpha _0=N_\alpha \pi `$ and $`\beta _{2\pi }\beta _0=N_\beta \pi `$.

The expectation value of the zero-field Hamiltonian in the set of trial wave-functions gives

$`E^{(0)}[\{u_k\}]={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}𝑑ku_k\left|H_k^{(0)}\right|u_k={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}𝑑k\left[{\displaystyle \frac{1}{2}}\mathrm{cos}\left(2\mathrm{\Theta }_k\right)+2t\mathrm{cos}\left({\displaystyle \frac{k}{2}}\right)\mathrm{sin}\left(2\mathrm{\Theta }_k\right)\mathrm{cos}\gamma _k\right],`$ (55)
where $`\gamma _k=\alpha _k\beta _k`$.

The polarization for the trial state is

$`P[\{u_k\}]={\displaystyle \frac{ie}{\pi }}{\displaystyle _0^{2\pi }}𝑑ku_k\left|{\displaystyle \frac{}{k}}\right|u_k={\displaystyle \frac{e}{2\pi }}\left\{\left[\alpha _k+\beta _k\right]_0^{2\pi }+{\displaystyle _0^{2\pi }}𝑑k\mathrm{cos}\left(2\mathrm{\Theta }_k\right){\displaystyle \frac{\gamma _k}{k}}\right\}.`$ (56)
Minimization of $`E^{(0)}[\{u_k\}]`$ with respect to $`\mathrm{\Theta }_k`$ and $`\gamma _k`$, by setting $`E_k^{(0)}/\mathrm{\Theta }_k=0`$ and $`E_k^{(0)}/\gamma _k=0`$, with $`E_k^{(0)}=u_k\left|H_k^{(0)}\right|u_k`$, leads to

$`\mathrm{tan}(2\mathrm{\Theta }_k)=4t\mathrm{cos}\left({\displaystyle \frac{k}{2}}\right)\mathrm{cos}\gamma _k;`$ (57)
$`\{\begin{array}{c}\mathrm{sin}\gamma _k=0\\ \mathrm{or}\\ 2t\mathrm{cos}\left(\frac{k}{2}\right)\mathrm{sin}\left(2\mathrm{\Theta }_k\right)=0.\end{array}`$ (61)
At $`\mathrm{cos}(k/2)=0`$, the solution of Eqs. 57 and 61 leads to $`\mathrm{sin}(2\mathrm{\Theta }_k)=0`$ and also implies that $`\gamma _k`$ is undefined. At $`\mathrm{cos}(k/2)0`$, a minimum solution is obtained by setting

$`\{\begin{array}{c}\mathrm{sin}\gamma _k=0\gamma _k=N_\gamma \times 2\pi ,\\ \mathrm{and}\\ \mathrm{tan}(2\mathrm{\Theta }_k)=4t\mathrm{cos}\left(\frac{k}{2}\right).\end{array}`$ (65)
The ground-state solution is given by $`\gamma _k=0`$ (i.e., $`\alpha _k=\beta _k`$ as in Eq. 44) for all values of $`k`$, with $`\mathrm{\Theta }_k`$ defined by Eq. 65. Note that a solution where $`\gamma _k`$ jumps by a multiple of $`2\pi `$ at $`k=_{}^+\pi `$ is also consistent with Eqs. 57-65, but not with the restriction that $`u_k`$ be differentiable with respect to $`k`$. Note also that, due to inversion symmetry, the zero-field ground-state polarization must vanish (modulo $`e`$). This is what is obtained from Eq. 56, by setting $`\gamma _k/k=0`$, $`\alpha _k=\beta _k`$, and $`\alpha _{2\pi }\alpha _0=N_\alpha \pi `$.

We consider now a trial wavefunction where $`\mathrm{\Theta }_k`$ is the same as in the ground-state solution, while $`\gamma _k`$ behaves as shown in Fig. 1, where it jumps by a value of $`2\pi `$ over an interval $`\mathrm{\Delta }k`$ centered at an arbitrary value of $`k`$. We show now that in the $`\mathrm{\Delta }k0`$ limit this function can be tailored to give an arbitrary value of the polarization, while its energy differs from the ground-state by an infinitesimal amount, of order $`\mathrm{\Delta }k`$.

The change in polarization for this state, with respect to the ground-state solution, is given by

$`\mathrm{\Delta }P={\displaystyle \frac{e}{2\pi }}{\displaystyle _0^{2\pi }}𝑑k\mathrm{cos}\left(2\mathrm{\Theta }_k\right){\displaystyle \frac{\gamma _k}{k}}e\mathrm{cos}\left(2\mathrm{\Theta }_k\right).`$ (66)
$`\mathrm{\Delta }P`$ in the above equation assumes values between $`e`$ and $`e/(1+16t^2)^{1/2}`$. By adding another kink in the definition of $`\gamma _k`$, where this function changes by $`2\pi `$, we can build a solution having any arbitrary value of $`P`$ in the interval $`e[0,1]`$.

Let us consider the change in energy of the trial state. The function $`\mathrm{cos}\gamma _k`$, as shown in Fig. 2, differs from one over a small interval of the order of $`\mathrm{\Delta }k`$. The change in energy with respect to the ground state is then

$`\mathrm{\Delta }E={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}𝑑k2t\mathrm{cos}\left({\displaystyle \frac{k}{2}}\right)\mathrm{sin}\left(2\mathrm{\Theta }_k\right)\left(\mathrm{cos}\gamma _k1\right){\displaystyle \frac{2t}{\pi }}\mathrm{cos}\left({\displaystyle \frac{k}{2}}\right)\mathrm{sin}\left(2\mathrm{\Theta }_k\right)\times \mathrm{\Delta }k.`$ (67)
So, for the trial state $`\mathrm{\Delta }E0`$ when $`\mathrm{\Delta }k0`$. The $`E^{(0)}(P)`$ curve becomes flat in this limit. This is the same pathology as the one discussed by NV in the real-space case. Discretization of the $`k`$-space mesh in the Brillouin zone is thus essential for the numerical stability of the energy functional.

Now, we show that for the discretized version of the formulation, a change in $`P`$ implies a finite change in the energy. The discretized polarization is written

$`P\left[u_{k_j}\right]`$ $`=`$ $`{\displaystyle \frac{e}{\pi }}{\displaystyle \underset{k=1}{\overset{N_k}{}}}\mathrm{Im}\left[\mathrm{ln}u_{k_j}|u_{k_{j+1}}\right]`$ (68)
$`=`$ $`{\displaystyle \frac{e}{\pi }}\left\{\mathrm{\Delta }N_\beta +\mathrm{Im}{\displaystyle \underset{k=1}{\overset{N_k}{}}}\mathrm{ln}\left[\mathrm{cos}\mathrm{\Theta }_{k_j}\mathrm{cos}\mathrm{\Theta }_{k_{j+1}}e^{i\left(\gamma _{k_{j+1}}\gamma _{k_j}\right)}+\mathrm{sin}\mathrm{\Theta }_{k_j}\mathrm{sin}\mathrm{\Theta }_{k_{j+1}}\right]\right\};`$ (69)
with the energy given by

$`E^{(0)}\left[u_{k_j}\right]={\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}u_{k_j}\left|H_{k_j}^{(0)}\right|u_{k_{j+1}}={\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\left[{\displaystyle \frac{1}{2}}\mathrm{cos}\left(2\mathrm{\Theta }_{k_j}\right)+2t\mathrm{cos}\left({\displaystyle \frac{k_j}{2}}\right)\mathrm{sin}\left(2\mathrm{\Theta }_{k_j}\right)\mathrm{cos}\gamma _{k_j}\right].`$ (70)
Again we consider the ground-state solution $`\gamma _{k_j}=0`$, with $`\mathrm{\Theta }_{k_j}`$ given by Eq. 65. An arbitrary change in polarization can be introduced by setting $`\gamma _{k_j}0`$ at a given $`k_j`$, while keeping the values of $`\mathrm{\Theta }`$ and $`\gamma `$ at all the other $`k`$-points unchanged. In this case, it can be immediately seen that a finite change $`\mathrm{\Delta }E_0=(4t/N_k)\mathrm{cos}\left(\frac{k_j}{2}\right)\mathrm{sin}\left(2\mathrm{\Theta }_{k_j}\right)(\mathrm{cos}\gamma _{k_j}1)`$ is introduced in the discretized energy.

### F Summary

The theoretical treatment of a periodic insulator placed in an homogeneous electric field is plagued by severe conceptual difficulties: (1) the potential associated with an electric field is non-periodic and unbounded; (2) for that reason the spectrum of electronic states changes non-analytically upon the application of a homogeneous electric field; (3) the quantity conjugated to the electric field, namely the polarization, cannot be computed as the expectation value of the position (or any other) operator; (4) local minima of the energy functional can be defined only in an infinitesimally small region as a function of the polarization, the energy functional being perfectly flat otherwise.

In order to address problems 1 and 2, following Nenciu, we restrict ourselves to periodic-polarized-insulating states, of which the lowest in energy is a long-lived resonance of the unrestricted system. Keeping this restriction in mind, we show that the position operator can be decomposed, in the crystal momentum representation, into a non-periodic part and a periodic part. The latter can be introduced in an ansatz Hamiltonian acting on the periodic part of the Bloch functions, from which the Berry phase formulation of the polarization is recovered, solving problem 3 as well.

We are aware that this line of thought does not yet justify rigorously the use of this Hamiltonian: a more careful derivation, in the spirit of the mathematical work of Nenciu, would be needed. However, this rather simple Hamiltonian allows to recover all the previously known lowest-order expressions for the polarization and its derivatives, and to derive other low-order expressions as well as generic expressions to all orders, as we shall see in the coming sections.

Problem 4 is solved by introducing a regularization procedure in reciprocal space, similar in spirit to the real-space cutoff radius introduced by NV. For the regularized energy functional, the local minima have a finite basin of attraction as a function of the polarization.

## III Perturbation theory applied to a variational total-energy functional

In view of the application of perturbation theory to Eq. 10, we summarize now the variational formulation of DFPT, as presented in Ref. . We consider the formalism at its non-self-consistent level, without including the Hartree and exchange-correlation terms of the perturbative expressions.

One considers a perturbative expansion of a variational principle applied to the electronic total-energy functional. In terms of the small parameter $`\lambda `$ associated with the perturbation, the perturbation series reads

$`𝒪(\lambda )`$ $`=`$ $`𝒪^{(0)}+\lambda 𝒪^{(1)}+\lambda ^2𝒪^{(2)}+\lambda ^3𝒪^{(3)}+\mathrm{},`$ (71)
$`𝒪^{(n)}`$ $`=`$ $`{\displaystyle \frac{1}{n!}}{\displaystyle \frac{d^n𝒪(\lambda )}{d\lambda ^n}}|_{\lambda =0}`$ (72)
for a generic observable $`𝒪`$. The system Hamiltonian is $`H=K+v_{ext}`$, and the total-energy functional is

$`E={\displaystyle \underset{\alpha =1}{\overset{N}{}}}\phi _\alpha \left|\left(K+v_{ext}\right)\right|\phi _\alpha ,`$ (73)
where $`K`$ and $`v_{ext}`$ are the kinetic-energy and external-potential operators. The total-energy functional is to be minimized under the orthonormality constraints for the occupied wavefunctions

$`\phi _\alpha |\phi _\beta =\delta _{\alpha \beta }.`$ (74)
Using the Lagrange-multiplier method, the functional

$`F={\displaystyle \underset{\alpha =1}{\overset{N}{}}}\phi _\alpha \left|\left(T+v_{ext}\right)\right|\phi _\alpha {\displaystyle \underset{\alpha ,\beta =1}{\overset{N}{}}}\mathrm{\Lambda }_{\beta \alpha }\left[\phi _\alpha |\phi _\beta \delta _{\alpha \beta }\right]`$ (75)
is minimized with respect to the wavefunctions. The minimum condition, $`\delta F/\delta \phi _\alpha ^{}=0`$, leads to the Euler-Lagrange equation

$`H|\phi _\alpha ={\displaystyle \underset{\beta =1}{\overset{N}{}}}\mathrm{\Lambda }_{\beta \alpha }|\phi _\beta .`$ (76)
Eq. 76 represents a set of generalized eigenvalue equations which assume the form of the usual eigenvalue equations when the so-called diagonal gauge is chosen to fix the phase arbitrariness of the wavefunctions.  Here, we keep the generalized form, as needed for the choice of gauge to be used in our theory. An expression for the Lagrange-multiplier matrix is obtained by multiplying Eq. 76 by an occupied wavefunction, leading to

$`\mathrm{\Lambda }_{\beta \alpha }=\phi _\beta \left|H\right|\phi _\alpha .`$ (77)
We consider now the perturbation expansion of Eqs. 74-77. The orthonormalization condition becomes

$`{\displaystyle \underset{j=0}{\overset{i}{}}}\phi _\alpha ^{(j)}|\phi _\beta ^{(ij)}=0\mathrm{for}i1.`$ (78)
The expansion of Eq. 76 gives the generalized Sternheimer equation

$`{\displaystyle \underset{j=0}{\overset{i}{}}}H^{(j)}|\phi _\alpha ^{(ij)}={\displaystyle \underset{j=0}{\overset{i}{}}}{\displaystyle \underset{\beta =1}{\overset{N}{}}}\mathrm{\Lambda }_{\beta \alpha }^{(j)}|\phi _\beta ^{(ij)},`$ (79)
where $`H^{(i)}=T^{(i)}+v_{ext}^{(i)}`$ is the $`ith`$-order term in the expansion of the Hamiltonian.

The expansion of the Lagrange-multiplier matrix is given by

$`\mathrm{\Lambda }_{\beta \alpha }^{(i)}={\displaystyle \underset{j=0}{\overset{i}{}}}{\displaystyle \underset{k=0}{\overset{i}{}}}\phi _\beta ^{(j)}\left|H^{(ijk)}\right|\phi _\alpha ^{(k)}.`$ (80)
Finally, a generic term in the perturbative expansion of the total-energy functional in Eq. 75 is written

$`E^{(i)}`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{N}{}}}{\displaystyle \underset{l=0}{\overset{j}{}}}{\displaystyle \underset{k=0}{\overset{i}{}}}{\displaystyle \underset{l^{}=0}{\overset{j}{}}}\delta (ilkl^{})\phi _\alpha ^{(l)}\left|H^{(k)}\right|\phi _\alpha ^{(l^{})}`$ (81)
$``$ $`{\displaystyle \underset{\alpha ,\beta =1}{\overset{N}{}}}{\displaystyle \underset{l=0}{\overset{j}{}}}{\displaystyle \underset{k=0}{\overset{ij1}{}}}{\displaystyle \underset{l^{}=0}{\overset{j}{}}}\delta (ilkl^{})\mathrm{\Lambda }_{\beta \alpha }^{(k)}\phi _\alpha ^{(l)}|\phi _\beta ^{(l^{})},`$ (82)
where $`i=2j`$ or $`i=2j+1`$. We remark that only wavefunctions derivatives up to order $`\lambda ^j`$ appear in the $`i`$th-order term of the energy, as a result of the $`2n+1`$-theorem. Moreover, a minimum principle holds for $`E^{(2j)}`$ with respect to the $`j`$th-order variations of the wave functions, i.e., $`\delta E^{(2j)}/\delta \phi _\alpha ^{(j)}=0`$.

A particularly useful result derived in Ref.  is a set of non-variational expressions for the second-order derivative of the energy. In the present work, the Hamiltonian is of first-order in the perturbation ($`v_{ext}^{(i)}=0\mathrm{for}i2`$), in which case the non-variational expressions are given by

$`E^{(2)}`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{N}{}}}\phi _\alpha ^{(1)}\left|v_{ext}^{(1)}\right|\phi _\alpha ^{(0)}={\displaystyle \underset{\alpha =1}{\overset{N}{}}}\phi _\alpha ^{(0)}\left|v_{ext}^{(1)}\right|\phi _\alpha ^{(1)}`$ (83)
$`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{N}{}}}{\displaystyle \frac{1}{2}}\phi _\alpha ^{(1)}\left|v_{ext}^{(1)}\right|\phi _\alpha ^{(0)}+{\displaystyle \frac{1}{2}}\phi _\alpha ^{(0)}\left|v_{ext}^{(1)}\right|\phi _\alpha ^{(1)}.`$ (84)
The zeroth-order wave functions are chosen to obey the unperturbed eigenvalue equation $`H^{(0)}|\phi _\alpha ^{(0)}=\epsilon _\alpha ^{(0)}|\phi _\alpha ^{(0)}`$. From Eq. 80, the zeroth-order Lagrange-multiplier matrix is given by

$`\mathrm{\Lambda }_{\beta \alpha }^{(0)}=\delta _{\beta \alpha }\epsilon _\alpha ^{(0)}.`$ (85)
In the present work, we use the so-called “parallel-transport” gauge, as discussed in Ref. . In this gauge, the following condition is imposed on the derivatives of the wave functions

$`\phi _\alpha ^{(0)}|\phi _\beta ^{(i)}\phi _\alpha ^{(i)}|\phi _\beta ^{(0)}=0;`$ (86)
which allows us to rewrite the expansion of the orthonormalization condition as

$`\phi _\alpha ^{(0)}|\phi _\beta ^{(i)}=\{\begin{array}{ccc}{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{i1}{}}}\phi _\alpha ^{(j)}|\phi _\beta ^{(ij)}\hfill & \mathrm{for}& i>1;\\ 0\hfill & \mathrm{for}& i=1.\end{array}`$ (89)
## IV Perturbation theory applied to the continuous form

### A Perturbation expansion and proof of gauge invariance

Following the discussion in Sec. II, we can develop a perturbation expansion for the electric-field problem. In this section, we discuss the continuous form of the theory. The cell-periodic Hamiltonian, including the perturbation term, is given in Eq. 13. We apply the machinery of the variational DFPT to this Hamiltonian, by postulating that the expression

$`E[\left\{u_{nk}\right\};]={\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k\left[{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}\left|H_k^{(0)}+ie{\displaystyle \frac{}{k}}\right|u_{nk}\right]`$ (90)
is to be minimized with respect to the $`\left\{u_{nk}\right\}`$, under the constraints

$$u_{mk}|u_{nk}=\delta _{mn}.$$
(91)
A local minimum will exist for the functional in Eq. 90 provided that a discretization of the $`k`$-space integrals is performed. The continuum formulation which is considered in this section is valid only at infinitesimal fields.

Applying Eq. 75, we write

$$F[\left\{u_{nk}\right\};]=\frac{a}{\pi }_0^{\frac{2\pi }{a}}𝑑k\left[\underset{n=1}{\overset{N}{}}u_{nk}\left|H_k^{(0)}+ie\frac{}{k}\right|u_{nk}\underset{m,n=1}{\overset{N}{}}\left\{u_{nk}|u_{mk}\delta _{nm}\right\}\mathrm{\Lambda }_{mn}(k)\right].$$
(92)
The unconstrained minimization of this functional is obtained by setting $`\delta F\left[\left\{u_{nk}\right\}\right]/\delta u_{nk}`$ = 0, leading to the corresponding Euler-Lagrange equation

$`\left(H_k^{(0)}+ie{\displaystyle \frac{}{k}}\right)|u_{nk}{\displaystyle \underset{m=1}{\overset{N}{}}}\mathrm{\Lambda }_{mn}(k)|u_{mk}=0.`$ (93)
Next, we consider separately the perturbation expansions of Eqs. 92 and 93. In both cases, we will demonstrate explicitly that the general expansion term transforms properly under a general unitary transformation of the occupied orbitals.

#### 1 Lagrange multipliers and orthonormalization constraints

In the present case, Eq. 89 for the orthonormalization constraints reads

$`u_{mk}^{(0)}|u_{nk}^{(i)}=\{\begin{array}{cc}{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{i1}{}}}u_{mk}^{(j)}|u_{nk}^{(ij)},\hfill & i>1;\\ 0,\hfill & i=1;\end{array}`$ (96)
giving the occupied-subspace projection of $`u_{nk}^{(i)}`$ in terms of the lower-order solutions for the periodic functions.

Since $`H_k^{(i)}0`$ for all $`i2`$, and $`H_k^{(1)}=ie\frac{}{k}`$, the expansion of the Lagrange multipliers becomes

$`\mathrm{\Lambda }_{mn}^{(i)}(k)={\displaystyle \underset{j=0}{\overset{i}{}}}u_{mk}^{(j)}\left|H_k^{(0)}\right|u_{nk}^{(ij)}+{\displaystyle \underset{j=0}{\overset{i1}{}}}u_{mk}^{(j)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(ij1)}.`$ (97)
In the following development, we will make explicit use of the expressions for $`\mathrm{\Lambda }_{mn}^{(0)}(k)`$, $`\mathrm{\Lambda }_{mn}^{(1)}(k)`$, and $`\mathrm{\Lambda }_{mn}^{(2)}(k)`$. From Eq. 85, $`\mathrm{\Lambda }_{mn}^{(0)}(k)=\epsilon _{nk}^{(0)}\delta _{mn}`$. Since $`u_{mk}^{(0)}|u_{nk}^{(1)}=0`$ from Eq. 96, and $`H_k^{(0)}|u_{nk}^{(0)}=\epsilon _{nk}^{(0)}|u_{nk}^{(0)}`$, $`\mathrm{\Lambda }_{mn}^{(1)}(k)`$ is given by

$`\mathrm{\Lambda }_{mn}^{(1)}(k)=u_{mk}^{(0)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)}.`$ (98)
The second-order term reads

$`\mathrm{\Lambda }_{mn}^{(2)}(k)`$ $`=`$ $`u_{mk}^{(0)}\left|H_k^{(0)}\right|u_{nk}^{(2)}+u_{mk}^{(2)}\left|H_k^{(0)}\right|u_{nk}^{(0)}+u_{mk}^{(1)}\left|H_k^{(0)}\right|u_{nk}^{(1)}`$ (100)
$`+u_{mk}^{(0)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(1)}+u_{mk}^{(1)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)}.`$
#### 2 Energy

The perturbation expansion for the energy is obtained from Eq. 81. We analyze even and odd terms separately. For the even-order terms we write

$`E^{(2i)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}dk[{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(i)}\left|H_k^{(0)}\right|u_{nk}^{(i)}+u_{nk}^{(i1)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(i)}+u_{nk}^{(i)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(i1)}`$ (102)
$`{\displaystyle \underset{m,n=1}{\overset{N}{}}}{\displaystyle \underset{j,j^{}=1}{\overset{i}{}}}{\displaystyle \underset{l=0}{\overset{i1}{}}}\delta (2ijj^{}l)\mathrm{\Lambda }_{mn}^{(l)}(k)u_{nk}^{(j)}|u_{mk}^{(j^{})}],`$
while the odd terms are given by

$$E^{(2i+1)}=\frac{a}{\pi }_0^{\frac{2\pi }{a}}𝑑k\left[\underset{n=1}{\overset{N}{}}u_{nk}^{(i)}\left|ie\frac{}{k}\right|u_{nk}^{(i)}\underset{n,m=1}{\overset{N}{}}\underset{j,j^{},l=1}{\overset{i}{}}\delta (2i+1jj^{}l)\mathrm{\Lambda }_{mn}^{(l)}(k)u_{nk}^{(j)}|u_{mk}^{(j^{})}\right].$$
(103)
An important aspect concerns the invariance of these expressions with respect to the choice of phases of the Bloch orbitals. More generally, we must consider unitary transformations that keep the subspace of occupied states invariant. We show in Appendix A that Eqs. 102 and 103 can be rewritten in such a way as to display the required gauge-invariance property explicitly. The lower-order derivatives are usually of more practical interest, and for that reason the invariant form of the energy terms up to fourth order are written explicitly here, along with the general expansion term.

The second-order energy derivative is obtained by setting $`i=1`$ in Eq. 102. After some manipulation, this quantity can be written in the following form:

$`E^{(2)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}dk{\displaystyle \underset{n=1}{\overset{N}{}}}[u_{nk}^{(1)}\left|(H_k^{(0)}\epsilon _{nk}^{(0)})\right|u_{nk}^{(1)}+u_{nk}^{(1)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(0)}u_{mk}^{(0)}\left|\right)\right|u_{nk}^{(0)}`$ (104)
$`u_{nk}^{(0)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(0)}u_{mk}^{(0)}\left|\right)\right|u_{nk}^{(1)}].`$ (105)
As expected, our formula for $`E^{(2)}`$ is identical to the linear-response expression.

The non-variational expression for $`E^{(2)}`$ (see Eq. 83) is given by

$`E^{(2)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}\left|\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)\right|u_{nk}^{(1)}`$ (106)
$`=`$ $`{\displaystyle \frac{a}{2\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(0)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)}`$ (107)
$`u_{nk}^{(0)}|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(0)}u_{mk}^{(0)}|\right)|u_{nk}^{(1)}.`$ (108)
The fourth-order energy term is written

$`E^{(4)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}dk{\displaystyle \underset{n=1}{\overset{N}{}}}[u_{nk}^{(2)}\left|(H_k^{(0)}\epsilon _{nk}^{(0)})\right|u_{nk}^{(2)}+u_{nk}^{(2)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(1)}u_{mk}^{(0)}\left|\right)\right|u_{nk}^{(0)}`$ (109)
$`u_{nk}^{(0)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(0)}u_{mk}^{(1)}\left|\right)\right|u_{nk}^{(2)}].`$ (110)
It is worth pointing out the simplicity of the expression for $`E^{(4)}`$, which mirrors that of $`E^{(2)}`$ almost exactly.

The general even-order energy term for $`i>2`$ is written

$`E^{(2i)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}dk{\displaystyle \underset{n=1}{\overset{N}{}}}\{u_{nk}^{(i)}\left|(H_k^{(0)}\epsilon _{nk}^{(0)})\right|u_{nk}^{(i)}+u_{nk}^{(i)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(i1)}u_{mk}^{(0)}\left|\right)\right|u_{nk}^{(0)}`$ (111)
$`u_{nk}^{(0)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(0)}u_{mk}^{(i1)}\left|\right)\right|u_{nk}^{(i)}+{\displaystyle \underset{j,j^{}=1}{\overset{i}{}}}{\displaystyle \underset{l=2}{\overset{i1}{}}}\delta (2ijj^{}l)\times `$ (112)
$`[{\displaystyle \underset{l^{}=0}{\overset{l}{}}}{\displaystyle \underset{m=1}{\overset{N}{}}}u_{mk}^{(l^{})}\left|H_k^{(0)}\right|u_{nk}^{(ll^{})}u_{nk}^{(j)}|u_{mk}^{(j^{})}{\displaystyle \underset{l^{}=0}{\overset{l1}{}}}u_{nk}^{(j)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(j^{})}u_{mk}^{(l^{})})\right|u_{nk}^{(ll^{}1)}]\}.`$ (113)
In the above and the following expressions, we use the notation $`\left(ie\frac{}{k}|u_{mk}u_{mk}|\right)`$ to indicate that $`/k`$ acts only on the quantities embraced in parenthesis. In order to demonstrate that the energy derivatives fulfill the gauge-invariance requirement, we consider a general gauge transformation  among the occupied states at each k-point:

$`|\stackrel{~}{u}_{nk}={\displaystyle \underset{m}{}}U_{mn}(k)|u_{nk},`$ (115)
where $`U`$ is an unitary transformation, i.e., $`UU^{}=1`$. It follows immediately that

$`{\displaystyle \frac{}{k}}\left({\displaystyle \underset{n}{}}|\stackrel{~}{u}_{nk}^{(i)}\stackrel{~}{u}_{nk}^{(j)}|\right)={\displaystyle \frac{}{k}}\left({\displaystyle \underset{lm}{}}{\displaystyle \underset{n}{}}U_{nl}(k)U_{nm}^{}(k)|u_{lk}^{(i)}u_{mk}^{(j)}|\right)={\displaystyle \frac{}{k}}\left({\displaystyle \underset{n}{}}|u_{nk}^{(i)}u_{nk}^{(j)}|\right).`$ (116)
Since $`\frac{}{k}`$ acts only on gauge-invariant quantities, Eqs. 104-111 are themselves gauge-invariant. The same argument holds for the odd-order derivatives we derive below.

Turning now to the odd-order derivatives of the energy, we set $`i=1`$ in Eq. 103 to write the third-order term as

$`E^{(3)}={\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(1)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)}.`$ (117)
This expression for $`E^{(3)}`$ is identical to the one previously derived by Dal Corso and Mauri.

The general odd-order term for $`i>1`$ is written

$`E^{(2i+1)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}dk{\displaystyle \underset{n=1}{\overset{N}{}}}\{u_{nk}^{(i)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(i)}u_{mk}^{(0)}\left|\right)\right|u_{nk}^{(0)}`$ (118)
$`+{\displaystyle \underset{j,j^{}=1}{\overset{i}{}}}{\displaystyle \underset{l=2}{\overset{i}{}}}\delta (2i+1jj^{}l)[{\displaystyle \underset{l^{}=0}{\overset{l}{}}}{\displaystyle \underset{m=1}{\overset{N}{}}}u_{mk}^{(l^{})}\left|H_k^{(0)}\right|u_{nk}^{(ll^{})}u_{nk}^{(j)}|u_{mk}^{(j^{})}`$ (119)
$`{\displaystyle \underset{l^{}=0}{\overset{l1}{}}}u_{nk}^{(j)}\left|\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(j^{})}u_{mk}^{(l^{})})\right|u_{nk}^{(ll^{}1)}]\}.`$ (120)
In appendix A, we demonstrate how Eqs. 104-111, 117, and 118 are obtained from Eqs. 102 and 103.

#### 3 Sternheimer equation

The projection of the wave functions on the subspace of occupied unperturbed states is given by Eq. 96. The projection onto the subspace of unoccupied states is given by the projection of the Sternheimer equation in that subspace. The perturbation series for the Sternheimer equation can be obtained either by expanding Eq. 93, or more directly from Eq. 102 above, by setting $`\delta E^{(2i)}/\delta u_{nk}^{(i)}=0`$. The general expansion term is written

$`P_{ck}\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)P_{ck}|u_{nk}^{(i)}=ieP_{ck}{\displaystyle \frac{}{k}}|u_{nk}^{(i1)}+{\displaystyle \underset{m=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{i1}{}}}\mathrm{\Lambda }_{mn}^{(j)}(k)P_{ck}|u_{mk}^{(ij)}.`$ (121)
where $`P_{ck}`$ is the projector onto the subspace of unoccupied unperturbed states. This equation can be solved for $`u_{nk}^{(i)}`$, once the lower-order derivatives of $`u_{nk}`$ and $`\mathrm{\Lambda }_{mn}(k)`$ have been obtained.

Using the invariant form for the even terms of the energy, Eqs. 104-111, we set $`\delta E^{(2i)}/\delta u_{nk}^{(i)}=0`$ to obtain the explicitly invariant form of the Sternheimer equation. For the $`i=1`$ and $`i=2`$ terms we obtain

$`P_{ck}\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)P_{ck}|u_{nk}^{(1)}=P_{ck}\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(0)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)},`$ (122)
(123)
$`P_{ck}\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)P_{ck}|u_{nk}^{(2)}=P_{ck}\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(1)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)}.`$ (124)
For the complete specification of $`u_{nk}^{(1)}`$ and $`u_{nk}^{(2)}`$, the projections onto the unperturbed occupied subspace are obtained from Eq. 96:

$`P_{vk}|u_{nk}^{(1)}`$ $`=`$ $`0,`$ (125)
$`P_{vk}|u_{nk}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(0)}u_{mk}^{(1)}|u_{nk}^{(1)};`$ (126)
where $`P_{vk}`$ is the projection operator for the occupied states.

The higher-order terms for the Sternheimer equation are given by

$`P_{ck}\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)P_{ck}|u_{nk}^{(i)}=`$ $``$ $`P_{ck}\left[{\displaystyle \underset{j=1}{\overset{i1}{}}}{\displaystyle \underset{l=0}{\overset{j1}{}}}\left(ie{\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}\right|u_{mk}^{(ij)}u_{mk}^{(l)}\left|\right)\right|u_{nk}^{(jl1)}`$ (127)
$``$ $`{\displaystyle \underset{j=2}{\overset{i1}{}}}{\displaystyle \underset{l=0}{\overset{j}{}}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(ij)}u_{mk}^{(l)}\left|H_k^{(0)}\right|u_{nk}^{(jl)}],i>2.`$ (128)
The valence-band component of $`|u_{nk}^{(i)}`$ is given by

$`P_{vk}|u_{nk}^{(i)}={\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{i1}{}}}|u_{nk}^{(0)}u_{mk}^{(j)}|u_{nk}^{(ij)}.`$ (129)
As in the case of the energy terms, in Eqs. 122, 124, and 127 the derivative $`/k`$ acts only on gauge-invariant quantities. This completes our development of the perturbation expansion, and the proof of gauge invariance of the continuous formulation.

### B Discretized form of lower-order expressions

We examine now the discretized form of the lower-order terms for the energy and the Sternheimer equation. In practical calculations, it is mandatory to use a discrete set of $`k`$-points to evaluate the Brillouin-zone integrals. However, when the focus is on $`E^{(2)}`$, the discretization of the $`\frac{}{k}`$ operation can be avoided, as the projection on the conduction bands of the derivative of the wavefunctions versus $`k`$ can be computed from a Sternheimer equation . This has the disadvantage to add a significant coding and computational step in the whole procedure.

We choose the following symmetric finite-difference expansion for the derivatives with respect to $`k`$:

$`{\displaystyle \frac{}{k}}|u_{nk}u_{nk}|{\displaystyle \frac{1}{2\mathrm{\Delta }k}}\left(|u_{nk_{j+1}}u_{nk_{j+1}}||u_{nk_{j1}}u_{nk_{j1}}|\right),`$ (130)
where $`\mathrm{\Delta }k=k_{j+1}k_j=(2\pi /aN_k)`$. Clearly, this expression retains the gauge invariance of the continuous form. Next, Eq. 130 is used in the derivation of explicit discretized expressions for the energy derivatives up to the fourth-order, and for the Sternheimer equation up to the second order.

#### 1 Energy

From Eqs. 104 and 130 we obtain the discretized formula for $`E^{(2)}`$:

$`E^{(2)}`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\{{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk_j}^{(1)}\left|(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)})\right|u_{nk_j}^{(1)}`$ (133)
$`+{\displaystyle \frac{ie}{2\mathrm{\Delta }k}}{\displaystyle \underset{n,m=1}{\overset{N}{}}}[u_{nk_j}^{(1)}\left|\left(\right|u_{mk_{j+1}}^{(0)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(0)}\left|\right)\right|u_{nk_j}^{(0)}`$
$`u_{nk_j}^{(0)}\left|\left(\right|u_{mk_{j+1}}^{(0)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(1)}\left|\right)\right|u_{nk_j}^{(0)}]\}.`$
The non-variational expression for $`E^{(2)}`$ is written

$`E^{(2)}`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk_j}^{(1)}\left|\left(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)}\right)\right|u_{nk_j}^{(1)}`$ (134)
$`=`$ $`{\displaystyle \frac{ie}{2\mathrm{\Delta }k}}{\displaystyle \underset{n,m=1}{\overset{N}{}}}[u_{nk_j}^{(1)}\left|\left(\right|u_{mk_{j+1}}^{(0)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(0)}\left|\right)\right|u_{nk_j}^{(0)}`$ (136)
$`u_{nk_j}^{(0)}\left|\left(\right|u_{mk_{j+1}}^{(0)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(1)}\left|\right)\right|u_{nk_j}^{(0)}].`$
The discretized versions of Eqs. 117 and 109 are

$`E^{(3)}`$ $`=`$ $`{\displaystyle \frac{ie}{N_k\mathrm{\Delta }k}}{\displaystyle \underset{n,m=1}{\overset{N}{}}}u_{nk_j}^{(1)}|\left(|u_{mk_{j+1}}^{(1)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(1)}u_{mk_{j1}}^{(0)}|\right)|u_{nk_j}^{(0)}.`$ (137)
$`E^{(4)}`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\{{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk_j}^{(2)}\left|(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)})\right|u_{nk_j}^{(2)}`$ (141)
$`+{\displaystyle \frac{ie}{2\mathrm{\Delta }k}}{\displaystyle \underset{n,m=1}{\overset{N}{}}}[u_{nk_j}^{(2)}\left|\left(\right|u_{mk_{j+1}}^{(1)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(1)}u_{mk_{j1}}^{(0)}\left|\right)\right|u_{nk_j}^{(0)}`$
$`u_{nk_j}^{(0)}\left|\left(\right|u_{mk_{j+1}}^{(0)}u_{mk_{j+1}}^{(1)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(1)}\left|\right)\right|u_{nk_j}^{(2)}]\}.`$
#### 2 Sternheimer equation

The discretized expressions for the $`i=1`$ and $`i=2`$ terms of the Sternheimer equation are given by

$`P_{ck_j}\left(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)}\right)P_{ck_j}|u_{nk_j}^{(1)}={\displaystyle \frac{ie}{2\mathrm{\Delta }k}}P_{ck_j}{\displaystyle \underset{m=1}{\overset{N}{}}}\left(|u_{mk_{j+1}}^{(0)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(0)}|\right)|u_{nk_j}^{(0)},`$ (142)
(143)
$`P_{ck_j}\left(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)}\right)P_{ck_j}|u_{nk_j}^{(2)}={\displaystyle \frac{ie}{2\mathrm{\Delta }k}}P_{ck_j}{\displaystyle \underset{m=1}{\overset{N}{}}}\left(|u_{mk_{j+1}}^{(1)}u_{mk_{j+1}}^{(0)}||u_{mk_{j1}}^{(0)}u_{mk_{j1}}^{(1)}|\right)|u_{nk_j}^{(0)}.`$ (144)
The required gauge invariance of the expansion terms is preserved in the discretized expressions. Note that the solutions at a given k-point are now coupled to the first-neighbor k-points in the reciprocal-space grid.

## V Perturbation theory applied to the discretized polarization

### A General perturbation expansion

Let us consider now the perturbation treatment of the problem on the basis of the energy functional given in Eq. 16, where the polarization is written in a discretized form. This formulation can be viewed as the reciprocal-space analog of the NV  real-space functional. In this approach the gauge-invariance of the energy is guaranteed by the fact that Eq. 3 is itself gauge invariant.

We seek a minimum for Eq. 16 with respect to the occupied orbitals $`\left\{u_{nk_j}\right\}`$, under the constraints $`u_{nk_j}|u_{mk_j}=\delta _{mn}`$. Lagrange multipliers are introduced to write the unconstrained functional

$`[\left\{u_{nk_j}\right\};]`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}[{\displaystyle \underset{n,m=1}{\overset{N}{}}}u_{mk_j}\left|H_{k_j}^{(0)}\right|u_{nk_j}\delta _{mn}(u_{nk_j}|u_{mk_j}\delta _{mn})\mathrm{\Lambda }_{mn}(k_j)`$ (145)
$``$ $`\left({\displaystyle \frac{e}{\mathrm{\Delta }k}}\right)\mathrm{Im}\{\mathrm{ln}det\left[S_{nm}(k_j,k_{j+1})\right]\}].`$ (146)
In Appendix B, we prove the following result:

$`{\displaystyle \frac{\delta _j\mathrm{Im}\left\{\mathrm{ln}det\left[S_{m^{}m}(k_j,k_{j+1})\right]\right\}}{\delta u_{nk_j}^{}}}={\displaystyle \frac{i}{2}}{\displaystyle \underset{m=1}{\overset{N}{}}}\left[|u_{mk_{j+1}}S_{mn}^1(k_j,k_{j+1})|u_{mk_{j1}}S_{mn}^1(k_j,k_{j1})\right],`$ (147)
which can be used to derive the Euler-Lagrange equation from Eq. 145, as follows:

$`{\displaystyle \frac{\delta }{\delta u_{nk_j}^{}}}`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}\left\{H_{k_j}^{(0)}\right|u_{nk_j}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk_j}\mathrm{\Lambda }_{mn}(k_j)`$ (148)
$`+`$ $`\left({\displaystyle \frac{ie}{2\mathrm{\Delta }k}}\right){\displaystyle \underset{m=1}{\overset{N}{}}}\left[\right|u_{mk_{j+1}}S_{mn}^1(k_j,k_{j+1})|u_{mk_{j1}}S_{mn}^1(k_j,k_{j1})]\}=0.`$ (149)
Below, we consider the perturbation expansions of the Lagrange multipliers, the energy, and the Sternheimer equation. The expansion of the orthonormalization condition was already developed in the previous section, and remains unaltered in the present case.

#### 1 Lagrange multipliers

We multiply Eq. 148 on the left by $`u_{mk_j}^{}`$ to write the Lagrange multipliers

$`\mathrm{\Lambda }_{mn}(k_j)=u_{mk_j}\left|H_{k_j}^{(0)}\right|u_{nk_j}.`$ (150)
It can be readily seen that the terms involving the overlap matrix cancel out, since

$`{\displaystyle \underset{l}{}}\left[u_{mk_j}|u_{lk_{j+1}}S_{ln}^1(k_j,k_{j+1})u_{mk_j}|u_{lk_{j1}}S_{ln}^1(k_j,k_{j1})\right]`$ (151)
$`={\displaystyle \underset{l}{}}\left[S_{ml}(k_j,k_{j+1})S_{ln}^1(k_j,k_{j+1})S_{ml}(k_j,k_{j1})S_{ln}^1(k_j,k_{j1})\right]=0.`$ (152)
The perturbation expansion of Eq. 150 takes the simple form

$`\mathrm{\Lambda }_{mn}^{(i)}(k_j)={\displaystyle \underset{j=0}{\overset{i}{}}}u_{mk_j}^{(ij)}\left|H_{k_j}^{(0)}\right|u_{nk_j}^{(j)}.`$ (153)
#### 2 Energy

In order to write the perturbation expansion of Eq. 145, we need the expansion of the polarization on the basis of the $`2n+1`$ theorem. The variation-perturbation framework allows us to focus on the part of $`E_{pol}^{(2i)}`$ or $`E_{pol}^{(2i+1)}`$ that comes from variation of the wavefunctions up to order $`i`$ only. For these quantities, we introduce the notation $`E_{pol}^{(2i,i)}`$ or $`E_{pol}^{(2i+1,i)}`$. The even terms are written

$`E_{pol}^{(2i,i)}`$ $`=`$ $`\left[P\left({\displaystyle \underset{j=0}{\overset{i}{}}}^ju_{nk_j}^{(j)}\right)\right]^{(2i)}={\displaystyle \frac{1}{(2i)!}}{\displaystyle \frac{^{2i}}{^{2i}}}\left[P\left({\displaystyle \underset{j=0}{\overset{i}{}}}^ju_{nk_j}^{(j)}\right)\right]|_{=0}`$ (154)
$`=`$ $`\left[P\left({\displaystyle \underset{j=0}{\overset{i}{}}}^ju_{nk_j}^{(j)}\right)\right]^{(2i1)}.`$ (155)
By the same token, for the odd terms we obtain

$`E_{pol}^{(2i+1,i)}=\left[P\left({\displaystyle \underset{j=0}{\overset{i}{}}}^ju_{nk_j}^{(j)}\right)\right]^{(2i+1)}=\left[P\left({\displaystyle \underset{j=0}{\overset{i}{}}}^ju_{nk_j}^{(j)}\right)\right]^{(2i)}.`$ (156)
In these expressions, the $`2n+1`$-theorem implies that only the contributions of order $`i`$ from the perturbed wave functions will appear, when we consider the contribution of the polarization term to the total-energy derivatives. More explicit formulas for computing these polarization derivatives will be given below.

With these results, we can expand Eq. 145. From Eqs. 81 and 154 we obtain the even terms

$`E^{(2i)}`$ $`=`$ $`E_{pol}^{(2i,i)}+{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}{\displaystyle \underset{m,n=1}{\overset{N}{}}}[u_{mk_j}^{(i)}\left|H_{k_j}^{(0)}\right|u_{nk_j}^{(i)}\delta _{mn}`$ (158)
$`{\displaystyle \underset{l^{},l^{\prime \prime }=1}{\overset{i}{}}}{\displaystyle \underset{l=0}{\overset{i1}{}}}\delta (2ill^{}l^{\prime \prime })\mathrm{\Lambda }_{mn}^{(l)}(k_j)u_{nk_j}^{(l^{})}|u_{mk_j}^{(l^)}],`$
while from Eqs. 81 and 156, the odd terms are written

$`E^{(2i+1)}=E_{pol}^{(2i+1,i)}+{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}{\displaystyle \underset{m,n=1}{\overset{N}{}}}{\displaystyle \underset{l^{},l^{\prime \prime }=1}{\overset{i}{}}}{\displaystyle \underset{l=1}{\overset{i}{}}}\delta (2i+1ll^{}l^{\prime \prime })\mathrm{\Lambda }_{mn}^{(l)}(k_j)u_{nk_j}^{(l^{})}|u_{mk_j}^{(l^)}.`$ (159)
#### 3 Sternheimer equation

The perturbation expansion of Eq. 148 yields the Sternheimer equation

$`P_{ck_j}(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)})P_{ck_j}|u_{nk_j}^{(i)}={\displaystyle }_{m=1}^N\{{\displaystyle \underset{l=1}{\overset{i1}{}}}P_{ck_j}|u_{mk_j}^{(il)}\mathrm{\Lambda }_{mn}^{(l)}(k_j)`$ (160)
$`{\displaystyle \frac{ie}{2\mathrm{\Delta }k}}{\displaystyle \underset{l=0}{\overset{i1}{}}}P_{ck_j}\left[\right|u_{mk_{j+1}}^{(il1)}S_{mn}^{1(l)}(k_j,k_{j+1})|u_{mk_{j1}}^{(il1)}S_{mn}^{1(l)}(k_j,k_{j1})]\}=0.`$ (161)
The conduction-band projector $`P_{ck}`$ in Eqs. 121 and 160 appears due the fact that the Lagrange multipliers in these two equations may be different, due to $`k`$-gauge freedom. Thus, only conduction-band contributions can be identified while comparing the two formulations.

By comparing Eqs. 121 and 160, we see that in the present formulation the term $`ieP_{ck}\frac{}{k}|u_{nk}`$ is approximated by the finite-difference formula

$`D(\mathrm{\Delta }k)={\displaystyle \frac{ie}{2\mathrm{\Delta }k}}{\displaystyle \underset{m=1}{\overset{N}{}}}P_{ck_j}\left[|u_{mk_{j+1}}S_{mn}^1(k_j,k_{j+1})|u_{mk_{j1}}S_{mn}^1(k_j,k_{j1})\right].`$ (162)
The theory of finite-difference approximations to derivatives could now be applied to Eq. 162, as a function of $`\mathrm{\Delta }k`$. Note that this expression is invariant under $`\mathrm{\Delta }k\mathrm{\Delta }k`$, as it induces simultaneous exchange of $`k_{j+1}`$ and $`k_{j1}`$. Thus

$`D(\mathrm{\Delta }k)=D(0)+𝒪(\mathrm{\Delta }k)^2.`$ (163)
One can now define

$`D(2\mathrm{\Delta }k)={\displaystyle \frac{ie}{2.2\mathrm{\Delta }k}}{\displaystyle \underset{m=1}{\overset{N}{}}}P_{ck_j}\left[|u_{mk_{j+2}}S_{mn}^1(k_j,k_{j+2})|u_{mk_{j2}}S_{mn}^1(k_j,k_{j2})\right],`$ (164)
giving a higher-order approximation of $`D(0)`$ as

$`D^{higherorder}=\left[4D(\mathrm{\Delta }k)D(2\mathrm{\Delta }k)\right]/3=D(0)+𝒪(\mathrm{\Delta }k)^4.`$ (165)
This improved expression also derives from a total energy functional. Instead of Eq. 16, one must start from

$`E^{higherorder}[\left\{u_{nk_j}\right\};]`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}u_{nk_j}\left|H_{k_j}^{(0)}\right|u_{nk_j}`$ (166)
$``$ $`{\displaystyle \frac{2e}{N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\left\{{\displaystyle \frac{4}{3}}\mathrm{ln}det\left[S_{mn}(k_j,k_{j+1})\right]{\displaystyle \frac{1}{6}}\mathrm{ln}det\left[S_{mn}(k_j,k_{j+2})\right]\right\}.`$ (167)
Despite its interest, we will not explore this topics further in the present paper.

### B Lower-order expressions

In Appendix B, we derive the Taylor expansion of Eq. 3 for the polarization, which allows us to obtain explicit expressions for Eqs. 154 and 156. Here, we look at the lower-order expressions for the energy and the Sternheimer equation.

#### 1 Energy

From Appendix B, the second-order polarization term is given by

$`E_{pol}^{(2,1)}`$ $`=`$ $`\left[P\left(\left\{u_{nk_j}^{(0)}+u_{nk_j}^{(1)}\right\}\right)\right]^{(1)}={\displaystyle \frac{2e}{N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\left\{\mathrm{Tr}\left[S^{(1)}(k_j,k_{j+1})Q(k_j,k_{j+1})\right]\right\}`$ (168)
$`=`$ $`{\displaystyle \frac{2e}{N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\left\{{\displaystyle \underset{m,n=1}{\overset{N}{}}}\left[u_{nk_j}^{(1)}|u_{mk_{j+1}}^{(0)}+u_{nk_j}^{(0)}|u_{mk_{j+1}}^{(1)}\right]Q_{mn}(k_j,k_{j+1})\right\},`$ (169)
where $`Q(k_j,k_{j+1})`$, obeying

$$\underset{l}{}Q_{ml}(k_j,k_{j+1})S_{ln}^{(0)}(k_j,k_{j+1})=\delta _{mn},$$
(170)
is the inverse of the zeroth-order overlap matrix $`S_{nm}^{(0)}(k_j,k_{j+1})=u_{nk_j}^{(0)}|u_{mk_{j+1}}^{(0)}.`$

The second-order expression for the energy is then given by the $`i=1`$ term in Eq. 158

$`E^{(2)}`$ $`=`$ $`{\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk_j}^{(1)}\left|\left(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)}\right)\right|u_{nk_j}^{(1)}`$ (172)
$`{\displaystyle \frac{e}{\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\left\{{\displaystyle \underset{m,n=1}{\overset{N}{}}}\left(u_{nk_j}^{(1)}|u_{mk_{j+1}}^{(0)}+u_{nk_j}^{(0)}|u_{mk_{j+1}}^{(1)}\right)Q_{mn}(k_j,k_{j+1})\right\}.`$
The third-order derivative of the energy is given by by the $`i=1`$ in Eq. 159. The first-order contribution to the Lagrange multipliers vanishes, since from Eqs. 96 and 153 we have $`\mathrm{\Lambda }_{mn}^{(1)}(k_j)=u_{mk_j}^{(1)}\left|H_{k_j}^{(0)}\right|u_{nk_j}^{(0)}+u_{mk_j}^{(0)}\left|H_{k_j}^{(0)}\right|u_{nk_j}^{(1)}=\epsilon _{nk_j}^{(0)}u_{mk_j}^{(1)}|u_{nk_j}^{(0)}+\epsilon _{mk_j}^{(0)}u_{mk_j}^{(0)}|u_{nk_j}^{(1)}=0`$, which leads to $`E^{(3)}=E_{pol}^{(3,1)}`$. From the results in Appendix B

$`E^{(3)}`$ $`=`$ $`\left[P\left(\left\{u_{nk_j}^{(0)}+Fu_{nk_j}^{(1)}\right\}\right)\right]^{(2)}`$ (173)
$`=`$ $`{\displaystyle \frac{e}{N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\{\mathrm{Tr}[2S^{(2)}(k_j,k_{j+1})Q(k_j,k_{j+1})`$ (175)
$`S^{(1)}(k_j,k_{j+1})Q(k_j,k_{j+1})S^{(1)}(k_j,k_{j+1})Q(k_j,k_{j+1})]\}`$
$`=`$ $`{\displaystyle \frac{e}{N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\{{\displaystyle \underset{m,n=1}{\overset{N}{}}}2u_{nk_j}^{(1)}|u_{mk_{j+1}}^{(1)}Q_{mn}(k_j,k_{j+1})`$ (178)
$`{\displaystyle \underset{m,n,l,l^{}=1}{\overset{N}{}}}\left[u_{mk_j}^{(1)}|u_{nk_{j+1}}^{(0)}+u_{mk_j}^{(0)}|u_{nk_{j+1}}^{(1)}\right]Q_{nl}(k_j,k_{j+1})`$
$`[u_{lk}^{(1)}|u_{l^{}k_{j+1}}^{(0)}+u_{lk}^{(0)}|u_{l^{}k_{j+1}}^{(1)}]Q_{l^{}m}(k_j,k_{j+1})\}.`$
For the fourth- and higher-order energy derivatives, the expansion yields very involved expressions. We end this section by considering the fourth-order term for the energy in a more compact notation \[we drop the $`(k_j,k_{j+1})`$ matrix arguments\]:

$`E_{pol}^{(4,2)}`$ $`=`$ $`\left[P\left(\left\{u_{nk_j}^{(0)}+Fu_{nk_j}^{(1)}+F^2u_{nk_j}^{(2)}\right\}\right)\right]^{(3)}`$ (179)
$`=`$ $`{\displaystyle \frac{2e}{3N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\left\{\mathrm{Tr}\left[3S^{(3)}Q3S^{(2)}QS^{(1)}Q+S^{(1)}QS^{(1)}QS^{(1)}Q\right]\right\}`$ (180)
$`=`$ $`{\displaystyle \frac{2e}{3N_k\mathrm{\Delta }k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}\mathrm{Im}\{3{\displaystyle \underset{m,n=1}{\overset{N}{}}}[u_{nk_j}^{(1)}|u_{mk_{j+1}}^{(2)}+u_{nk_j}^{(2)}|u_{mk_{j+1}}^{(1)}]Q_{mn}`$ (181)
$``$ $`3{\displaystyle \underset{m,n=1}{\overset{N}{}}}\left[u_{nk_j}^{(0)}|u_{nk_{j+1}}^{(2)}+u_{nk_j}^{(1)}|u_{mk_{j+1}}^{(1)}+u_{nk_j}^{(2)}|u_{mk_{j+1}}^{(0)}\right]\left(QS^{(1)}Q\right)_{mn}`$ (182)
$`+`$ $`S^{(1)}QS^{(1)}QS^{(1)}Q\}`$ (183)
In this expression, we write explicitly only the terms containing $`u_{nk_j}^{(2)}`$, which will determine the second-order term of the Sternheimer equation. The corresponding fourth-order energy is given by

$`E^{(4)}={\displaystyle \frac{2}{N_k}}{\displaystyle \underset{k=1}{\overset{N_k}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk_j}^{(2)}\left|\left(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)}\right)\right|u_{nk_j}^{(2)}+E_{pol}^{(4,2)}.`$ (184)
#### 2 Sternheimer equation

From Eq. 160, the first-order term for the Sternheimer equation reads

$`P_{ck_j}\left(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)}\right)P_{ck_j}|u_{nk_j}^{(1)}={\displaystyle \frac{ie}{2\mathrm{\Delta }k}}P_{ck}{\displaystyle \underset{m=1}{\overset{N}{}}}\left[|u_{mk_{j+1}}^{(0)}Q_{mn}(k_j,k_{j+1})|u_{mk_{j1}}^{(0)}Q_{mn}(k_j,k_{j1})\right];`$ (185)
while the second-order derivative is written

$`P_{ck_j}(H_{k_j}^{(0)}\epsilon _{nk_j}^{(0)})P_{ck_j}|u_{nk_j}^{(2)}={\displaystyle \frac{ie}{2\mathrm{\Delta }k}}P_{ck}\{{\displaystyle \underset{m=1}{\overset{N}{}}}\left[\right|u_{mk_{j+1}}^{(1)}Q_{mn}(k_j,k_{j+1})|u_{mk_{j1}}^{(1)}Q_{mn}(k_j,k_{j1})]`$ (187)
$`{\displaystyle \underset{m,l,l^{}=1}{\overset{N}{}}}\left[\right|u_{mk_{j+1}}^{(0)}Q_{ml}(k_j,k_{j+1})S_{l,l^{}}^{(1)}(k_j,k_{j+1})Q_{l^{}n}(k_j,k_{j+1})`$ (188)
$`|u_{mk_{j1}}^{(0)}Q_{ml}(k_j,k_{j1})S_{l,l^{}}^{(1)}(k_j,k_{j1})Q_{l^{}n}(k_j,k_{j1})]\}.`$ (189)
These two expressions can also be consistently obtained by taking the conduction-band projection of the corresponding terms in the expansion of $`\delta E^{(2i)}/\delta u_{nk_j}^{(i)}=0`$.

## VI Numerical tests

In this section, we illustrate the present theory by applying it the two-band 1D-TB Hamiltonian introduced in Sec. II. For this model, exact analytical expressions can be written for the continuous formulation. Our purpose in this simple application is to demonstrate consistency between the continuous and the discretized versions of the theory, and also to make a preliminary assessment of the convergence of the energy derivatives obtained in the two discretized formulations, with respect to $`k`$-point sampling in the Brillouin zone. Appendix C contains the detailed derivations of the results presented here. As such, it presents a step-by-step example of the use of the formalism developed in the present paper.

### A Response to a homogeneous electric field

In the continuous formulation, the first-order change of the valence state is obtained from the corresponding term in the Sternheimer equation. By setting $`i=1`$ in Eq. 121, this is given by

$`P_{ck}\left(H_k^{(0)}\epsilon _{vk}^{(0)}\right)P_{ck}|u_{vk}^{(1)}=P_{ck}i{\displaystyle \frac{}{k}}|u_{vk}^{(0)}.`$ (190)
The solution to this equation is given in Appendix C.

Having solved Eq. 190, $`E^{(2)}`$ can then be obtained from the simplest non-variational expression Eq. 106 as

$`E^{(2)}={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}𝑑k{\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}\left({\displaystyle \frac{\mathrm{\Theta }_k}{k}}\right)^2,`$ (191)
where $`\mathrm{\Delta }\epsilon _k=\epsilon _{ck}^{(0)}\epsilon _{vk}^{(0)}`$, and $`\mathrm{\Theta }_k/k`$ is obtained from Eqs. 41 and 48 (see Eq. C5).

The second-order change of the valence state is obtained from the $`i=2`$ term in Eq. 121 and from Eq. 126. In our 1D model, the Sternheimer equation for the conduction-band projection reads

$`P_{ck}\left(H_k^{(0)}\epsilon _{vk}^{(0)}\right)P_{ck}|u_{vk}^{(2)}=P_{ck}\left[i{\displaystyle \frac{}{k}}|u_{vk}^{(1)}\mathrm{\Lambda }_k^{(1)}|u_{vk}^{(1)}\right].`$ (192)
Using the solution for $`|u_{vk}^{(2)}`$ given in Appendix C, we can now arrive at the expression for the fourth-order energy. Applying Eq. 109, we write

$`E^{(4)}={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{2\pi }}𝑑k\left[\left({\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}\right)^3\left({\displaystyle \frac{\mathrm{\Theta }_k}{k}}\right)^4\left({\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}\right)\left\{{\displaystyle \frac{}{k}}\left({\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}{\displaystyle \frac{\mathrm{\Theta }_k}{k}}\right)\right\}^2\right].`$ (193)
(Expressions for $`E^{(2)}`$ and $`E^{(4)}`$ in terms of elliptic integrals of the second kind are given in Appendix C.)

Turning now to the discretized expressions, using the non-variational Eq. 106 we compute the DAPE expression for the second-order energy:

$`E^{(2)}={\displaystyle \frac{2}{N_k}}{\displaystyle \underset{j=1}{\overset{N_k}{}}}{\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _{k_j}}}\left[{\displaystyle \frac{1}{2\mathrm{\Delta }k}}\mathrm{sin}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j1}\right)\mathrm{cos}\left(2\mathrm{\Theta }_j\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j1}\right)\right]^2.`$ (194)
The PEAD second-order energy expression takes a different form:

$`E^{(2)}={\displaystyle \frac{2}{N_k}}{\displaystyle \underset{k=1}{\overset{N_k}{}}}{\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _j}}\left\{{\displaystyle \frac{1}{2\mathrm{\Delta }k}}\left[\mathrm{tan}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\mathrm{tan}\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_j\right)\right]\right\}^2.`$ (195)
With the results in Appendix C, we have all the ingredients to write analytical expressions for the discretized versions of $`E^{(4)}`$, but these are quite cumbersome and will not be reproduced here. The numerical results obtained using these expressions are discussed below.

### B Numerical results

In order to test the consistency between the three formulations, we checked that by sufficiently increasing the number of $`k`$-points used in the evaluation of the discretized energy derivatives, we obtain an agreement between continuous and discretized expressions within stringent degrees of accuracy. For example, for $`t=1`$, 80 $`k`$-points in the full Brillouin zone are needed for the three expressions for $`E^{(2)}`$ (Eqs. 191, 194, and 195) to agree within $`1`$%, while 240 $`k`$-points are needed to get the same level of agreement for the $`E^{(4)}`$ expressions.

By decreasing the number of $`k`$-points, thus worsening the level of accuracy of the discretized expressions, the differences between them become more apparent. In Fig. 3, we show the quantity $`\left[E_{discr.}^{(2)}E_{exact}^{(2)}\right]/E_{exact}^{(2)}`$ giving the percentual error in the evaluation of $`E^{(2)}`$ for the two discretized formulations, using 20 $`k`$-points, with the hoping parameter varying over the  interval. In Fig. 4 the corresponding quantity for $`E^{(4)}`$, $`\left[E_{discr}^{(4)}E_{exact}^{(4)}\right]/E_{exact}^{(4)}`$, is shown for a sampling of 80 $`k`$-points.

It is clear from Figs. 3 and 4 that the energy derivatives obtained from the PEAD converge faster with respect to the number of $`k`$-points, at least for the present 1D model. However, it must be borne in mind that this formulation involves the calculation of the inverse of the zero-field overlap matrix, and in practical applications the additional cost associated with this operation could offset the gain in $`k`$-point convergence, specially when the two formulations are applied to systems with large numbers of atoms in the unit cell. This point remains to be further addressed when the theory is applied in the context of realistic tight-binding and ab initio calculations.

We also computed the norm of the first- and second-order wavefunction derivatives $`\left[u_{vk}^{(i)}|u_{vk}^{(i)}\right]^{1/2}`$ as a function of $`k`$, for the continuous solutions and the two discretized forms, with a value of $`t=1`$ for the hoping parameter and samplings of 20 and 80 $`k`$-points. As expected from the above results for $`E^{(2)}`$ and $`E^{(4)}`$, we observe that the wavefunctions in the PEAD are better approximations to the exact ones from the continuous formulation.

## VII Summary

The goal of this work was to obtain second- and higher-order derivatives of the total energy of periodic insulators, with respect to an applied homogeneous electric field. Related physical properties are the linear and non-linear dielectric susceptibilities (connected to linear and non-linear optical constants).

Although a variation-perturbation framework had been formulated earlier for the computation of derivatives of the total energy with respect to many different perturbations, several formal and technical difficulties must be addressed when considering the specific case of a homogeneous-electric-field perturbation.

At the level of the electric-field-dependent energy functional, we proposed a basic expression, Eq. 10, that we argue to be valid in the space of states possessing a periodic density. It is directly linked to the modern theory of polarization, proposed by King-Smith and Vanderbilt nearly a decade ago, and allows to recover easily the Berry phase polarization formula, central to this theory.

Unfortunately, when the polarization is varied, this energy functional leads to local minima with a basin of attraction of infinitesimal extent. A regularization procedure, based on a $`k`$-point discretization of the reciprocal-space integrals, must be used in order to lead to a finite size basin. This is the reciprocal-space analog of the real-space cut-off introduced by Nunes and Vanderbilt in their treatment of polarized Wannier functions.

Having thus defined a suitable energy functional, that depends on the applied homogeneous electric field, we were allowed to proceed with the application of the variation-perturbation machinery. Interestingly, the derivation of the canonical formulas at all orders of perturbation can be done either on the basis of the energy functional already regularized, or on the basis of the unregularized one, followed by regularization at each order. The formulas derived in the two cases differ from each other. Working with the regularized energy functional gives more cumbersome expressions, however perfectly consistent with an energy functional, while the a posteriori application of the regularization at each order is not consistent with the regularized energy functional. The two procedures will tend to the same limit when the discretization is refined further and further. The expression for the third-order derivative of the total energy, previously proposed by Mauri and Dal Corso, is recovered, as an instance of the ”a posteriori” application of the regularization technique.

We applied this formalism to a model one-dimensional two-band Hamiltonian, showing explicitly the pathology of a non-regularized energy functional and its cure, as well as the differences related to the order in which the perturbation expansion and the discretization procedure are applied. The two discretized formulations are shown to agree in the continuum limit, although the “perturbation expansion after discretization” (PEAD) formulation seems to be closer to the exact answer than the “discretization after perturbation expansion” (DAPE) formulation, for an equivalent grid.

## A

In this appendix, we describe the algebraic manipulations needed to transform Eqs 102, 103, and 121 into their gauge-invariant forms presented in Secs. IV A 2 and IV A 3. We will use the following result:

$`x\left|{\displaystyle \frac{}{k}}\right|y|z=x|y{\displaystyle \frac{}{k}}|z+|z{\displaystyle \frac{}{k}}x|y\left({\displaystyle \frac{}{k}}\right|zx|)|y.`$ (A1)
Furthermore, from Eq. 78 we have

$`{\displaystyle \underset{l=0}{\overset{i}{}}}u_{mk}^{(l)}|u_{nk}^{(il)}=0;i1.`$ (A2)
With the help of Eq. A1 and $`\delta _{mn}=u_{mk}^{(0)}|u_{nk}^{(0)}`$ , it is straightforward to show that

$`\left[{\displaystyle \frac{}{k}}\delta _{mn}u_{mk}^{(0)}\left|{\displaystyle \frac{}{_k}}\right|u_{nk}^{(0)}\right]|u_{mk}^{(i)}=\left({\displaystyle \frac{}{k}}|u_{mk}^{(i)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)}.`$ (A3)
We recall that the notation “$`\left(_k|uu|\right)`$” indicates that $`_k`$ acts only on the expression in parenthesis.

Now, for reference we repeat Eqs. 102 and 103 here:

$`E^{(2i)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}dk[{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(i)}\left|H_k^{(0)}\right|u_{nk}^{(i)}+u_{nk}^{(i1)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(i)}+u_{nk}^{(i)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(i1)}`$ (A5)
$`{\displaystyle \underset{m,n=1}{\overset{N}{}}}{\displaystyle \underset{j,j^{}=1}{\overset{i}{}}}{\displaystyle \underset{l=0}{\overset{i1}{}}}\delta (2ijj^{}l)\mathrm{\Lambda }_{mn}^{(l)}(k_j)u_{nk}^{(j)}|u_{mk}^{(j^{})}].`$
$`E^{(2i+1)}={\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k\left[{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(i)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(i)}{\displaystyle \underset{n,m=1}{\overset{N}{}}}{\displaystyle \underset{j,j^{},l=1}{\overset{i}{}}}\delta (2i+1jj^{}l)\mathrm{\Lambda }_{mn}^{(l)}(k_j)u_{nk}^{(j)}|u_{mk}^{(j^{})}\right].`$ (A6)
Let us examine the $`i=1`$ and $`i=2`$ terms in Eq. A5, and the $`i=1`$ term in Eq. A6. We use Eq. 98 to write

$`E^{(2)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}\left|\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)\right|u_{nk}^{(1)}+u_{nk}^{(0)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(1)}+u_{nk}^{(1)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)}.`$ (A8)
Eq. 96 and the orthonormality relation $`u_{mk}^{(0)}|u_{nk}^{(0)}=\delta _{mn}`$ allow us to rewrite the last two terms in this equation in the explicit gauge-invariant form of Eq. 104, as follows:

$`u_{nk}^{(1)}\left|{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)}`$ $`=`$ $`{\displaystyle \underset{m}{}}u_{nk}^{(1)}|_ku_{mk}^{(0)}u_{mk}^{(0)}|u_{nk}^{(0)}+u_{nk}^{(1)}|u_{mk}^{(0)}_ku_{mk}^{(0)}|u_{nk}^{(0)}`$ (A9)
$`=`$ $`u_{nk}^{(1)}|\left({\displaystyle \frac{}{k}}{\displaystyle \underset{m}{}}|u_{mk}^{(0)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)},`$ (A10)
where we use the notation $`|_ku_{mk}\frac{}{k}|u_{mk}`$. The second term in Eq. A8 is obtained simply as the hermitian conjugate of this latter equation.

For the third and fourth-order terms, we have

$`E^{(3)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(1)}{\displaystyle \underset{m,n=1}{\overset{N}{}}}u_{nk}^{(1)}|u_{mk}^{(1)}u_{mk}^{(0)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)},`$ (A11)
and

$`E^{(4)}`$ $`=`$ $`{\displaystyle \frac{a}{\pi }}{\displaystyle _0^{\frac{2\pi }{a}}}𝑑k{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(2)}\left|\left(H_k^{(0)}\epsilon _{nk}^{(0)}\right)\right|u_{nk}^{(2)}+u_{nk}^{(1)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(2)}+u_{nk}^{(2)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(1)}`$ (A13)
$`{\displaystyle \underset{m,n=1}{\overset{N}{}}}\left(u_{nk}^{(1)}|u_{mk}^{(2)}+u_{nk}^{(2)}|u_{mk}^{(1)}\right)u_{mk}^{(0)}\left|ie{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)}.`$
We now apply Eq. A3 to recombine terms in these two expressions as follows:

$`{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}\left|{\displaystyle \frac{}{k}}\right|u_{nk}^{(j)}{\displaystyle \underset{m,n=1}{\overset{N}{}}}u_{nk}^{(1)}|u_{mk}^{(j)}u_{mk}^{(0)}\left|{\displaystyle \frac{}{k}}\right|u_{nk}^{(0)}`$ (A14)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{N}{}}}u_{nk}^{(1)}|\left({\displaystyle \frac{}{k}}{\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk}^{(j)}u_{mk}^{(0)}|\right)|u_{nk}^{(0)};`$ (A15)
where $`j=1`$ applies to $`E^{(3)}`$ and $`j=2`$ to $`E^{(4)}`$. This leads to the gauge-invariant expressions for these quantities, Eqs. 117 and 109, respectively.

For the general energy derivative, besides the terms corresponding to those above for $`E^{(3)}`$ and $`E^{(4)}`$, we must also consider terms of the form

$`\mathrm{\Lambda }_{mn}^{(l)}(k)u_{nk}^{(j)}|u_{mk}^{(j^{})}`$ $`=`$ $`{\displaystyle \underset{l^{}=0}{\overset{l}{}}}u_{mk}^{(l^{})}\left|H_k^{(0)}\right|u_{nk}^{(ll^{})}u_{nk}^{(j)}|u_{mk}^{(j^{})}+{\displaystyle \underset{l^{}=0}{\overset{l1}{}}}u_{mk}^{(l^{})}\left|ie{\displaystyle \frac{}{_k}}\right|u_{nk}^{(ll^{}1)}u_{nk}^{(j)}|u_{mk}^{(j^{})}.`$ (A16)
For the second term on the right, we use Eqs. A1 and A2 to write

$`u_{nk}^{(j)}\left|\left({\displaystyle \underset{l^{}=0}{\overset{l1}{}}}u_{mk}^{(l^{})}\left|{\displaystyle \frac{}{_k}}\right|u_{nk}^{(ll^{}1)}\right)\right|u_{mk}^{(j^{})}={\displaystyle \underset{l^{}=0}{\overset{l1}{}}}u_{nk}^{(j)}|\left({\displaystyle \frac{}{_k}}|u_{mk}^{(j^{})}u_{mk}^{(l^{})}|\right)|u_{nk}^{(ll^{}1)}.`$ (A18)
These results demonstrate the gauge invariance of the general expansion terms for the energy. The proper gauge-transformation properties of the Sternheimer equation can also be proven explicitly along the same lines, but it follows more simply from the invariance of the even-order terms for the energy.

## B

In order to develop the PEAD formulation, we examine the following expression appearing in Eqs. 3 and 16:

$`lndet\left[S_{nm}(k_j,k_{j+1})()\right]`$ $`=`$ $`lndet\left[S_{nm}^{(0)}(k_j,k_{j+1})+S_{nm}^{(1)}(k_j,k_{j+1})+\mathrm{}\right]`$ (B1)
$`=`$ $`lndet\left[S_{nm}^{(0)}(k_j,k_{j+1})\right]+{\displaystyle _0^{}}𝑑{\displaystyle \frac{}{}}ln𝑑et\left[S_{nm}(k_j,k_{j+1})()\right],`$ (B2)
where the perturbation expansion of $`S_{nm}(k_j,k_{j+1})()`$ is defined according to Eq. 71. With the help of the “magic” formula

$`{\displaystyle \frac{}{\lambda }}\left\{lndet\left[A(\lambda )\right]\right\}=tr\left[A^1{\displaystyle \frac{A}{\lambda }}\right]={\displaystyle \underset{m,n}{}}A_{mn}^1{\displaystyle \frac{A_{nm}}{\lambda }},`$ (B3)
we can rewrite Eq. B1 (to simplify the notation, in the following equations we drop the $`k`$-point arguments) as follows:

$`lndet\left[S_{nm}()\right]=lndet\left[S_{nm}^{(0)}\right]+{\displaystyle _0^{}}𝑑{\displaystyle \underset{mn}{}}{\displaystyle \frac{S_{nm}()}{}}S_{mn}^1().`$ (B4)
To proceed further, we let $`\mathrm{\Delta }S()=S()S^{(0)}=S^{(1)}()+^2S^{(2)}()+^3S^{(3)}()+\mathrm{}`$, and the inverse overlap matrix can be written as

$`S^1()=\left\{S^{(0)}\left[I+Q\mathrm{\Delta }S()\right]\right\}^1=\left\{I+{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}(1)^i\left[Q\mathrm{\Delta }S()\right]^i\right\}Q,`$ (B5)
where, following the notation introduced in Sec. V B, we use $`Q=[S^{(0)}]^1`$ for the inverse of the zeroth-order overlap matrix.

To obtain the lower-order terms in the expansion of Eq. B1, we write the terms up to third-order explicitly:

$`S^1()`$ $`=`$ $`QQS^{(1)}Q+^2\left[QS^{(2)}Q+QS^{(1)}QS^{(1)}Q\right]`$ (B6)
$`+`$ $`^3[QS^{(3)}Q+QS^{(2)}QS^{(1)}Q+QS^{(1)}QS^{(2)}Q`$ (B7)
$``$ $`QS^{(1)}QS^{(1)}QS^{(1)}Q]+𝒪(^4).`$ (B8)
Combining Eqs. B4, B6 and $`S()/=S^{(1)}+2S^{(2)}+3^2S^{(3)}+4^3S^{(4)}+𝒪(^4)`$, we arrive at

$`lndet\left[S()\right]`$ $`=`$ $`lndet\left[S^{(0)}\right]+tr\left[S^{(1)}Q\right]+{\displaystyle \frac{^2}{2}}\mathrm{Tr}\left[2S^{(2)}QS^{(1)}QS^{(1)}Q\right]`$ (B9)
$`+`$ $`{\displaystyle \frac{^3}{3}}tr\left[3S^{(3)}Q3S^{(2)}QS^{(1)}Q+S^{(1)}QS^{(1)}QS^{(1)}Q\right]`$ (B10)
$`+`$ $`{\displaystyle \frac{^4}{4}}tr[4S^{(4)}Q4S^{(3)}QS^{(1)}Q2S^{(2)}QS^{(2)}Q`$ (B11)
$`+`$ $`4S^{(2)}QS^{(1)}QS^{(1)}QS^{(1)}QS^{(1)}QS^{(1)}QS^{(1)}Q]+𝒪(^5).`$ (B12)
Finally, we note that the magic formula may also be used to derive the Euler-Lagrange equation given in Eq. 148, as follows:

$`{\displaystyle \frac{\delta \left\{\mathrm{ln}det\left[\frac{}{}S_{nm}(k_j,k_{j+1})\right]\right\}}{\delta u_{nk}^{}}}=\mathrm{Tr}\left[S^1{\displaystyle \frac{\delta S}{\delta u_{nk}^{}}}\right]={\displaystyle \underset{m=1}{\overset{N}{}}}|u_{mk_{j+1}}S_{mn}^1(k_j,k_{j+1}).`$ (B13)
## C

In this appendix, we provide all the steps for the derivation of the energy derivatives for the 1D-TB model, given in Sec. VI. The expressions for the discretized versions of $`E^{(4)}`$ are not written explicitly.

Eq. 49 for the unperturbed total energy is a complete elliptic integral of the second kind , which is given in its general form by

$`_{\frac{n}{2}}={\displaystyle _0^{\frac{\pi }{2}}}𝑑y\left(1+A^2\mathrm{cos}^2y\right)^{\frac{n}{2}},`$ (C1)
where $`n`$ is a positive or negative integer, and $`A=4t`$. Several such integrals will be encountered in the course of our derivation.

### 1 Continuous formulation

We begin by developing some useful preliminary results. For our 1D-TB model, the derivatives with respect to $`k`$ of the unperturbed valence and conduction states (Eq. 44) are written analytically as

$`i{\displaystyle \frac{}{k}}|u_{vk}^{(0)}={\displaystyle \frac{\alpha _{vk}}{k}}|u_{vk}^{(0)}i{\displaystyle \frac{\mathrm{\Theta }_k}{k}}e^{i\mathrm{\Delta }\alpha _k}|u_{ck}^{(0)}`$ (C2)
$`i{\displaystyle \frac{}{k}}|u_{ck}^{(0)}={\displaystyle \frac{\alpha _{ck}}{k}}|u_{ck}^{(0)}+i{\displaystyle \frac{\mathrm{\Theta }_k}{k}}e^{i\mathrm{\Delta }\alpha _k}|u_{vk}^{(0)},`$ (C3)
where $`\mathrm{\Delta }\alpha _k=\alpha _{ck}\alpha _{vk}`$.

Recalling that $`u_{vk}^{(0)}|u_{vk}^{(1)}=0`$, hence $`P_{vk}|u_{vk}^{(1)}=0`$, the solution of Eq. 190 gives the first-order change of the valence state

$`|u_{vk}^{(1)}=P_{ck}|u_{vk}^{(1)}={\displaystyle \frac{i}{\mathrm{\Delta }\epsilon _k}}{\displaystyle \frac{\mathrm{\Theta }_k}{k}}e^{i\mathrm{\Delta }\alpha _k}|u_{ck}^{(0)},`$ (C4)
where $`\mathrm{\Delta }\epsilon _k=\epsilon _{ck}^{(0)}\epsilon _{vk}^{(0)}=\left[1+A^2\mathrm{cos}^2\left(\frac{k}{2}\right)\right]^2`$. From Eqs. 41 and 48, the partial derivative appearing in this expression is simply

$`{\displaystyle \frac{\mathrm{\Theta }_k}{k}}={\displaystyle \frac{tsin\left(\frac{k}{2}\right)}{\left[1+A^2\mathrm{cos}^2\left(\frac{k}{2}\right)\right]}}.`$ (C5)
Using Eq. C5, the second-order energy term (Eq. 191) is written

$`E^{(2)}={\displaystyle \frac{1}{4\pi }}\left[_{\frac{3}{2}}\left(1+A^2\right)_{\frac{5}{2}}\right].`$ (C6)
The second-order change of the valence state is obtained from Eqs. 124 and 126. The first-order Lagrange multiplier in Eq. 124 is obtained from Eqs. 44, 98, and C2:

$`\mathrm{\Lambda }_k^{(1)}=u_{vk}^{(0)}\left|i{\displaystyle \frac{}{k}}\right|u_{vk}^{(0)}={\displaystyle \frac{\alpha _{vk}}{k}}.`$ (C7)
Simple algebraic manipulations involving Eqs. C2C4, and C7, combined with Eqs. 124 and 126, yield the second-order wave-function derivatives

$`P_{ck}|u_{vk}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}\left[{\displaystyle \frac{}{k}}\left({\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}{\displaystyle \frac{\mathrm{\Theta }_k}{k}}\right)\right]e^{i\mathrm{\Delta }\alpha _k}|u_{ck}^{(0)},`$ (C8)
$`P_{vk}|u_{vk}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}|u_{vk}^{(0)}u_{vk}^{(1)}|u_{vk}^{(1)}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _k}}{\displaystyle \frac{\mathrm{\Theta }_k}{k}}\right)^2|u_{vk}^{(0)}.`$ (C9)
These results lead to Eq. 193 for $`E^{(4)}`$. In terms of elliptic integrals, this quantity is written

$`E^{(4)}`$ $`=`$ $`{\displaystyle \frac{37\left(1+A^2\right)^2}{64\pi }}_{\frac{11}{2}}{\displaystyle \frac{\left(1+A^2\right)}{32\pi }}\left[18\left(1+A^2\right)+25\right]_{\frac{9}{2}}+{\displaystyle \frac{1}{64\pi }}\left[48\left(1+A^2\right)+17\right]_{\frac{7}{2}}`$ (C11)
$`{\displaystyle \frac{1}{4\pi }}_{\frac{5}{2}}.`$
### 2 DAPE formulation

We now apply the DAPE expressions obtained in Sec. IV B to our 1D-TB model. The first-order valence state is given by

$`P_{ck_j}\left[H_{k_j}^{(0)}\epsilon _{vk_j}^{(0)}\right]P_{ck_j}|u_{vk_j}^{(1)}={\displaystyle \frac{i}{2\mathrm{\Delta }k}}P_{ck_j}\left[|u_{vk_{j+1}}^{(0)}u_{vk_{j+1}}^{(0)}|u_{vk_j}^{(0)}|u_{vk_{j1}}^{(0)}u_{vk_{j1}}^{(0)}|u_{vk_j}^{(0)}\right]`$ (C12)
$`={\displaystyle \frac{i}{2\mathrm{\Delta }k}}|u_{ck_j}^{(0)}\left[u_{ck_j}^{(0)}|u_{vk_{j+1}}^{(0)}u_{vk_{j+1}}^{(0)}|u_{vk_j}^{(0)}u_{ck_j}^{(0)}|u_{vk_{j1}}^{(0)}u_{vk_{j1}}^{(0)}|u_{vk_j}^{(0)}\right].`$ (C13)
Using Eq. 44, this expression becomes

$`|u_{vk_j}^{(1)}={\displaystyle \frac{i}{2\mathrm{\Delta }k\mathrm{\Delta }\epsilon _j}}\mathrm{sin}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j1}\right)\mathrm{cos}\left(2\mathrm{\Theta }_j\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j1}\right)e^{i\mathrm{\Delta }\alpha _j}|u_{ck_j}^{(0)};`$ (C14)
where we use a simplified notation $`\mathrm{\Theta }_{k_j}\mathrm{\Theta }_j`$ (likewise for $`\mathrm{\Delta }\alpha `$ and $`\mathrm{\Delta }\epsilon `$). From Eqs. C14 and 134, we obtain Eq. 194.

The second-order wavefunction is obtained from Eqs. 126 and 144. The Sternheimer equation is written

$`P_{ck_j}\left[H_{k_j}^{(0)}\epsilon _{vk_j}^{(0)}\right]P_{ck_j}|u_{vk_j}^{(2)}={\displaystyle \frac{i}{2\mathrm{\Delta }k}}P_{ck_j}\left[|u_{vk_{j+1}}^{(1)}u_{vk_{j+1}}^{(0)}|u_{vk_j}^{(0)}|u_{vk_{j1}}^{(1)}u_{vk_{j1}}^{(0)}|u_{vk_j}^{(0)}\right]`$ (C15)
$`={\displaystyle \frac{i}{2\mathrm{\Delta }k}}|u_{ck_j}^{(0)}\left[u_{ck_j}^{(0)}|u_{vk_{j+1}}^{(1)}u_{vk_{j+1}}^{(0)}|u_{vk_j}^{(0)}u_{ck_j}^{(0)}|u_{vk_{j1}}^{(1)}u_{vk_{j1}}^{(0)}|u_{vk_j}^{(0)}\right].`$ (C16)
By combining Eqs. 44, 126, and C14, we arrive at

$`P_{vk_j}|u_{vk_j}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{8\mathrm{\Delta }k^2\mathrm{\Delta }\epsilon _j^2}}\mathrm{sin}^2\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j1}\right)\mathrm{cos}^2\left(2\mathrm{\Theta }_j\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j1}\right)|u_{vk_j}^{(0)}`$ (C17)
$`P_{ck_j}|u_{vk_j}^{(2)}`$ $`=`$ $`{\displaystyle \frac{e^{i\mathrm{\Delta }\alpha _j}}{4\mathrm{\Delta }k^2\mathrm{\Delta }\epsilon _j}}[{\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _{j+1}}}\mathrm{sin}(\mathrm{\Theta }_{j+2}\mathrm{\Theta }_j)\mathrm{cos}(2\mathrm{\Theta }_{j+1}\mathrm{\Theta }_{j+2}\mathrm{\Theta }_j)\mathrm{cos}^2(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j)`$ (C19)
$`{\displaystyle \frac{1}{\mathrm{\Delta }\epsilon _{j1}}}\mathrm{sin}(\mathrm{\Theta }_j\mathrm{\Theta }_{j2})\mathrm{cos}(2\mathrm{\Theta }_{j1}\mathrm{\Theta }_{j2}\mathrm{\Theta }_j)\mathrm{cos}^2(\mathrm{\Theta }_j\mathrm{\Theta }_{j1})]|u_{ck_j}^{(0)}.`$
Eqs. C14 and C17, combined with Eq. 141, lead to an analytic expression for $`E^{(4)}`$.

### 3 PEAD formulation

Here, we apply the PEAD expressions from Sec. V to the model system. The zero-field overlap matrix, $`S^{(0)}(k_j,k_{j+1})=u_{vk_j}^{(0)}|u_{vk_{j+1}}^{(0)}`$ is

$`S^{(0)}(k_j,k_{j+1})=e^{i\left(\alpha _{v_{j+1}}\alpha _{v_j}\right)}\mathrm{cos}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right).`$ (C20)
Note that for this model, $`S`$ is a 1$`\times `$1 matrix, and hence the inverse overlap matrix is simply

$`Q(k_j,k_{j+1})={\displaystyle \frac{1}{S^{(0)}(k_j,k_{j+1})}}={\displaystyle \frac{e^{i\left(\alpha _{v_{j+1}}\alpha _{v_j}\right)}}{\mathrm{cos}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)}}.`$ (C21)
From Eq. 185, the first-order Sternheimer equation reads

$`P_{ck_j}\left[H_{k_j}^{(0)}\epsilon _{vk_j}^{(0)}\right]P_{ck_j}|u_{vk_j}^{(1)}=P_{ck_j}{\displaystyle \frac{i}{2\mathrm{\Delta }k}}\left[|u_{vk_{j+1}}^{(0)}Q(k,k_{j+1})|u_{vk_{j1}}^{(0)}Q(k,k_{j1})\right].`$ (C22)
Using Eq. 44, we obtain

$`|u_{vk_j}^{(1)}={\displaystyle \frac{ie^{i\mathrm{\Delta }\alpha _j}}{2\mathrm{\Delta }k\mathrm{\Delta }\epsilon _k}}|u_{ck}^{(0)}\left[\mathrm{tan}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\mathrm{tan}\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_j\right)\right].`$ (C23)
To obtain the second-order wavefunctions, we need the first-order term for the overlap matrix $`S^{(1)}(k_j,k_{j+1})=u_{vk_j}^{(0)}|u_{vk_{j+1}}^{(1)}+u_{vk_j}^{(1)}|u_{vk_{j+1}}^{(0)}`$. Thanks to Eqs. 44 and C23, we have

$`S^{(1)}(k_j,k_{j+1})`$ $`=`$ $`{\displaystyle \frac{ie^{i\left(\alpha _{v_{j+1}}\alpha _j\right)}}{2\mathrm{\Delta }k}}\mathrm{sin}(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j)[{\displaystyle \frac{\mathrm{tan}\left(\mathrm{\Theta }_{j+2}\mathrm{\Theta }_{j+1}\right)\mathrm{tan}\left(\mathrm{\Theta }_j\mathrm{\Theta }_{j+1}\right)}{\mathrm{\Delta }\epsilon _{j+1}}}+`$ (C25)
$`{\displaystyle \frac{\mathrm{tan}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\mathrm{tan}\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_j\right)}{\mathrm{\Delta }\epsilon _j}}].`$
Plugging the above results for $`Q`$, $`S^0`$, $`S^{(1)}`$, $`u_{vk}^{(0)}`$, and $`u_{vk}^{(1)}`$ in Eqs. 126 and 187 we get

$`P_{vk_j}|u_{vk_j}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{8\mathrm{\Delta }k^2\mathrm{\Delta }\epsilon _j^2}}|u_{vk}^{(0)}\left[\mathrm{tan}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\mathrm{tan}\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_j\right)\right]^2`$ (C26)
$`P_{ck_j}|u_{vk_j}^{(2)}`$ $`=`$ $`|u_{ck_j}^{(0)}{\displaystyle \frac{e^{i\mathrm{\Delta }\alpha _j}}{4\mathrm{\Delta }^2k\mathrm{\Delta }\epsilon _j}}\{{\displaystyle \frac{\left[1+\mathrm{tan}^2\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\right]\left[\mathrm{tan}\left(\mathrm{\Theta }_{j+2}\mathrm{\Theta }_{j+1}\right)\mathrm{tan}\left(\mathrm{\Theta }_j\mathrm{\Theta }_{j+1}\right)\right]}{\mathrm{\Delta }\epsilon _{j+1}}}+`$ (C29)
$`{\displaystyle \frac{\left[1+\mathrm{tan}^2\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_j\right)\right]\left[\mathrm{tan}\left(\mathrm{\Theta }_{j2}\mathrm{\Theta }_{j1}\right)\mathrm{tan}\left(\mathrm{\Theta }_j\mathrm{\Theta }_{j1}\right)\right]}{\mathrm{\Delta }\epsilon _{j1}}}+`$
$`{\displaystyle \frac{\left[\mathrm{tan}^2\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\mathrm{tan}^2\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_j\right)\right]\left[\mathrm{tan}\left(\mathrm{\Theta }_{j+1}\mathrm{\Theta }_j\right)\mathrm{tan}\left(\mathrm{\Theta }_{j1}\mathrm{\Theta }_{j1}\right)\right]}{\mathrm{\Delta }\epsilon _j}}\}.`$
As in the DAPE case, from these expressions for $`u_{vk_j}^{(1)}`$ and $`u_{vk_j}^{(2)}`$ follows the PEAD analytic form for $`E^{(4)}`$.