11institutetext: UniversitÃ degli studi di Milano, Italy 22institutetext: INFN
sezione di Milano, Milano, Italy 33institutetext: UniversitÃ Paris CitÃ, CNRS,
Astroparticule et Cosmologie, Paris, France 44institutetext: UniversitÃ PSL,
Observatoire de Paris, AstroParticule et Cosmologie, Paris, France
55institutetext: Waterloo Centre for Astrophysics, University of Waterloo,
Waterloo, Canada 66institutetext: Department of Physics and Astronomy,
University of Waterloo, Waterloo, Canada 77institutetext: UniversitÃ di Roma -
La Sapienza, Italy 88institutetext: UniversitÃ di Milano â Bicocca, Italy
99institutetext: Instituto Argentino de RadioastronomÃa, Buenos Aires,
Argentina 1010institutetext: Consejo Nacional de Investigaciones CientÃficas y
TÃcnicas, Ciudad de Buenos Aires, Argentina 1111institutetext: Facultad de
Ciencias AstronÃmicas y GeofÃsicas, Universidad Nacional de La Plata,
Argentina

# The advantage of Bolometric Interferometry for controlling Galactic
foreground contamination in CMB primordial B-modes measurements

E. Manzan 1122<EMAIL_ADDRESS>M. Regnier 33 J-Ch. Hamilton 33 A.
Mennella 1122 J. Errard 33 L. Zapelli 1122 S.A. Torchinsky 3344 S. Paradiso
5566 E. Battistelli 77 M. Bersanelli 1122 P. De Bernardis 77 M. De Petris 77
G. DâAlessandro 77 M. Gervasi 88 S. Masi 77 M. Piat 33 E. Rasztocky 99 G.E
Romero 99 C.G. Scoccola 1010 11 11 M. Zannoni 88 the QUBIC Collaboration

###### Abstract

In the quest for the faint primordial B-mode polarization of the Cosmic
Microwave Background, three are the key requirements for any present or future
experiment: an utmost sensitivity, excellent control over instrumental
systematic effects and over Galactic foreground contamination.

Bolometric Interferometry (BI) is a novel technique that matches them all by
combining the sensitivity of bolometric detectors, the control of instrumental
systematics from interferometry and a software-based, tunable, in-band
spectral resolution due to its ability to perform band-splitting during data
analysis (spectral imaging).

In this paper, we investigate how the spectral imaging capability of BI can
help in detecting residual contamination in case an over-simplified model of
foreground emission is assumed in the analysis. To mimic this situation, we
focus on the next generation of ground-based CMB experiment, CMB-S4, and
compare its anticipated sensitivities, frequency and sky coverage with a
hypothetical version of the same experiment based on BI, CMB-S4/BI, assuming
that line-of-sight (LOS) frequency decorrelation is present in dust emission
but is not accounted for during component separation.

We show results from a Monte-Carlo analysis based on a parametric component
separation method (FGBuster), highlighting how BI has the potential to
diagnose the presence of foreground residuals in estimates of the tensor-to-
scalar ratio $r$ in the case of unaccounted Galactic dust LOS frequency
decorrelation.

##  1 Introduction

As the current best upper limit on the tensor-to-scalar ratio $r$ of the CMB
primordial B-modes tightens to $r<0.032$ Tristram_2022 , measuring such a
faint signal requires instruments with utmost sensitivity, excellent control
over instrumental systematic effects and over foreground contamination. The
latter is achieved through multifrequency observations, to exploit the
different spectral behaviour of the blackbody CMB emission with respect to the
Galactic synchrotron and dust grains emissions, which are conventionally
assumed to scale over frequency as a power-law and a modified blackbody (MBB),
respectively. However, since both synchrotron and thermal dust are brighter
than the CMB emission in polarization, modern CMB experiments are also
relying, now more than ever, on improved foreground models to mitigate as much
as possible the residuals due to improper or overly simplified modeling.

Indeed, over the years many models have been implemented in PySM Thorne_2017
(Python Sky Model) to account for additional complexities in the Galactic dust
emission, such as different dust grain compositions Hensley_2017 , the
superimposition of coherently aligned molecular clouds with different
temperatures and spectral indices Finkbeiner_1999 ; Martinez_Solaeche_2018  or
line-of-sight (LOS) frequency decorrelations caused by a frequency-dependent
polarization angle due to misalignments of various molecular clouds along the
LOS, following the statistical approach described in Vansyngel_2018 .

In this paper, we investigate how the spectral imaging capability of an
unconventional technique for CMB polarimetry, Bolometric Interferometry (BI)
2020.QUBIC.PAPER1 , can help in diagnosing the presence of foreground
residuals in case of improper or overly simplified foreground modeling.
Precisely, we assume that LOS frequency decorrelations are present in the dust
emission but are not accounted for during component separation and we compare
the estimated tensor-to-scalar ratio, $r$, after a Monte-Carlo analysis based
on a parametric component separation, in the case of the next generation of
ground-based CMB experiment, CMB-S4 Abazajian_2022 , and of a hypothetical
version of the same experiment based on BI, CMB-S4/BI.

##  2 Bolometric Interferometry in a nutshell

Bolometric Interferometry combines the use of bolometers, which are state-of-
the-art cryogenic broadband detectors providing high sensitivity, with the
control of instrumental systematic effects typical of interferometry
1981MNRAS.196.1067C ; Bigot-Sazy2013 . A schematic of BI is shown in the left
panel of Fig. 1: an array of antennas receives the sky signal and re-emits it
through a back array of twin antennas towards an optical combiner that focuses
the radiation onto the focal plane, conversely to a traditional imager where
each antenna is coupled to a single detector instead. Therefore, the signal
received by each detector is the sky signal convolved with the joint angular
response of all the antennas simultaneously, which we call synthesized beam,
and exhibits a multiple-peaked shape as shown in the right panel of Fig. 1.

The interested reader can find more details on BI and on its current state-of-
the-art, represented by the QUBIC experiment, in 2020.QUBIC.PAPER1 ;
2020.QUBIC.PAPER2 ; 2020.QUBIC.PAPER3 ; 2020.QUBIC.PAPER4 ; 2020.QUBIC.PAPER5
; 2020.QUBIC.PAPER6 ; 2020.QUBIC.PAPER7 ; 2020.QUBIC.PAPER8 .

Figure 1: Left panel: schematic of a BI instrument. The sky signal is received
by an array of back-to-back antennas and re-imaged onto the bolometric focal
planes where the field interferes additively. A polarizer and a rotating half-
wave plate make the instrument sensitive to linear polarization. Right panel:
azimuth cut of the synthesized beam (solid lines) at 131 GHz (red line) and at
169 GHz (blue line) for a detector at the center of the focal plane. Dashed
lines represent the beam pattern of a single feed-horn antenna. One can
appreciate the frequency-dependent position of the secondary peaks.

###  2.1 Spectral Imaging

The width of the peaks in the synthesized beam, $\theta_{\mathrm{FWHM}}$, and
their angular distance, $\Theta$, depend on the signal wavelength, $\lambda$,
on the number of antennas along the maximum axis of the antenna array, $P$,
and on the separation between two consecutive antennas, $\Delta h$, as
2020.QUBIC.PAPER2 :

$\theta_{\mathrm{FWHM}}=\frac{\lambda}{(P-1)\Delta
h},\,\,\,\,\,\,\,\,\,\,\Theta=\frac{\lambda}{\Delta h}.$ (1)

The linear dependence of the angular separation $\Theta$ on the wavelength
$\lambda$ causes the position of the secondary peaks to significantly move as
a function of the frequency, as shown in the right panel of Fig. 1. Therefore,
during data analysis one can recover the sky signal for two
frequencies$,\nu_{1}$ and $\nu_{2}$, within the instrument frequency band as
long as the secondary peaks are well-resolved. This occurs if
$\Theta(\nu_{2})-\Theta(\nu_{1})>\theta_{\mathrm{FWHM}}(\sqrt{\nu_{1}\nu_{2}})$,
namely if $\frac{\Delta\nu}{\nu}\geq\frac{1}{P-1}$.

We call this technique spectral imaging and, since it happens at the data
analysis level, it allows us to tune the spectral resolution and also to re-
analize the data with various spectral configurations to search for biases.

##  3 Simulation set-up

###  3.1 Sky models

In our simulation we consider three cases of dust emissions: i. a MBB with
spatially varying temperature, $T_{d}(\hat{n})$, and spectral index,
$\beta_{d}(\hat{n})$, modelled using the PySM d1 model as

$I_{d}(\hat{n},\nu)\ \bigg{\rvert}_{MBB}=A_{d,\nu_{0}}(\hat{n})\
\frac{B_{\nu}\left(T_{d}(\hat{n})\right)}{B_{\nu_{0}}\left(T_{d}(\hat{n})\right)}\left(\frac{\nu}{\nu_{0}}\right)^{\beta_{d}(\hat{n})},$
(2)

ii. a MBB with constant temperature, $T_{d}$, and spectral index, $\beta_{d}$,
derived from Equation 2 using the d0 model and iii. a deviation from a single
MBB that accounts for LOS frequency decorrelation, computed using the d6 model
Vansyngel_2018  by multiplying the d1 emission in Equation 2 for a
decorrelation factor, $D$, as

$I_{d}(\hat{n},\nu)=D(\nu,\nu_{0},\ell_{corr})\ I_{d}(\hat{n},\nu)\
\bigg{\rvert}_{MBB}$ (3)

The decorrelation factor, $D(\nu,\nu_{0},\ell_{corr})$, is a function of the
simulated frequency $\nu$, a reference frequency $\nu_{0}$ for which
$D(\nu=\nu_{0})=1$, and of the correlation length $\ell_{corr}$, a parameter
that quantifies the deviation from a single MBB. For each simulated frequency
$\nu$, the decorrelation factor $D$ is randomly sampled from a Gaussian
distribution with center $\mu=1$ and standard deviation $\sigma$ that scales
as the inverse of the correlation length, $\sigma=\sigma(1/\ell_{corr})$.
Therefore, a smaller correlation length causes larger decorrelation, and
viceversa, as it is shown in Fig. 2. This model mimics the presence of a
frequency-varying polarization angle caused by some misalignments of the
magnetic field (and therefore of the molecular clouds) along the LOS, without
making any physical assumption on the misalignment itself.

For the synchrotron, we use the s0 and s1 models, which simulate a power-law
with spatially constant or varying spectral index, $\beta_{s}$, respectively.

Figure 2: Dispersion of the dust SED for different correlation lengths of the
PySM d6 model normalized by the single MBB emission (d1 model). The colored
areas represent the statistical deviation from a MBB for a given correlation
length, evaluated over 500 realizations. As one can see, smaller correlation
lengths $\ell_{corr}$ cause larger decorrelation. Moreover, for a given
$\ell_{corr}$, the decorrelation decreases as $\nu$ gets closer to $\nu_{0}$,
where $\nu_{0}=353$ GHz for polarization, as described in Vansyngel_2018 .

###  3.2 Instrumental configurations

To simulate the next generation of ground-based CMB experiment, CMB-S4, we
consider its anticipated nine frequency channels, $\Delta\nu_{i}$,
sensitivities, $\sigma_{i}$, and 3$\%$ circular sky patch following
Abazajian_2022 . To simulate a BI version of CMB-S4, we assume to perform
spectral imaging during data analysis on the six highest frequency channels
where the dust emission dominates over the synchrotron. Precisely, we split
each frequency band $\Delta\nu_{i}$ into a number of sub-bands
$n_{\mathrm{sub}}$ and correspondingly increase the white noise level,
$\sigma^{\mathrm{BI}}_{i}$, as

$\Delta\nu^{\mathrm{BI}}_{i}=\frac{\Delta\nu_{i}}{n_{\mathrm{sub}}},\,\,\,\,\,\,\,\,\,\,\sigma^{\mathrm{BI}}_{i}=\sigma_{i}\times\sqrt{n_{\mathrm{sub}}}\times\varepsilon\quad\mathrm{for}\
i=4,...,9.$ (4)

