Thermodynamic perturbation theory for non-interacting quantum particles with application to spin-spin interactions in solids

Abstract

The determination of the Landau free energy (the grand thermodynamic potential) by a perturbation theory is advanced to arbitrary order for the specific case of non-interacting fermionic systems perturbed by a one-particle potential. Peculiar features of the formalism are highlighted, and its applicability for bosons is indicated. The results are employed to develop a more explicit approach describing exchange interactions between spins of Anderson’s magnetic impurities in metals, semiconductors, and insulators. Within the fourth order our theory provides on the equal footing formulae for the Ruderman-Kittel-Kasuya-Yosida, Bloembergen-Rowland, superexchange, and two-electron exchange integrals at non-zero temperature.

I Introduction

Significant effort has been devoted since the 1950s to the development of a quantum perturbation theory of thermodynamic potentials for gases of interacting fermions. The theory, as presented in Refs. 1, 2, 3, 4, provides exact formulae for diagrammatic expansion of thermodynamic potentials in powers of a two-particle interaction potential. The starting point of such a theory is one-particle Hamiltonian, whose eigenstates and eigenergies are known, for example, having been obtained by exact diagonalization of the one-particle Hamiltonian.

In general, however, the diagonal form of one-particle Hamiltonian is unknown, and one has to resort to some perturbation theory to (approximately) diagonalize the Hamiltonian even in non-interacting cases. The standard tool to perform such calculations, the Rayleigh-Schrödinger (RS) quantum perturbation theory, requires significant amounts of algebra at each consecutive order of the expansion. Moreover, in practice, application of the standard RS perturbation theory may be complicated by possible degeneracy in the set of zeroth-order eigenstates, which requires prescribing a gauge.

As we show in this Article, such a procedure is not necessary within the quantum perturbation theory applied to the Landau free energy. The equations we obtain for successive orders of expansion in powers of the one-particle potential are relatively simple, and have a surprising feature which manifests beyond the second order.

We demonstrate how this new formalism can serve to study exchange interactions in solids containing magnetic impurities described by the Anderson Hamiltonian. In particular, we provide general formulae for exchange integrals between two Anderson’s magnetic impurities at non-zero temperature and show that in the leading (fourth) order they describe on the equal footing the superexchange, two-electron, and electron-hole (Bloembergen-Rowland) mechanisms in insulators and, additionally, the Ruderman-Kittel-Kasuya-Yosida coupling in metals and extrinsic semiconductors.

Our paper is organized as follows. We first recall, in Sec. II, the standard 4th order formula for the energy and discuss its shortcomings. Our alternative approach is discussed qualitatively in Sec. III, whereas the formal perturbation theory for the Landau free energy is presented in Sec. IV. Its application for the case of exchange interactions in solids containing Anderson’s magnetic impurities is exposed in Sec. V.

II Standard approach

The standard theory of exchange interactions between magnetic impurity spins in solids, as developed by Larson et al. [5] and Savoyant et al. [6], following the work of Anderson, Falicov and others [7, 8, 9], is based on the following formulation of the fourth-order quantum perturbation theory for fermions: the matrix element of an effective Hamiltonian $H_{\text{eff}}$ between the initial state $\left|i\right>$ and the final state $\left|f\right>$, for a particle subject to a perturbing interaction described by a one-particle operator $V$, is given by the sum of paths over intermediate states $I_{1}$, $I_{2}$, $I_{3}$ (see Eq. 4.3 of Ref. 5),

$$\left<f\middle|H_{\text{eff}}\middle|i\right>=\sum_{I_{1},I_{2},I_{3}}\frac{\left<f\middle|V\middle|I_{1}\right>\left<I_{1}\middle|V\middle|I_{2}\right>\left<I_{2}\middle|V\middle|I_{3}\right>\left<I_{3}\middle|V\middle|i\right>}{(E_{0}-E_{1})(E_{0}-E_{2})(E_{0}-E_{3})}.$$ (1)

Now, to calculate the exchange interaction energy, one assumes as $i$ the state with the first spin up and the second spin down, and as $f$ the state with the first spin down and the second spin up. This approach is valid at zero temperature and breaks down at energy crossings, $E_{0}=E_{i}$. Moreover, it has been noted [10, 11] that taking in Eq. 1 the unperturbed ground state energy as $E_{0}$ is (in general) an inappropriate approximation.

III Alternative approach

Here we develop an alternative more formal and strict approach in which the Landau free energy is determined by the perturbation theory [12]. Our approach is valid at non-zero temperature (as long as the perturbation is smaller than $k_{B}T$) and handles properly divergences associated with energy crossings, $E_{0}=E_{i}$ in denominators of Eq. (1).

Our model system consists of two magnetic impurities with singly occupied $d$ orbitals hybridizing with extended band states via the hybridization operator $V=V_{\text{hyb}}$. We assume the local spins are in spin coherent states[13], parameterized by the spin direction, and calculate the free energy of the electronic subsystem (the occupied band states) as the function of the directions of the localized spins. We need to calculate the fourth order perturbation to the energy of a band state $I_{0}$

$$\sum_{I_{1},I_{2},I_{3}}\frac{\left<I_{0}\middle|V\middle|I_{1}\right>\left<I_{1}\middle|V\middle|I_{2}\right>\left<I_{2}\middle|V\middle|I_{3}\right>\left<I_{3}\middle|V\middle|I_{0}\right>}{(E_{0}-E_{1})(E_{0}-E_{2})(E_{0}-E_{3})},$$ (2)

and to sum over all the band states with the Fermi-Dirac occupation function factor $f(E)$ as follows

$$E^{(4)}=\sum_{I_{0},I_{1},I_{2},I_{3}}\frac{f(E_{0})}{(E_{0}-E_{1})(E_{0}-E_{2})(E_{0}-E_{3})}{}$$ (3)

$$\qquad\qquad\times\left<I_{0}\middle|V\middle|I_{1}\right>\left<I_{1}\middle|V\middle|I_{2}\right>\left<I_{2}\middle|V\middle|I_{3}\right>\left<I_{3}\middle|V\middle|I_{0}\right>,$$

where now the final state coincides with the initial one, and all energies are the unperturbed ones. We have verified directly using a computer algebra system that the last equation is rigorous, with the requirement that the sum over the states $(I_{0},I_{1},I_{2},I_{3})$ has to be extended over all states, also including those combinations of $(I_{0},I_{1},I_{2},I_{3})$ for which the denominator vanishes, either because of the degeneracy of the unperturbed Hamiltonian or (more importantly) due to repetitions in the sequence $(I_{0},I_{1},I_{2},I_{3})$. Indeed, although only $f(E_{0})$ appears explicitly in the equation, if l’Hôpital’s rule is applied in a straight-forward manner [14], the derivatives $f^{\prime}(E_{0})$, $f^{\prime\prime}(E_{0})$ and $f^{\prime\prime\prime}(E_{0})$ appear, as required by the chain rule for the derivatives of $f(E(\lambda))$. We detail the regularization prescription below.

IV Thermodynamic perturbation theory for systems of non-interacting particles

