# Thermalization of a Disordered Interacting System under an Interaction
Quench

Eric Dohner  Department of Physics, University at Albany (SUNY), Albany, New
York 12222, USA  Hanna Terletska  Department of Physics and Astronomy, Middle
Tennessee State University, Murfreesboro, TN 37132, USA  Herbert F Fotso
Department of Physics, University at Buffalo SUNY, Buffalo, New York 14260,
USA

###### Abstract

Although most studies of strongly correlated systems away from equilibrium
have focused on clean systems, it is well known that disorder may
significantly modify observed properties in various nontrivial ways. The
nonequilibrium interplay of interaction and disorder in these systems thus
requires further investigation. In the present paper, we use the recently
developed nonequilibrium DMFT+CPA embedding scheme, that combines both the
dynamical mean field theory (DMFT) and the coherent potential approximation
(CPA) nonequilibrium extensions, to characterize the relaxation and the
thermalization of a disordered interacting system described by the Anderson-
Hubbard model under an interaction quench. The system, initially in
equilibrium at a given temperature, has the interaction abruptly switched from
zero to a finite value at a given time. To investigate the role of disorder,
we use our effective medium approach to calculate, for different values of the
final interaction and of the disorder strength, the distribution functions as
the system evolves in time. This allows us to determine the effective
temperature after the quench and to analyze the effects of disorder on the
thermalization for various interaction strengths. We find that, for moderate
interactions after the interaction quench, disorder can tune the final
temperature of the system across a broad range of values with increased
disorder strength leading to lower effective temperature.

##  I Introduction

The dynamics of quantum systems away from equilibrium has been the subject of
increased interest as a result of the recent experimental advances extending
from quantum information processing platforms to time-resolved spectroscopies.
A salient question that has garnered a great deal of attention is that of how
quantum systems thermalize (or not) when they are abruptly driven out of
equilibrium. Beyond the theoretical question of how thermalization arises in
quantum systems that are supposed to be governed in their dynamics by unitary
time evolution operatorsDeutsch_PRA1994 ; Srednicki_PRE1994 ;
RigolEtAl_Nature2008 , these research questions have important experimental
consequences. For instance, it is often typical in the analysis of pump-probe
spectroscopy experiments to use a so-called “hot” electrons model whereby
electrons are driven by the pump pulse into an equilibrated state that is
thermalized at a higher temperature than that of the initial
systemPerfettiEtAl_PRL2006 ; PerfettiEtAl_PRL2007 . This brings into focus the
importance of the relevant relaxation scenarios and the associated timescales.
Also, experiments simulating various lattice models in optical lattices are
either intrinsically out of equilibrium or can be used to simulate, through
their high degree of tunability, the dynamics of nonequilibrium quantum
systemsBlochDaliwerger_RMP2008 ; Greinerbloch_Nat2002 ; BlochNatPhys2005 ;
BakrGreiner_Nat2009 . This further highlights the need for accurate modeling
and benchmarking.

While numerous efforts have been dedicated to the investigation of the
thermalization of correlated quantum systems away from equilibriumDMFT_noneq ;
FK_NonEq_DMFT08 ; DMFT_noneq_Aoki , little has been done to explore the effect
of disorder which we can anticipate, in some circumstances, to have
significant impacts on the dynamicsNandkishore_Huse_AnnRevCondMatPhys2015 ;
Kondov_DeMarco_PRL2015  and which we know to be ubiquitous in most systems of
interest. In particular, nonequilibrium dynamical mean field theory (DMFT) was
used to investigate the thermalization of correlated systems in a variety of
nonequilibrium scenarios extending from interaction
quenchesEcksteinKollarWerner_PRL2009 ; EcksteinKollar_PRL2008 , to DC field-
driven systemsthermalizationSciRep2014 ; NoneqFDT_Frontiers ; FreericksPRB2004
; FreericksTurkowskiZlatic_PRL2006 ; Freericks_PRB2008 , to simulations of
time-resolved spectroscopiesEcksteinWerner_PRB2011 ; MoritzEtAl_PRL2013 .
However, the effect of disorder in the thermalization of these nonequilibrium
systems remains generally understudied.

In this paper we use the recently developed nonequilibrium DMFT+CPA embedding
scheme NEDMFTCPA_PRB2022  that combines the nonequilibrium extensions of both
DMFTDMFT ; DMFT_2 ; DMFT_3 ; DMFT_4 ; DMFT_FK ; DMFT_noneq ; FK_NonEq_DMFT08 ;
DMFT_noneq_Aoki  and CPA (coherent potential approximation)CPA_Soven_1967 ;
CPA_Kirkpatrick ; CPA_Velicky_1969 ; CPA_Yonezawa_1973 ; NonEqCPA_1 ;
NonEqCPA_2 , to investigate the thermalization dynamics of a correlated
disordered system modeled by the Anderson-Hubbard model under an interaction
quench. In this way, we are able to assess the impact of the disorder on the
relaxation of the system and, specifically, to evaluate the temperature of the
system once it has settled into its long-time thermal state. We analyze the
nonequilibrium distribution functions calculated after the quench for various
values of the final interaction strengths and as a function of disorder
strength. We find that, for moderate interactions after the interaction
quench, disorder can tune the final temperature of the system across a broad
range of values with increased disorder strength leading to lower effective
temperature.

The rest of the paper is structured as follows: In Section II, we briefly
discuss the model and review the nonequilibrium DMFT+CPA formalism and its
numerical implementation. In Section III, we present the results that describe
the thermalization of the system after relaxation of the system following the
interaction quench. We end the paper with our conclusion in Section IV.

##  II Model and Methods

###  II.1 Model

We consider a correlated disordered system described by the Anderson-Hubbard
model initially in equilibrium at temperature $1/\beta$. The Hamiltonian is
given by Eq.(1). Where $t_{ij}=t_{hop}$ is the hopping amplitude between
nearest-neighbor sites (denoted by $\langle ij\rangle$), $U(t)$ is the Coulomb
interaction strength, and $V_{i}$ is the random onsite disorder for site $i$.
$c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ are respectively the creation and
the annihilation operators for a particle of spin $\sigma=\uparrow,\downarrow$
at site $i$. $n_{i\sigma}$ is the number of particles of spin
$\sigma=\uparrow,\downarrow$ at site $i$ and $\mu$ is the chemical potential.
We study the system at half-filling, such that $\mu=U/2$.

$\displaystyle H=-\sum_{\langle ij\rangle\sigma}$ $\displaystyle t_{ij}$
$\displaystyle(c^{\dagger}_{i\sigma}c_{j\sigma}+h.c.)+\sum_{i}U(t)n_{i\uparrow}n_{i\downarrow}$
(1) $\displaystyle+$
$\displaystyle\sum_{i\sigma}\left(V_{i}-\mu\right)n_{i\sigma},$

In equilibrium, the Coulomb interaction is constant $U(t)=U$. In the
nonequilibrium scenario of interest in this work, it is given by a step
function $U(t)=\Theta(t-t_{quench})U_{2}$ with $t_{quench}=0$, such that the
interaction is $U_{1}=0$ for negative times and some constant $U_{2}\neq 0$
for positive times. The onsite disorder $V_{i}$ is constant in time and
follows a uniform distribution such that
$P(V_{i})=\frac{1}{2W}\Theta(W-|V_{i}|)$, where $W$ is the disorder strength.
We use the notation $\langle...\rangle_{\\{V\\}}$ to indicate averaging over
all disorder values in the angle brackets. Here, we focus on the model for the
Bethe lattice in the limit of infinite coordination number.

###  II.2 Nonequilibrium DMFT+CPA

The nonequilibrium many-body formalism can be formulated on the Keldysh
contour whereby the system is evolved forward in time from an early
$t=t_{min}$ to times of physical interest up to a maximum value $t_{max}$ and
then back backward to the early times againKeldysh64_65 ;
StefanucciLeeuwen_CUP2013 ; rammer_2007 . The formalism involves several types
of two-time Green’s functions among which $G^{<}(t,t^{\prime})$ (the lesser),
$G^{>}(t,t^{\prime})$ (the greater), and $G^{R}(t,t^{\prime})$ (the retarded)
Green’s functions. In the context of a system initially in equilibrium at an
initial temperature $T=1/\beta$, a vertical spur of imaginary times of length
$-i\beta$ is added to the Keldysh contour resulting in the so called Kadanoff-
Baym-Keldysh contourKeldysh64_65 ; BaymKadanoff62 . In this situation, one
should add to the previous types of Green’s functions in the formalism, the
Matsubara Green’s function $G^{\tau}$, and the mixed time Green’s functions,
where one of the times is on either one of the horizontal branches of real
times, while the other is on the vertical branch of imaginary times. The
solution for a given problem can be either formulated in terms of the
different Green’s functions $G^{<}$, $G^{>}$ $G^{R}$, $G^{\tau}$, etc.
Alternatively, it can be formulated in terms of the contour-ordered Green’s
function $G_{c}(t,t^{\prime})$ from which all the others can be extracted. It
is this latter approach that we use in this work. The contour-ordered
quantities have time ordering performed with respect to time advance along the
entire contour. Hereafter we drop the subscript $c$ from the contour-ordered
quantities for convenience.

Our solution for the above described Anderson-Hubbard model under an
interaction quench is performed within the recently developed nonequilibrium
DMFT+CPA formalism which builds on the equilibrium formalismDMFT_CPA_1 ;
DMFT_CPA_2 ; DMFT_CPA_3 ; DMFT_CPA_4 ; DMFT_CPA_5  and is described
extensively in Ref.[NEDMFTCPA_PRB2022, ]. Here, for the sake of completeness,
we briefly summarise the algorithm. The method maps the lattice problem onto
that of an impurity embedded in a self-consistently determined medium
characterized by the hybridization $\Delta(t,t^{\prime})$ that is consistent
with that of DMFT for the clean system and with that of the disordered non-
interacting system for CPA.

In practice, the algorithm consists of the following self-consistency
procedure. From an initial guess of the hybridization function
$\Delta(t,t^{\prime})$, one obtains the noninteracting Green’s function for
each disorder configuration given by:

$\mathcal{G}_{V_{i}}(t,t^{\prime})=\left(\left(i\partial_{t}+\mu-
V_{i}\right)\delta_{c}-\Delta)\right)^{-1}(t,t^{\prime})$ (2)

From this, one obtains the Coulomb interaction self-energy. Here, similar to
Ref[NEDMFTCPA_PRB2022, ], we focus on the weak-to-moderate interaction and
disorder strengths regime, and we use second order perturbation theory so that
the self-energy is given by:

$\Sigma_{V_{i}}(t,t^{\prime})=-U(t)U(t^{\prime})\mathcal{G}_{V_{i}}(t,t^{\prime})^{2}\mathcal{G}_{V_{i}}(t^{\prime},t).$
(3)

After obtaining the self-energy for all disorder configurations, we evaluate
the disorder-averaged Green’s function:

$G_{ave}(t,t^{\prime})=\langle\left(G_{V_{i}}\right)\rangle_{\\{V\\}}$ (4)

where $G_{V_{i}}(t,t^{\prime})$ is the Green’s function for the disorder
configuration $\\{V_{i}\\}$:

$G_{V_{i}}(t,t^{\prime})=\left[\mathcal{G}_{V_{i}}^{-1}-\Sigma_{V_{i}}\right]^{-1}(t,t^{\prime}).$
(5)

This is followed by the evaluation of the updated hybridization function which
in the present case of the Bethe lattice with infinite coordination is given
by $\Delta(t,t^{\prime})={t^{*}}^{2}G_{ave}(t,t^{\prime})$ and the self-
consistency loop is repeated starting from the calculation of the new Coulomb
interaction self-energies and proceeds until convergence of the self-energy
within a desired criterion. $t^{*}$ is the hopping amplitude rescaled with the
coordination number $z$ so that $t_{hop}=\frac{t^{*}}{\sqrt{z}}$. We use
$t^{*}=0.25$ and thus set the bandwidth to be our energy unit and its inverse
to be the time unit.

Figure 1: Imaginary part of the retarded Green’s function as a function of
relative time for a range of average times, for $U_{2}=3t^{*},\;W=t^{*}$. The
long-time (relaxed) retarded Green’s function is represented by the dashed red
lines, and the causality time is marked for each average time with a black
dot. Note that the $t_{rel}$ at which $G^{R}$ begins to diverge from its
relaxed form is greater and greater with increasing $T_{ave}$. Inset:
illustration of the relationship between $(T_{ave}$, $t_{rel})$ and
$(t,t^{\prime})$, with the blue vertical and horizontal lines indicating the
time at which the quench occurs in $t$ and $t^{\prime}$ and black dots
indicating $t_{causality}$.

###  II.3 Numerical Implementation

Our implementation of the nonequilibrium DMFT+CPA follows the discrete time
construction of Refs.[NEDMFTCPA_PRB2022, ; Freericks_PRB2008, ]. The Kadanoff-
Baym-Keldysh contour is discretized into $(2N_{t}+N_{\tau})$ time steps, where
$N_{t}$ is the number of time steps on each leg of the horizontal real-time
branch of the contour and $N_{\tau}$ is the number of time steps on the
vertical branch of imaginary time. The step sizes are $\Delta
t=\left(t_{max}-t_{min}\right)/N_{t}$ for real time and
$\Delta\tau=\beta/N_{\tau}$ for imaginary time. In this paper, $t_{min}=-5$
and $t_{max}=20$ while the initial temperature of the system is such that
$\beta_{initial}=15$.

In this context, the contour-ordered quantities such as $G(t,t^{\prime})$
become square complex matrices $G_{ij}$ of size
$(2N_{t}+N_{\tau})\times(2N_{t}+N_{\tau})$. Convolutions of contour-ordered
quantities becomes matrix multiplications, and the continuous matrix inverse
becomes a discrete matrix inverse. The analysis is often performed by
switching from the $(t,t^{\prime})$ time coordinates to the Wigner coordinates
$(T_{ave},t_{rel})$ where $T_{ave}$ can be viewed as the effective time of the
system while frequency domain information is obtained by Fourier transforming
with respect to $t_{rel}$. Observables calculated from the discretized
contour, such as the distribution functions and the energy, are often obtained
for multiple step sizes then extrapolated to the continuum limit $\Delta t\to
0$. We use standard Lagrange interpolating polynomials to quadratic order.

Figure 2: Extrapolated potential, kinetic, and total energies for
$U_{2}=2t^{*}$. The vertical black line shows the time at which we evaluate
the relaxed distribution function. This time is well after the relaxation of
the system.

##  III Results

The system is initially in equilibrium at temperature $T=1/\beta$ with
$\beta_{initial}=15$. While keeping the disorder strength $W$ constant, the
interaction quench is applied at time $t=t_{quench}=0$ with the interaction
abruptly changing from an initial value $U_{1}=0$ to a final value $U=U_{2}$.
We are interested in tracking the thermalization of the system at long times.
Our analysis is guided by two fundamental quantities: the density of states
and the distribution function. For a thermalized system, the former is given
by the retarded Green’s function, while the latter is given by the lesser
Green’s function. Namely:

$\rho(\omega)=-i\mathrm{Im}G^{R}(\omega)/\pi$ (6)

and according to the fluctuation dissipation theorem, for a thermalized
system,

$G^{<}(\omega)=-2iF(\omega)\mathrm{Im}G^{R}(\omega).$ (7)

Where $F(\omega)$ is the distribution function. In the nonequilibrium
formalism, we can track these quantities as a function of average time.

###  III.1 Density of states

For the system at half-filling, we know that the real part, in the time
domain, of the retarded Green’s function vanishes for all average times
thermalizationSciRep2014 . Thus, the density of states is fully defined by the
imaginary part of the retarded Green’s function in the time domain. For this
reason, we can track the dynamics of the density of states through the
imaginary part of the retarded Green’s function in the time domain. Fig. 1
presents the typical behavior of $G^{R}(T_{ave},t_{rel})$ as a function of
$t_{rel}$ for a series of $T_{ave}$ values. Note that the relative time axis
(represented by the green lines in the insert), for earlier values $T_{ave}$
(magenta lines in the insert), has segments of time coordinates
$(t,t^{\prime})$ for which one (or both) of the times is (are) before the
interaction quench leading to a mixed character of the corresponding $t_{rel}$
coordinates. The blue lines in the main figure correspond to successive
$T_{ave}$ values after the quench, while the dashed red line corresponds to an
average time value after the quench for which all $t_{rel}$ involves both $t$
and $t^{\prime}$ that have the new interaction strength $U_{2}$. The black
circles correspond to the causality time beyond which $t_{rel}$ has mixed
character. One can see on this figure that the solid blue curves overlap with
the dashed red curve up to the causality time and that the retarded Green’s
function is only constrained by causality. So, the density of states of states
is immediately established after the quench. The relaxation of the system can
thus be tracked through the distribution function.

Figure 3: Imaginary parts of lesser and retarded Green’s functions as a
function of frequency for $U_{2}=2t^{*},\;W=2t^{*}$ at the analysis time. The
shaded box shows the region over which we evaluate the distribution function
$F(\omega)=-Im[G^{<}]/(2Im[G^{R}])$. Outside of this region, the ratio is
prone to numerical instabilities due the Gibbs phenomenon in the frequency
data obtained and to the division by small numbers.

###  III.2 Distribution function $F(\omega)$

In the present study, we are interested in the thermalization of the system
after it has undergone its early transient following the quench. Fig. 2 shows,
for different disorder strengths and for $U_{2}=2t^{*}$, the evolution in time
of the kinetic, potential and total energies of the system evaluated following
Refs.[ NEDMFTCPA_PRB2022, ; HubbardQuenchEckstein, ]. The quench is performed
at time $t=0$. After an initial nontrivial response to the quench, the
observables settle into a constant value for the remaining duration of the
simulation. The vertical black line indicates the time $t=7.5$ at which the
long-time analysis is performed.

Given that the density of states is established immediately after the quench
and is only constrained by causality, this analysis time is chosen so as to
allow a range of $t_{rel}$ values that enables a reliable Fourier transform.
To obtain the distribution function, we will use the fluctuation-dissipation
theorem as expressed by Eq.(7). To this end, we first Fourier transform the
lesser and retarded Green’s functions $G^{R/<}(T_{ave},t_{rel})$ in relative
time to yield $G^{R/<}(T_{ave},\omega)$. The result of this operation is
illustrated for $U_{2}=2t^{*}$ and $W=2t^{*}$ in Fig. 3. To avoid numerical
instabilities, the distribution function is only evaluated in a frequency
range around $\omega=0$ for which both $G^{R}(T_{ave},\omega)$ and
$G^{<}(T_{ave},\omega)$ remain finite as illustrated by the shaded box in Fig.
3.

Figure 4: Relaxation of $F(\omega)$ at $U_{2}=3t^{*},W=t^{*}$ soon after the
quench, but before thermalization, demonstrating the non-thermal form of the
distribution function at the early stages of the relaxation. Figure 5: Post-
relaxation $F(\omega)$ for times slightly before and slightly after our
analysis time for $U_{2}=3t^{*},W=t^{*}$, demonstrating that the distribution
function changes minimally around this analysis time at which we evaluate the
relaxed $F(\omega)$.

Figs. 4 and (5) show the extracted distribution function for $U_{2}=3t^{*}$
and $W=t^{*}$ for different average times. One can readily observe that
following the interaction quench at time $t=0$, the distribution function
initially changes in a highly nontrivial way and may in fact clearly
correspond to a non-thermal system (Fig. 4). However, around our analysis
time, corresponding to $T_{ave}=7.5$, the distribution function is seen to
change very little for different values of the average time and the different
curves essentially overlap (Fig. 5). For this reason, the system can be
assumed to have settled into its long-time state at time $T_{ave}=7.5$. It is
in this regime that we evaluate a long time effective temperature of the
system after the quench.

Figure 6: Distribution function for the equilibrium system with
$U=2t^{*},\;W=2t^{*},\;\beta=15$, and after relaxation for the quenched system
with $U_{2}=2t^{*},\;W=2t^{*},\;\beta_{initial}=15$. The dashed line shows the
fit to the quenched system distribution function after the transient. The
shaded box indicates the region over which the fit is performed. Here we fit
the Fermi function,
$F_{Fit}(\omega)=1/(1+\mathrm{exp}\left(\beta\omega)\right)$ with $\beta$ as a
free parameter, to the calculated $F(\omega)$, and this allows us to extract
an effective temperature. Figure 7: Inverse effective temperature as a
function of the final interaction strength for different disorder strengths.
The systems is initially at a temperature such that $\beta_{initial}=15$.
Inset: Effective inverse temperature $\beta$ as a function of the disorder
strength for different interaction strengths. Increased disorder strength for
moderate interaction strengths leads to a lower long-time temperature.