The sub-optimality factor, $\varepsilon\in[1.2,1.6]$, acts as a proxy that
increases the white noise level to account for the BI correlated noise, as
described in 2020.QUBIC.PAPER2 , and $n_{\mathrm{sub}}$ ranges from 2 to 8,
for a total of seven simulated CMB-S4/BI cases. The comparison between the
CMB-S4 sensitivity and three cases of CMB-S4/BI is shown in Fig. 3, along with
the simulated sky patch.

Figure 3: Left Panel: Polarization sensitivity of CMB-S4 and three examples of
CMB-S4/BI, with $n_{\mathrm{sub}}=3,5,7$ respectively. The three lowest
frequency bands in grey are identical for all the instruments because we
choose not to split the synchrotron dominated frequency channels. Right Panel:
The anticipated CMB-S4 sky patch (white patch), compared to the dust emission
map at 150 GHz plotted underneath.

###  3.3 Simulation pipeline

Our simulation pipeline consists of 500 iterations of a Monte-Carlo chain. For
each iteration, we generate a CMB realization, a noise realization for each
frequency channel of each instrumental configuration and the foreground maps
in one of the three cases: d0s0, d1s1 or d6s1 model. We apply a parametric
component separation method to the input maps, assuming the dust model to be a
MBB independently of the true input dust emission. In the d0s0 and d1s1 case,
this corresponds to fitting with the correct model, whereas in the d6s1 it
mimics the situation in which our model is overly simplified due to our lack
of knowledge of the details of the dust emission. We compute the cross-power
spectra from each of the 500 reconstructed CMB maps and the corresponding
(Gaussian) likelihood on $r$, from which we extract the maximum-likelihood
value of $r$.

In this paper we show the results obtained using the FGBuster component
separation code Stompor . However, we repeated the study with a different
parametric code, Commander eriksen_2006 , and found coherent results. The
interested reader can find more details in regnier2023identifying .

##  4 Results

The average maximum likelihood value of $r$ and standard deviation as a
function of the number of sub-bands is shown in Fig. 4. A classical imager,
like CMB-S4, would only recover the case $n_{\mathrm{sub}}=1$. In the left
panel, our results show that when we fit for the correct model (d0 or d1), we
expect the result from BI to be independent of $n_{\mathrm{sub}}$111The small
bias is due to an $E\rightarrow B$ modes leakage caused by the power spectra
computation on a sky patch, where the spherical harmonics are no longer
orthogonal. This bias could be mitigated by increasing the apodization radius
of the mask at the expense of a smaller effective sky fraction ($<3\%$), but
this optimization is outside the scope of the study.. Instead, when dust LOS
frequency decorrelations are present but not accounted for during component
separation, a bolometric interferometer would recover a decreasing value of
$r$ as a function of $n_{\mathrm{sub}}$: increasing the number of sub-bands
therefore allows us to reduce the bias with respect to a traditional imager.
Moreover, the decreasing trend as a function of $n_{\mathrm{sub}}$, conversely
to the expected constant trend, is itself a direct hint of the contamination
from foreground residuals. This advantage of BI over a classical imager is
mantained also for smaller levels of decorrelation, as shown in the right
panel of Fig. 4.

Figure 4: Left Panel: Average maximum-likelihood value of $r$ and standard
deviation as a function of the number of sub-bands, computed from the
distribution shown vertically as a violin plot. The cases where we fit for the
correct model are shown in purple, orange and green (d0 or d1, with input
$r=0$ or $r=0.006$). The cyan curve shows the result when LOS frequency
decorrelations are present in the dust emission (model d6s1 with
$\ell_{\mathrm{corr}}=10$) but not accounted for in the analysis. Right Panel:
Summary of the result on $r$ for all the simulated foreground models (d0s0,
d1s1 and several $\ell_{\mathrm{corr}}$ cases of d6s1) with input $r=0$. For
ease of reading, only CMB-S4 and CMB-S4/BI with $n_{\mathrm{sub}}=3,5,7$ are
shown.

##  5 Conclusions

In this paper we have shown that the spectral imaging capability of Bolometric
Interferometry can help diagnosing the presence of foreground residuals in
estimates of the tensor-to-scalar ratio $r$ when LOS frequency decorrelations
are present in the dust emission but not properly accounted for during
component separation, whereas an imager of similar performance could be
biased.

Although neglecting the impact of instrumental systematic effects, this study
illustrates the potential of BI in the context of future CMB polarization
experiments that will be challenged by complex foregrounds in their quest for
B-modes detection.

## References

  * (1) M. Tristram, et al., Physical Review D 105, 083524 (2022)
  * (2) B. Thorne, et al., MNRAS 469, 2821 (2017)
  * (3) B.S. Hensley, B.T. Draine, APJ 836, 179 (2017)
  * (4) D.P. Finkbeiner, M. Davis, D.J. Schlegel, APJ 524, 867 (1999)
  * (5) G. Martínez-Solaeche, et al., MNRAS 476, 1310 (2018)
  * (6) F. Vansyngel, et al., A&A 618, C4 (2018)
  * (7) J.C. Hamilton, et al., JCAP 2022, 034 (2022)
  * (8) K. Abazajian, et al., APJ 926, 54 (2022)
  * (9) T.J. Cornwell, P.N. Wilkinson, MNRAS 196, 1067 (1981)
  * (10) M.A. Bigot-Sazy, et al., A&A 59, 1 (2013)
  * (11) L. Mousset, et al., JCAP 2022, 035 (2022)
  * (12) S.A. Torchinsky, et al., JCAP 2022, 036 (2022)
  * (13) M. Piat, et al., JCAP 2022, 037 (2022)
  * (14) S. Masi, et al., JCAP 2022, 038 (2022)
  * (15) G. D’Alessandro, et al., JCAP 2022, 039 (2022)
  * (16) F. Cavaliere, et al., JCAP 2022, 040 (2022)
  * (17) C. O’Sullivan, et al., JCAP 2022, 041 (2022)
  * (18) R. Stompor, et al., MNRAS 392, 216 (2008)
  * (19) H.K. Eriksen, et al., APJ 641, 665 (2006)
  * (20) M. Regnier, et al. (2023), 2309.02957


# Level occupation switching with density functional theory

Nahual Sobrino<EMAIL_ADDRESS>Donostia International Physics Center
(DIPC), Paseo Manuel de Lardizabal 4, E-20018 San Sebastián, Spain
Departamento de Polímeros y Materiales Avanzados: Física, Química y
Tecnología, Universidad del País Vasco UPV/EHU, Av. Tolosa 72, E-20018 San
Sebastián, Spain  David Jacob<EMAIL_ADDRESS>Departamento de Polímeros y
Materiales Avanzados: Física, Química y Tecnología, Universidad del País Vasco
UPV/EHU, Av. Tolosa 72, E-20018 San Sebastián, Spain  IKERBASQUE, Basque
Foundation for Science, Plaza Euskadi 5, E-48009 Bilbao, Spain  Stefan Kurth
Departamento de Polímeros y Materiales Avanzados: Física, Química y
Tecnología, Universidad del País Vasco UPV/EHU, Av. Tolosa 72, E-20018 San
Sebastián, Spain  IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5,
E-48009 Bilbao, Spain  Donostia International Physics Center (DIPC), Paseo
Manuel de Lardizabal 4, E-20018 San Sebastián, Spain

###### Abstract

The charge transport properties of zero-temperature multi-orbital quantum dot
systems with one dot coupled to leads and the other dots coupled only
capacitatively are studied within density functional theory. It is shown that
the setup is equivalent to an effective single impurity Anderson model. This
allows to understand the level occupation switching effect as transitions
between ground states of different integer occupations in the uncoupled dots.
Level occupation switching is very sensitive to small energy differences and
therefore also to the details of the parametrized exchange-correlation
functionals. An existing functional already captures the effect on a
qualitative level but we also provide an improved parametrization which is
very accurate when compared to reference numerical renormalization group
results.

##  I Introduction

Quantum dots (QDs) are an ideal testbed to investigate the interplay between
quantum many-body physics and transport phenomena. They can be fabricated in
the lab from a large variety of materials and techniques, such as metallic
nanoparticles[1], lateral confinement of a two-dimensional electron gas (2DEG)
in GaAs/AlGaAs heterostructures (for a review see Ref. 2 and references
therein), carbon nanotubes (CNTs)[3, 4], and molecular junctions[5]. Indeed
important many-body phenomena such as the Kondo effect[6] and Coulomb blockade
(CB)[7], characteristic for so-called strongly correlated electrons, where
electronic interactions dominate over the kinetic energy, have been measured
in transport setups of QDs.[8, 9, 5, 10, 11]

The physics of QDs can be further enriched by the existence of multiple
electronic levels (or orbitals), or by coupling of QDs. The interplay between
strong electronic correlations and the spin and orbital degrees of freedom in
multi-orbital QDs, may lead to new physical phenomena, such as the SU(4) and
underscreened Kondo effects, which have both been measured in CNTs.[12, 13,
14] An interesting effect may occur in multi-orbital QDs when one of the QD
levels couples more strongly to the leads than the other levels. In this case
abrupt changes in the conductance and transmission phases between Coulomb
blockade peaks have been observed.[15, 16, 17, 18, 19, 20] These may be
attributed to the so-called level occupation switching (LOS), where the
strongly coupled level is abruptly emptied, while the weakly coupled level(s)
are abruptly filled, or vice versa. [21, 22, 23, 24, 25, 26] In essence this
phenomenon is a result of the competition between the kinetic energy of the
strongly coupled level and its electrostatic repulsion with the weakly coupled
levels.[23]

Density functional theory (DFT) is one of the most successful and popular
approaches for computing the electronic structure of molecules and solids
owing to its relative simplicity and computational efficiency.[27, 28, 29]
Since DFT is in principle an exact (many-body) theory for the ground-state
energy and density of a many-electron system, it should also be capable of
describing strong electronic correlation phenomena such as Kondo effect,
Coulomb blockade, and ultimately the LOS effect. However, in practice
approximations need to be made in DFT for the exchange-correlation (xc) part
of the total energy functional. And unfortunately, the most popular
approximations to DFT, such as the local-density[28] and generalized-gradient
approximations[30, 31, 32] in condensed-matter physics, and the so-called
hybrid functionals in chemistry,[33] are known to fail for strongly correlated
systems.

Nevertheless, if equipped with proper approximations for the xc part of the
functional, DFT is indeed capable of describing strongly correlated phenomena
such as Coulomb blockade and Kondo effect in transport through nanoscale
devices,[34, 35, 36] and the Mott-Hubbard gap in solids.[37] More recently, it
has also been shown that the actual many-body spectral function may be
extracted from a DFT calculation by making use of an extension of DFT called
i-DFT.[38, 39] This DFT framework also allows to describe the Mott metal-
insulator transition, one of the hallmarks of strong electronic
correlations.[40] The crucial ingredient for the description of these
phenomena within DFT are steps at integer occupations in the xc
potentials.[41] These steps are related to the famous derivative discontinuity
of exact DFT [42], which is missing in the standard approximations. In the
context of the Anderson impurity model, the step feature in the xc potential
gives rise to the pinning of the Kohn-Sham (KS) impurity level to the Fermi
energy, which results in a plateau for the zero-bias conductance as a function
of the applied gate, in accordance with the Kondo effect.[29, 30, 31]

Here we show how DFT can be used to study the LOS phenomenon which occurs in
multi-level quantum dot (MQD) systems in the asymmetric situation where only
one of the levels is connected to leads. In Section II we introduce the model
and show how it can exactly be mapped onto an effective single impurity
Anderson model (SIAM). In Section III we use lattice DFT for the effective
SIAM and an energy minimization argument to decide which configuration of
integer occupations of the disconnected dots is realized. In Section IV we
show that an existing parametrization of the SIAM Hxc potential already
qualitatively captures the LOS effect although not always at the correct value
of the gate. The origin of these deviations is investigated and remedied by a
re-parametrization of the Hxc functional. Finally, we present our conclusions
in Section V.

##  II Model

We consider a multi-orbital quantum dot consisting of $M$ impurities which can
hold up to two electrons and which are all capacitively coupled among each
other. We consider the situation when only one of the impurities is connected
to two leads and we also restrict ourselves to the zero-temperature limit. The
total Hamiltonian of the system is given as the sum of the Hamiltonians of the
isolated dot and leads, as well as the coupling between them and reads