Consider now a system of $n_{f}$ non-interacting fermions numbered $n=1,2,\ldots,n_{f}$, with the corresponding annihilation and creation operators $({\hat{a}}_{n},{\hat{a}}_{n}^{\dagger})$ and eigenenergies $E_{n}^{(0)}$, subject to a one-particle perturbation $\hat{V}=\sum_{mn}V_{mn}\,{\hat{a}}_{m}^{\dagger}{\hat{a}}_{n}$. The Hamiltonian of the system reads

$$\hat{H}(\lambda)=\sum_{n}E_{n}^{(0)}{\hat{a}}_{n}^{\dagger}{\hat{a}}_{n}+\lambda\hat{V}.$$ (4)

The grand thermodynamic potential, called also the Landau free energy, assumes the form

$$\Omega(\lambda)=-k_{B}T\,\mathop{\mathrm{Tr}}\Biggl(\ln\Biggl\{\hat{\mathbbm{1}}{}$$ (5)

$$\qquad\qquad+\exp\left[-\frac{1}{k_{B}T}\left(\hat{H}_{1}(\lambda)-\mu\hat{N}_{1}\right)\right]\Biggr\}\Biggr),$$

where the lower subscript $1$ indicates a one-particle operator, $\mu$ is the chemical potential, $\hat{N}=\sum_{n}{\hat{a}}_{n}^{\dagger}{\hat{a}}_{n}$ is the number of particles operator ($\hat{N}_{1}$ is the identity operator). We expand $\Omega(\lambda)$ in powers of $\lambda$,

$$\Omega(\lambda)=\Omega^{(0)}+\lambda\,\Omega^{(1)}+\lambda^{2}\,\Omega^{(2)}\ldots,$$ (6)

where

$$\Omega^{(0)} = -k_{B}T\sum_{n}\ln\left\{1+\exp\left[-\frac{1}{k_{B}T}\left(E_{n}^{(0)}-\mu\right)\right]\right\},\quad$$ (7)

$$\Omega^{(1)} = \sum_{n_{1}}f\left(E_{n_{1}}^{(0)}\right)V_{n_{1}n_{1}},$$ (8)

$$\Omega^{(2)} = \sum_{n_{1},n_{2}}\frac{f\left(E_{n_{1}}^{(0)}\right)}{E_{n_{1}}^{(0)}-E_{n_{2}}^{(0)}}V_{n_{1}n_{2}}V_{n_{2}n_{1}},$$ (9)

$$\Omega^{(3)} = \sum_{n_{1},n_{2},n_{3}}\frac{f\left(E_{n_{1}}^{(0)}\right)}{\left(E_{n_{1}}^{(0)}-E_{n_{2}}^{(0)}\right)\left(E_{n_{1}}^{(0)}-E_{n_{3}}^{(0)}\right)}$$ (10)

$$\qquad\qquad{}\times V_{n_{1}n_{2}}V_{n_{2}n_{3}}V_{n_{3}n_{1}},$$

$$\Omega^{(4)} = \sum_{n_{1},n_{2},n_{3},n_{4}}$$ (11)

$$\frac{f\left(E_{n_{1}}^{(0)}\right)}{\left(E_{n_{1}}^{(0)}-E_{n_{2}}^{(0)}\right)\left(E_{n_{1}}^{(0)}-E_{n_{3}}^{(0)}\right)\left(E_{n_{1}}^{(0)}-E_{n_{4}}^{(0)}\right)}$$

$$\qquad\qquad{}\times V_{n_{1}n_{2}}V_{n_{2}n_{3}}V_{n_{3}n_{4}}V_{n_{4}n_{1}},$$

with the Fermi-Dirac distribution function

$$f(E)=\frac{1}{1+\exp\left(\frac{E-\mu}{k_{B}T}\right)}.$$ (12)

This result is partly formal, because the denominators may vanish, e.g. they do if there are repetitions in the sequence of summation indices. We will explain the exact meaning of these expressions in such cases below, and give a more explicit formula in terms of a ratio of simple determinants in Appendix B.

By symmetrizing with respect to the summation indices we obtain the same expression in a more standard notation,

$$\Omega^{(2)}=\frac{1}{2}\sum_{m,n}\frac{f\left(E_{m}^{(0)}\right)-f\left(E_{n}^{(0)}\right)}{E_{m}^{(0)}-E_{n}^{(0)}}\left|V_{mn}\right|^{2},$$ (13)

sometimes also written as[15]

$$\Omega^{(2)}=\sum_{m,n}\frac{f\left(E_{m}^{(0)}\right)\left[1-f\left(E_{n}^{(0)}\right)\right]}{E_{m}^{(0)}-E_{n}^{(0)}}\left|V_{mn}\right|^{2}.$$ (14)

However, surprisingly, the higher orders [Eqs. (10) and (11)] lack the “expected” factors of $\prod_{i>1}\left[1-f\left(E_{n_{i}}^{(0)}\right)\right]$ [cf. Eq. (83)], with $i$ corresponding to intermediate states other than the state being perturbed. We interpret this fact as follows: the expectation behind writing $\Omega^{(2)}$ as in Eq. (14) involves the Pauli principle, which would require the intermediate states to be empty to allow virtual transitions to those states. However, since the eigenstates of the full Hamiltonian remain orthogonal on perturbation, an admixture of an occupied eigenstate (which is just other wording for a virtual transition to that state) does not violate the Pauli principle. Moreover, e.g. $\Omega^{(3)}$ given in Eq. (10) obeys a particle-hole symmetry, not shared by its modification [Eq. (83)].

Most importantly, however, various literature theoretical prescriptions (see, e.g., Lewiner, Gaj, and Bastard [15]) do not specify a procedure to handle singular denominators, like those which occur identically for $m=n$ [$(\nu,k)=(\nu^{\prime},k^{\prime})$ in Lewiner’s et al. Eq. 2]. Our observation is that properly regularizing and including such terms in the total r.h.s. is crucial to be able to interpret the result as a perturbation to the free energy of the system. Indeed, the latter includes — via the chain rule — the derivatives $f^{\prime}(E_{n})$, which are properly accounted for by applying l’Hôpital’s rule to the terms in Eq. (13) with $m=n$.

We stress that the above formulae are valid for any degeneracy in the set of energies. However, since the denominators are singular if degeneracies are present, the expressions may need a regularization. For example

$$E_{1},E_{2}\to E\Rightarrow\frac{f\left(E_{1}\right)}{E_{1}-E_{2}}+\frac{f\left(E_{2}\right)}{E_{2}-E_{1}}{}$$ (15)

$$\qquad\qquad=\frac{f\left(E_{1}\right)-f\left(E_{2}\right)}{E_{1}-E_{2}}\to f^{\prime}(E).$$

At any order, prior to regularization the iterated sums in Eqs. (9) to (11) have to be symmetrized with respect to summation indices. That is, one calculates the total of all terms corresponding to combinations $(n_{1},n_{2},\ldots)$ which differ by a permutation. Any permutation is allowed. Indeed, the fraction consisting of the Fermi-Dirac numerator and the energetic denominators is preserved by any permutation which keeps the first index in place, while the remaining product of matrix elements is preserved by a cycle of all the indices (and these two kinds of permutations generate the full symmetric group). Below we show that despite singular denominators in the individual terms, each total is a well defined expression. In the following we assume that $x_{1},x_{2},\ldots\to x$, $y_{1},y_{2},\ldots\to y$, $z_{1},z_{2},\ldots\to z$. Possible combinations of the degeneracy in the set of up to 4 energies are summarized in Table 1, which serves as an index to the relevant equations.