###  III.3 Effective temperature

The effective temperature is obtained by fitting a Fermi-Dirac distribution
function ($F(\omega)=1/(1+\mathrm{exp}\left(\beta\omega)\right)$ with $\beta$
as a free parameter) to the extracted distribution function over a frequency
window around $\omega=0$ as illustrated in Fig. 6. As indicated above, after
the quench but before relaxation, the distribution function can take non-
thermal forms (Fig.4). Consequently, an effective temperature cannot be traced
over the entire time evolution of the system. However, this procedure is well-
defined for the chosen analysis time for the long-time behavior.

Fig. 7 shows the long time effective temperature of the system as a function
of the final interaction strength $U_{2}$ with different solid lines
corresponding to different values of the disorder strength $W$. The inset
shows the same data but with the disorder strength on the $x$-axis and
different dashed lines corresponding to different values of the final
interaction strength. The figures show the significant dependence of the final
inverse temperature on disorder strength. For weak $U_{2}$ values, increased
disorder strength leads to small increase in the long time temperature.
However, as the interaction strength $U_{2}$ is increased, we observe that
increasing the disorder strength leads lower long-time effective temperatures.
This shows that under an interaction quench, the long-time temperature can
vary over a broad range of values depending on the disorder strength, with
increased disorder strength leading to lower final temperature.

##  IV Conclusion

We have analyzed the relaxation of a disordered interacting system after an
interaction quench where, with the disorder strength held constant, the
interaction strength is abruptly switched from zero to a finite value $U_{2}$
at which it is subsequently kept. We have used the recently developed
nonequilibrium DMFT+CPA formalism that maps the lattice problem onto an
effective mean field that is equivalent to that of the dynamical mean field
theory (DMFT) for the clean system and to that of the coherent potential
approximation (CPA) for the disordered noninteracting system. By extracting
the distribution function from the Green’s function using the fluctuation-
dissipation theorem, we showed that while the early transient does not follow
the fluctuation dissipation theorem, at longer times, the system settles into
a thermal state at a constant temperature. This long time temperature is
lowered by increased disorder strengths at moderate values of the interaction.
Altogether our studies demonstrate that after the interaction quench, disorder
can tune the long-time temperature of the system over a broad range of values.

## Acknowledgments

HFF is supported by the National Science Foundation under Grant No.
PHY-2014023. HT has been supported by NSF DMR-1944974 grant.

## References

  * (1) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * (2) M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * (3) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854–858 (2008).
  * (4) L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and M. Wolf, Phys. Rev. Lett., 97, 067402 (2006).
  * (5) L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Eisaki, and M. Wolf, Phys. Rev. Lett., 99, 197001 (2007).
  * (6) I Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys., 80, 885 (2008).
  * (7) Markus Greiner, Olaf Mandel, Tilman Esslinger, Theodor W. Hänsch, and Immanuel Bloch, Nature 415, 39 (2002).
  * (8) Immanuel Bloch, Nature Physics 1, 23 (2005).
  * (9) Waseem S. Bakr, Jonathon I. Gillen, Amy Peng, Simon Fölling, and Markus Greiner, Nature 462, 74 (2009).
  * (10) J. K. Freericks, V. M. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006).
  * (11) J. K. Freericks, Phys. Rev. B 77, 075109 (2008).
  * (12) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * (13) Rahul Nandkishore, and David A. Huse, Annual Review of Condensed Matter Physics, 6, 15, (2015).
  * (14) S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco, Phys. Rev. Lett. 114, 083002 (2015).
  * (15) M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. 103, 056403 (2009).
  * (16) M. Eckstein, and M. Kollar, Phys. Rev. Lett. 100, 120404 (2008).
  * (17) H. F. Fotso, K. Mikelsons and J. K. Freericks, Scientific Reports 4, 4699 (2014).
  * (18) H. F. Fotso and J. K. Freericks, Frontiers in Physics, 8, 324 (2020).
  * (19) J. K. Freericks, Phys. Rev. B 70, 195342 (2004).
  * (20) J. K. Freericks, V. M. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006).
  * (21) J. K. Freericks, Phys. Rev. B 77, 075109 (2008).
  * (22) M. Eckstein, and P. Werner, Phys. Rev. B 84, 035122 (2011).
  * (23) B. Moritz, A. F. Kemper, M. Sentef, T. P. Devereaux, and J. K. Freericks, Phys. Rev. Lett. 111, 077401 (2013).
  * (24) E. Dohner, H. Terletska, K.-M. Tam, J. Moreno, and H. F. Fotso, Phys. Rev. B 106, 195156 (2022).
  * (25) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * (26) Y. Kuramoto, Springer Series in Solid State Science Vol. 62, edited by T. Kasuya and T. Sao (Springer, 1985), p. 152.
  * (27) E. Müller-Hartmann, Z. Phys. B 74, 507 (1989).
  * (28) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg Rev. Mod. Phys. 68, 13 (1996).
  * (29) J. K. Freericks and V. Zlatić, Rev. Mod. Phys. 75, 1333 (2003).
  * (30) P. Soven, Phys. Rev.156, 809 (1967).
  * (31) S. Kirkpatrick, B. Velický, H. Ehrenreich, Phys. Rev. B 1, 3250 (1970).
  * (32) P. Soven, Phys. Rev.184, 614 (1969).
  * (33) F. Yonezawa and K. Morigaki, Prog. Theor. Phys. Supp. 53, 1 (1973).
  * (34) Y. Zhu, L. Liu, and H. Guo, Phys. Rev. B 88, 205415 (2013).
  * (35) A. V. Kalitsov, M. G. Chshiev, and J. P. Velev, Phys. Rev. B 85, 235111 (2012).
  * (36) L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1945 (1964) [Sov. Phys. JETP 20, 1018 (1964)].
  * (37) G. Stefanucci and R. van Leeuwen, Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction, Cambridge University Press (2013).
  * (38) Jörgen Rammer, Quantum Field Theory of Non-equilibrium States, Cambridge University Press (2007).
  * (39) L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics ( Benjamin, New York, 1962).
  * (40) D. Semmler, K. Byczuk, and W. Hofstetter, Phys. Rev. B 84, 115113 (2011).
  * (41) E. Miranda, and V. Dobrosavljevic, ArXiv:1112.6184 (2011).
  * (42) E. Z. Kuchinskii, N. A. Kuleeva, I. A. Nekrasov, and M. V. Sadovskii, JETP 110 (2), 325 (2010).
  * (43) A. Weh, Y. Zhang, A. Östlin, H. Terletska, D. Bauernfeind, K.-M. Tam, H. G. Evertz, K. Byczuk, D. Vollhardt, and L. Chioncel, Phys. Rev. B 104, 045127 (2021).
  * (44) V. Janiš and D. Vollhart, Phys. Rev. B 46, 15712 (1992).
  * (45) M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. B 81, 115131 (2010).


Binary Cross Entropy loss. The labels for benign and adversarial examples are
$\boldsymbol{1}$ and $\boldsymbol{0}$ respectively, where $\boldsymbol{1}$
($\boldsymbol{0}$) is a vector with all the elements are $1$ ($0$).

When applying SAM to the benign term loss, our optimization procedure is as
follows: First, compute the gradient of other terms in the total loss with
respect to the current model parameters and save it as $g_{other}$. Second,
zero out current gradient then follow the procedure of SAM to compute the SAM
gradient of only the benign term, $g_{benign}$. Then average $g_{other}$ and
$g_{benign}$ to get the final gradient update.

###  D.2 Experimental Settings

##### General Settings:

Cifar-100 consists of 50000 training images and 10000 test images over 100
classes. All the images are in the shape of $3\times 32\times 32$. We apply on
this dataset the same list of augmentations as for Cifar-10 [2], except the
pixel values are are scaled to the range $[0,1]$ without normalization. We
also use the same backbone, projection head, classification head and
discriminator as for Cifar-10 to run on Cifar-100.

##### SSL Settings:

For two datasets, the default training configuration of the discriminator is
as follows: The optimizer is plain SGD with the learning rate of $0.5$,
$\lambda_{global}$ is $0.5$. For sharpness-aware minimization, we adopt A-SAM
[57] with $\rho=1.0$.

###  D.3 Experimental Results

In this part, for the ease of readability, we repeat the results on Cifar-10
in Table 5, and present additional results on Cifar-100 in Table 4.

From Table 4, first, consider settings without D and without SAM, it can be
seen that as we increase $\lambda_{benign}$, clean and robust accuracy under
LE and AT-LE both increase. Furthermore, when the benign term loss is
incorporated with SAM, i.e. settings with SAM and without D, the model
performance can be further improved.

Second, consider settings without SAM, with the introduction of the
discriminator, model achieves better robust accuracy in all settings except
the one under AT-LE with $\lambda_{benign}=2.0$ in which comparable value is
obtained. Interestingly, clean accuracy is also improved consistently over all
settings.

Finally, when combining SAM and D, full method, we can generally improve clean
and robust accuracy. In particular, when $\lambda_{benign}=1.0$, our full
method obtains the highest scores; when $\lambda_{benign}=0.5$, under AT-LE,
it outperforms other variants; and when $\lambda_{benign}=2.0$, it obtains the
highest robust performance and comparable clean accuracy to its variants.