$\displaystyle\hat{H}=\hat{H}_{{\rm MQD}}+\sum_{\alpha
k\sigma}\varepsilon_{\alpha k}\hat{c}^{\dagger}_{\alpha
k\sigma}\hat{c}_{\alpha k\sigma}$ $\displaystyle+\sum_{\alpha
k\sigma}\left(t_{\alpha k}\hat{d}^{\dagger}_{1\sigma}\hat{c}_{\alpha
k\sigma}+{\rm H.c.}\right)$ (1)

where

$\displaystyle\hat{H}_{{\rm
MQD}}=\sum_{i=1}^{M}v_{i}\hat{n}_{i}+\sum_{i=1}^{M}U_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}+\sum_{i<j}^{M}U_{ij}\hat{n}_{i}\hat{n}_{j}$
(2)

describes the capacitatively coupled MQD where
$\hat{n}_{i}=\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}$ is the number operator
for level $i$ with
$\hat{n}_{i\sigma}=\hat{d}_{i\sigma}^{\dagger}\hat{d}_{i\sigma}$ and
$\hat{d}_{i\sigma}$ ($\hat{d}_{i\sigma}^{\dagger}$) are the annihilation
(creation) operators for electrons with spin $\sigma$ in orbital $i$. In Eq.
2, $v_{i}$ and $U_{i}$ are the on-site energy (also referred to as gate) and
the intra-Coulomb repulsion of level $i$ while $U_{ij}$ is the inter-Coulomb
repulsion between levels $i$ and $j$. The second term in Section II describes
the non-interacting electrons in left (L) and right (R) leads ($\alpha=L,R$)
while the third term accounts for the (symmetric) coupling of the first
impurity to left and right leads. The resulting broadening functions
$\Gamma_{\alpha}({\omega})=2\pi\sum_{k}|t_{\alpha
k}|^{2}\delta({\omega}-\varepsilon_{\alpha k})$ are assumed to describe
featureless leads and therefore become independent of frequency, i.e., we take
the wide band limit (WBL) with $\Gamma_{\alpha}({\omega})=\gamma_{\alpha}$. In
Fig. 1 a schematic representation of the MQD setup is shown.

Figure 1: Schematic representation of the transport setup for the MQD. The
first impurity is connected to the electron reservoirs while the rest only
interact through electrostatic repulsions with each other.

Since the impurities $j=2,\ldots,M$ are not connected to the reservoirs, the
multi-orbital quantum dot system can be mapped exactly onto an effective
SIAM.[23] This follows from the fact that the operators $\hat{n}_{j\sigma}$
for $j>1$ all commute with the Hamiltonian of Section II. Therefore, all the
many-body eigenstates of $\hat{H}$ can be chosen to be eigenstates of
$\hat{n}_{j\sigma}$ (for $j>1$) and the corresponding eigenvalues take values
$n_{j\sigma}=0,1$. Therefore, the total occupations $n_{j}$ of the
disconnected dots can only take the integer values $n_{j}=0,1,2$. As a
consequence of this commutation property, when the Hamiltonian is applied to
an eigenstate of $\hat{n}_{j\sigma}$ ($j>1$), the problem is seen to be
equivalent to an effective SIAM related to the first impurity with an
effective potential

$\displaystyle v_{1}^{\rm eff}=v_{1}+\sum_{j=2}^{M}U_{1j}n_{j}$ (3)

and an additional constant contribution to the total energy given by
$\sum_{j\neq 1}U_{j}\delta_{n_{j},2}$. Therefore, the only effect of charging
and discharging the levels $j=2,\ldots,M$ is to modify the average Hartree
potential felt by the first level.

##  III Lattice density functional theory

With the observations of the previous Section it is clear that any many-body
method which can accurately treat the SIAM can be employed to obtain the
ground state energy and density of $\hat{H}$. Here our method of choice is
lattice DFT which in the past has successfully been used for the SIAM. [34,
35, 36, 43].

###  III.1 Energy functional

In lattice DFT, for a given set of gates $\mathbf{v}=(v_{1},\ldots,v_{M})$ the
total energy functional for the interacting system described by $\hat{H}$
reads

$\displaystyle E[\mathbf{n}]=T_{s}[n_{1}]+E_{\rm
Hxc}[\mathbf{n}]+\sum_{i=1}^{M}n_{i}v_{i},$ (4)

where $T_{s}[n_{1}]$ is the non-interacting kinetic energy of the only
impurity connected to leads (the kinetic energy of the disconnected dots
vanishes). $E_{\rm Hxc}[\mathbf{n}]$ is the Hartree-exchange-correlation (Hxc)
energy which is a functional of all densities and we used the notation
$\mathbf{n}=(n_{1},\ldots,n_{M})$. The non-interacting kinetic energy can be
expressed in terms of the (local) Green function $G^{R}_{d}(\omega)$ of the
connected dot as
$T_{s}=\frac{\gamma}{\pi}\int_{-\varepsilon_{C}}^{0}d\omega\text{Re}\left[G_{d}^{R}(\omega)\right]$.
Here $\varepsilon_{C}$ is an energy cutoff which ensures the convergence of
the integral. In the WBL the non-interacting kinetic energy contribution can
be expressed in a closed form (see Appendix for the detailed derivation) as

$\displaystyle
T_{s}[n]=\frac{\gamma}{2\pi}\log\left[\frac{\frac{\gamma^{2}}{4}\left(\tan{\left(\frac{\pi}{2}(1-n)\right)}^{2}+1\right)}{\varepsilon_{C}^{2}+\frac{\gamma^{2}}{4}}\right].$
(5)

Note that the energy cutoff only acts as a constant shift to the total energy
but does not change its shape. Therefore both the energy minimum and any
energy difference between different configurations are independent of the
value of $\varepsilon_{C}$.

The Hxc energy functional $E_{\rm Hxc}[\mathbf{n}]$ of the total system can be
simplified using the considerations of the previous Section. Since for the
disconnected dots $i\in\\{2,\ldots,M\\}$, the only possible occupations are
$n_{i}=0,1,2$ (and the KS system has to reproduce these occupations), $E_{\rm
Hxc}$ may be written as

$\displaystyle E_{\rm Hxc}[\mathbf{n}]=E_{\rm Hxc}^{\rm
SIAM}[n_{1}]+\sum_{i<j}^{M}U_{ij}n_{i}n_{j}+\sum_{j=2}^{M}U_{j}\delta_{n_{j},2}$
(6)

where $E_{\rm Hxc}^{\rm SIAM}[n]$ is the Hxc functional for a simple SIAM.
Inserting Eq. 6 into Eq. 4 it becomes immediately clear that the resulting
total energy functional has the form of the energy of a single Anderson
impurity but with effective on-site potential $v_{1}^{\rm eff}$ given by Eq.
3.

For a given set of gate levels ($v_{1},\ldots,v_{M}$), there are $3^{M-1}$
different configurations of occupations ($n_{2},\ldots,n_{M}$) of the
disconnected dots. The ground state energy and the resulting set of
occupancies of the system can then be obtained by comparing the energies of
the different configurations of available states and taking the one
corresponding to the minimum of energy. In the degenerate case $v_{i}=v$, the
ground state energy of the system can be found by minimizing the universal
functional $F[n]=T_{s}[n]+E_{Hxc}[n]$, since the last term of Eq. 4 is
constant (at given fixed total occupation) due to the one-to-one
correspondence between total occupation $N$ and the external potential $v$.

A LOS event exactly corresponds to a change in the ground state energy of the
system, and provided an accurate parametrization for the SIAM Hxc energy
$E_{\rm Hxc}^{\rm SIAM}$ is used, this event can be completely captured within
DFT.

###  III.2 Kohn-Sham equation

For a given configuration of (integer) occupations ($n_{2},\ldots,n_{M}$) of
the disconnected dots, the (non-integer) density on the connected dot can now
be found by the Hohenberg-Kohn variational principle, i.e., by searching for
the value $n_{1}$ which minimizes the total energy (4). Therefore we need to
solve $\frac{\partial E[\mathbf{n}]}{\partial n_{1}}=0$ which is easily shown
to be equivalent to solving the KS equation

$\displaystyle n_{1}=1-\frac{2}{\pi}\arctan{\left(2\frac{v_{1}^{\rm
eff}+v_{\rm Hxc}^{\rm SIAM}[n_{1}]}{\gamma_{1}}\right)},$ (7)

where we have used the definition (3) of the effective potential $v_{1}^{\rm
eff}$ and defined the SIAM Hxc potential as

$\displaystyle v_{\rm Hxc}^{\rm SIAM}[n]=\frac{{\rm d}E_{\rm Hxc}^{\rm
SIAM}[n]}{{\rm d}n}\;.$ (8)

We note in passing that Eq. 7 is equivalent to expressing the density as[44]

$\displaystyle n_{1}=2\int_{-\infty}^{0}\frac{{\rm
d}{\omega}}{2\pi}\;A_{s,1}({\omega})$ (9)

where, without loss of generality, we assumed a vanishing equilibrium chemical
potential in the leads. The KS spectral function for the connected dot is

$\displaystyle
A_{s,1}(\omega)=\frac{{\gamma}_{1}}{\frac{{\gamma}_{1}^{2}}{4}+(\omega-
v_{s,1})^{2}}\;,$ (10)

consistent with WBL approximation used to derive Eq. 5, and
$v_{s,1}=v_{1}^{\rm eff}+v_{\rm Hxc}^{\rm SIAM}$.

According to Friedel sum rule, in the zero-temperature limit, the impurity
spectral function at the Fermi energy is completely determined by the impurity
density.111See, e.g., Ch. 5.2 in the book by Hewson[6] Since in the setup
considered here only one impurity level is connected to the leads, the zero-
bias conductance is directly given by the spectral function at the Fermi
level. In this case the Friedel sum rule implies that the electrical
conductance of the system is fully determined by the equilibrium density at
the impurity and reads[46, 47, 48]

$\displaystyle G=\sin{\left(\frac{\pi n_{1}}{2}\right)}^{2}.$ (11)

On the other hand, the correct description of the electrical conductance
requires the access to the actual (many-body) electrical current of the
system, and therefore one possibility is the inclusion of the xc corrections
to the bias of the system. [49, 50, 51] However, since here the exact zero-
bias conductance can be expressed explicitly in terms of the ground state
density alone (which is a quantity accessible to standard ground state DFT),
already the KS zero-bias conductance $G_{s}$ is exact [34, 35, 36, 52, 53] and
can be expressed as

$\displaystyle
G=G_{s}=\frac{\frac{\gamma^{2}}{4}}{v_{s,1}^{2}+\frac{\gamma^{2}}{4}}.$ (12)

Eqs.(11) and (12) provide two equivalent expressions in the zero temperature
regime for the zero-bias conductance. [46, 47, 48]

Figure 2: Comparison of the local occupancies obtained with the exchange-
correlation functional of Ref. 35 (DFT 1), and those obtained with the new
parametrization of Eq. 13 (DFT 2) to the NRG results of Ref. 22. Both $n_{i}$
are shown as function of the gate level for the degenerate case
$v=v_{1}=v_{2}$ and for strong correlations $\gamma/U=0.1$.

##  IV Results

As discussed in Section II, the problem of the capacitively coupled MQD with
only one impurity coupled to leads can exactly be mapped onto a SIAM. In the
context of lattice DFT (Sec. III) this implies that the only quantity to be
approximated is the Hxc functional $E_{\rm Hxc}^{\rm SIAM}$ from which $v_{\rm
Hxc}^{\rm SIAM}$ follows by differentiation. Once such a parametrization is
given, the MQD DFT problem is solved as follows: For a given set of gate
potentials $\mathbf{v}$, for each configuration of (integer) occupations of
the disconnected impurities (covering the $3^{M-1}$ possible configurations)
we solve the KS equation (7) for $n_{1}$ with the effective potential of Eq. 3
and calculate the corresponding total energies. The lowest of these energies
then corresponds to the configuration of the ground state.