$$\frac{f(x_{1})}{x_{1}-x_{2}}+\frac{f(x_{2})}{x_{2}-x_{1}}\to f^{\prime}(x),$$ (16)

$$\frac{f(y_{1})}{(y_{1}-x_{1})(y_{1}-x_{2})}+\frac{f(x_{1})}{(x_{1}-y_{1})(x_{1}-x_{2})}{}$$ (17)

$$\qquad\qquad\qquad\qquad{}+\frac{f(x_{2})}{(x_{2}-y_{1})(x_{2}-x_{1})}\to\frac{\frac{f(x)-f(y)}{x-y}-f^{\prime}(x)}{y-x},$$

$$\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}+{}$$ (18)

$$\qquad\qquad\qquad\qquad{}+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}\to\frac{f^{\prime\prime}(x)}{2},$$

$$\frac{f(x_{1})}{(x_{1}-y_{1})(x_{1}-x_{2})(x_{1}-z_{1})}+\frac{f(y_{1})}{(y_{1}-x_{1})(y_{1}-x_{2})(y_{1}-z_{1})}+{}$$

$$\qquad\qquad{}+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-y_{1})(x_{2}-z_{1})}+\frac{f(z_{1})}{(z_{1}-x_{1})(z_{1}-y_{1})(z_{1}-x_{2})}\to{}$$ (19)

$$\frac{1}{2}\frac{f(y)-f(z)}{y-z}\left[\frac{1}{(y-x)^{2}}+\frac{1}{(z-x)^{2}}\right]+\frac{1}{(y-x)(z-x)}\Biggl[f^{\prime}(x)+{}\Biggr.$$

$$\qquad\qquad\qquad\qquad\qquad\qquad\Biggl.{}+\left(f(x)-\frac{f(y)+f(z)}{2}\right)\left(\frac{1}{y-x}+\frac{1}{z-x}\right)\Biggr],$$

$$\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-y_{1})(x_{1}-y_{2})}+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-y_{1})(x_{2}-y_{2})}+{}$$

$$\qquad{}+\frac{f(y_{1})}{(y_{1}-x_{1})(y_{1}-x_{2})(y_{1}-y_{2})}+\frac{f(y_{2})}{(y_{2}-x_{1})(y_{2}-x_{2})(y_{2}-y_{1})}\to{}$$ (20)

$$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\frac{[f^{\prime}(x)+f^{\prime}(y)]-2\frac{f(x)-f(y)}{x-y}}{(x-y)^{2}},$$

$$\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-y_{1})}+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})(x_{2}-y_{1})}+{}$$

$$\qquad\qquad{}+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})(x_{3}-y_{1})}+\frac{f(y_{1})}{(y_{1}-x_{1})(y_{1}-x_{2})(y_{1}-x_{3})}\to{}$$ (21)

$$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\frac{\frac{\frac{f(x)-f(y)}{x-y}-f^{\prime}(x)}{y-x}-\frac{f^{\prime\prime}(x)}{2}}{y-x},$$

$$\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})}+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})(x_{2}-x_{4})}+{}$$ (22)

$${}+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})(x_{3}-x_{4})}+\frac{f(x_{4})}{(x_{4}-x_{1})(x_{4}-x_{2})(x_{4}-x_{3})}\to\frac{1}{6}f^{\prime\prime\prime}(x).$$

order	degeneracy	equation number
2	$xx$	(16)
3	$xxy$	(17)
3	$xxx$	(18)
4	$xxyz$	(19)
4	$xxyy$	(20)
4	$xxxy$	(21)
4	$xxxx$	(22)

Table 1: Summary of possible degeneracies in a set of up to four eigenenergies.

To prove Eqs. (8)–(11) we make use of the relation given in Appendix A [Eq. (47)],

$$\frac{1}{2}\mathop{\mathrm{Tr}}(H_{0}+\lambda V)+\Omega(\lambda)=-\frac{1}{\beta}\ln\det\left[-\frac{\partial}{\partial\tau}-(H_{0}+\lambda V)\right],$$ (23)

where $H_{0}$ is the unperturbed Hamiltonian and $V$ is the perturbation (as usual, $\beta=1/k_{B}T$ and we drop the constant factor $\det\left[-\frac{\partial}{\partial\tau}\right]$ which affects only the normalization of the partition function). Let

$$\mathcal{G}_{0}^{-1}=-\frac{\partial}{\partial\tau}-H_{0},\qquad\delta\mathcal{G}^{-1}=-\lambda V.$$ (24)

We rewrite the above relation as

$$\frac{1}{2}\mathop{\mathrm{Tr}}(H_{0}+\lambda V)+\Omega(\lambda)=-\frac{1}{\beta}\mathop{\mathrm{Tr}}\left\{\ln\left[\mathcal{G}_{0}^{-1}+\mathcal{G}^{-1}\right]\right\},$$ (25)

which can be easily expanded into power series [16]:

$$\frac{1}{2}\mathop{\mathrm{Tr}}(H_{0}+\lambda V)+\Omega(\lambda)$$ (26)

$$= \left[\frac{1}{2}\mathop{\mathrm{Tr}}(H_{0})+\Omega(0)\right]+\mathop{\mathrm{Tr}}\left[\frac{(-1)^{n}}{n}\left(\mathcal{G}_{0}\mathcal{G}^{-1}\right)^{n}\right].$$

In the spectral representation,

$$\mathop{\mathrm{Tr}}\left[\left(\mathcal{G}_{0}\mathcal{G}^{-1}\right)^{n}\right]=\sum_{m=-\infty}^{\infty}\sum_{i_{1}}\sum_{i_{2}}\ldots\sum_{i_{n}}$$ (27)

$$(-\lambda)^{n}\left[\prod_{a=1}^{n}\frac{1}{i\omega_{m}-(E_{i_{a}}-\mu)}\right]V_{i_{1}i_{2}}V_{i_{2}i_{3}}\ldots V_{i_{n}i_{1}},$$

where $\omega_{m}=(2m+1)\pi/\beta$ is the fermionic Matsubara frequency, and $i_{1},i_{2},\ldots,i_{n}$ label eigenstates of the unperturbed Hamiltonian. The quantity (double braces indicate a multiset)

$$s_{n}(\left\{\left\{E_{i_{1}},E_{i_{2}},\ldots,E_{i_{n}}\right\}\right\})$$ (28)

$$= \frac{1}{\beta}\sum_{m=-\infty}^{\infty}\prod_{a=1}^{n}\frac{1}{i\omega_{m}-(E_{i_{a}}-\mu)},$$