Table 4: Clean and robust accuracy of a ResNet-18 model trained with and without our method’s components, across different weighting for the benign loss term. Dataset: Cifar-100. $\lambda_{\mathrm{benign}}$ |  SAM |  D | LE | AT-LE
---|---|---|---|---
|  |  Clean |  Adv. |  Clean |  Adv.
0 |  ✗ |  ✗ |  26.23 |  13.86 |  23.47 |  14.95
✗ |  ✓ |  27.45 |  14.44 |  24.94 |  15.65
0.5 |  ✗ |  ✗ |  31.70 |  15.06 |  27.53 |  16.25
✗ |  ✓ |  32.19 |  15.44 |  27.88 |  16.40
✓ |  ✗ |  31.79 |  15.14 |  27.67 |  16.28
✓ |  ✓ |  32.94 |  15.31 |  28.76 |  16.79
1.0 |  ✗ |  ✗ |  36.35 |  16.11 |  31.22 |  17.32
✗ |  ✓ |  36.81 |  16.63 |  31.77 |  17.78
✓ |  ✗ |  36.58 |  16.37 |  31.32 |  17.79
✓ |  ✓ |  37.72 |  17.02 |  32.34 |  17.99
2.0 |  ✗ |  ✗ |  43.36 |  17.06 |  35.65 |  18.88
✗ |  ✓ |  43.86 |  17.43 |  35.90 |  18.83
✓ |  ✗ |  43.58 |  17.29 |  35.94 |  18.94
✓ |  ✓ |  43.26 |  17.62 |  35.72 |  18.98
Table 5: Clean and robust accuracy of a ResNet-18 model trained with and without our method’s components, across different weighting for the benign loss term. Dataset: Cifar-10. $\lambda_{\mathrm{benign}}$ |  SAM |  D | LE | AT-LE
---|---|---|---|---
|  |  Clean |  Adv. |  Clean |  Adv.
0 |  ✗ |  ✗ |  70.54 |  38.97 |  66.28 |  42.73
✗ |  ✓ |  71.04 |  39.08 |  66.15 |  42.57
0.5 |  ✗ |  ✗ |  77.54 |  40.72 |  72.84 |  44.23
✗ |  ✓ |  77.66 |  40.83 |  73.32 |  44.46
✓ |  ✗ |  77.83 |  40.09 |  73.33 |  44.23
✓ |  ✓ |  77.76 |  40.05 |  73.35 |  44.17
1.0 |  ✗ |  ✗ |  81.14 |  41.08 |  76.36 |  44.69
✗ |  ✓ |  81.02 |  41.34 |  76.95 |  44.90
✓ |  ✗ |  81.31 |  41.14 |  76.98 |  44.95
✓ |  ✓ |  81.64 |  41.23 |  77.33 |  45.28
2.0 |  ✗ |  ✗ |  84.80 |  41.78 |  81.35 |  45.71
✗ |  ✓ |  84.12 |  41.93 |  80.92 |  45.89
✓ |  ✗ |  85.09 |  41.97 |  81.64 |  46.30
✓ |  ✓ |  85.01 |  41.82 |  81.86 |  46.45


ULB-TH/21-06

# Dark matter from dark photons

L. VANDERHEYDEN  Service de Physique Théorique, Université Libre de Bruxelles,
Boulevard du Triomphe, CP225, 1050 Brussels, Belgium

###### Abstract

In this talk, based on [1], I explained how the observed dark matter (DM)
relic abundance can be accounted for in models composed of three sectors (the
DM, the Standard Model (SM) and a light mediator) connected to each other.
This scenario is explored in the context of the massive dark photon model in
which the DM is connected to the SM through the kinetic mixing between the
hypercharge gauge boson and the vector gauge boson associated to a new
$U(1)^{\prime}$ gauge symmetry. In such portal models, the DM relic abundance
can be produced in no less than nine regimes along five dynamical mechanisms.
In particular, I emphasize the sequential freeze-in dynamical mechanism which
can occurs when the dark photon is massive and consists in two successive
freeze-in mechanisms.

##  1 Introduction

The idea that the DM could belong to a hidden sector (HS) is an interesting
possibility, especially given the fact that so far no experiments has come
with a clear evidence for the nature of DM as a particle, even if a HS could
be connected to the visible sector (VS) made of SM particles through a portal
[2]. How to account for the observed DM relic abundance in portal models has
already been studied in various works, see in particular the detailed work of
for the kinetic mixing portal case with a massless dark photon [3]. We review
here the massive dark photon case, as studied at length in [1]. The kinetic
mixing model introduces a new $U(1)^{\prime}$ gauge symmetry under which the
DM candidate, a Dirac fermion $\chi$, is charged. The gauge vector boson
$B^{\prime\mu}$ associated to this new symmetry is called dark photon and
kinetically mixes with the SM massless boson $B^{\mu}$ through the parameter
$\hat{\epsilon}$. The Lagrangian of this portal model can be written as
following,

$\displaystyle\mathcal{L}^{\prime}=-\frac{1}{4}B^{\prime\mu\nu}B^{\prime}_{\mu\nu}-\frac{\hat{\epsilon}}{2}B^{\mu\nu}B^{\prime}_{\mu\nu}+\frac{1}{2}m_{\gamma^{\prime}}^{2}B^{\prime\mu}B^{\prime}_{\mu}+i\bar{\chi}\not{D}\chi-
m_{\rm DM}\bar{\chi}\chi,$ (1)

where the covariant derivative of $\chi$ is expressed in terms of the dark
coupling $e^{\prime}$: $D^{\mu}=\partial^{\mu}+ie^{\prime}B^{\mu}$.

The massive dark photon case brings an all range of new phenomenons which do
not appear in the massless case. In order to see that, let us extract, from
the Lagrangian given in Eq. 1, the actual relevant interactions for DM
production.

##  2 Three sectors - three connectors

It is necessary to re-express the Lagrangian of Eq. 1 in the basis of physical
states (i.e. kinetic and mass eigenstates). To do so, one needs to perform
first a non-orthogonal transformation followed by the Weinberg rotation
($\theta_{W}$) and, finally, an orthogonal diagonalisation of the remaining
mass terms, see [1] for more details. At leading order in the mixing parameter
$\epsilon\equiv\hat{\epsilon}\cos\theta_{W}$, one has

$\displaystyle\mathcal{L}_{I}=-\frac{e}{\sin\theta_{W}\cos\theta_{W}}J^{\mu}_{Z}Z_{\mu}-e^{\prime}J^{\mu}_{\chi}\left(\epsilon\sin\theta_{W}Z_{\mu}-\gamma^{\prime}_{\mu}\right)+eJ^{\mu}_{EM}\left(\gamma_{\mu}-\epsilon\cos\theta_{W}\gamma^{\prime}_{\mu}\right),$
(2)

for the interaction part of the Lagrangian. One sees from Eq. 2 that new
connections arise in the massive dark photon scenario with respect to the
massless case, see Figure 1. Indeed, when the dark photon is massless, one has
two massless gauge bosons which can be rotated to go in a basis where the two
last diagrams of Figure 1 do not exist and in which one cannot produce
directly the dark photon from the SM bath as well known, see e.g. [4].
Processes depicted in Figure 1 are all driven by one of the three connections
of the theory: $\alpha^{\prime}\equiv\frac{e^{\prime 2}}{4\pi}$ rules the
interaction between the DM and dark photon particles, $\epsilon$ rules
interactions between the dark photon and SM particles and
$\kappa^{\prime}\equiv\epsilon\sqrt{\frac{\alpha^{\prime}}{\alpha}}$ rule
interactions between the DM and SM particles. These three populations (the DM,
mediator and SM particles) are then connected to each other through three
connections instead of two as in the massless case, see Figure 2. Notice that
taking into account finite temperature corrections the $\epsilon$ coupling is
effectively replaced by $\epsilon_{\rm eff}\equiv\epsilon\times
m_{\gamma^{\prime}}^{2}/(m_{\gamma^{\prime}}^{2}-\Pi_{\gamma,T})$ where
$\Pi_{\gamma,T}$ is the self-energy of the transverse propagating photons
$\Pi_{\gamma,T}$,[4, 5, 6] so that taking the limit
$m_{\gamma^{\prime}}\rightarrow 0$ leads well to the massless case with only 2
physical connections.

Figure 1: Feynmann diagrams of the relevant processes for DM and dark photon
production in the kinetic mixing portal model for the massive dark photon
case.

Depending on the strength of these three couplings, the three populations (the
DM, mediator and SM particles) may or may not thermalise with each others,
leading to DM relic abundance production mechanisms which are
phenomenologically very different. The DM relic abundance depends on four
independent parameters which are the DM and dark photon masses ($m_{\rm DM}$
and $m_{\gamma^{\prime}}$) as well as two of the three couplings (we choose
$\alpha^{\prime}$ and $\kappa^{\prime}$). Then, fixing the masses, we can plot
contours of the DM relic abundance in the $\kappa^{\prime}-\alpha^{\prime}$
plane starting with an empty dark sector (i.e. assuming that no DM and no dark
photon are produced at the end of inflation) then going through all production
regimes all the way from phases where the DM particles never thermalise with
any other particles to phases where all particles in the early Universe
thermal bath thermalise with each other. This has been done in Figure 3 for
$\left(m_{\rm DM},m_{\gamma^{\prime}}\right)=\left(3\,{\rm GeV},1\,{\rm
GeV}\right)$ in the left panel and $\left(m_{\rm
DM},m_{\gamma^{\prime}}\right)=\left(100\,{\rm GeV},10\,{\rm GeV}\right)$ in
the right panel.

In order to get the phase diagrams shown in left and right panels of Figure 3,
one has to integrate a set of Boltzmann equations for the DM and the mediator
yields as functions of $x\equiv m_{\rm DM}/T$:

Figure 2: Three sectors and their connections in the massive (left) and
massless (right) dark photon cases.

$\displaystyle xHs\frac{\mathrm{d}Y_{\rm DM}}{\mathrm{d}x}$ $\displaystyle=$
$\displaystyle
2\langle\sigma_{\gamma^{\prime}\gamma^{\prime}\rightarrow\chi\bar{\chi}}v\rangle
n^{2}_{\gamma^{\prime}}-2\langle\sigma_{\chi\bar{\chi}\rightarrow\gamma^{\prime}\gamma^{\prime}}v\rangle
n^{2}_{\chi}$ (3) $\displaystyle+2\langle\sigma_{\chi\bar{\chi}\rightarrow
f\bar{f}}v\rangle\left[({n_{\chi}^{\rm
eq}})^{2}-n^{2}_{\chi}\right]+2\langle\Gamma^{D}_{Z\rightarrow\chi\bar{\chi}}\rangle\frac{n_{Z}^{\rm
eq}}{(n_{\chi}^{\rm eq})^{2}}\left[{(n_{\chi}^{\rm
eq}})^{2}-n^{2}_{\chi}\right]$ $\displaystyle
xHs\frac{\mathrm{d}Y_{\gamma^{\prime}}}{\mathrm{d}x}$ $\displaystyle=$
$\displaystyle
2\langle\sigma_{\chi\bar{\chi}\rightarrow\gamma^{\prime}\gamma^{\prime}}v\rangle
n^{2}_{\chi}-2\langle\sigma_{\gamma^{\prime}\gamma^{\prime}\rightarrow\chi\bar{\chi}}v\rangle
n^{2}_{\gamma^{\prime}}$ (4)
$\displaystyle+\left(\langle\sigma_{f\gamma\rightarrow
f\gamma^{\prime}}v\rangle+\langle\sigma_{f\bar{f}\rightarrow\gamma\gamma^{\prime}}v\rangle\right)\left[({n_{\gamma^{\prime}}^{\rm
eq}})-n_{\gamma^{\prime}}\right]+2\langle\Gamma^{D}_{{\rm
H}\rightarrow\gamma^{\prime}\gamma^{\prime}}\rangle\frac{n_{\rm H}^{\rm
eq}}{(n_{\gamma^{\prime}}^{\rm eq})^{2}}\left[{(n_{\gamma^{\prime}}^{\rm
eq}})^{2}-n^{2}_{\gamma^{\prime}}\right].$