For the SIAM functional, we start by considering the new parametrization at
zero temperature suggested in Ref. 35. In Fig. 2 we present the local
occupancies for the case of the double quantum dot $M=2$ as a function of the
gate level in the degenerate case $v=v_{1}=v_{2}$ and strong correlations
$\gamma/U=0.1$. The results labeled DFT 1 correspond to the self-consistent
densities obtained with the parametrization of Ref. 35. The first thing to
note is that our approach does capture the LOS transitions and gives densities
which qualitatively agree with the reference NRG results of Ref. 22. However,
we also notice that the LOS events take place at values of $v$ considerably
different from the many-body results. Since, as mentioned above, the only
possible source of error in our approach is the parametrization of $E_{\rm
Hxc}^{\rm SIAM}$, below we propose a reparametrization of the functional (DFT
2 in Fig. 2) which correctly and accurately captures the LOS transitions.

Figure 3: Self-consistent densities of the simple SIAM obtained with the
parametrization of $v_{\rm Hxc}^{\rm SIAM}$ given in Ref. 35 (DFT 1) and our
new parametrization (DFT 2) as function of gate for different values of
$\gamma/U$. The upper insets highlights small differences in the mixed-valence
regime while the lower inset shows the Hxc potentials for the two
parametrizations. Figure 4: Energy difference between the states with
$n_{2}=0$ and $n_{2}=1$ of the double quantum dot for different coupling
strengths in the case $v=v_{1}=v_{2}$. The LOS events exactly correspond to
the degeneracy of the states $E_{n_{2}=1}=E_{n_{2}=0}$. Note how the new
parametrization shifts the crossings of zero of the energy differences, i.e.,
the gate values at which the LOS event occurs.

The parametrization of Ref. 35 depends on two parameters which are both
functions of $\gamma/U$. In order to correctly capture the LOS events we here
propose to keep the same functional form but reparametrize the parameter
$\sigma$ of Eq. (16) of Ref. 35. Our fit to those numerical values of $\sigma$
which best reproduces the positions of the LOS events for the parameters of
Ref. 22 is given as

$\displaystyle\sigma=0.07\arctan\left(\frac{171.358(\gamma/U)^{2}}{2+\gamma/U}\right)\;.$
(13)

This newly parametrized Hxc functional for the SIAM will be denoted as DFT 2
in the following. It now accurately captures the gate positions of the LOS
events, see Fig. 2.

The gate positions of the LOS events are highly sensitive to the
parametrization of the Hxc functional, especially in the mixed valence regime.
This can be appreciated in Fig. 3, where we show the self-consistent SIAM
densities produced with the DFT 1 (solid line) and the DFT 2 (dashed line)
functionals for different coupling strengths. The discrepancies between the
corresponding densities are almost negligible, with the maximum disagreement
in the mixed valence regime, i.e., in the transition from empty to half
occupation and from half occupation to full occupation. In the inset of Fig.
3, both parametrizations of $v_{\rm Hxc}^{\rm SIAM}[n]$ are compared. Note
that although the corrections are very small (of the order of $\sim
10^{-3}U$), they are crucial in order to accurately capture the LOS events.

Figure 5: Comparison of the zero-bias conductance obtained with the exchange-
correlation functional of Ref. 35 (DFT 1) with the new parametrization (DFT 2)
and NRG calculations. $G$ is shown as function of the gate level in the
degenerate case $v=v_{1}=v_{2}$ and for different coupling strengths
$\gamma/U=0.1,0.2,0.25,0.3$ from a) to d), respectively. All energies in units
of U.

In Fig. 4, we show as an illustrative example the difference between the
computed energies for two different (integer) $n_{2}$, $E_{n_{2}=1}$ and
$E_{n_{2}=0}$, obtained self-consistently with the different parametrizations.
The effect of the new parametrization is considerably larger for some values
of the coupling strength ($\gamma/U=0.1,0.2,0.3$), while for others it is not
appreciable ($\gamma/U=0.25$). In particular, this difference is relevant at
$E_{n_{2}=1}-E_{n_{2}=0}=0$, which exactly corresponds to the gate at which
the the LOS event occurs. Although not shown in the present paper, the new
parametrization does not introduce any change in the $E_{n_{2}=2}=E_{n_{2}=0}$
crossing (only found at $v/U=-1,-2$) while the transition from $E_{n_{2}=2}$
to $E_{n_{2}=1}$ follows exactly the same correction as the one illustrated
for gates centered at $v/U=-2$.

###  IV.1 Results for the double quantum dot

Figure 6: Zero-bias conductance as function of the gate voltage
$v_{g}=v-U_{12}-U/2$ in units of $U$ for a) $U_{12}=0.99U$ and b)
$U_{12}=0.5U$. In both panels $\gamma/U=0.1$.

In Fig. 5 we present the conductance as function of the gate level in the
degenerate case $v=v_{1}=v_{2}$ for different coupling strengths
$\gamma/U=0.1,0.2,0.25,0.3$ from a) to d), respectively. For small coupling
strength, six LOS events occur, two of them pinned at the gate values
$v/U=-2,-1$, related to two peaks usually referred to as CB peaks,[54] since
similar structures are present in the CB regime. The other four LOS
transitions lead to a widening of the three Kondo plateaus. We observe that
the DFT 1 results are qualitatively correct, while still failing to predict
the evolution of the LOS events with increasing coupling strength, except for
the ones related to the two CB peaks. On the other hand, the DFT 2 results
accurately reproduce the NRG calculations, showing the correct evolution of
the three Kondo plateaus to wider structures and the two CB peaks into dips
for higher values of $\gamma/U=0.3$.

Figure 7: Conductance as function of the gate level in units of $U$ for a)
$\Delta v=0.005$ and b) $\Delta v=0.0075$. In both panels $\gamma/U=0.1$.

The results shown so far are related to the fully degenerate case, i.e,
$v=v_{1}=v_{2}$ and $U=U_{1}=U_{2}=U_{12}$, but our DFT study can go further.
Since the dependence on the inter-Coulomb repulsions only enters into the
effective potential felt by the first impurity, we can easily investigate the
effects of changing it. In Fig. 6 we show the conductances for a)
$U_{12}=0.99U$ and b) $U_{12}=0.5U$ as function of the gate voltage
$v_{g}=v-U_{12}-U/2$ in units of $U$. Again, we find that DFT 2 gives an
excellent agreement with NRG except for the central plateau in panel a), where
the NRG calculations predict a slightly wider central structure. When the
inter-Coulomb repulsion is decreased to be half of the intra-Coulomb
repulsion, the two CB peaks and the central Kondo plateau merge into a smooth
plateau leading to the disappearance of four of the LOS events. On the other
hand, in Fig. 7 we explore the effect of considering a finite difference
between the gate levels $\Delta v=v_{1}-v_{2}$. In Fig. 7 a) the DFT 2
parametrization correctly captures the evolution of the two CB peaks for
$\Delta v=0.005$, while for $\Delta v=0.0075$ (Fig. 7 b)) the DFT 2 results
present a finite difference with the reference NRG results. The origin of this
discrepancy is not completely clear: it could either be due to subtle details
of the Hxc potentials which our parametrization does not capture. However, it
could also be due to small numerical effects due to the sensitivity of the LOS
events to small energy differences. However, we do find that our results for
$\Delta v=0.01$ completely agree with the NRG ones for $\Delta v=0.0075$.

###  IV.2 Results for more than two dots

We further apply our new parametrization of $v_{\rm Hxc}^{\rm SIAM}$ to the
situation of more than two orbitals with only one of them being connected to
leads. The generalization is straightforward and only requires to find the
ground state energy between the different $3^{M-1}$ states corresponding to
the different configurations of integer occupancies of the disconnected dots.

Some results for the fully degenerate case of the triple quantum dot
($v=v_{1}=v_{2}=v_{3}$) for $\gamma/U=0.2$ are presented in Fig. 8. In panel
a) the densities predict a total of twelve LOS events. Since both the second
and third impurity levels are completely degenerate, a swap in their local
occupancies leaves the problem invariant and the corresponding many-body
eigenstates are degenerate. Therefore, the (average) occupancies may take
semi-integer values for some values of the gate. Following the same reasoning
we observe that for the fully degenerate MQD, the $M-1$ degenerate levels of
disconnected dots can only reach occupancies of integer multiples of
$(M-1)^{-1}$. In panel b) the conductance present four CB peaks and five Kondo
plateaus, which evolve with the coupling strength in an analogous manner but
reaching the inversion of the CB peaks into valleys around $\gamma/U\sim 0.5$.

Finally, in Fig. 9 we present the densities and the related conductances for
the triple (panels a) and c)) and quadruple (panels b) and d)) quantum dots in
the nondegenerate case. For the triple quantum dot we choose the values
$U=\frac{5}{4}U_{12}=\frac{5}{3}U_{13}=2U_{23}$ and
$v=v_{1}=v_{2}-0.1=v_{3}-0.4$ and for the quadruple quantum dot
$U=U_{i}=2U_{ij}$ for all $i\neq j$ and
$v=v_{1}=v_{2}-0.02=v_{3}-0.02=v_{4}-0.04$. In both cases we choose
$U=10\gamma$. For the selected set of parameters, the LOS events always
correspond to an abrupt filling of a weakly coupled impurity and an abrupt
emptying of the strongly coupled one. In both systems we observe that as the
gate is decreased, a LOS event for $n_{1}<1$ ($n_{1}>1$) implies a sudden
decreasing (increasing) of the conductance.

Figure 8: Densities and electrical conductance of the triple quantum dot as
function of the gate level for the fully degenerate case and $\gamma/U=0.2$.
Figure 9: DFT densities and zero-bias conductance as function of the gate
level. Panels a) and c): Triple quantum dot with
$U=\frac{5}{4}U_{12}=\frac{5}{3}U_{13}=2U_{23}$ and
$v=v_{1}=v_{2}-0.1=v_{3}-0.4$. Panels b) and d): Quadruple quantum dot with
$U=U_{i}=2U_{ij}$ for all $i\neq j$ and
$v=v_{1}=v_{2}-0.02=v_{3}-0.02=v_{4}-0.04$. In both cases we chose
$U=10\gamma$.

##  V Conclusions

In this work we have studied from a DFT perspective the LOS effect which
occurs in multi-orbital quantum dots subject to inter and intra Coulomb
repulsions when only one of the dots is coupled to leads. The system can be
mapped into an effective SIAM problem for the coupled impurity which
experiences an effective gate potential due to electrostatic interactions with
the other impurities. The density of the system is obtained by choosing the
minimum total energy (expressed in terms of DFT quantities) among the
available $3^{M-1}$ different configurations of integer occupation of the
uncontacted levels. A LOS event occurs at that gate for which the energy
minimum passes from one configuration of integer occupations to another one.
We have modified an already quite accurate parametrization of the SIAM Hxc
functional in order to correctly capture the coupling strength dependence of
the LOS events which is quite sensitive to details of the functional. Our new
parametrization yields very small energy differences (of the order of $\sim
10^{-3}U$) as compared to the previous one. This produces almost no effect on
the selfconsistent densities of SIAM, but is essential to shift the gates at
which the LOS events occur to the correct positions. DFT calculations
employing the parametrized Hxc functionals for the double quantum dot show
excellent agreement with many-body NRG calculations. We have further presented
results for the triple and quadruple quantum dot.

##  VI Acknowledgments

We acknowledge funding through the grant “Grupos Consolidados UPV/EHU del
Gobierno Vasco” (IT1453-22). We also acknowledges funding through a grant of
the ”Ministerio de Ciencia y Innovación (MCIN)” (Grant No.
PID2020-112811GB-I00).

##  Appendix A Kinetic energy of a non-interacting impurity coupled to two
leads

The purpose of this Appendix is to derive the explicit functional of the non-
interacting kinetic energy of a non-interacting impurity connected to left (L)
and right (R) leads. We start by writing the Hamiltonian of this non-
interacting impurity in second quantized form as

$\displaystyle\hat{H}=v_{s}\hat{d}^{\dagger}\hat{d}+\sum_{\alpha=L,R}\sum_{k\sigma}\epsilon_{\alpha
k}\hat{c}^{\dagger}_{\alpha k\sigma}\hat{c}_{\alpha k\sigma}+\hat{H}_{\rm
kin}$ (14)