can be calculated in the standard manner as a contour integral [17] over $z=i\omega_{m}$, where the integrand is the summand of Eq. (28) multiplied by the occupation function $f(z)$. Assuming no degeneracy, each residue at $E_{i_{a}}\in\left\{\left\{E_{i_{1}},E_{i_{2}},\ldots,E_{i_{n}}\right\}\right\}$ contributes to the value of $s_{n}$ [Eq. (28)] the term

$$f(E_{i_{a}})\mathop{\prod\limits_{b=1}^{n}}\limits_{b\neq a}\frac{1}{E_{i_{a}}-E_{i_{b}}},$$ (29)

[cf. Eqs. (8)–(11)], while for $n=1$ the contour at infinity yields a contribution of $1/2$ which cancels the term $\frac{1}{2}\mathop{\mathrm{Tr}}(\lambda V)$ on the l.h.s. of Eq. (26).

If the set of energies is degenerate, e.g. $E_{i_{a}}=E_{i_{b}}$, higher order singularities are present in the integrand. Either the residue at a higher order singularity can be obtained as a derivative of the regular part of the integrand at $z=E_{i_{a}}=E_{i_{b}}$, or the limit of $s_{n}$ can be taken as $E_{i_{b}}\to E_{i_{a}}$. If a given eigenstate of the unperturbed Hamiltonian appears multiple times, e.g. $i_{a}=i_{b}$, then $E_{i_{a}}=E_{i_{b}}$, as if a degeneracy was present. In such a case one can either use the former prescription, or take the limit of $s_{n}$ as $E_{i_{b}}\to E_{i_{a}}$, paying careful attention to consider the symbols $E_{i_{a}}$ and $E_{i_{b}}$ in Eq. (28) as independent variables. This ends the proof.

Equations (16)–(22) can be generalized to higher orders by using the recursion relation given in Appendix B [Eq. (54)].

We have verified in the second order of perturbation that our result is valid for bosons, with the Fermi-Dirac occupation function replaced by its bosonic counterpart.

A relation of our results to the Rayleigh-Schrödinger perturbation theory and to the path integral formulation is discussed in Appendices C and D, respectively.

V Spin interactions in solids

Let us now return to the specific case of a spin-spin interaction in a solid containing magnetic impurities described by the Anderson Hamiltonian. The Hilbert space of our model system is a direct sum of the sector of localized electronic states (for example $3d$ states of a magnetic dopant) and the continuum (band states). The localized states are spin-degenerate. This describes the situation where a magnetic impurity (such as Mn) is present in a solid matrix (specifically a semiconductor such as CdTe or HgTe). Hence, the Hamiltonian (the Anderson Hamiltonian for multiple impurities),

$$H=\sum_{k}E_{k}\,a_{k}^{\dagger}a_{k}+\sum_{l}H_{l}+\lambda(V_{hyb}+V_{\text{hyb}}^{\dagger}),$$ (30)

where $\{E_{k}\}$ is the continuous spectrum, and the Hamiltonian of an impurity with label $l\in L$, $H_{l}$, is of the form

$$H_{l}=E_{d}\,(a_{l\uparrow}^{\dagger}a_{l\uparrow}+a_{l\downarrow}^{\dagger}a_{l\downarrow})+U\,a_{l\uparrow}^{\dagger}a_{l\uparrow}a_{l\downarrow}^{\dagger}a_{l\downarrow}.$$ (31)

The one-particle operator $V_{\text{hyb}}$ describes virtual transitions from the extended to the localized states. The term $V_{\text{hyb}}+V_{\text{hyb}}^{\dagger}$ is considered as a perturbation to the electronic Hamiltonian. Hence, our perturbation theory formalism applies. One easily observes that a process corresponding to a given set of states $(n_{1},n_{2},\ldots)$ can contribute to a free energy perturbation if the states $(n_{1},n_{2},\ldots)$ are alternately extended and localized, and that the leading order in which a spin-spin interaction occurs is the fourth one. Moreover, in the leading contributing processes, the labels $(n_{1},n_{2},n_{3},n_{4})$ correspond to two $d$ states $(d_{1},d_{2})$ associated with two different impurities and two band states $k$ and $k^{\prime}$, thus the only possible repetition is $k=k^{\prime}$. (This observation stays behind the success of the standard approach to spin-spin interactions in solids and demonstrates that the standard approach may fail at higher orders). We denote by $E_{k}$ and $E_{k^{\prime}}$ the energies of the band states, by $E_{d}$ the energy of the singly-occupied localized $d$ state, and by $U$ the Coulomb energy of a doubly-occupied $d$ state (this energy adds to $E_{d})$. Our goal is calculate the exchange integral describing the interaction between two localized spins. Therefore, we assume the localized states are occupied by electrons with spin either up or down, and consider the spins of those electrons as the interacting objects, and the remaining electrons as a non-interacting gas. Thus, the present approach neglects the potential exchange within the carriers, which augments the ferromagnetic portion of RKKY coupling [18], as well as between carriers and localized spins, in particular the intra-ion $sp-d(f)$ exchange that may control the strength of spin-spin interactions between rare earth impurities [19], for which $V_{\text{hyb}}$ is small.

We calculate the free energy of the non-interacting electronic gas as the total of the free energy of the electronic states with spin up and those with spin down (this a further simplification which is valid if no spin-orbit interaction is present). If the directions of the local spins agree, and the direction of the electronic spin is the same, we use Eq. (19) with $x=E_{d}$ (the singly-occupied $d$ state — double occupation is prohibited by the Pauli principle), $y=E_{k}$ and $z=E_{k^{\prime}}$ (the two intermediate band states). We do so, because since the dopants are identical, the $d$ states are degenerate. If the spin of the electron is opposite to the direction of the local spins, $x=E_{d}+U$ (the unoccupied $d$ state energy, including the Coulomb interaction with the electron in the $d$ state which is considered as the local spin). However, if the directions of the local spins are opposite, the degeneracy of the two intermediate states is lifted, and Eq. (11) can be used directly with the energies $E_{k}$, $E_{k^{\prime}}$, $E_{d}$ and $E_{d}+U$. We calculate the exchange interaction energy as the difference of the electronic free energies for agreeing and opposite directions of the local spins. Thus the Fermi-Dirac occupation function and the energetic denominators contribute the factor (let us denote it $A_{k,k^{\prime}}$):

$$A_{k,k^{\prime}} = \frac{1}{U}\left[\frac{2E_{d}^{2}+2E_{k}E_{k^{\prime}}+(E_{k}+E_{k^{\prime}})U-2E_{d}(E_{k}+E_{k^{\prime}}+U)}{(E_{d}-E_{k})^{2}(E_{d}-E_{k^{\prime}})^{2}}f\left(\frac{E_{d}-\mu}{k_{B}T}\right)+{}\right.$$ (32)

$$\left.{}+\frac{-2(E_{d}-E_{k})(E_{d}-E_{k^{\prime}})+3(E_{k}+E_{k^{\prime}}-2E_{d})U-4U^{2}}{(E_{d}+U-E_{k})^{2}(E_{d}+U-E_{k^{\prime}})^{2}}f\left(\frac{E_{d}+U-\mu}{k_{B}T}\right)\right]+{}$$