where $\Gamma^{D}$ refers to the decay rate of the SM $Z$ boson into two DM
particles or of the SM $H$ boson into two $\gamma^{\prime}$. We kept a sum
over all SM fermions channels implicit for simplicity.

Integrating Eqs. 3 and 4 gives rise to no less than nine regimes along five
distinct dynamical mechanisms111In Fig. 3 each corner, that the dashed line
does, indicates a change of regime. In the massless mediator case, five
regimes and four dynamical mechanisms were already found and discussed in [3]
for the Kinetic Mixing model and a scalar portal.. This five dynamical
production mechanisms which are the freeze-in (I), the sequential freeze-in
(II), the reannihilation (III), the secluded freeze-out (IV) and, finally, the
freeze-out (V) are distinguishable looking at which of the three connections
of left panel of Figure 2 thermalises or not as well as looking at which
connection dominates over the others.

Figure 3: Contours of the DM relic density in the
$\alpha^{\prime}-\kappa^{\prime}$ plane for $\left(m_{\rm
DM},m_{\gamma^{\prime}}\right)=\left(3\,{\rm GeV},1\,{\rm GeV}\right)$ (left
panel) and $\left(m_{\rm DM},m_{\gamma^{\prime}}\right)=\left(100\,{\rm
GeV},10\,{\rm GeV}\right)$ (right panel) in the kinetic mixing portal model
for the massive dark photon case. The dashed contour indicates $\Omega_{\rm
DM}h^{2}=0.1188$.

As an example, let us explain in more details the Sequential Freeze-In regime
(II) which corresponds to the vertical regime on the extreme left of each
panel of Figure 3. In this part of the parameter space, one can check that
none of the three sectors thermalises with any of the two others. That is to
say that,

$\displaystyle\Gamma_{\rm SM\leftrightarrow DM}<H\hskip
14.22636pt{\rm\&}\hskip 14.22636pt\Gamma_{{\rm
SM}\leftrightarrow\gamma^{\prime}}<H\hskip 14.22636pt{\rm\&}\hskip
14.22636pt\Gamma_{{\rm DM}\leftrightarrow\gamma^{\prime}}<H,$ (5)

where $\Gamma_{\rm X\leftrightarrow Y}$ stems for the total interaction rate
between sector $X$ and sector $Y$ and with $H$ the Hubble expansion rate.
Thus, as no sector thermalises with any another one and starting with an empty
dark sector,222The opposite option (more initial condition dependent and which
we do not study here) would be that the final amount of DM would be mainly due
to the presence of many HS particles at the end of inflation with little role
of the portal. the DM relic abundance is slowly produced by slow out-of-
equilibrium processes from SM. Along this regime the value of the DM-to-SM
coupling, $\kappa^{\prime}$, turns out to be too low to directly play this
role. Thus DM cannot be produced from an usual $SM\rightarrow DM$ freeze-in
regime (first diagram of Figure 1). Nevertheless, it turns out along this
regime that the SM can slowly produces enough dark photon through out-of-
equilibrium ${\rm SM}\rightarrow\gamma^{\prime}$ processes, for these dark
photons to slowly produce subsequently enough DM through
$\gamma^{\prime}\rightarrow$ DM processes, see the two last diagrams depicted
in Figure 1. In other words these dark photons never enter thermal
equilibrium, but can still slowly produce DM through out-of-equilibrium
processes $\gamma^{\prime}\rightarrow{\rm DM}$ (second diagram of Figure 1).
Therefore, it is possible that the whole DM observed today can be produced
through these two successive freeze-in mechanisms, a ”Sequential Freeze-In”.

##  3 Summary

The DM relic abundance in a three sectors-three connections scenario can be
set through five different dynamical mechanisms, and this along nine distinct
regimes. We illustrated these results through the massive dark photon model of
Eq. 1, generalizing the massless dark photon case [3]. Among the four new
dynamical mechanisms, a new DM production regime which we dubbed ”Sequential
Freeze-In” (II) consists first of a slow out-of-equilibrium production of dark
photons from SM particles followed by a slow out-of-equilibrium production of
DM particles from dark photons particles. The other new regimes are DM freeze-
in production from a thermalised population of dark photons (Ib),
reannihilation (IIIa) and secluded freeze-out (IVa) through the freeze-in
production of dark photons, see [1]. Although we considered the kinetic mixing
portal model to illustrate our findings, we argue that the behaviours we
presented are actually more generics. They are characteristics of DM scenarios
in which the DM is connected to the SM bath through a portal carried by a new
particle which would give rise to three sectors connected to each others by
three connections due to two different couplings, here $\alpha^{\prime}$ and
$\epsilon$. This also holds in particular for Higgs portal models, see [1] and
also [7].

## Acknowledgments

I thank the organizers for this special edition of the Rencontres de Moriond.
I also thank T. Hambye, M.H.G. Tytgat and J. Vandecasteele for collaboration
on the work discussed in this note. My work is supported by the FRIA grant
from FRS-FNRS.

## References

## References

  * [1] T. Hambye, M. H. G. Tytgat, J. Vandecasteele, and L. Vanderheyden, Phys. Rev. D 100, 095018 (2019)
  * [2] B. Patt and F. Wilczek, arXiv:hep-ph/0605188.
  * [3] X. Chu, T. Hambye, and M. H. G. Tytgat, J. Cosmol. Astropart. Phys. 05 (2012) 034
  * [4] J. Redondo, J. Cosmol. Astropart. Phys. 07 (2008) 008
  * [5] H. An, M. Pospelov, and J. Pradler, Phys. Lett. B 725, 190 (2013)
  * [6] J. Redondo and G. Raffelt, J. Cosmol. Astropart. Phys. 08 (2013) 034
  * [7] G. Bélanger, C. Delaunay, A. Pukhov and B. Zaldivar, Phys. Rev. D 102, 035017 (2020)


# Higher Derivative Muffin Tin Orbitals (HDMTO) and Higher Derivative Koringa
Khon and Rostoker (HDKKR) methods

Garry Goldstein1<EMAIL_ADDRESS>
###### Abstract

In this work we introduce a Linearized version of the Koringa Khon and
Rostoker method (LKKR) and show it to be equivalent to the Linearized Muffin
Tin Orbitals method (LMTO). We then present higher derivative versions of both
methods, e.g. HDKKR and HDMTO and show them to be partially distinct (not
equivalent). In particular HDKKR basis set does not have an equivalent ground
state for the Khon Sham (KS) Hamiltonian as the HDKKR basis set and has
greater variational power then the HDMTO one. Because the KS method, for
Density Functional Theory (DFT), is variational HDKKR will give better ground
state energies then HDMTO. However HDKKR is much harder to work with then
HDMTO requiring much greater computer resources so HDMTO can often be
preferred.

##  I Introduction

The choice of basis set is fundamental for an efficient solution of a Density
Functional Theory problem (DFT) - the diagonalization of the Khon Sham (KS)
Hamiltonian. The basis set must be efficient, in that a small number of basis
elements well approximate the exact Khon Sham (KS) wavefunctions at least in
the middle of the valence band. The basis must also be simple to manipulate -
with various practical calculations associated with finding matrix element of
the KS Hamiltonian with respect to the basis easily implementable on computer.
The basis set must be transferable, that is one does not want to use a
significantly different basis for every single many body problem or
equivalently for diagonalizing every KS Hamiltonian. The Koringa Khon Rostoker
(KKR) and the Augmented plane Wave (APW) basis sets are such basis sets,
however they are deficient in that the basis elements, of both basis sets,
depend explicitly on the energy of the eigenstate, to be computed, which leads
to a self consistent calculation of the basis set thereby increasing
computational costs by easily an order of magnitude [1, 2]. The key idea to
overcome this difficulty is due to O. K. Andersen who linearized the basis set
(with the linearization energy being chosen in the middle of the valance band)
thereby obtaining a fixed basis for each iteration of the self consistency
loop for the solution of the self consistent KS problem [3, 1, 2, 4], each
step requires simply the diagonalization of a matrix. This lead to the
introduction of the Linearized Augmented Plane Wave (LAPW) [3, 1, 2] and
Linearized Muffin Tin Orbitals (LMTO) methods [6, 5]. In both these methods
the solid is divided into a Muffin Tin (MT) sphere part (with spheres centered
around atomic nuclei) and an interstitial part [3, 1, 2]). In the interstitial
part the KS potential is assumed to be sufficiently smooth that plane waves or
spherical Bessel functions form a good basis for the region while the basis
set is adapted to the MT part to be the solution of the spherically averaged
KS Hamiltonian. The solution to the spherically averaged KS Hamiltonian is
chosen at some linearization energy (typically in the middle of the valence
band) and the solution $\psi_{l\mu}^{E}\left(r\right)$ as well as
$\dot{\psi}_{l\mu}^{E}\left(r\right)=\frac{\partial}{\partial
E}\psi_{l\mu}^{E}\left(r\right)$ \- its derivative with respect to energy -
are used as a basis set inside the MT spheres. Further accuracy may be
obtained for the HDLO method where
$\ddot{\psi}_{l\mu}^{E}\left(r\right)=\frac{\partial^{2}}{\partial
E^{2}}\psi_{l\mu}^{E}\left(r\right)$ is used to augment the basis set [7, 2].
In this work after a brief review of the usual KKR method we linearize the KKR
method to LKKR and show it to be equivalent to LMTO. We then extend these
ideas to higher derivative versions HDKKR and HDMTO and show that HDKKR basis
set has greater variational power then the HDMTO basis set. However HDKKR is
much harder to implement then HDLMTO and requires much greater computer
resources so often HDMTO is preferred.

##  II Regular KKR method (review)

We now write the KKR basis set wavefunctions:

$\displaystyle
a_{l}^{\mu}\Phi_{lm}^{1\mu}\left(E,\mathbf{r}\right)=Y_{lm}\left(\widehat{\mathbf{r-\mathbf{r}}_{\mu}}\right)\times$
$\displaystyle\times\left\\{\begin{array}[]{cc}\psi_{l\mu}^{E}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)-b_{l}^{\mu}J_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\leq
S^{\mu}\\\
a_{l}^{\mu}K_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|>S^{\mu}\end{array}\right.$
(3)