where

$\displaystyle\hat{H}_{\rm kin}=\sum_{\alpha=L,R}\sum_{k\sigma}\left(V_{\alpha
k,d}\hat{d}^{\dagger}_{\sigma}\hat{c}_{\alpha k\sigma}+{\rm H.c.}\right)\;.$
(15)

The (non-interacting) kinetic energy is then given as

$\displaystyle T_{s}=\langle\hat{H}_{\rm kin}\rangle=2\sum_{\alpha
k}\left(-iV_{\alpha k,d}G_{\alpha k,d}^{<}(t^{+},t)+{\rm c.c.}\right)$ (16)

where the prefactor two comes from spin

$\displaystyle G_{\alpha
k,d}^{<}(t,t^{\prime})=i\langle\hat{d}^{\dagger}_{\sigma,H}(t^{\prime})\hat{c}_{\alpha
k\sigma,H}(t)\rangle$ (17)

is the (spin-independent) matrix element of the lesser Green function between
single-particle basis states $|\alpha k\sigma\rangle$ and $|d\sigma\rangle$.
By standard Green function techniques [55] Eq. 16 can be written as

$\displaystyle T_{s}=2\sum_{\alpha k}\left(V_{\alpha
k,d}\int_{-\varepsilon_{C}}^{0}\frac{{\rm d}\omega}{2\pi}A_{\alpha
k,d}(\omega)+{\rm c.c.}\right)$ $\displaystyle=2\sum_{\alpha k}\left(V_{\alpha
k,d}\int_{-\varepsilon_{C}}^{0}\frac{{\rm d}\omega}{2\pi}i(G_{\alpha
k,d}^{R}(\omega)-G_{\alpha k,d}^{A}(\omega))+{\rm c.c.}\right)$

where $G^{R}$ ($G^{A}$) are the retarded (advanced) Green functions and we
introduced an energy cutoff $\varepsilon_{C}$ in order for the integral to
converge. The retarded Green function $\hat{G}^{R}(\omega)$ at energy $\omega$
is defined through [55, 44]

$\displaystyle\left((\omega+i\eta)\hat{\mathbb{I}}-\hat{H}\right)\hat{G}^{R}(\omega)=\hat{\mathbb{I}}.$
(19)

In the single-particle basis, all Hamiltonian matrix elements directly
connecting left and right leads vanish, i.e., $\langle
Lk\sigma|\hat{H}|Rk^{\prime}\sigma^{\prime}\rangle=\langle
Rk\sigma|\hat{H}|Lk^{\prime}\sigma^{\prime}\rangle=0$. Then, from Eq. 19 one
can derive

$\displaystyle G_{\alpha k,d}^{R}(\omega)=\frac{V_{\alpha
k,d}^{*}}{\omega-\epsilon_{k}+i\eta}G_{d}(\omega)$ (20) $\displaystyle
G_{d}^{R}(\omega)=\frac{1}{\omega-v_{s}-\Delta^{R}(\omega)}$ (21)

with $\eta\to 0^{+}$ and
$\Delta^{R}(\omega)=\sum_{\alpha}\Delta_{\alpha}^{R}(\omega)$ is the total
embedding self energy with

$\displaystyle\Delta_{\alpha}^{R}(\omega)=\sum_{k}\frac{|V_{\alpha
k,d}|^{2}}{\omega-\varepsilon_{\alpha k}+i\eta}$ (22)

In the wide-band limit we have
$\Delta_{\alpha}^{R}(\omega)=-i\frac{\gamma_{\alpha}}{2}$, independent of
$\omega$. Inserting Eqs. (20) and (22) into LABEL:Ts_kin2 we arrive at

$\displaystyle T_{s}=\frac{\gamma}{\pi}\int_{-\varepsilon_{C}}^{0}{\rm
d}\omega\;{\rm
Re}[G_{d}^{R}(\omega)]=\frac{\gamma}{2\pi}\log\left[\frac{v_{s}^{2}+\frac{\gamma^{2}}{4}}{\varepsilon_{C}^{2}+\frac{\gamma^{2}}{4}}\right]$
(23)

The non-interacting density-potential relation is

$\displaystyle n=1-\frac{2}{\pi}\arctan\left(\frac{2v_{s}}{\gamma}\right)$
(24)

which can easily be inverted to give

$\displaystyle v_{s}=\frac{\gamma}{2}\tan\left(\frac{\pi}{2}(1-n)\right)\;.$
(25)

Inserting Eq. 25 into Eq. 23 then gives the final result for the non-
interacting kinetic energy functional

$\displaystyle
T_{s}[n]=\frac{\gamma}{2\pi}\log\left[\frac{\frac{\gamma^{2}}{4}\left(\tan{\left(\frac{\pi}{2}(1-n)\right)}^{2}+1\right)}{\varepsilon_{C}^{2}+\frac{\gamma^{2}}{4}}\right].$
(26)

## References

  * Petta and Ralph [2001] J. R. Petta and D. C. Ralph, Phys. Rev. Lett. 87, 266801 (2001).
  * Hanson  _et al._ [2007] R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007).
  * Nygård  _et al._ [2000] J. Nygård, D. H. Cobden, and P. E. Lindelof, Nature 408, 342 (2000).
  * Buitelaar  _et al._ [2002] M. R. Buitelaar, A. Bachtold, T. Nussbaumer, M. Iqbal, and C. Schönenberger, Phys. Rev. Lett. 88, 156801 (2002).
  * Park  _et al._ [2002] J. Park, A. N. Pasupathy, J. I. Goldsmith, C. Chang, Y. Yaish, J. R. Petta, M. Rinkoski, J. P. Sethna, H. D. Abruña, P. L. McEuen, and D. C. Ralph, Nature 417, 722 (2002).
  * Hewson [1997] A. C. Hewson,  _The Kondo problem to heavy fermions_ (Cambr. Univ. Press, Cambridge, 1997).
  * Averin and Likharev [1986] D. V. Averin and K. K. Likharev, Journal of Low Temperature Physics 62, 345 (1986).
  * Shtrikman  _et al._ [1998] H. Shtrikman, D. Mahalu, D. Abusch-Magder, U. Meirav, M. A. Kastner, and D. Goldhaber-Gordon, Nature 391, 156 (1998).
  * Goldhaber-Gordon  _et al._ [1998] D. Goldhaber-Gordon, J. Göres, M. A. Kastner, H. Shtrikman, D. Mahalu, and U. Meirav, Phys. Rev. Lett. 81, 5225 (1998).
  * Champagne  _et al._ [2005] A. R. Champagne, A. N. Pasupathy, and D. C. Ralph, Nano Lett. 5, 305 (2005).
  * Parks  _et al._ [2007] J. J. Parks, A. R. Champagne, G. R. Hutchison, S. Flores-Torres, H. D. Abruña, and D. C. Ralph, Phys. Rev. Lett. 99, 026601 (2007).
  * Jarillo-Herrero  _et al._ [2005] P. Jarillo-Herrero, J. Kong, H. S. J. van der Zant, C. Dekker, L. P. Kouwenhoven, and S. DeFranceschi, Nature 434, 484 (2005).
  * Roch  _et al._ [2009] N. Roch, S. Florens, T. A. Costi, W. Wernsdorfer, and F. Balestro, Phys. Rev. Lett. 103, 197202 (2009).
  * Parks  _et al._ [2010] J. J. Parks, A. R. Champagne, T. A. Costi, W. W. Shum, A. N. Pasupathy, E. Neuscamman, S. Flores-Torres, P. S. Cornaglia, A. A. Aligia, C. A. Balseiro, G. K.-L. Chan, H. D. Abruña, and D. C. Ralph, Science 328, 1370 (2010).
  * Avinun-Kalish  _et al._ [2005] M. Avinun-Kalish, M. Heiblum, O. Zarchin, D. Mahalu, and V. Umansky, Nature 436, 529 (2005).
  * Yacoby  _et al._ [1995] A. Yacoby, M. Heiblum, D. Mahalu, and H. Shtrikman, Phys. Rev. Lett. 74, 4047 (1995).
  * Schuster  _et al._ [1997] R. Schuster, E. Buks, M. Heiblum, D. Mahalu, V. Umansky, and H. Shtrikman, Nature 385, 417 (1997).
  * Aikawa  _et al._ [2004] H. Aikawa, K. Kobayashi, A. Sano, S. Katsumoto, and Y. Iye, J. Phys. Soc. Jpn. 73, 3235 (2004).
  * Ji  _et al._ [2002] Y. Ji, M. Heiblum, and H. Shtrikman, Phys. Rev. Lett. 88, 076601 (2002).
  * Zaffalon  _et al._ [2008] M. Zaffalon, A. Bid, M. Heiblum, D. Mahalu, and V. Umansky, Phys. Rev. Lett. 100, 226601 (2008).
  * Lim  _et al._ [2006] J. S. Lim, M.-S. Choi, M. Choi, R. López, and R. Aguado, Phys. Rev. B 74, 205119 (2006).
  * Kleeorin and Meir [2017a] Y. Kleeorin and Y. Meir, Phys. Rev. B 96, 045118 (2017a).
  * Büsser  _et al._ [2011] C. Büsser, E. Vernek, P. Orellana, G. Lara, E. Kim, A. Feiguin, E. Anda, and G. Martins, Phys. Rev. B 83, 125404 (2011).
  * Roura-Bas  _et al._ [2011] P. Roura-Bas, L. Tosi, A. Aligia, and K. Hallberg, Phys. Rev. B 84, 073406 (2011).
  * Weymann  _et al._ [2018] I. Weymann, R. Chirla, P. Trocha, and C. P. Moca, Phys. Rev. B 97, 085404 (2018).
  * Silvestrov and Imry [2007] P. Silvestrov and Y. Imry, New J. Phys. 9, 125 (2007).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Dreizler and Gross [1990] R. M. Dreizler and E. K. U. Gross,  _Density Functional Theory_ (Springer, Berlin, 1990).
  * J.P. Perdew [1985] J.P. Perdew, Phys. Rev. Lett. 55, 1665 (1985).
  * Becke [1988] A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * A.D. Becke [1993] A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
  * Stefanucci and Kurth [2011] G. Stefanucci and S. Kurth, Phys. Rev. Lett. 107, 216401 (2011).
  * Bergfield  _et al._ [2012] J. P. Bergfield, Z.-F. Liu, K. Burke, and C. A. Stafford, Phys. Rev. Lett. 108, 066801 (2012).
  * Tröster  _et al._ [2012] P. Tröster, P. Schmitteckert, and F. Evers, Phys. Rev. B 85, 115409 (2012).
  * Lima  _et al._ [2002] N. A. Lima, L. N. Oliveira, and K. Capelle, Europhys. Lett. 60, 601 (2002).
  * Jacob and Kurth [2018] D. Jacob and S. Kurth, Nano Lett. 18, 2086 (2018).
  * Kurth  _et al._ [2019] S. Kurth, D. Jacob, N. Sobrino, and G. Stefanucci, Phys. Rev. B 100, 085114 (2019).
  * Jacob  _et al._ [2020] D. Jacob, G. Stefanucci, and S. Kurth, Phys. Rev. Lett. 125, 216401 (2020).
  * Sobrino  _et al._ [2020] N. Sobrino, S. Kurth, and D. Jacob, Phys. Rev. B 102, 035159 (2020).
  * Perdew  _et al._ [1982] J. P. Perdew, R. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
  * Kurth and Stefanucci [2016a] S. Kurth and G. Stefanucci, Phys. Rev. B 94, 241103(R) (2016a).
  * Kurth and Stefanucci [2017] S. Kurth and G. Stefanucci, J. Phys.: Condens. Matter 29, 413002 (2017).
  * Note [1] See, e.g., Ch. 5.2 in the book by Hewson[6].
  * Mera  _et al._ [2010] H. Mera, K. Kaasbjerg, Y. Niquet, and G. Stefanucci, Phys. Rev. B 81, 035110 (2010).
  * Langreth [1966] D. C. Langreth, Phys. Rev. 150, 516 (1966).
  * Mera and Niquet [2010] H. Mera and Y. Niquet, Phys. Rev. Lett. 105, 216408 (2010).
  * Stefanucci and Kurth [2015] G. Stefanucci and S. Kurth, Nano Lett. 15, 8020 (2015).
  * Sobrino  _et al._ [2021] N. Sobrino, F. Eich, G. Stefanucci, R. D’Agosta, and S. Kurth, Phys. Rev. B 104, 125115 (2021).
  * Sobrino  _et al._ [2019] N. Sobrino, R. D’Agosta, and S. Kurth, Phys. Rev. B 100, 195142 (2019).
  * Kurth and Stefanucci [2016b] S. Kurth and G. Stefanucci, Phys. Rev. B 94, 241103(R) (2016b).
  * Stefanucci and Kurth [2013] G. Stefanucci and S. Kurth, phys. stat. sol. (b) 250, 2378 (2013).
  * Kleeorin and Meir [2017b] Y. Kleeorin and Y. Meir, Phys. Rev. B 96, 045118 (2017b).
  * Stefanucci and van Leeuwen [2013] G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ (Cambridge University Press, 2013).