$${}+U^{2}\left[\frac{1}{(E_{d}-E_{k})^{2}(E_{d}-E_{k}+U)^{2}}\frac{f\left(\frac{E_{k}-\mu}{k_{B}T}\right)}{E_{k}-E_{k^{\prime}}}+\frac{1}{(E_{d}-E_{k^{\prime}})^{2}(E_{d}-E_{k^{\prime}}+U)^{2}}\frac{f\left(\frac{E_{k^{\prime}}-\mu}{k_{B}T}\right)}{E_{k^{\prime}}-E_{k}}\right]+{}$$

$${}+\frac{1}{k_{B}T}\left[\frac{1}{(E_{d}-E_{k})(E_{d}-E_{k^{\prime}})}f^{\prime}\left(\frac{E_{d}-\mu}{k_{B}T}\right)+\frac{1}{(E_{d}+U-E_{k})(E_{d}+U-E_{k^{\prime}})}f^{\prime}\left(\frac{E_{d}+U-\mu}{k_{B}T}\right)\right].$$

The final expression for the difference between the free energy of the carriers for parallel local spins [$\Omega_{c}(\uparrow\uparrow)$] and for antiparallel local spins [$\Omega_{c}(\uparrow\downarrow)$] is proportional to the matrix element of $V_{\text{hyb}}$ in the fourth power:

$$\Omega_{c}(\uparrow\uparrow)-\Omega_{c}(\uparrow\downarrow) = \sum_{k,k^{\prime}}A_{k,k^{\prime}}\left<k\middle|V_{\text{hyb}}^{\dagger}\middle|d_{1}\right>\left<d_{1}\middle|V_{\text{hyb}}\middle|k^{\prime}\right>\left<k^{\prime}\middle|V_{\text{hyb}}^{\dagger}\middle|d_{2}\right>\left<d_{2}\middle|V_{\text{hyb}}\middle|k\right>$$ (33)

$$= \frac{1}{2}\sum_{k,k^{\prime}}A_{k,k^{\prime}}\biggl(\left<k\uparrow\middle|V_{\text{hyb}}^{\dagger}\middle|d_{1}\uparrow\right>\left<d_{1}\uparrow\middle|V_{\text{hyb}}\middle|k^{\prime}\uparrow\right>\left<k^{\prime}\uparrow\middle|V_{\text{hyb}}^{\dagger}\middle|d_{2}\uparrow\right>\left<d_{2}\uparrow\middle|V_{\text{hyb}}\middle|k\uparrow\right>+{}$$

$$\qquad{}+\left<k\downarrow\middle|V_{\text{hyb}}^{\dagger}\middle|d_{1}\downarrow\right>\left<d_{1}\downarrow\middle|V_{\text{hyb}}\middle|k^{\prime}\downarrow\right>\left<k^{\prime}\downarrow\middle|V_{\text{hyb}}^{\dagger}\middle|d_{2}\downarrow\right>\left<d_{2}\downarrow\middle|V_{\text{hyb}}\middle|k\downarrow\right>\biggr).$$ (34)

We have used the fact that in the absence of a spin-orbit interaction the matrix elements factorize into the product of the orbital and the spin part,

$$\left<d_{i}\uparrow\middle|V_{\text{hyb}}\middle|k\uparrow\right>=\left<d_{i}\downarrow\middle|V_{\text{hyb}}\middle|k\downarrow\right>=\left<d_{i}\middle|V_{\text{hyb}}\middle|k\right>,$$ (35)

$$\left<d_{i}\uparrow\middle|V_{\text{hyb}}\middle|k\downarrow\right>=\left<d_{i}\downarrow\middle|V_{\text{hyb}}\middle|k\uparrow\right>=0,$$ (36)

and only the orbital parts $\left<d_{i}\middle|V_{\text{hyb}}\middle|k\right>$ enter the equation for $\Omega_{c}(\uparrow\uparrow)-\Omega_{c}(\uparrow\downarrow)$.

Hence, the non-interacting approximation to the classical exchange integral between the two spins is given by Eqs. (32) and (34) via

$$J_{12}\approx-\Bigl[\Omega_{c}(\uparrow\uparrow)-\Omega_{c}(\uparrow\downarrow)\Bigr].$$ (37)

Then, the spin-dependent part of the effective Hamiltonian for two spins can be reconstructed from its spin coherent state[13] representation as

$$\hat{H}^{(4)}_{\text{eff}}=-2J_{12}\,\hat{S}_{1}\cdot\hat{S}_{2},$$ (38)

with spin $S=1/2$ angular momentum operators $\hat{S}_{1}$ and $\hat{S}_{2}$ representing the degrees of freedom of the two spins.

In a crystalline solid, the matrix elements include the plane-wave phase factor:

$$\left<d_{i}\middle|V_{\text{hyb}}\middle|k\right>=\exp(-i\kappa\cdot R_{i})\tilde{V}_{k},$$ (39)

with $R_{i}$ being the real-space position of the impurity, and $\kappa$ the wave vector corresponding to the band state $k$. Therefore, the product of matrix elements is proportional to $\exp[-i(\kappa-\kappa^{\prime})\cdot(R_{2}-R_{1})]$, and we obtain

$$\Omega_{c}(\uparrow\uparrow)-\Omega_{c}(\uparrow\downarrow)=\sum_{k,k^{\prime}}{}$$ (40)

$$\quad{}\exp[-i(\kappa-\kappa^{\prime})\cdot(R_{2}-R_{1})]\left|\tilde{V}_{k}\right|^{2}\left|\tilde{V}_{k^{\prime}}\right|^{2}A_{k,k^{\prime}}.$$

Here, $\exp[-i(\kappa-\kappa^{\prime})\cdot(R_{2}-R_{1})]$ can be replaced with $\cos[(\kappa-\kappa^{\prime})\cdot(R_{2}-R_{1})]$ by the possibility of symmetrization with respect to the interchange $k\leftrightarrow k^{\prime}$. In general, if a spin orbit interaction is present, and the spin directions are not collinear, this symmetrization must be suppressed.

Additional degeneracies may be present. In particular, in the limit $E_{k^{\prime}}\to E_{k}$ we obtain:

$$A_{k,k}=\frac{2U^{2}(2E_{d}-2E_{k}+U)}{(E_{d}-E_{k})^{3}(E_{d}+U-E_{k})^{3}}f\left(\frac{E_{k}-\mu}{k_{B}T}\right){}$$ (41)

$${}+\frac{1}{k_{B}T}\frac{U^{2}}{(E_{d}-E_{k})^{2}(E_{d}+U-E_{k})^{2}}f^{\prime}\left(\frac{E_{k}-\mu}{k_{B}T}\right)$$

$${}+\ldots$$

The second term is proportional to $\frac{1}{k_{B}T}f^{\prime}$, i.e. to the Dirac delta at the Fermi level in the zero-temperature limit. In insulators (no density of states at the Fermi level), its contribution vanishes as $T\to 0$. However, this is not the case in metals, where this term represents the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [20]. In contrast, Eq. 4.5 of Ref. 5 fails even to produce a meaningful (finite) result for $\varepsilon_{n}(k)=\varepsilon_{n^{\prime}}(k^{\prime})$ ($E_{k}=E_{k^{\prime}}$).