Here then $1$ in $\Phi_{lm}^{1\mu}\left(E,\mathbf{r}\right)$ is for later
notational use. Here:

$\left[-\frac{d^{2}}{dr^{2}}+\frac{l\left(l+1\right)}{r^{2}}+\bar{V}_{KS}\left(r\right)\right]r\psi_{l\mu}^{E}\left(r\right)=Er\psi_{l\mu}^{E}\left(r\right)$
(4)

and $\bar{V}_{KS}\left(r\right)$ is the spherically average Khon Sham (KS)
potential. Where furthermore:

$\displaystyle\kappa^{2}$ $\displaystyle=E$ $\displaystyle
K_{l}^{\kappa}\left(r\right)$
$\displaystyle=-\kappa^{l+1}\left\\{\begin{array}[]{cc}n_{l}^{\kappa}\left(r\right)&\kappa^{2}>0\\\
n_{l}^{\kappa}\left(r\right)-ij_{l}^{\kappa}\left(r\right)&\kappa^{2}<0\end{array}\right.$
(7) $\displaystyle J_{l}^{\kappa}\left(r\right)$
$\displaystyle=\kappa^{-l}j_{l}^{\kappa}\left(r\right)$ (8)

Here $n_{l}$ and $j_{l}$ are the spherical Neumann and Bessel functions
respectively. Now we look for eigenfunctions of the form:

$\chi\left(E,\mathbf{k},\mathbf{r}\right)=\sum_{\mathbf{R}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\sum_{lm}A_{lm}^{\mu}\left(\mathbf{k},E\right)\Phi_{lm}^{1\mu}\left(E,\mathbf{r}-\mathbf{R}\right)$
(9)

Now we write:

$\displaystyle\sum_{\mathbf{R}+\mathbf{r}_{\nu}\neq\mathbf{r}_{\mu}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)K_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{R}-\mathbf{r}_{\nu}\right|\right)Y_{lm}\left(\widehat{\mathbf{r}-\mathbf{R}-\mathbf{r}_{\nu}}\right)$
$\displaystyle=\sum_{l^{\prime},m^{\prime}}B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)J_{l^{\prime}}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)Y_{l^{\prime}m^{\prime}}\left(\widehat{\mathbf{r}-\mathbf{r}_{\mu}}\right)$
(10)

For the some structure constants
$B_{l^{\prime}m^{\prime};lm}^{\left\\{\mathbf{r}_{\mu}\right\\}}\left(\mathbf{k}\right)$
5 [5]. Now we wish to eliminate all tails in all basis functions that is
Bessel like wavefunction components inside all MT spheres. As such we want for
$\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\leq S^{\mu}$:

$\displaystyle\chi\left(E,\mathbf{k},\mathbf{r}\right)$
$\displaystyle=\sum_{lm}A_{lm}^{\mu}Y_{lm}\left(\widehat{\mathbf{r}-\mathbf{r}_{\mu}}\right)\psi_{l}^{E}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)$
(11)

However we have that:

$\displaystyle\chi\left(E,\mathbf{k},\mathbf{r}\right)$
$\displaystyle=\sum_{l^{\prime}m^{\prime}}A_{l^{\prime}m^{\prime}}^{\mu}\left(\mathbf{k},E\right)Y_{l^{\prime}m^{\prime}}\left(\widehat{\mathbf{r}-\mathbf{r}_{\mu}}\right)a_{l^{\prime}}^{\mu-1}\times$
$\displaystyle\times\left[\psi_{l^{\prime}}^{E}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)-b_{l^{\prime}}^{\mu}J_{l^{\prime}}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)\right]+$
$\displaystyle+\sum_{\nu
lm}A_{lm}^{\nu}\left(\mathbf{k},E\right)\sum_{l^{\prime},m^{\prime}}B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)\times$ $\displaystyle\times
J_{l}^{\kappa}\left(\mathbf{r}-\mathbf{r}_{\mu}\right)Y_{l^{\prime}m^{\prime}}\left(\widehat{\mathbf{r}-\mathbf{r}_{\mu}}\right)$
(12)

Equating the Right Hand Sides of Eqs. (11) and (12) we obtain the relationship

$\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)-\frac{b_{l}^{\mu}}{a_{l}^{\mu}}\delta_{\mu\nu;ll^{\prime};mm^{\prime}}\right]A_{lm}^{\mu}\left(\mathbf{k},E\right)=0$
(13)

Furthermore we wish for the wavefunction to be continuous and continuously
differentiable MT radius, this implies that we want that:

$\displaystyle\left(\psi_{l}\left(S^{\mu}\right),S^{\mu}\frac{\partial}{\partial
S^{\mu}}\psi_{l}\left(S^{\mu}\right)\right)$
$\displaystyle=\left(a_{l}^{\mu},b_{l}^{\mu}\right)\left(\begin{array}[]{cc}K_{l}^{\kappa}\left(S^{\mu}\right)&S^{\mu}\frac{\partial}{\partial
S^{\mu}}K_{l}^{\kappa}\left(S^{\mu}\right)\\\
J_{l}^{\kappa}\left(S^{\mu}\right)&S^{\mu}\frac{\partial}{\partial
S^{\mu}}J_{l}^{\kappa}\left(S^{\mu}\right)\end{array}\right)$ (16)

Solving we get that :

$\frac{b_{l}^{\mu}}{a_{l}^{\mu}}=-\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\psi}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\psi}\left(S^{\mu}\right)}$
(17)

Where

$\mathcal{D}_{f}\left(S^{\mu}\right)=S^{\mu}\frac{f^{\prime}\left(S^{\mu}\right)}{f\left(S^{\mu}\right)}$
(18)

As such we obtain the KKR equation (see Eq. (13)):

$\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)+\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\psi_{\mu}}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\psi_{\mu}}\left(S^{\mu}\right)}\delta_{\nu\mu;ll^{\prime};mm^{\prime}}\right]A_{lm}^{\mu}\left(\mathbf{k},E\right)=0$
(19)

##  III LKKR

We now linearize KKR (LKKR). We write a new set of basis functions:

$\displaystyle\dot{a}_{l}^{\mu}\Phi_{lm}^{2\mu}\left(E,\mathbf{r},S^{\mu}\right)=Y_{lm}\left(\widehat{\mathbf{r-\mathbf{r}}_{\mu}}\right)$
$\displaystyle\left\\{\begin{array}[]{cc}\dot{\psi}_{l\mu}^{E}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)-\dot{b}_{l}^{\mu}J_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\leq
S^{\mu}\\\
\dot{a}_{l}^{\mu}K_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|>S^{\mu}\end{array}\right.$
(22)

For some linearization energy $E$. Furthermore we will assume that

$\displaystyle\Phi_{lm}^{1\mu}\left(E,\mathbf{r},S^{\mu}\right),\>\frac{\partial}{\partial\mathbf{r}}\Phi_{lm}^{1\mu}\left(E,\mathbf{r},S^{\mu}\right),$
$\displaystyle\dot{\Phi}_{lm}^{2\mu}\left(E,\mathbf{r},S^{\mu}\right),\>\frac{\partial}{\partial\mathbf{r}}\dot{\Phi}_{lm}^{2\mu}\left(E,\mathbf{r},S^{\mu}\right)$
(23)

are continuous. We now look for wavefunctions of the form:

$\displaystyle\chi_{lm}^{\mu}\left(E,\mathbf{k},\mathbf{r}\right)$
$\displaystyle=\sum_{\mathbf{R}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\Phi_{lm}^{\mu}\left(E,\mathbf{r}-\mathbf{R},S^{\mu}\right)$
$\displaystyle\dot{\chi}_{lm}^{\mu}\left(E,\mathbf{k},\mathbf{r}\right)$
$\displaystyle=\sum_{\mathbf{R}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\dot{\Phi}_{lm}^{\mu}\left(E,\mathbf{r}-\mathbf{R},S^{\mu}\right)$
(24)

We now define the wavefunctions:

$\displaystyle\chi\left(E,\mathbf{k},\mathbf{r}\right)$
$\displaystyle=\sum_{\mu
lm}A_{lm}^{\mu}\left(\mathbf{k},E\right)\chi_{lm}^{\mu}\left(E_{l}^{\mu},\mathbf{k},\mathbf{r}\right)$
$\displaystyle+\sum_{\mu
lm}B_{lm}^{\mu}\left(\mathbf{k},E\right)\dot{\chi}_{lm}^{\mu}\left(E_{l}^{\mu},\mathbf{k},\mathbf{r}\right)$
(25)

Where furthermore demand:

$\displaystyle\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)+\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\psi_{\mu}}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\psi_{\mu}}\left(S^{\mu}\right)}\delta_{\nu\mu;ll^{\prime};mm^{\prime}}\right]A_{lm}^{\mu}\left(\mathbf{k}\right)$
$\displaystyle+\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)+\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\dot{\psi}_{\mu}}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\dot{\psi}_{\mu}}\left(S^{\mu}\right)}\delta_{\nu\mu;ll^{\prime};mm^{\prime}}\right]B_{lm}^{\mu}\left(\mathbf{k}\right)$
$\displaystyle=0$ (26)