# Mechanism of magnetic phase transition in correlated magnetic metal: insight
into itinerant ferromagnet Fe3-δGeTe2

Yuanji Xu  Institute for Applied Physics, University of Science and Technology
Beijing, Beijing 100083, China  Yue-Chao Wang  Laboratory of Computational
Physics, Institute of Applied Physics and Computational Mathematics, Beijing
100088, China  Xintao Jin  Institute for Applied Physics, University of
Science and Technology Beijing, Beijing 100083, China  Haifeng Liu  Laboratory
of Computational Physics, Institute of Applied Physics and Computational
Mathematics, Beijing 100088, China  Yu Liu  Laboratory of Computational
Physics, Institute of Applied Physics and Computational Mathematics, Beijing
100088, China  Haifeng Song<EMAIL_ADDRESS>Laboratory of
Computational Physics, Institute of Applied Physics and Computational
Mathematics, Beijing 100088, China  Fuyang Tian<EMAIL_ADDRESS>Institute
for Applied Physics, University of Science and Technology Beijing, Beijing
100083, China

###### Abstract

Developing a comprehensive magnetic theory of correlated itinerant magnets is
a challenging task due to the difficulty in reconciling both local moments and
itinerant electrons. In this work, we investigate the microscopic process of
magnetic phase transition in ferromagnet metal Fe3-δGeTe2. A new paradigm is
proposed to describe the magnetic phase transition in correlated metallic
ferromagnets, where Hund’s coupling dominates the spectral weight transfer
between different spin channels, rather than spin-splitting as described by
the Stoner model. We recognize that our theory should be universal for
itinerant magnets. Additionally, we reveal an efficient way to achieve novel
quantum states from various competing orders in multi-site crystal structures.
Our research shows that Fe1 are proximate to Mott physics, while Fe2 exhibit
Hund physics due to their distinct atomic environments. These competing orders
work together to produce heavy fermion behavior within ferromagnetic long-
range order through well-defined quasiparticle bands, which are promoted by
Hund’s coupling and further hybridized with relative itinerant bands. The
complex interactions of competing orders drive correlated magnetic metal to a
new frontier for discovering outstanding quantum states and exotic phenomena
in condensed matter physics.

Magnetism plays a crucial role in condensed matter physics, especially in
strongly correlated systems. It is closely linked to various emergent
phenomena such as heavy fermion behavior [1, 2], unconventional
superconductivity [3, 4], Hund metal [5, 6], and quantum phase transition [7,
8]. These fascinating behaviors arise due to the dual nature of correlated
electrons, which creates challenges in understanding the role of magnetism.
One essential scenario is the delocalization of $f$-electrons at low
temperatures in heavy fermion systems, achieved through the hybridization of
local moments with conductive electrons via RKKY interaction [9]. Another
critical instance is the $d$-electron Hund metals, which appear even more
mysterious as the systems do not have explicit two fluids [10]. The odd
coexistence of itinerant and localized behaviors in Hund metals, dominated by
Hund’s coupling, has scored significant victories in iron-based
superconductors [5, 11, 12, 13, 14]. However, from a theoretical perspective,
the two commonly used models, one starting from the Heisenberg local moment
model and the other from the Stoner itinerant electron model [15, 16], have
their limitations in treating these systems. The microscopy mechanism of
considering both local moments and itinerant electrons in correlated magnetic
metal still needs to be explored.

Recently, van der Walls (vdW) itinerant ferromagnet Fe3GeTe2 (FGT) has
garnered significant attention. It is a potential candidate for voltage-
controlled spintronic devices due to its exceptional exfoliation properties
and robust high Curie temperature [17, 18]. In bulk material, the
ferromagnetic (FM) Curie temperature $T_{c}$ is around 220 K, and it can even
reach room temperature with an ionic gate in its thin flakes [17, 19, 20].
Furthermore, experiments have observed significant electron mass enhancement
[20, 21, 22, 23], which is also supported by previous theoretical calculations
[24]. Considering its high ferromagnetic transition temperature and strongly
correlated effects, the itinerant ferromagnet FGT presents an ideal platform
for investigating the dual nature of magnetism and its interplay with other
exotic phenomena [25].

Here, we highlight some mysteries that may be related to the issue of
magnetism in FGT. First and foremost is the question of how to describe the
Fe-$d$ electrons accurately. Although the itinerant Stoner model is commonly
accepted [20, 26], recent experiments suggest that the Heisenberg model may be
necessary to describe the ferromagnetism in FGT correctly [27, 28, 29, 30].
For instance, the angel-resolved photoemission spectroscopy (APRES)
experiments find that the temperature-dependent evolution towards the
ferromagnetic transition is strikingly weak [27, 30]. Secondly, as is commonly
discussed in iron-based superconductors [5, 31, 32], it is unclear whether the
irons in FGT exhibit Hundness or Mottness. Kim et al. suggest that it is a
“site-differentiated” Hund metal [33], while another study suggests that it is
close to an orbital selective Mott transition (OSMT) [34]. Finally, previous
studies have reported heavy fermion behavior in FGT with a large Sommerfeild
coefficient ($\gamma_{0}=135\,$mJ/mol K2) [20, 21, 22]. This is striking since
$d$-electron heavy fermion systems are rare [35, 36, 37]. While in FGT, heavy
fermion behavior appears to be promoted within ferromagnetic order [22]. These
controversies encourage us to explore the complex interaction of magnetism
with possible Hund, Mott, or Kondo physics in this system.

In this work, we systematically investigate the correlation effects in
itinerant ferromagnet Fe3-δGeTe2 using density functional theory combined with
dynamical mean-field theory (DFT+DMFT). Our findings reveal a new paradigm of
magnetic transition in correlated metallic ferromagnets, characterized by the
spectral weight transfer between two spin channels caused by the correlation
effect, rather than spin splitting in the mean-field Stoner model. Similarly
to Hund metals, the quasiparticles form in itinerant ferromagnet as
temperature decreases. Interestingly, the Fe atoms in FGT exhibit two distinct
behaviors: one follows Mott physics, while the other exhibits Hundness. This
suggests that FGT may be the first material to exhibit both Hundness and
Mottness, stemming from the multi-site nature of correlated atoms in its
crystal structure. Furthermore, our results demonstrate that Hund’s coupling
localizes quasiparticles and promotes their hybridization with other itinerant
bands, leading to heavy fermion behavior. Our findings not only naturally
reconcile the conflicts between heavy fermion behavior and long-range magnetic
order, but also open up a new frontier for discovering new quantum states in
itinerant magnets.

Figure 1:  Crystal structure and DFT+U electronic structures in FGT. a,
Illustration of the crystal structure of Fe3GeTe2. b, The electronic
structures and partial density of states of non-spin polarized calculations of
Fe3GeTe2 in DFT+U calculations. c, The electronic structures and partial
density of states of ferromagnetic order in DFT+U calculations.

Results
DFT+U calculations
Fe3-δGeTe2 is a ternary multi-site ferromagnet, as illustrated in Fig. 1a,
with layered Fe3Ge substructures sandwiched by two layers of Te atoms and a
vdW gap between adjacent Te layers. FGT crystallizes into a hexagonal
structure with the space group of $P6_{3}/mmc$ (No. 194), and Fe atoms occupy
two inequivalent Wyckoff sites, labeled Fe1 and Fe2. Notably, the occupational
deficiency ($0\textless\delta\textless 0.3$) occurs only at the Fe2 site [19,
21]. Our DFT results are consistent with available calculations (Supplementary
Information Fig.1) [21, 34]. Considering the strong correlation effects in
FGT, we perform DFT+U calculations using the correlation values suggested by
previous work [24]. In non-magnetic (NM) calculations, the Fe-$d$ electrons
dominate near the Fermi level in Fig. 1b. However, in ferromagnetic order, the
correlation effects strongly localize Fe-$d$ electrons, causing the Fe-bands
to be pushed away from the Fermi level. As a result, the Fe-bands in the spin-
up channel are almost entirely filled and pushed down to 5 eV below the Fermi
level, leaving only itinerant bands with Te character in Fig 1c. In contrast,
the bands in the spin-down channel move up slightly, leaving a small peak of
Fe-$d$ density of states (DOS) at the Fermi level. Besides, the partial DOS of
Fe1 and Fe2 have similar peak structures, indicating strong hybridization
between them. We note that the DOS at the Fermi level in non-magnetic GGA+U
calculations is reduced to 2.9 states/eV f.u., much smaller than 13.6 in GGA
calculations [21]. The much smaller DOS in DFT+U calculations, combined with
the weak temperature-dependent evolution of Fe-$d$ bands in APRES
measurements, suggests to reconsider the application of Stoner model to FGT
[27, 30].

Figure 2: The results of DFT+DMFT calculations of FGT at low and high
temperatures. a,b, The calculated density of states and spectral functions in
high-temperature non-magnetic order and low-temperature ferromagnetic order.
In ferromagnetic order, the pink arrow in DOS shows the region of the spectral
weight increasing, and the green arrow shows the spectral weight decreasing.
c, The real and imaginary part of self-energy of Fe1 and Fe2 atoms in no-
magnetic order. d, The several highest probabilities of the Fe atomic states
from the impurity solver in non-magnetic calculations. The $\uparrow$ arrow
means $S_{z}$ is positive and the $\downarrow$ arrow means $S_{z}$ is
negative. e, The real part of self-energy of Fe1 and Fe2 atoms in
ferromagnetic order. f, The several highest probabilities of the Fe atomic
states from the impurity solver in ferromagnetic order.

Mechanism of magnetic phase transition using DFT+DMFT calculations
To further explore the electronic structures and magnetism with dynamical
correlation effects, we systematically perform DFT+DMFT calculations on FGT
[38, 39]. This method has been successfully applied to correlated materials
such as Hund metals and heavy fermion materials [9, 5]. For the magnetic
calculations using DFT+DMFT, we consider ten Fe-$d$ orbitals and introduce
symmetry breaking of the self-energy at the initial stage. As expected, our
DFT+DMFT simulations show that the system converges to non-magnetic order at
high temperature in Fig. 2a. The spectral function exhibits large blurred
regions, and the DOS shows very broad peaks, mainly from the incoherent states
of Fe atoms around the Fermi level (from -4 eV to 4 eV). This incoherence
arises from the Hubbard bands of Fe atoms, originating from the strong Coulomb
interaction. Moreover, the system converges to the ferromagnetic order at low
temperature in Fig. 2b. Compared with high-temperature results, the bands
above the Fermi level become sharp in the spin-up channel, while the bands
remain blurred below the Fermi level. Conversely, in the spin-down channel,
the upper Hubbard bands are retained, and the bands below the Fermi level
become sharp. Additionally, the bands close to the Fermi level become more
distinct in both spin channels, thus forming quasiparticles, which result in
multiple peaks in DOS around the Fermi level.

Figure 3: The results of DFT+DMFT calculations of nonstoichiometric
Fe2.82GeTe2 with temperatures. a,b, The spectral functions of the spin-up and
down channels with temperature from 1000 K to 100 K. c, The DFT+DMFT Fermi
surfaces at different temperatures. d, The total density of states of spin-up
and spin-down channels with varying temperatures. The colored arrows indicate
the regions of spectral weight transfer. e, The calculated DOS is compared
with STM measuerments [22, 23].