Although the resulting expression for $A_{k,k^{\prime}}$ is valid for arbitrary temperature, we take the zero temperature limit in order to single out various contributions discussed in the literature. Furthermore, we assume the case of an insulator, which allows us to omit terms proportional to the derivatives of the Fermi-Dirac occupation function.

Following literature conventions, we decompose our zero-temperature exchange integrals into three terms, $J_{dd}=J_{hh}+J_{eh}+J_{ee}$. The first term is called superexchange and corresponds to occupied intermediate band states. The last term, called the two-electron term, corresponds to empty intermediate band states. The remaining term, called the electron-hole term, turns out in fact to comprise two contributions ($eh$ and $he$), and describes the Bloembergen-Rowland interaction [21]. The total numerical factor resulting from the energy denominators is

$$A_{k,k^{\prime}} \to \Theta(\mu-E_{k})\Theta(\mu-E_{k^{\prime}})A^{hh}$$ (42)

$${}+\Theta(\mu-E_{k})\Theta(E_{k^{\prime}}-\mu)A^{he}$$

$${}+\Theta(E_{k}-\mu)\Theta(\mu-E_{k^{\prime}})A^{eh}$$

$${}+\Theta(E_{k}-\mu)\Theta(E_{k^{\prime}}-\mu)A^{ee}.$$

where the individual terms are given below ($A^{he}$ corresponds to an occupied state with energy $E_{k}$, while $A^{eh}$ to an occupied state with energy $E_{k^{\prime}}$).

$$A^{hh}=\frac{2}{(E_{d}+U-E_{k})(E_{d}+U-E_{k^{\prime}})U}+\frac{1}{(E_{d}+U-E_{k})(E_{d}+U-E_{k^{\prime}})}\left(\frac{1}{E_{d}+U-E_{k}}+\frac{1}{E_{d}+U-E_{k^{\prime}}}\right),$$ (43)

$$A^{he}=\frac{2}{(E_{d}+U-E_{k})(E_{d}-E_{k^{\prime}})U}+\frac{1}{E_{k}-E_{k^{\prime}}}\left(\frac{1}{E_{d}-E_{k^{\prime}}}-\frac{1}{E_{d}+U-E_{k}}\right)^{2},$$ (44)

$$A^{eh}=\frac{2}{(E_{d}-E_{k})(E_{d}+U-E_{k^{\prime}})U}+\frac{1}{E_{k^{\prime}}-E_{k}}\left(\frac{1}{E_{d}-E_{k}}-\frac{1}{E_{d}+U-E_{k^{\prime}}}\right)^{2},$$ (45)

$$A^{ee}=\frac{2}{(E_{d}-E_{k})(E_{d}-E_{k^{\prime}})U}-\frac{1}{(E_{d}-E_{k})(E_{d}-E_{k^{\prime}})}\left(\frac{1}{E_{d}-E_{k}}+\frac{1}{E_{d}-E_{k^{\prime}}}\right).$$ (46)

We underline that Eqs. (43)–(42) do not take into account the RKKY term [see Eq. (41)] included in Eq. (32).

The series of Eqs. (43)–(46) lead to exchange integrals in agreement with formulae obtained by Blinowski, Kacman, and Majewski for $d^{5}$ transition-metal ions in insulators [22]. This agreement reconfirms the fact that luckily the regularization procedure is not necessary in this case. These formulae were later modified for the $d^{4}$ configuration [23]. The theoretical model for $d^{4}$ ions was successfully employed to describe ferromagnetic ordering of Mn spins observed experimentally in insulating $(\mathrm{Ga},\mathrm{Mn})\mathrm{N}$ with $\mathrm{Mn}^{3+}$ ions [24, 25]. This agreement between theoretical and experimental results indicates that inclusion of higher order terms is not necessary, at least when the short range superexchange dominates, the case of $(\mathrm{Ga},\mathrm{Mn})\mathrm{N}$.

The key accomplishment of this section is presented in Eq. 32 that allows describing various contributions to the exchange coupling of magnetic impurity pairs on an equal footing and at arbitrary temperature. However, in the case of insulators the low temperature approximation usually holds, since the energy distance of the Fermi level to both band edges and $d$ states is significantly larger than $k_{B}T$. A similar approximation is metals is valid (as long as the RKKY term is being kept) as the condition of a strong degeneracy of the carrier gas is fulfilled, $\mu\gg k_{B}T$. In contrast to metals, however, this condition is often not satisfied in extrinsic semiconductors [26], and indeed non-standard temperature effects in the latter have been noted [27].

VI Conclusions

Landau’s thermodynamic perturbation theory [12] is a simple-to-use, general, strict, formal, systematic and elegant tool. Here, it has been advanced to arbitrary order in the specific case of non-interacting fermionic systems (its extension to non-interacting bosonic systems appears as straightforward). The developed approach has been applied to the case of exchange interactions between Anderson’s magnetic impurities in solids. In the fourth order our results generalize literature formulae to non-zero temperatures and handle properly zeros of the energetic denominators due to repeated level indices and possible overlaps or crossings of energy levels. In particular, the present results apply rigorously to situations, in which band carriers mediate a long-ranged interaction between localized spins, the case of dilute magnetic metals and extrinsic dilute magnetic semiconductors.

VII Acknowledgments

The work has been supported within the Master project of National Center of Science in Poland (2011/02/A/ST3/00125). The International Centre for Interfacing Magnetism and Superconductivity with Topological Matter project is carried out within the International Research Agendas programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund.

Appendix A Partition function as a functional determinant

The partition function for fermions can be represented as a functional determinant in the imaginary time formalism,

$$1+\exp(-\beta E)=\exp(-\beta E/2)\frac{\det\left(-\frac{\partial}{\partial\tau}-E\right)}{\det\left(-\frac{\partial}{\partial\tau}\right)}.$$ (47)

Indeed, the ratio on the r.h.s. of this relation can be represented as follows and making use of the Euler representation for the meromorphic function $\sin(z)$ we obtain:

$$\frac{\det\left(-\frac{\partial}{\partial\tau}-E\right)}{\det\left(-\frac{\partial}{\partial\tau}\right)}=\frac{\prod_{n=-\infty}^{\infty}\frac{(2n+1)i\pi}{\beta}-E}{\prod_{n=-\infty}^{\infty}\frac{(2n+1)i\pi}{\beta}}$$ (48)

$$= \prod_{n=-\infty}^{\infty}\left(1+i\frac{\beta E}{(2n+1)\pi}\right)$$ (49)

$$= \prod_{n=0}^{\infty}\left[1+\left(\frac{\beta E}{(2n+1)\pi}\right)^{2}\right]$$ (50)

$$= \frac{\prod_{n=1}^{\infty}\left[1+\left(\frac{\beta E}{n\pi}\right)^{2}\right]}{\prod_{n=1}^{\infty}\left[1+\left(\frac{\beta E}{2n\pi}\right)^{2}\right]}$$ (51)

$$= \frac{\sinh(\beta E)}{\sinh(\beta E/2)}=2\cosh(\beta E/2)$$ (52)