That is all the terms proportional to Bessel and Neumann Functions vanish
inside the spheres $\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\leq S^{\mu}$ (see
Eq. (19). Let us pick a basis of solutions of Eq. (26)

$\mathcal{A}_{1,\mu
lm}^{n}\left(\mathbf{k}\right)=A_{lm}^{\mu}\left(\mathbf{k}\right),\>\mathcal{A}_{2,\mu
lm}^{n}\left(\mathbf{k}\right)=B_{lm}^{\mu}\left(\mathbf{k}\right)$ (27)

We now define the energy matrices to be:

$\displaystyle\bar{H}_{\mathbf{k}}^{m,n}$ $\displaystyle=\sum_{\mu lm,\nu
l^{\prime}m^{\prime}}\mathcal{A}_{i,\mu
lm}^{m*}\left(\mathbf{k}\right)\left[-\bar{\Delta}^{\mathbf{k}}+\bar{V}_{KS}^{\mathbf{k}}\right]_{\mu
lm,\nu l^{\prime}m^{\prime}}^{i,j}\mathcal{A}_{j,\nu
l^{\prime}m^{\prime}}^{n}\left(\mathbf{k}\right)$
$\displaystyle\bar{O}_{\mathbf{k}}^{m,n}$ $\displaystyle=\sum_{\mu lm,\nu
l^{\prime}m^{\prime}}\mathcal{A}_{i,\mu
lm}^{m*}\left(\mathbf{k}\right)\left[\bar{O}^{\mathbf{k}}\right]_{\mu lm,\nu
l^{\prime}m^{\prime}}^{i,j}\mathcal{A}_{j,\nu
l^{\prime}m^{\prime}}^{n}\left(\mathbf{k}\right)$ (28)

We now obtain the KS equation:

$\sum_{n}\bar{H}_{\mathbf{k}}^{m,n}V_{n}^{\mathfrak{A}}=\varepsilon^{\mathfrak{A}}\sum_{n}\bar{O}_{\mathbf{k}}^{m,n}V_{n}^{\mathfrak{A}}$
(29)

We have as such setup a LKKR Khon Sham like system for solving many body DFT
problems. Here $O$ is the overlap matrix.

##  IV Equivalence of LKKR to LMTO

###  IV.1 LMTO review

We note that the regular LMTO method is based on the following wavefunction \-
which is a sum of three different wavefunctions given by:

$\displaystyle\Theta_{lm}^{\mu}\left(E_{l}^{\mu},\mathbf{r}\right)=\Theta_{lm}^{1\mu}\left(E,\mathbf{r}\right)+\Theta_{lm}^{2\mu}\left(E,\mathbf{r}\right)+\Theta_{lm}^{3\mu}\left(E,\mathbf{r}\right)$
$\displaystyle\Theta_{lm}^{1\mu}\left(E,\mathbf{r}\right)=\left\\{\begin{array}[]{cc}Y_{lm}\left(\widehat{\mathbf{r-\mathbf{r}}_{\mu}}\right)\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}\Psi_{l}^{\mu}\left(\mathcal{D}_{K},\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\leq
S^{\mu}\\\ 0&otherwise\end{array}\right.$ (32)
$\displaystyle\Theta_{lm}^{2\mu}\left(E,\mathbf{r}\right)=\left\\{\begin{array}[]{cc}\sum_{l^{\prime}m^{\prime}}Y_{l^{\prime}m^{\prime}}\left(\widehat{\mathbf{r-\mathbf{r}}_{\nu}}\right)\frac{J_{l^{\prime}}^{\kappa}\left(S^{\nu}\right)}{\Psi_{l^{\prime}}\left(\mathcal{D}_{J},S^{\nu}\right)}B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)\Psi_{l^{\prime}}^{\nu}\left(\mathcal{D}_{J},\left|\mathbf{r-\mathbf{r}}_{\nu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\nu}\right|\leq
S^{\nu}\\\ 0&otherwise\end{array}\right.$ (35)
$\displaystyle\Theta_{lm}^{3\mu}\left(E,\mathbf{r}\right)=\left\\{\begin{array}[]{cc}Y_{lm}\left(\widehat{\mathbf{r-\mathbf{r}}_{\mu}}\right)K_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\nu}-\mathbf{R}\right|\geq
S^{\nu},\forall\nu,\,\mathbf{R}\\\ 0&otherwise\end{array}\right.$ (38)

Where

$\displaystyle\Psi_{l}^{\mu}\left(\mathcal{D},\mathbf{r}\right)$
$\displaystyle=\psi_{l\mu}^{E}\left(r\right)+\omega_{l}\left(\mathcal{D}\right)\dot{\psi_{l\mu}^{E}}$
$\displaystyle\omega_{l}\left(\mathcal{D}\right)$
$\displaystyle=-\frac{\psi_{l\mu}\left(S^{\mu}\right)}{\dot{\psi_{l\mu}}\left(S^{\mu}\right)}\cdot\frac{\mathcal{D-\mathcal{D}_{\psi}}}{\mathcal{D}-\mathcal{D}_{\dot{\psi}}}$
(39)

and

$\chi_{\mu,lm}^{MTO}\left(\mathbf{k}\right)=\sum_{\mathbf{R}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\Theta_{lm}^{\mu}\left(E_{l}^{\mu},\mathbf{r}-\mathbf{R},\mathbf{k}\right)$
(40)

is the LMTO basis set

###  IV.2 Equivalence between LKKR and LMTO calculation

We wish to show that LKKR and LMTO are equivalent. To do se we now need to
check if:

$\displaystyle\chi_{\mu,lm}^{MTO}\left(\mathbf{k}\right)$ $\displaystyle\in
Span\left\\{\chi_{l^{\prime}m^{\prime}}^{\mu}\left(E,\mathbf{r}\right),\dot{\chi}_{l^{\prime}m^{\prime}}^{\mu}\left(E,\mathbf{r}\right)\right\\}$
$\displaystyle\Rightarrow\chi_{\mu,lm}^{MTO}\left(\mathbf{k}\right)$
$\displaystyle=\sum_{\nu l^{\prime}m^{\prime}}A_{\nu
l^{\prime}m^{\prime}}^{lm\Theta}\left(\mathbf{k}\right)\chi_{l^{\prime}m^{\prime}}^{\nu}\left(E,\mathbf{r}\right)+$
$\displaystyle+\sum_{l^{\prime}m^{\prime}}B_{\nu
l^{\prime}m^{\prime}}^{lm\Theta}\left(\mathbf{k}\right)\dot{\chi}_{l^{\prime}m^{\prime}}^{\nu}\left(E,\mathbf{r}\right)$
(41)

for some $A_{\nu l^{\prime}m^{\prime}}^{lm\Theta}\left(\mathbf{k}\right)$ and
$B_{\nu l^{\prime}m^{\prime}}^{lm\Theta}\left(\mathbf{k}\right)$. As such we
would see that the LMTO method is spanned by the LKKR method, however by
counting the total number of basis set elements we would then see the two are
equivalent. However Eq. (41) is a direct check - see Appendix A.

##  V HDKKR

We now write the HDKKR wavefunctions. We write the following basis set
wavefunctions:

$\displaystyle\ddot{a}_{l}^{\mu}\Phi_{lm}^{3\mu}\left(E,\mathbf{r},S^{\mu}\right)=Y_{lm}\left(\widehat{\mathbf{r-\mathbf{r}}_{\mu}}\right)\times$
$\displaystyle\times\left\\{\begin{array}[]{cc}\ddot{\psi}_{l\mu}^{E}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)-\ddot{b}_{l}^{\mu}J_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\leq
S^{\mu}\\\
\ddot{a}_{l}^{\mu}K_{l}^{\kappa}\left(\left|\mathbf{r}-\mathbf{r}_{\mu}\right|\right)&\left|\mathbf{r}-\mathbf{r}_{\mu}\right|>S^{\mu}\end{array}\right.$
(44)

For some linearization energy $E$. and

$\Phi_{lm}^{i\mu}\left(E,\mathbf{r},S^{\mu}\right),\>\frac{\partial}{\partial\mathbf{r}}\Phi_{lm}^{i\mu}\left(E,\mathbf{r},S^{\mu}\right);\>i=1,2,3$
(45)

Are continuous. We now define the wavefunctions:

$\chi_{lm}^{i\mu}\left(E,\mathbf{k},\mathbf{r}\right)=\sum_{\mathbf{R}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\Phi_{lm}^{i\mu}\left(E_{l}^{\mu},\mathbf{r}-\mathbf{R},S^{\mu}\right)$
(46)

for $i=1,2,3$. We now look for wavefunctions of the form:

$\chi\left(E,\mathbf{k},\mathbf{r}\right)=\sum_{i\mu
lm}A_{lm}^{i\mu}\left(\mathbf{k}\right)\chi_{lm}^{i\mu}\left(E_{l}^{\mu},\mathbf{k},\mathbf{r}\right)$
(47)

Where furthermore demand:

$\displaystyle\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)+\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\psi_{\mu}}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\psi_{\mu}}\left(S^{\mu}\right)}\delta_{\nu\mu;ll^{\prime};mm^{\prime}}\right]A_{lm}^{1\mu}\left(\mathbf{k}\right)$
$\displaystyle+\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)+\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\dot{\psi}_{\mu}}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\dot{\psi}_{\mu}}\left(S^{\mu}\right)}\delta_{\nu\mu;ll^{\prime};mm^{\prime}}\right]A_{lm}^{2\mu}\left(\mathbf{k}\right)$
$\displaystyle+\sum_{\mu lm}\left[B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)+\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{J_{l}^{\kappa}\left(S^{\mu}\right)}\times\frac{\mathcal{D}_{K}\left(S^{\mu}\right)-\mathcal{D}_{\ddot{\psi}_{\mu}}\left(S^{\mu}\right)}{\mathcal{D}_{J}\left(S^{\mu}\right)-\mathcal{D}_{\ddot{\psi}_{\mu}}\left(S^{\mu}\right)}\delta_{\nu\mu;ll^{\prime};mm^{\prime}}\right]A_{lm}^{3\mu}\left(\mathbf{k}\right)$
$\displaystyle=0$ (48)

So that Bessel like terms vanish inside the MT spheres (see Eqs. (26) and
(19)). Let us pick a basis of solutions of Eq. (48):

$\mathcal{A}_{i,\mu
lm}^{n}\left(\mathbf{k}\right)=A_{lm}^{i\mu}\left(\mathbf{k}\right)$ (49)

We now use Eqs. (28) and (29) to complete calculations.

##  VI HDMTO

The HDMTO basis set is based on the same wavefunctions as in Eq. (38) where
however:

$\Psi_{l}^{\mu}\left(\mathcal{D},\mathbf{r}\right)=\psi_{l\mu}^{E}\left(r\right)+\dot{\omega}_{l}^{\mu}\left(\mathcal{D},\mathcal{D}^{2}\right)\dot{\psi_{l\mu}^{E}}+\ddot{\omega}_{l}^{\mu}\left(\mathcal{D},\mathcal{D}^{2}\right)\ddot{\psi}_{l\mu}^{E}$
(50)

Where

$\mathcal{D}_{f}\left(S\right)=S^{2}\frac{f"\left(S\right)}{f\left(S\right)}$
(51)

We want that:

$\displaystyle\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}\left(\begin{array}[]{c}\Psi_{l}\left(S^{\mu}\right)\\\
S^{\mu}\frac{\partial}{\partial S^{\mu}}\Psi_{l}\left(S^{\mu}\right)\\\
\left(S^{\mu}\right)^{2}\frac{\partial^{2}}{\partial\left(S^{\mu}\right)^{2}}\Psi_{l}\left(S^{\mu}\right)\end{array}\right)$
(55)
$\displaystyle=\left(\begin{array}[]{c}K_{l}^{\kappa}\left(S^{\mu}\right)\\\
S^{\mu}\frac{\partial}{\partial S^{\mu}}K_{l}^{\kappa}\left(S^{\mu}\right)\\\
\left(S^{\mu}\right)^{2}\frac{\partial^{2}}{\partial\left(S^{\mu}\right)^{2}}K_{l}^{\kappa}\left(S^{\mu}\right)\end{array}\right)$
(59)