Here, we conclude several significant findings from further inspecting the
microscopic processes in DFT+DMFT calculations. The most important observation
is the spectral weight transfer in Hubbard bands, rather than the spin
splitting suggested by the Stoner model. The colored arrows in Fig. 2b
indicate the primary areas where spectral weight transfer occurs. This
behavior is significantly different from those observed in DFT and DFT+U
calculations. Previous ARPES experiments have strongly indicated this spectral
weight transfer behavior [26]. They found that the spectral weight transfer is
independent from the details of band structures, with the enhanced DOS at 500
meV and suppressed DOS at 200 meV below the Fermi energy in ferromagnetic
state. Furthermore, recent ARPES measurements also show a weak temperature-
dependent evolution of band shift, suggesting that local moments may play a
crucial role in ferromagnetic order [27, 30]. In short, our proposed mechanism
of spectral weight transfer agrees well with ARPES measurements.

Secondly, a crossover is observed from high-temperature incoherence states to
low-temperature coherence states in FGT, similar to what occurs in Hund metals
[5, 10]. It is supported by recent optical experiments, indicating that Hund’s
coupling is essential in FGT [28]. The incoherent Hubbard bands originate from
the significant spin scattering of Fe-$d$ electrons, which can be clearly seen
from the self-energy. The imaginary parts of the self-energy in Fig. 2c have
significant values at zero frequency, approximately 2 eV for Fe1 and 0.1 eV
for Fe2. However, at low temperatures, the imaginary parts of the self-energy
approach zero, leading to the quasiparticle bands becoming distinct near the
Fermi level, as shown in Fig. 2b. Furthermore, correlation effects result in
the renormalization of the Fe-$d$ bands near the Fermi level, which can be
estimated from the real parts of self-energy with quasiparticle weight,
$Z=(1-\partial_{\omega}\rm{Re}\Sigma|_{\omega=0})^{-1}$. Table 1 lists the
mass enhancements ($m^{*}/m_{\rm LDA}$) of each orbital, calculated as
$Z^{-1}$ in DMFT. Our results demonstrate that all the Fe-$d$ orbitals are
renormalized in both non-magnetic and ferromagnetic phases.

Thirdly, it is vital to shed light on the differences between FGT and other
paramagnetic Hund metals [12, 40]. We focus on the detailed evolution of
electronic structures during the magnetic transition. The highest atomic
probabilities from the DMFT impurity solver is shown in Fig. 2d. Despite the
fact that the average total spin of Fe1 and Fe2 are zero in paramagnetic
order, the most significant cases are $S_{z}=\pm 5/2$, and they have the same
probabilities. This implies that the Fe atoms have large local moments and
fluctuate significantly over time. In the ferromagnetic state, the highest
probabilities of Fe1 and Fe2 are primarily composed of only positive values of
$S_{z}$ in Fig. 2f, indicating the formation of magnetic order. Compared with
paramagnetic Hund metals, we recognize that the long-range magnetic order will
rapidly suppress the spin fluctuations, making it easier to form
quasiparticles in the magnetic phase.

Finally, we recognize that Fe1 and Fe2 atoms exhibit different physical
behavior, with Fe1 being described by Mottness and Fe2 characterized by
Hundness. This difference is evident in the histogram associated with each
atomic configuration (Supplementary Information Fig.2). In Fig. 2d and f, Fe1
appear concentrated in a single charge state, $N=5$, typical of the Mott type
[10]. However, the histogram of Fe2 extend over many charge states, a direct
signature of Hundness [5]. Another piece of evidence is the feature of DOS.
The partial DOS of Fe1 has a pseudogap at the Fermi level between two broad
humps near $\pm$1 eV in Fig. 2a, which is very similar to the Mott system V2O3
[32]. Conversely, the partial DOS of Fe2 is characterized by a single broad
feature, similar to the Hund metal Sr2RuO4 [32]. We realize that the different
physics of Fe1 and Fe2 may be due to their local environments, with the
occupation of Fe1 and Fe2 being 5.0 and 5.8 in our DFT+DMFT calculations,
corresponding to Fe13+ and Fe22+. In multi-orbital correlated systems, Hund’s
coupling has a Janus-faced effect depending on the electron occupation [41].
As Fe1 are closed to half-filling, the Mott mechanism is promoted. In
contrast, Fe2 are away from half-filling, so the Hund interaction is dominant.
The valences of Fe atoms in our calculations are obtained using an exact
double-counting scheme in DFT+DMFT calculations [42] and are consistent with
previous predictions [17, 19].

Table 1:  The mass enhancements of Fe-orbitals at paramagnetic and
ferromagnetic order in FGT. $m_{{}_{\rm DFT}}$ represents the mass from DFT
calculations. The respective Fe-3$d$ states in FGT split into three classes: a
$d_{z^{2}}$-like ${\rm a_{1g}}$ orbital, two degenerate ${\rm e_{g}^{,}}$
orbitals and two degenerate ${\rm e_{g}}$ orbitals. The $\uparrow$ arrow means
spin-up channel and the $\downarrow$ arrow means spin-down channel.
| $m^{*}_{\rm a_{1g}}/m_{{}_{\rm DFT}}$ | $m^{*}_{\rm e_{g}^{,}}/m_{{}_{\rm DFT}}$ | $m^{*}_{\rm e_{g}}/m_{{}_{\rm DFT}}$
---|---|---|---
PM | Fe1 | 9.1 | 4.5 | 7.6
Fe2 | 1.5 | 1.8 | 2.3
FM | Fe1 | $\uparrow$ | 1.4 | 1.4 | 1.3
$\downarrow$ | 1.3 | 1.3 | 1.2
Fe2 |  $\uparrow$ | 1.4 | 1.4 | 1.3
$\downarrow$ | 1.6 | 1.4 | 2.0

Effect of Fe vacancy on electronic structures
We have noticed a systematic discrepancy between our electronic structures of
FGT and ARPES measurements [22, 27, 30]. Although surface effects could
partially explain the discrepancy [43], the nonstoichiometry in FGT is
significant due to the nonnegligible vacancy of Fe2 atoms [19, 21, 44]. Fig. 3
shows the temperature-dependent spectral functions in nonstoichiometric
Fe3-δGeTe2 ($\delta\sim 0.18$). As expected, the main findings from
stoichiometric FGT calculations are retained, including the spectral weight
transfer in Hubbard bands, low-temperature quasiparticles with ferromagnetic
order, and the coexistence of Hundness and Mottness. We observe that the band
location remains relatively constant with temperature in Fig. 3a and b, which
agrees well with recent APRES measurements [27, 30]. Furthermore, the Fermi
surfaces in Fig. 3c consist of several circular and hexagonal-shaped Fermi
pockets near the $\Gamma$ point and elliptical-shaped pockets at the $K$
point. Moreover, there is a complex-shaped Fermi pocket near the $M$ point in
the spin-down channel, which is consistent with experiments [22, 27, 30]. At
low temperatures, some quasiparticle bands emerge in the spin-up channel near
-0.5 eV, also observed in the ARPES experiments at the same binding energy
[22, 27]. Fig. 3d shows that spectral weight transfer occurs as temperatures
decrease. Meanwhile, the DOS has many small sharp peaks due to the complex
flat quasiparticle bands formed at low temperatures. The agreement between our
hole-doping calculations and ARPES measurements suggests the vital role of Fe2
vacancies in FGT.

Origin of heavy fermion behavior in ferromagnets
Another critical issue is the emergence of heavy fermion behavior within
ferromagnetic order [22], as the two often contradict each other. DFT+DMFT
calculations in Fig. 3e show similar peaks near the Fermi level compared to
STM measurements [22, 23]. As shown in Fig. 3d and e, DFT+DMFT calculations
illustrate some peaks just above the Fermi level, indicated by a red circle.
As experiments suggest, these peaks possibly originate from the Kondo
resonance [22, 23]. According to our calculations, these peaks come from the
hole-like flat quasiparticle bands located at the $\Gamma$ point, promoted by
the rapid decrease of spin fluctuations in ferromagnetic order. Furthermore,
as temperature decreases, the quasiparticle bands at the $\Gamma$ point show
strong hybridization with other relatively itinerant bands, where the partial
DOS of Fe1, Fe2, and Ge atoms have similar shapes and peaks positions over a
wide range (Supplementary Information Fig.3). These $\Gamma$-centered hole-
like bands have a small Fermi velocity in Fig. 3c. Based on this evidence, we
conclude that Hund’s coupling can promote the formation of ferromagnetic order
and heavy fermion behavior simultaneously. This mechanism demonstrates, for
the first time, the complex competition and cooperation of magnetic order,
Hundness, Mottness, and heavy fermion behavior in a detailed manner.

Figure 4: The schematic diagram of electronic structure evolution in itinerant
ferromagnets. The left and right sides show the local Heisenberg model and the
itinerant Stoner model, respectively. In the Heisenberg model, the electrons
are localized and cannot move through the crystal. In the Stoner model,
electrons can move in the crystal, and the spin splitting comes from the
different spin channels that moves in the opposite direction. In correlated
itinerant magnets, dynamical correlation effects will lead to collective
excitation, making spectral weight transfer and quasiparticle forming through
magnetic phase transition.

Discussion
In strongly correlated physics, various exotic phenomena are closely related
to the local moments of correlated atoms. One such outstanding material is the
strongly correlated ferromagnetic metal FGT. However, the origin of magnetism
in itinerant ferromagnets has been a controversial topic for several decades.
As shown in Fig. 4, the Heisenberg model is used to describe the origin of
magnetism in insulators [15]. In contrast, the Stoner model, based on the
mean-field approximation, provides a phenomenological description of spin
splitting in magnetic metals [16]. However, the Stoner model has been proven
inadequate for these correlated systems [27, 30]. Even though Moriya and
Kawabata attempted to consider spin fluctuations and correlation effects
through self-consistent renormalization (SCR) theory, a detailed physical
picture still needs to be provided [45, 20]. Our magnetic DFT+DMFT
calculations aim to resolve this problem from the many-body dynamical
correlation point of view. In ferromagnetic metals, Hund’s coupling can
produce collective excitation modes at low temperatures and create
renormalized flat quasiparticle bands through dynamical correlation processes.
Meanwhile, the spectral weight transfer occurs during magnetic phase
transition. Additionally, the formation of ferromagnetic long-range order, in
turn, promotes quasiparticle formation through positive feedback by rapidly
reducing spin fluctuations in the microscopic mechanism. Ultimately, this
leads to the simultaneous occurrence of spectral weight transfer, long-range
magnetic order, and even other competing orders, such as heavy fermion
behavior.

In summary, through systematic magnetic DFT+DMFT calculations on FGT, we have
uncovered a new microscopic description of electronic structures in correlated
ferromagnetic metals based on dynamical correlation effects. The Hund’s
coupling results in long-range ferromagnetic order accompanied by spectral
weight transfer, rather than the spin splitting described by the Stoner model.
This long-range magnetic order suppress spin fluctuation rapidly, leading to
the formation of renormalized quasiparticle bands. Our work also show that
heavy fermion behavior can be promoted by long-range magnetic order in
itinerant magnets, through flat quasiparticle bands around the Fermi level
further hybridizing with other itinerant bands at lower temperatures. In
addition, the two Fe sites exhibit different physics, with Fe1 characterized
by Mottness and Fe2 by Hundness, induced by the diversity of their atomic
environment. Overall, our research indicates that multi-site itinerant magnets
present a new frontier for discovering new quantum states and their
competitions.

Methods
DFT calculations
The first-principles DFT calculations were performed using the projector-
augmented wave (PAW) method [46], as implemented in the Vienna ab initio
simulation package (VASP) [47] with the Perdew–Burke–Ernzerhof (PBE)
functional [48]. The strongly correlated correction was considered with the
GGA+U method, taking the suggested correlated values as $U_{\rm eff}=$ 4 eV
from previous work [24]. We also used doubly screened Coulomb correction
(DSCC) approach to check the correlated values and got similar values [49].
The cut-off energy of the plane wave basis-set was set to 500 eV. The
Monkhorst-Pack gird was chosen as 15$\times$15$\times$3 to ensure convergence.
The total energies were converged to 10-8 eV. In our calculations, we adopted
experimental lattice parameters with $a=b=3.99\,$Å and $c=16.33\,$Å [19].