$$= \exp(\beta E/2)[1+\exp(-\beta E)].$$ (53)

Appendix B Recursion relations

In this Appendix we show that the quantity $s_{n}(\left\{\left\{E_{i_{1}},E_{i_{2}},\ldots,E_{i_{n}}\right\}\right\})$ given by Eq. (28) can be obtained recursively. The recursion has the following desirable feature: at each step of recursion, at most one singular denominator occurs, which allows to apply directly l’Hôpital’s rule. The recursion provides an expression for $s_{n}$ for the multiset of energies $\epsilon=\left\{\left\{E_{i_{1}},E_{i_{2}},\ldots,E_{i_{n}}\right\}\right\}$ in terms of $s_{n-1}$:

$$s_{n}(\epsilon)=\frac{1}{n(n-1)}\sum_{a,b=1,2,\ldots,n\atop a\neq b}{}$$ (54)

$$\qquad\frac{s_{n-1}(\epsilon\setminus\{\{E_{i_{a}}\}\})-s_{n-1}(\epsilon\setminus\{\{E_{i_{b}}\}\})}{E_{i_{b}}-E_{i_{a}}}.$$

The summand has the form $\frac{g(E_{i_{b}})-g(E_{i_{a}})}{E_{i_{b}}-E_{i_{a}}}$, with $g(y)=s_{n-1}(\epsilon\setminus\{\{E_{i_{a}},E_{i_{b}}\}\}\cup\{\{y\}\})$.

Alternatively, the following representation in terms of the Vandermonde determinants is possible:

$$s_{n}(\left\{\left\{E_{1},E_{2},\ldots,E_{n}\right\}\right\})={}$$

$$\left|\begin{array}[]{cccc}1&1&\ldots&1\\ E_{1}&E_{2}&\ldots&E_{n}\\ (E_{1})^{2}&(E_{2})^{2}&\ldots&(E_{n})^{2}\\ \vdots&\vdots&\cdots&\vdots\\ (E_{1})^{n-2}&(E_{2})^{n-2}&\ldots&(E_{n})^{n-2}\\ f(E_{1})&f(E_{2})&\ldots&f(E_{n})\end{array}\right|\Biggr/$$

$$\qquad{}\left|\begin{array}[]{cccc}1&1&\ldots&1\\ E_{1}&E_{2}&\ldots&E_{n}\\ (E_{1})^{2}&(E_{2})^{2}&\ldots&(E_{n})^{2}\\ \vdots&\vdots&\cdots&\vdots\\ (E_{1})^{n-1}&(E_{2})^{n-1}&\ldots&(E_{n})^{n-1}\end{array}\right|.$$

If a degeneracy is present in the set of energies, a repeated column should be replaced by its consecutive derivatives with respect to the corresponding energy, identically in the numerator and the denominator. For example, the expression (19) can be written as

$$\left.\left|\begin{array}[]{cccc}1&0&1&1\\ x&1&y&z\\ x^{2}&2x&y^{2}&z^{2}\\ f(x)&f^{\prime}(x)&f(y)&f(z)\end{array}\right|\Biggm/\left|\begin{array}[]{cccc}1&0&1&1\\ x&1&y&z\\ x^{2}&2x&y^{2}&z^{2}\\ x^{3}&3x^{2}&y^{3}&z^{3}\end{array}\right|\right.,$$ (68)

and (20) as

$$\left.\left|\begin{array}[]{cccc}1&0&1&0\\ x&1&y&1\\ x^{2}&2x&y^{2}&2y\\ f(x)&f^{\prime}(x)&f(y)&f^{\prime}(y)\end{array}\right|\Biggm/\left|\begin{array}[]{cccc}1&0&1&0\\ x&1&y&1\\ x^{2}&2x&y^{2}&2y\\ x^{3}&3x^{2}&y^{3}&3y^{2}\end{array}\right|\right..$$ (69)

Appendix C Relation to Rayleigh-Schrödinger perturbation theory

We will now demonstrate how the third order perturbation can be derived using the standard Rayleigh-Schrödinger (RS) perturbation theory. The consecutive orders of expansion are,

$$E^{(1)}_{n}=\left<n\middle|V\middle|n\right>,\quad E^{(2)}_{n}=\sum_{k\neq n}\frac{\left|\left<k\middle|V\middle|n\right>\right|^{2}}{E_{n}-E_{k}},$$ (70)

$$E^{(3)}_{n}=\sum_{k\neq n}\sum_{m\neq n}\frac{\left<n\middle|V\middle|m\right>\left<m\middle|V\middle|k\right>\left<k\middle|V\middle|n\right>}{(E_{n}-E_{m})(E_{n}-E_{k})}-\left<n\middle|V\middle|n\right>\sum_{m\neq n}\frac{\left|\left<n\middle|V\middle|m\right>\right|^{2}}{(E_{n}-E_{m})^{2}}.$$ (71)

We see a subtracted term on the r.h.s. of the equation for $E^{(3)}_{n}$. To expand $\Omega(\lambda)$ we need only the energies:

$$\Omega(\lambda)=-\frac{1}{\beta}\sum_{n}\ln\left[1+\exp(-\beta E_{n}(\lambda))\right].$$ (72)

Indeed,

$$\Omega^{(3)}=\sum_{n}\left[f(E_{n})E^{(3)}_{n}+f^{\prime}(E_{n})E^{(1)}_{n}E^{(2)}_{n}+\frac{1}{6}f^{\prime\prime}(E_{n})\left(E^{(1)}_{n}\right)^{3}\right].$$ (73)

Direct substitution yields:

$$\Omega^{(3)}=\left[\sum_{k\neq n\atop m\neq n}\frac{f(E_{n})}{(E_{n}-E_{m})(E_{n}-E_{k})}\left<n\middle|V\middle|m\right>\left<m\middle|V\middle|k\right>\left<k\middle|V\middle|n\right>+{}\right.$$ (74)

$${}\left.-\sum_{m\neq n}\frac{f(E_{n})}{(E_{n}-E_{m})^{2}}\left<n\middle|V\middle|n\right>\left|\left<n\middle|V\middle|m\right>\right|^{2}+{}\right.$$

$$\left.{}+\sum_{m\neq n}\frac{f^{\prime}(E_{n})}{E_{n}-E_{m}}\left<n\middle|V\middle|n\right>\left|\left<m\middle|V\middle|n\right>\right|^{2}+\sum_{n}\frac{1}{6}f^{\prime\prime}(E_{n})\left(\left<n\middle|V\middle|n\right>\right)^{3}\right].$$

We regroup the terms as follows:

$$\Omega^{(3)}=\frac{1}{3}\left\{\sum_{k\neq m\neq n}\left[\frac{f(E_{k})}{(E_{k}-E_{m})(E_{k}-E_{n})}+\frac{f(E_{m})}{(E_{m}-E_{k})(E_{m}-E_{n})}+\frac{f(E_{n})}{(E_{n}-E_{k})(E_{n}-E_{m})}\right]\times\right.$$ (75)

$${}\qquad\qquad\qquad\qquad\qquad\qquad\times\left<n\middle|V\middle|m\right>\left<m\middle|V\middle|k\right>\left<k\middle|V\middle|n\right>+{}$$

$${}+\sum_{m\neq n}\left[\frac{f(E_{m})-f(E_{n})}{(E_{n}-E_{m})^{2}}+\frac{f^{\prime}(E_{n})}{E_{n}-E_{m}}\right]\left[\left<n\middle|V\middle|n\right>\left<n\middle|V\middle|m\right>\left<m\middle|V\middle|n\right>+{}\right.$$

$$\qquad\qquad\qquad\qquad\qquad\left.{}+\left<m\middle|V\middle|n\right>\left<n\middle|V\middle|n\right>\left<n\middle|V\middle|m\right>+\left<n\middle|V\middle|m\right>\left<m\middle|V\middle|n\right>\left<n\middle|V\middle|n\right>\right]+{}$$

$$\left.{}+\sum_{n}\frac{f^{\prime\prime}(E_{n})}{2}\left<n\middle|V\middle|n\right>\left<n\middle|V\middle|n\right>\left<n\middle|V\middle|n\right>\right\}.$$

This is the required form compatible with Eq. (10). Indeed, the first term is the direct symmetrization of Eq. (10). However, if $k$ or $m$ becomes equal $n$, one uses Eq. (17). Finally, if $k=m=n$, Eq. (18) applies. The purpose of repeating some terms is to make the combinatorial weights evident. We note that one can symmetrize either the part comprising the occupation function and the energetic denominators, or the product of matrix elements (or both).

Appendix D Relation to path integrals

In this Appendix we relate our observation to the properties of path integrals. The partition function $Z(\lambda)$ for a system of fermions can be written as a coherent state path integral,

$$Z(\lambda)=\mathop{\mathrm{Tr}}e^{-\beta H}=\int\mathcal{D}[\xi_{\alpha}(\tau),\overline{\xi}_{\alpha}(\tau)]\,e^{\mathcal{S}[\xi_{\alpha}(\tau),\overline{\xi}_{\alpha}(\tau)]},$$ (76)

where $\alpha$ labels fermionic (Grassmannian) degrees of freedom $\xi_{\alpha}$, $0<\tau<\beta$ is the imaginary time, and $\mathcal{S}$ is the usual action

$$\mathcal{S}[\xi_{\alpha}(\tau),\overline{\xi}_{\alpha}(\tau)]{}$$ (77)

$$=\int_{0}^{\beta}d\tau\,\left(-\overline{\xi}_{\alpha}(\tau)\frac{\partial\xi_{\alpha}(\tau)}{\partial\tau}-H\left(\xi_{\alpha}(\tau),\overline{\xi}_{\alpha}(\tau)\right)\right).\qquad$$

In simple terms, this is to be understood as follows: the integration domain $[0,\beta]$ is divided into $N$ equal intervals, each $\epsilon=\beta/N$ long. Then $\mathop{\mathrm{Tr}}e^{-\beta H}=\mathop{\mathrm{Tr}}\left[e^{-\epsilon H}e^{-\epsilon H}\ldots e^{-\epsilon H}\right]=\mathop{\mathrm{Tr}}\left[e^{-\epsilon H}\mathbbm{1}e^{-\epsilon H}\mathbbm{1}\ldots e^{-\epsilon H}\mathbbm{1}\right]$, where $\mathbbm{1}$ is a resolution of unity, $\mathbbm{1}=\int d\mu(z)\left|z\right>\left<z\right|$, and the trace turns out to be an $N$-dimensional integral (in the limit $N\to\infty$, a path integral). Assume $H=H_{0}+\lambda V$. We use the Trotter formula, $e^{-\epsilon H}\approx e^{-\epsilon H_{0}}e^{-\epsilon\lambda V}$, and approximate $e^{-\epsilon\lambda V}\approx 1-\epsilon\lambda V$. This leads to an integral representation for the terms of the expansion $Z(\lambda)=Z^{(0)}+\lambda Z^{(1)}+\lambda^{2}Z^{(2)}+\ldots$ (from which an expansion of the free energy $F=-\frac{1}{\beta}\ln Z$ can be directly obtained through series composition):

$$Z^{(0)} = \mathop{\mathrm{Tr}}e^{-\beta H_{0}};$$ (78)

$$Z^{(1)} = \int_{\tau_{1}=0}^{\beta}d\tau_{1}\,\mathop{\mathrm{Tr}}\left[e^{-\tau_{1}H_{0}}(-V)e^{-(\beta-\tau_{1})H_{0}}\right];$$ (79)

$$Z^{(2)} = \int_{\tau_{1}=0}^{\beta}d\tau_{1}\int_{\tau_{2}=\tau_{1}}^{\beta}d\tau_{2}\mathop{\mathrm{Tr}}\Bigl[e^{-\tau_{1}H_{0}}(-V)e^{-(\tau_{2}-\tau_{1})H_{0}}{}$$ (80)

$$\qquad\times(-V)e^{-(\beta-\tau_{2})H_{0}}\Bigr],$$

and $Z^{(3)}$ is given by a triple integral

$$Z^{(3)}=\int_{\tau_{1}=0}^{\beta}d\tau_{1}\int_{\tau_{2}=\tau_{1}}^{\beta}d\tau_{2}\int_{\tau_{3}=\tau_{2}}^{\beta}d\tau_{3}\,\mathop{\mathrm{Tr}}\left[e^{-\tau_{1}H_{0}}(-V){}\right.$$ (81)

$$\times\left.e^{-(\tau_{2}-\tau_{1})H_{0}}(-V)e^{-(\tau_{3}-\tau_{2})H_{0}}(-V)e^{-(\beta-\tau_{3})H_{0}}\right].\qquad$$

At this point we can substitute a diagonal non-interacting $H_{0}$ and a one-particle fermionic operator $V$ into these expressions. We obtain (using a computer algebra system), in the generic case $k\neq m\neq n$, terms in $\Omega^{(3)}$ of the form

$$\left[\frac{f(E_{k})}{(E_{k}-E_{m})(E_{k}-E_{n})}+\frac{f(E_{m})}{(E_{m}-E_{k})(E_{m}-E_{n})}+{}\right.$$ (82)

$$\left.{}+\frac{f(E_{n})}{(E_{n}-E_{k})(E_{n}-E_{m})}\right]{}$$

$$\qquad\qquad\qquad\qquad\times\left<k\middle|V\middle|m\right>\left<m\middle|V\middle|n\right>\left<n\middle|V\middle|k\right>,\quad$$

rather than

$$\left\{\frac{f(E_{k})[1-f(E_{m})][1-f(E_{n})]}{(E_{k}-E_{m})(E_{k}-E_{n})}+c.p.\right\}{}\times$$ (83)

$$\qquad\qquad{}\times\left<k\middle|V\middle|m\right>\left<m\middle|V\middle|n\right>\left<n\middle|V\middle|k\right>.$$

Accidentally, Eq. (82) is equivalent to the particle-hole symmetrized form of Eq. (83) obtained by taking a difference of Eq. (83) and the same expression with negated energies (this observation is specific to the third order perturbation).