so that the wavefunction along with its first and second derivatives are
continuous. The derivation is identical to the one in Appendix A. Now we
write:

$\chi_{\mu,lm}^{MTO}\left(\mathbf{k}\right)=\sum_{\mathbf{R}}\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\Theta_{lm}^{\mu}\left(E,\mathbf{r}-\mathbf{R},\mathbf{k}\right)$
(60)

And study this basis set.

###  VI.1 Differences and similarities between HDMTO and HDKKR

We note that HDKKR basis set includes HDMTO as HDKKR only demands continuity
of the wavefunction and its derivative while HDMTO also demands continuity of
the second derivative. Furthermore because of the similarities of the form of
the HDMTO and LMTO basis wavefunctions the HDMTO basis set eliminates all
Bessel like wavefunction components inside all MT spheres so HDMTO basis
wavefunctions are a type of HDKKR basis wavefunctions. As such HDKKR has
greater variational power then HDMTO. We note that this observation is true,
despite the fact that the exact solution of the KS Hamiltonian have all
derivatives continuous, as wavefunctions with discontinuous second derivatives
may help approximate ones with continuous second derivatives better. As the KS
Equations are variational HDKKR gives a better estimate of the ground state
energy then HDMTO. However the HDMTO method is much more efficient as the
basis is smaller and it does not require solving Eq. (48) on computer as the
basis set automatically satisfies it as there are no Bessel like functions
inside any MT spheres so in many case HDMTO is preferable to HDKKR.

##  VII Conclusion

In this work we have studied a linearized version of KKR (LKKR) and shown it
to be exactly equivalent to LMTO. We have extended these ideas to higher
derivative HDKKR and HDMTO basis sets. These two basis sets are not
equivalent. The HDKKR basis set has greater variational power then HDMTO basis
set as it includes it as a subset. However HDMTO is much easier to work with,
requires significantly less computational power and no auxiliary equations to
solve. As such in many cases, since HDMTO is easier to implement, it is likely
the preferred basis set method. In the future it would be of interest to set
up practical HDMTO calculations for real solid state crystal systems.

##  Appendix A Technical Calculations

As such for the interstitial region to match between LKKR and LMTO we must
have that:

$A_{\nu l^{\prime}m^{\prime}}^{lm\Theta}\left(\mathbf{k}\right)+B_{\nu
l^{\prime}m^{\prime}}^{lm\Theta}\left(\mathbf{k}\right)=\delta_{\mu\nu}\delta_{l,l^{\prime};m,m^{\prime}}$
(61)

Now we match the MT region. From which we read of:

$\displaystyle A_{\nu
l^{\prime}m^{\prime}}^{lm\mu\Theta}\left(\mathbf{k}\right)$
$\displaystyle=\delta_{\mu,\nu;l,l^{\prime};m,m^{\prime}}\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}a_{l}^{\mu}+\frac{J_{l^{\prime}}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l^{\prime}}\left(\mathcal{D}_{J},S^{\mu}\right)}B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)a_{l^{\prime}}^{\nu}$
$\displaystyle B_{\nu
l^{\prime}m^{\prime}}^{lm\mu\Theta}\left(\mathbf{k}\right)$
$\displaystyle=\delta_{\mu,\nu;l,l^{\prime};m,m^{\prime}}\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}\omega_{l}\left(\mathcal{D}_{K}\right)\dot{a}_{l}^{\mu}+\frac{J_{l^{\prime}}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l^{\prime}}\left(\mathcal{D}_{J},S^{\mu}\right)}B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)\omega_{l^{\prime}}\left(\mathcal{D}_{J}\right)\dot{a}_{l^{\prime}}^{\nu}$
(62)

So we want:

$\delta_{\mu,\nu;l,l^{\prime};m,m^{\prime}}\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}\left[a_{l}^{\mu}+\omega_{l}\left(\mathcal{D}_{K}\right)\dot{a}_{l}^{\mu}\right]+B_{\nu\mu
l^{\prime}m^{\prime};lm}\left(\mathbf{k}\right)\frac{J_{l^{\prime}}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l^{\prime}}\left(\mathcal{D}_{J},S^{\mu}\right)}\left[a_{l^{\prime}}^{\nu}+\omega_{l^{\prime}}\left(\mathcal{D}_{J}\right)\dot{a}_{l^{\prime}}^{\nu}\right]=\delta_{\mu\nu}\delta_{l,l^{\prime};m,m^{\prime}}$

Now we have that:

$\displaystyle
a_{l}^{\mu}+\omega_{l}\left(\mathcal{D}_{J}\right)\dot{a}_{l}^{\mu}$
$\displaystyle=a_{l}^{\mu}-\frac{\psi_{l\mu}\left(S^{\mu}\right)}{\dot{\psi_{l\mu}}\left(S^{\mu}\right)}\cdot\frac{\mathcal{D}_{J}\mathcal{-\mathcal{D}_{\psi}}}{\mathcal{D}_{J}-\mathcal{D}_{\dot{\psi}}}\dot{a}_{l}^{\mu}$
$\displaystyle=S^{\mu}\psi_{l\mu}\left(S^{\mu}\right)J_{l}^{\kappa}\left(S^{\mu}\right)\left(\mathcal{D}_{J}-\mathcal{D}_{\psi}\right)-\frac{\psi_{l\mu}\left(S^{\mu}\right)}{\dot{\psi_{l\mu}}\left(S^{\mu}\right)}\cdot\frac{\mathcal{D}_{J}\mathcal{-\mathcal{D}_{\psi}}}{\mathcal{D}_{J}-\mathcal{D}_{\dot{\psi}}}S^{\mu}\dot{\psi}_{l\mu}\left(S^{\mu}\right)J_{l}^{\kappa}\left(S^{\mu}\right)\left(\mathcal{D}_{J}-\mathcal{D}_{\dot{\psi}}\right)=0$
(63)

As such Eq. (61) is satisfied. Next we wish to show that the continuity
equation is satisfied for LMTO basis wavefunctions. Indeed:

$\displaystyle\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}\left[a_{l}^{\mu}+\omega_{l}\left(\mathcal{D}_{K}\right)\dot{a}_{l}^{\mu}\right]$
$\displaystyle=\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\Psi_{l}\left(\mathcal{D}_{K},S^{\mu}\right)}\left[a_{l}^{\mu}+\frac{\psi_{l\mu}\left(S^{\mu}\right)}{\dot{\psi_{l\mu}}\left(S^{\mu}\right)}\cdot\frac{\mathcal{D}_{K}\mathcal{-\mathcal{D}_{\psi}}}{\mathcal{D}_{K}-\mathcal{D}_{\dot{\psi}}}\dot{a}_{l}^{\mu}\right]$
$\displaystyle=\frac{K_{l}^{\kappa}\left(S^{\mu}\right)}{\psi_{l\mu}\left(S^{\mu}\right)\left[1-\frac{\mathcal{D}_{K}-\mathcal{D}_{\psi}}{\mathcal{D}_{K}-\mathcal{D}_{\dot{\psi}}}\right]}\left[S^{\mu}\psi_{l\mu}\left(S^{\mu}\right)J_{l}^{\kappa}\left(S^{\mu}\right)\left(\mathcal{D}_{J}-\mathcal{D}_{\psi}\right)+\frac{\psi_{l\mu}\left(S^{\mu}\right)}{\dot{\psi_{l\mu}}\left(S^{\mu}\right)}\cdot\frac{\mathcal{D}_{K}\mathcal{-\mathcal{D}_{\psi}}}{\mathcal{D}_{K}-\mathcal{D}_{\dot{\psi}}}\left[S^{\mu}\dot{\psi}_{l}\left(S^{\mu}\right)J_{l}^{\kappa}\left(S^{\mu}\right)\left(\mathcal{D}_{J}-\mathcal{D}_{\cdot\psi}\right)\right]\right]$
$\displaystyle=\frac{S^{\mu}K_{l}^{\kappa}\left(S^{\mu}\right)J_{l}^{\kappa}\left(S^{\mu}\right)}{\mathcal{D}_{\psi}-\mathcal{D}_{\dot{\psi}}}\left[\left[\mathcal{D}_{K}-\mathcal{D}_{\dot{\psi}}\right]\left(\mathcal{D}_{J}-\mathcal{D}_{\psi}\right)-\left[\mathcal{D}_{K}\mathcal{-\mathcal{D}_{\psi}}\right]\left(\mathcal{D}_{J}-\mathcal{D}_{\dot{\psi}}\right)\right]$
$\displaystyle=S^{\mu}K_{l}^{\kappa}\left(S^{\mu}\right)J_{l}^{\kappa}\left(S^{\mu}\right)\left[\mathcal{D}_{K}-\mathcal{D}_{J}\right]=1$
(64)

As such Eq. (26) is verified for the LMTO basis set as it is a continuous
continuously differentiable wavefunction in the span of the KKR basis set with
no components of the form of Bessel functions inside the MT spheres.

## References

  * 1 [2005] D. J. Singh and L. Nordstrom, Planewaves, Pesudopotentials and the LAPW Method (Springer, New York, 2005).
  * 2 [2020] R. M. Martin, Electronic Structures Basic Theory and Practical Methods (Cambridge University Press, Cambridge, 2020).
  * 3 [1975] O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
  * 4 [2009] D. Marx and J. Hutter, Ab Initio Molecular Dynamics Basic Theory and Advanced Methods (Cambridge University Press, Cambridge, 2009).
  * 5 [2010] J. M. Wills, M. Alouani, P. Anderson, A. Dellin, O. Eriksson, and O. Grechnyev, Full-Potential Electronic Structure Method Energy and Force Calculations with Density Functional Theory and Dynamical Mean Field Theory (Springer, New York, 2010).
  * 6 [1984] H. L. Skriver, The LMTO method Muffin-Tin Orbitals and Electronic Structure (Springer, New York, 1984).
  * 7 [2013] G. Michalicek, M. Betzinger, C. Freidrich and S. Blugel, Comp. Phys. Comm. 184, 2670 (2013).