DFT+DMFT calculations
The DFT+DMFT calculations were employed to treat the electronic correlations
of Fe-$d$ orbitals. Full-potential linearized augmented plane-wave method
implemented in the WIEN2k package was used in DFT+DMFT [50]. We employed the
hybridization expansion continuous-time quantum Monte Carlo method as the
impurity solver [51, 52] and applied the exact double-counting scheme [42].
The correlation parameters were set $U=$ 5 eV and $J=$ 0.8 eV, as suggested
from earlier work [24]. We used the exact double-counting scheme, and it can
get the proper occupied number of Fe atoms, consistent with experiments [17,
19]. To account for the effect of Fe vacancy, we modified our DFT+DMFT
calculations by changing the number of electrons. Previous research has
confirmed that changing electron numbers leads to similar results in
constructing supercells with Fe2 vacancies [44, 53]. In order to perform
magnetic DFT+DMFT calculations, we consider 10 Fe-$d$ orbitals and break the
symmetry with a change in the real-part energy of the initial self-energy of
both two types of Fe atoms.

## References

  * [1] Stewart, G. R. Heavy-fermion systems. _Rev. Mod. Phys._ 56, 755-787 (1984).
  * [2] Yang, Y.-F., Fisk, Z., Lee, H. O., Thompson, J. D. & Pines, D. Scaling the Kondo lattice. _Nature_ 454, 611-613 (2008).
  * [3] Dagotto, E. Correlated electrons in high-temperature superconductors. _Rev. Mod. Phys._ 66, 763-840 (1994).
  * [4] Scalapino, D. J. A common thread: The pairing interaction for unconventional superconductors. _Rev. Mod. Phys._ 84, 1383-1417 (2012).
  * [5] Yin, Z. P., Haule, K. & Kotliar, G. Kinetic frustration and the nature of the magnetic and paramagnetic states in iron pnictides and iron chalcogenides. _Nat. Mater._ 10, 932-935 (2011).
  * [6] Haule, K. & Kotliar, G. Coherence–incoherence crossover in the normal state of iron oxypnictides and importance of Hund’s rule coupling. _New J. Phys._ 11, 025021 (2009).
  * [7] Si, Q., Rabello, S., Ingersent, K. & Smith, J. L. Locally critical quantum phase transitions in strongly correlated metals. _Nature_ 413, 804-808 (2001).
  * [8] Shen, B. et al. Strange-metal behaviour in a pure ferromagnetic Kondo lattice. _Nature_ 579, 51-55 (2020).
  * [9] Shim, J. H., Haule, K. & Kotliar, G. Modeling the localized-to-itinerant electronic transition in the heavy fermion system CeIrIn5. _Science_ 318, 1615-1617 (2007).
  * [10] Georges, A. & Kotliar, G. The Hund-metal path to strong electronic correlations. _Physics Today_ 77, 46-53 (2024).
  * [11] Georges, A., Medici, L. d. & Mravlje, J. Strong correlations from Hund’s coupling. _Annu. Rev. Condens. Matter Phys._ 4, 137-178 (2013).
  * [12] Haule, K., Shim, J. H. & Kotliar, G. Correlated electronic structure of LaO1-xFxFeAs. _Phys. Rev. Lett._ 100, 226402 (2008).
  * [13] Yin, Z. P., Haule, K.& Kotliar, G. Magnetism and charge dynamics in iron pnictides. _Nat. Phys._ 7, 294-297 (2011).
  * [14] Liu, M. et al. Nature of magnetic excitations in superconducting BaFe1.9Ni0.1As2. _Nat. Phys._ 8, 376-381 (2012).
  * [15] Heisenberg, W. Zur theorie des Ferromagnetismus. _Z. Phys._ 49, 619 (1928).
  * [16] Stoner, E. C. Ferromagnetism. _Rep. Prog. Phys._ 11, 43 (1947).
  * [17] Deng, Y. et al. Gate-tunable room-temperature ferromagnetism in two-dimensional Fe3GeTe2. _Nature_ 563, 94-99 (2018).
  * [18] Fei, Z. et al. Two-dimensional itinerant ferromagnetism in atomically thin Fe3GeTe2. _Nat. Mater._ 17, 778-782 (2018).
  * [19] Deiseroth, H. J., Aleksandrov, Reiner, K., Kienle, C., L. & Kremer, R. K. Fe3GeTe2 and Ni3GeTe2 \- Two new layered transition‐metal compounds: crystal Structures, HRTEM investigations, and magnetic and electrical properties. _Eur. J. Inorg. Chem._ 2006, 1561-1567 (2006).
  * [20] Chen, B. et al. Magnetic properties of layered itinerant electron ferromagnet Fe3GeTe2. _J. Phys. Soc. Jpn._ 82, 124711 (2013).
  * [21] Verchenko, V. Y., Tsirlin, A. A., Sobolev, A. V., Presniakov, I. A. & Shevelkov, A. V. Ferromagnetic order, strong magnetocrystalline anisotropy, and magnetocaloric effect in the layered telluride Fe3-δGeTe2. _Inorg. Chem._ 54, 8598-8607 (2015).
  * [22] Zhang, Y. et al. Emergence of Kondo lattice behavior in a van der Waals itinerant ferromagnet, Fe3GeTe2. _Sci. Adv._ 4, eaao6791 (2018).
  * [23] Zhao, M. et al. Kondo holes in the two-dimensional itinerant Ising ferromagnet Fe3GeTe2. _Nano Lett._ 21, 6117-6123 (2021).
  * [24] Zhu, J.-X. et al. Electronic correlation and magnetism in the ferromagnetic metal Fe3GeTe2. _Phys. Rev. B_ 93, 144404 (2016).
  * [25] Kim, K. et al. Large anomalous Hall current induced by topological nodal lines in a ferromagnetic van der Waals semimetal. _Nat. Mater._ 17, 794-799 (2018).
  * [26] Zhuang, H. L., Kent, P. R. C. & Hennig, R. G. Strong anisotropy and magnetostriction in the two-dimensional Stoner ferromagnet Fe3GeTe2. _Phys. Rev. B_ 93, 134407 (2016).
  * [27] Xu, X. et al. Signature for non-Stoner ferromagnetism in the van der Waals ferromagnet Fe3GeTe2. _Phys. Rev. B_ 101, 201104(R) (2020).
  * [28] Corasaniti, M. et al. Electronic correlations in the van der Waals ferromagnet Fe3GeTe2 revealed by its charge dynamics. _Phys. Rev. B_ 102, 161109(R) (2020).
  * [29] Bao, S. et al. Neutron spectroscopy evidence on the dual nature of magnetic excitations in a van der Waals metallic ferromagnet Fe2.72GeTe2. _Phys. Rev. X_ 12, 011022 (2022).
  * [30] Wu, H. et al. Spectral evidence for local-moment ferromagnetism in the van der Waals metals Fe3GaTe2 and Fe3GeTe2. _Phys. Rev. B_ 109, 104410 (2024).
  * [31] Song, Y. et al. A Mott insulator continuously connected to iron pnictide superconductors. _Nat. Commun._ 7, 13879 (2016).
  * [32] Deng, X. et al. Signatures of Mottness and Hundness in archetypal correlated metals. _Nat. Commun._ 10, 2721 (2019).
  * [33] Kim, T. J., Ryee, S. & Han, M. J. Fe3GeTe2: a site-differentiated Hund metal. _npj Comput. Mater._ 8, 245 (2022).
  * [34] Bai, X. et al. Antiferromagnetic fluctuations and orbital-selective Mott transition in the van der Waals ferromagnet Fe3-xGeTe2. _Phys. Rev. B_ 106, L180409 (2022).
  * [35] Kondo, S. et al. LiV2O4: a heavy fermion trasition metal oxide. _Phys. Rev. Lett._ 78, 3729-3732 (1997).
  * [36] Cheng, J.-G. et al. Possible Kondo physics near a metal-insulator crossover in the a-site ordered perovskite CaCu3Ir4O12. _Phys. Rev. Lett._ 111, 176403 (2013).
  * [37] Takegami, D. et al. CaCu3Ru4O12: A high-Kondo-temperature transition-metal oxide. _Phys. Rev. X_ 12, 011017 (2022).
  * [38] Kotliar, G. et al. Electronic structure calculations with dynamical mean-field theory. _Rev. Mod. Phys._ 78, 865-951 (2006).
  * [39] Haule, K., Yee, C. & Kim, K. Dynamical mean-field theory within the full-potential methods: Electronic structure of CeIrIn5,CeCoIn5, and CeRhIn5. _Phys. Rev. B_ 81, 195107 (2010).
  * [40] Mravlje, J. et al. Coherence-incoherence crossover and the mass-renormalization puzzles in Sr2RuO4. _Phys. Rev. Lett._ 106, 096401 (2011).
  * [41] Medici, L. de’, Mravlje, J. & Georges, A. Janus-faced influence of Hund’s rule coupling in strongly correlated materials. _Phys. Rev. Lett._ 107, 256401 (2011).
  * [42] Haule, K. Exact double counting in combining the dynamical mean field theory and the density functional theory. _Phys. Rev. B_ 115, 196403 (2007).
  * [43] Wang, Y.-C. et al. First-principles study of the role of surface in the heavy-fermion compound CeRh2Si2. _Phys. Rev. B_ 103, 165140 (2021).
  * [44] Jang, S. W. et al. Origin of ferromagnetism and the effect of doping on Fe3GeTe2. _Nanoscale_ 12, 13501 (2020).
  * [45] Hatatani, M. & Moriya, T. Ferromagnetic spin fluctuations in two-dimensional metals. _J. Phys. Soc. Jpn._ 64, 3434 (1995).
  * [46] Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 50, 17953–17979 (1994).
  * [47] Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 54, 11169-11186 (1996).
  * [48] Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. _Phys. Rev. Lett._ 77, 3865-3868 (1996).
  * [49] Liu, B., Wang, Y., Liu, Y., Liu, H. & Song, H. Doubly screened Coulomb correction approach for strongly correlated systems. _J. Phys. Chem. Lett._ 14, 8930-8939 (2023).
  * [50] Blaha, P. et al. _WIEN2K, An augmented plane wave + local orbitals program for calculating crystal properties_ (Karlheinz Schwarz, Techn. Universität Wien, Austria, 2023).
  * [51] Werner, Comanac, P. A., Medici, L. de’, Troyer, M. & Millis, A. J. Continuous-time solver for quantum impurity models. _Phys. Rev. Lett._ 97, 076405 (2006).
  * [52] Haule, K. Quantum Monte Carlo impurity solver for cluster dynamical mean-field theory and electronic structure calculations with adjustable cluster base. _Phys. Rev. B_ 75, 155113 (2007).
  * [53] Lee, J.-E. et al. Electronic structure of above-room-temperature van der Waals ferromagnet Fe3GaTe2. _Nano Lett._ 23, 11526-11532 (2023).

Acknowledgements
The authors thank Yingying Cao, Beilei Liu, Huican Mao and Kun Zhai for the
fruitful discussions. This work is supported by the National Natural Science
Foundation of China (Grant Nos. U2230401, 12204033, 52371174), the National
Key Research and Development Program of China (Grant No. 2021YFB3501503), the
Foundation of LCP, the Fundamental Research Funds for the Central Universities
(Grant No. FRF-TP-22-097A1) and the Young Elite Scientist Sponsorship Program
by BAST (Grant No. BYESS2023301). We thank the Tianhe platforms at the
National Supercomputer Center in Tianjin.

Author contributions
F. T. and H. S. conceived the idea and supervised the project; Y. X. and Y. W.
performed the calculations; X. J., Y. L. and H. L. participated in the
discussions. Y. X. wrote the paper with input from all authors. All authors
contributed to the interpretation of the data and revised the paper. Y. X. and
Y. W. contributed equally to this work.