# Highly nonperturbative nature of the Mott metal-insulator transition:
Two-particle vertex divergences in the coexistence region

M. Pelza,b S. Adlera M. Reitnera A. Toschia aInstitute of Solid State Physics,
TU Wien, 1040 Vienna, Austria  bArnold Sommerfeld Center for Theoretical
Physics, Center for NanoScience, and Munich Center for Quantum Science and
Technology, Ludwig-Maximilians-Universität München, 80333 Munich, Germany

###### Abstract

We thoroughly analyze the divergences of the irreducible vertex functions
occurring in the charge channel of the half-filled Hubbard model in close
proximity to the Mott metal-insulator transition (MIT). In particular, by
systematically performing dynamical mean-field theory (DMFT) calculations on
the two-particle level, we determine the location and the number of the vertex
divergences across the whole coexistence region adjacent to the first-order
metal-to-insulator transition. We find that the lines in the parameter space,
along which the vertex divergences occur, display a qualitatively different
shape in the coexisting metallic and insulating phase, which is also
associated to an abrupt jump of the number of divergences across the MIT.
Physically, the systematically larger number of divergences on the insulating
side of the transition reflects the sudden suppression of local charge
fluctuation at the MIT. Further, a systematic analysis of the results
demonstrates that the number of divergence lines increases as a function of
the inverse temperature ${\beta\\!=\\!(k_{\mathrm{B}}T)^{-1}}$ by approaching
the Mott transition in the zero temperature limit. This makes it possible to
identify the zero-temperature MIT as an accumulation point of an infinite
number of vertex divergence lines, unveiling the highly nonperturbative nature
of the underlying transition.

††preprint:

##  I Introduction

The multifaceted manifestations [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] of the
breakdown of the self-consistent perturbation theory for the many-electron
problem have been recently in focus of several studies [11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. In particular, it has
been demonstrated [10] how the breakdown of the perturbative expansion
corresponds to the crossings of different solutions [3, 4] of the Luttinger-
Ward functional or, equivalently, to multiple divergences [1, 9, 13] of the
two-particle vertex functions irreducible in the charge channel, i.e. the
kernel of the Bethe-Salpeter equation describing the charge response of the
many-electron system under consideration.

On a less formal perspective, the crucial role [21, 23, 25, 27] of the local
moment formation in triggering the perturbation-theory breakdown as well as
the contra-intuitive physical consequences associated to nonperturbative
scattering processes [20] have been extensively investigated in the most
recent literature, e.g. for the Anderson impurity and the Hubbard model.

Hitherto, however, the analysis of an equally important aspect of this
problem, namely the precise relation linking the above mentioned
manifestations of the perturbative breakdown to the occurrence of Mott-Hubbard
metal-to-insulator transitions (MITs) [29], has been put, to some extent,
aside. In fact, in some of the earliest works [1, 9] on this subject, it was
suggested that the divergences of the irreducible vertex could be viewed as
“precursors” of the Mott-Hubbard MIT. In later studies, however, it was shown
[13] that multiple irreducible vertex divergences were also occurring in
cases, such as the Anderson Impurity model (AIM), where no Mott-Hubbard
transition takes place. The contradiction of this observation with the
previously proposed interpretation has left the full understanding of this
aspect of the perturbative breakdown unsolved.

In this work, we will hence address the still outstanding question of how the
Mott MIT, which represents an intrinsically nonperturbative phenomenon, is
actually related to the breakdown of the perturbation expansion by analyzing
one of its characterizing manifestations, the divergences of the irreducible
vertex functions.

To this aim, we will perform systematic dynamical mean-field theory (DMFT)
[30] calculations of the two-particle Green’s/vertex functions [31] of the
Hubbard model in one of its most delicate parameter regimes, the coexistence
region of the Mott MIT, which was not considered in preceding studies on this
specific topic.

In particular, we will determine the location, the number, and the properties
of the vertex divergences occurring in both the (coexisting) metallic and
insulating DMFT solution of the half-filled Hubbard model in the proximity of
the Mott transition. Hereby, we focus on the changes taking place across the
finite-$T$ first-order transition and on the extrapolation of the
corresponding results towards the ${T\\!=\\!0}$ second-order (quantum)
critical endpoint of the MIT. This procedure will eventually allow to draw
rigorous conclusions about the link between vertex divergences and the Mott-
Hubbard MIT.

The plan of the paper is the following: In Sec. II, we introduce the model and
the formalism necessary for our analysis and briefly recapitulate results of
previous studies relevant for our scopes; in Sec. III we illustrate our DMFT
results for the divergences of the irreducible vertex functions systematically
obtained in the coexistence region of the Mott MIT and the corresponding
extrapolation performed down to zero-temperature; then in Sec. IV we will
discuss the overall scenario emerging from our study, and, eventually, in Sec.
V we will present the conclusions and the outlook of our work.

##  II Model, Formalism and Methods

###  II.1 Model

In this study we consider a single-band Hubbard model (HM) [32] on the Bethe
lattice with infinite connectivity, whose density of state is semicircular
with half-bandwidth ${D\\!=\\!1}$, which serves as unit of energy throughout
the work.

The Hamiltonian reads

$\mathcal{H}=-t\sum_{\left<i\,j\right>\,\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}\,,$
(1)

where ${t\\!=\\!\frac{1}{2D}}$ is the (spin-independent) nearest neighbor
hopping between neighboring lattice sites $i$ and $j$, and
${c^{\dagger}_{i\sigma}}$, (${c_{i\sigma}}$) the fermionic creation
(annihilation) operator with spin ${\sigma\\!=\,\uparrow,\downarrow}$ at site
$i$, and $U$ is the local (Hubbard) repulsive interaction between two
electrons on the same lattice site
(${n_{i\sigma}\\!=\\!c^{\dagger}_{i\sigma}c_{i\sigma}}$ denoting the particle-
number operator at site $i$ for spin $\sigma$).

We set the density to half-filling
${\big{(}\big{<}n_{\uparrow}\big{>}\\!=\\!\big{<}n_{\downarrow}\big{>}\\!=\\!\frac{1}{2}\big{)}}$,
where the HM we consider, which can be exactly solved by means of dynamical
mean-field theory (DMFT), is known to feature a paradigmatic realization of
the Mott-Hubbard MIT transition. Specifically, we briefly recall here [30]
that, by neglecting possible SU(2) symmetry-broken phases, the DMFT solution
of Eq. (1) yields a first-order transition between a paramagnetic metallic
(PM) and a paramagnetic insulating phase (PI) along a $U_{c}(T)$ transition
line, ending with second-order critical points at finite
${T\\!=\\!T_{c}\\!\approx\\!\frac{1}{39}}$ and at ${T\\!=\\!0}$, as
schematically shown in the leftmost panel of Fig. 1. The first-order nature of
the transition is witnessed by the presence of a broad hysteresis region,
delimited by the lines $U_{c_{1}}(T)$ (${<\\!U_{c}(T)}$ on the left side) and
$U_{c_{2}}(T)$ (${>\\!U_{c}(T)}$ on the right side) where a PM and a PI
numerical DMFT solution coexists [33], with distinct physical properties (as
exemplified by the qualitatively different shape of the corresponding one
particle Green’s functions, shown in the inset of the leftmost panel in Fig.
1).

Figure 1: Location of the irreducible vertex divergence lines in the phase-
diagram of Hubbard model on a square lattice solved by dynamical mean field
theory (DMFT). Left panel: Schematic illustration of the MIT and its
coexistence region in DMFT. The insets show examples of the imaginary part of
the one-particle Green’s function on the Matsubara frequency axis for the PM
(left inset) and the PI (right inset) DMFT solution, featuring completely
different low-energy behaviors. The arrows represent the scanning direction
for the respective convergent solution of the two phases in the coexistence
region. Middle panel: Reproduced from Ref. [9], solid red and orange lines
[34] mark the first five $\bm{\Gamma}_{c}^{\infty}$-lines of the HM bending
around the metal-insulator transition (blue solid line). Red and orange dashed
lines show the first two $\bm{\Gamma}_{c}^{\infty}$-lines of the Hubbard atom
(HA) for comparison. Right panel: Close-up of the region near the MIT marked
in the middle panel. The blue shaded area indicates the coexistence region of
the MIT. The thermodynamic phase transition ($U_{c}$) is marked in blue and
the borders of the coexistence region ($U_{c1}$ and $U_{c2}$) are displayed as
dark blue lines. To the left in red and orange the closest vertex divergence
lines from [9] are visible.

###  II.2 Formalism

As mentioned in the Introduction, an evident manifestation of the breakdown of
the self-consistent perturbation expansion in many-electron problems is the
divergence of the kernel of the Bethe-Salpeter equation (BSE) for the system
response in the charge sector. Hence, the central quantity for our DMFT
investigation will be the on-site generalized susceptibility, which, in the
imaginary time-domain, is defined as:

$\displaystyle\bm{\chi}_{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}}(\tau_{1},\tau_{2},\tau_{3},\tau_{4}):=$
$\displaystyle\hskip
28.45274pt\bigg{[}\big{<}\mathcal{T}_{\tau}c^{\dagger}_{\sigma_{1}}(\tau_{1})c_{\sigma_{2}}(\tau_{2})c^{\dagger}_{\sigma_{3}}(\tau_{3})c_{\sigma_{4}}(\tau_{4})\big{>}\,-$
$\displaystyle\hskip
56.9055pt\big{<}\mathcal{T}_{\tau}c^{\dagger}_{\sigma_{1}}(\tau_{1})c_{\sigma_{2}}(\tau_{2})\big{>}\big{<}\mathcal{T}_{\tau}c^{\dagger}_{\sigma_{3}}(\tau_{3})c_{\sigma_{4}}(\tau_{4})\big{>}\bigg{]}$
(2) $\displaystyle\hskip
28.45274pt:=G^{(2)}_{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})\,-$
$\displaystyle\hskip
56.9055ptG_{\sigma_{1},\sigma_{2}}(\tau_{1},\tau_{2})\,G_{\sigma_{3},\sigma_{4}}(\tau_{3},\tau_{4}),$
(3)

in terms of the one ${G_{\sigma_{1},\sigma_{2}}(\tau_{1},\tau_{2})}$ and two-
particle
${G^{(2)}_{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})}$
local Green’s functions of the DMFT solution, where $\sigma_{i}$ denotes the
spins of the incoming/outgoing particles and $\tau_{i}$ their imaginary time
arguments [31]. By taking the Fourier transform into Matsubara frequencies,
and exploiting the time-translational invariance as well as the SU(2) symmetry
of the problem, one then obtains [31]

$\displaystyle\bm{\chi}_{ph,\sigma\sigma^{\prime}}^{\nu\nu^{\prime}\Omega}=\int_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}\,e^{-i\nu\tau_{1}}e^{i(\nu+\Omega)\tau_{2}}e^{-i(\nu^{\prime}+\Omega)\tau_{3}}\,\cdot$
$\displaystyle\hskip
42.67912pt\bm{\chi}_{\sigma,\sigma,\sigma^{\prime},\sigma^{\prime}}(\tau_{1},\tau_{2},\tau_{3},0),$
(4)

where, in the particle-hole convention we adopted, $\nu$, $\nu^{\prime}$
($\Omega$) denote fermionic (bosonic) Matsubara frequencies, respectively. We
recall that the expression in Eq. (4) can be directly connected to the
physical charge response $\chi_{c}(\Omega)$ of the system by performing a
summation over both fermionic Matsubara frequencies $\nu$, $\nu^{\prime}$ and
spin indices:

$\chi_{c}(\Omega)=\frac{1}{\beta^{2}}\,\sum_{\nu,\nu^{\prime}}\left(\bm{\chi}_{ph,\uparrow\uparrow}^{\nu\nu^{\prime}\Omega}+\bm{\chi}_{ph,\uparrow\downarrow}^{\nu\nu^{\prime}\Omega}\right).$
(5)

Consistent with the SU(2) symmetry of the problem, the BSE of the generalized
response in the charge sector
${\bm{\chi}_{c}^{\nu\nu^{\prime}\Omega}\\!=\\!\bm{\chi}_{ph,\uparrow\uparrow}^{\nu\nu^{\prime}\Omega}\\!+\\!\bm{\chi}_{ph,\uparrow\downarrow}^{\nu\nu^{\prime}\Omega}}$
can be recasted, for each value of the bosonic frequency $\Omega$, in a closed
matrix equation for the fermionic Matsubara frequency space, as
diagrammatically shown in Fig. 2.

Figure 2:  Diagrammatic representation of the Bethe-Salpeter equation (BSE)
${\bm{\chi}\\!=\\!\bm{\chi}_{0}\\!-\\!\bm{\chi}_{0}\bm{\Gamma}\bm{\chi}}$,
where $\bm{\chi}$ is the generalized susceptibility, $\bm{\chi}_{0}$ the
bubble term ${\bigl{(}\propto\\!-\\!\beta
G(\nu)\,G(\nu\\!+\\!\Omega)\,\delta_{\nu\nu^{\prime}}\bigr{)}}$, and
$\bm{\Gamma}$ contains all irreducible vertex diagrams [9, 35].

The equation describes the infinitely repeated insertions of (two-particle)
irreducible vertex corrections (${\bm{\Gamma}_{c}}$) to the independent
propagation of a particle and a hole (i.e., the so-called bubble term:
${\bm{\chi}_{c,0}^{\nu\nu^{\prime},\Omega}\\!=\\!-\beta\,G(\nu)\,G(\nu\\!+\\!\Omega)\,\delta_{\nu\nu^{\prime}}}$,
where $G(\nu)$ denotes the one-particle local DMFT Green’s function on the
Matsubara axis). Hence, the irreducible vertex function, which represents the
kernel of the BSE, is defined (and can be explicitly computed) as [9]

$\bm{\Gamma}_{c}^{\nu\nu^{\prime}\Omega}=\beta^{2}\bigg{(}\big{[}\bm{\chi}_{c}^{\Omega}\big{]}_{\nu\nu^{\prime}}^{-1}-\big{[}\bm{\chi}_{c,0}^{\Omega}\big{]}_{\nu\nu^{\prime}}^{-1}\bigg{)}.$
(6)

In this work, we consider the case of zero bosonic frequency, which is linked,
according to Eq. (5), to the isothermal (or static) charge response [36, 37]
of the system, as the corresponding vertex divergences are the first to occur
by increasing values of $U$ and are those directly related [10] to the
crossing of solutions of the Luttinger-Ward functional. Hence, for the sake of
readability, we will denote
${\bm{\Gamma}_{c}^{\nu\nu^{\prime}\Omega=0}\equiv\bm{\Gamma}_{c}^{\nu\nu^{\prime}}}$
and
${\bm{\chi}_{c}^{\nu\nu^{\prime}\Omega=0}\equiv\bm{\chi}_{c}^{\nu\nu^{\prime}}}$.

By a quick inspection of Eq. (6), one understands that, at any finite
temperature, no divergence of $\bm{\Gamma}_{c}$ can be caused by the inversion
of the (frequency diagonal) bubble term, as ${G(\nu)\\!\neq\\!0}$ for all
$\nu$. As a consequence, all divergences of $\bm{\Gamma}_{c}$ must originate
from the non-invertibility of the generalized susceptibility matrix or, more
formally, from the vanishing of one of its eigenvalues $\lambda_{i}$ [1, 9].
Hence, all points $(T,U)$, where a specific eigenvalue $\lambda_{i}$ of
$\bm{\chi}_{c}$ crosses zero, define a line in the phase-diagram of the half-
filled HM, where the corresponding irreducible charge vertex diverges
($\bm{\Gamma}_{c}^{\infty}$-line). It is important to stress, here, that since
for small $U$, ${\bm{\chi}_{c}\\!\approx\\!\bm{\chi}_{c,0}}$ , and the latter
is a positive definite (diagonal) matrix, the number of negative eigenvalues
($N_{\lambda<0}$) of the generalized charge susceptibility matrix at a
specific parameter set corresponds to the number of crossed
$\bm{\Gamma}_{c}^{\infty}$-lines coming from $U\\!=\\!0$ and can be also used
to approximate the shape of the $\bm{\Gamma}_{c}^{\infty}$-lines in close
proximity to the MIT (where they will become particularly dense) and analyze
their behavior towards ${T\\!\rightarrow\\!0}$.

###  II.3  Limitations of existing results

The occurrence of divergences of the irreducible vertex functions in several
fundamental many-electron problems has been recently investigated in different
publications [9, 12, 13, 14, 16, 17]. In particular, one of the most
systematic analysis made for the case of the half-filled HM on a square
lattice solved in DMFT (essentially equivalent [38] to the case considered
here) has been reported in [9], whose results are summarized in the
central/right most panels of Fig. 1. Specifically, the right panel of the
figure, a zoom of the parameter region marked by a black box in the DMFT
phase-diagram of the central panel, shows the first few
$\bm{\Gamma}_{c}^{\infty}$-lines of the DMFT solution of the HM, marked in red
or orange, depending whether the corresponding vertex divergence occurs in the
charge sector only or, simultaneously, in the charge and in the particle-
particle channels [39]. As it was already noted [9], all the
$\bm{\Gamma}_{c}^{\infty}$-lines tend to approach the corresponding vertex
divergences lines of the Hubbard atom (straight dashed lines starting from
${U\\!=\\!0}$) for large values of $T$ and $U$, while displaying a clear
bending around the Mott MIT.

Evidently the clear bending of the first divergences lines, and their somewhat
similar shape as the rightmost border of the MIT hysteresis $U_{c_{2}}(T)$,
may suggest the originally proposed interpretation of the occurrence of the
vertex divergences as precursors of the Mott MIT itself. However, as already
mentioned in the Introduction, subsequent studies have demonstrated the
occurrence of similar vertex divergences in models, such as the AIM [13],
where no Mott MIT takes place, and ascribed [8, 13, 21, 27] them to
suppressive effects of the on-site/impurity charge fluctuation triggered by
the formation of a local magnetic moment.

Hence, in order to clarify the nature of the relation linking the Mott MIT and
the occurrence of divergences in the irreducible vertex functions of the
charge sector, it is necessary to extend the DMFT studies performed hitherto
to the most challenging parameter regime, namely, the coexistence region
across the MIT, which represents the central goal of the present work.

###  II.4 Methods

The DMFT calculations of the generalized local charge susceptibility have been
performed by using a continuous-time quantum Monte Carlo (CT-QMC) solver [40]
to sample the one and two-particle Green’s functions in Eq. (3) for the
auxiliary AIM of the corresponding self-consistent DMFT solutions.
Specifically, we used the CT-QMC solver of the w2dynamics package [41, 42].
Further technical details about the numerical calculations are shortly
reported in Appendix A. Here, we want to concisely recall, how the PM and PI
DMFT-solutions in the coexistence region are obtained: Starting from outside
of the coexistence region, the interaction $U$ is changed step-by-step for a
fixed temperature $T$, whereby the previously converged DMFT calculations are
used as a starting point of the new self-consistent DMFT cycle (as
schematically illustrated by the two arrows in the leftmost panel of Fig. 1).
By entering the coexistence/hysteresis region, the variation-steps in $U$ must
be small (e.g. $\mathcal{O}(0.1)-\mathcal{O}(0.01)$), in order to allow for
convergence to different meta-stable DMFT solutions, depending on the initial
condition used. In this way two different solutions can be stabilized, at a
given temperature $T$, in the interval ${U_{c1}(T)\\!<\\!U\\!<\\!U_{c2}(T)}$,
the PM (PI) one being obtained along the path from the left to the right
${U_{c1}(T)\\!\to\\!U_{c2}(T)}$ (the right to the left
${U_{c2}(T)\\!\to\\!U_{c1}(T)}$). For each PM or PI converged DMFT-solution
obtained in the coexistence region, the corresponding on-site generalized
charge susceptibility is then computed as explained above, and Fourier
transformed in Matsubara frequencies. The diagonalization of their
corresponding matrix representation in the fermionic Matsubara frequencies
allows to determine the number of negative eigenvalues, $N_{\lambda<0}$,
which, as discussed in Sec. II B, corresponds to the number of crossed
$\bm{\Gamma}_{c}^{\infty}$-lines and can then be used to approximate the
$\bm{\Gamma}_{c}^{\infty}$-lines (see Appendix B for further details) in the
region of the phase-diagram close to the Mott MIT.

##  III Results

###  III.1 Metallic coexistence region

Figure 3: Phase diagrams of the MIT with PM solution in the coexistence region
(blue-shaded area) for the Hubbard model (HM) on the Bethe lattice. $U_{c}$
(blue), taken from Ref. [43], denotes the thermodynamic transition. Upper
panel: Coexistence region with phase points of performed DMFT calculations,
where green diamonds correspond to a metallic solution and red squares to an
insulating one. The numbers next to markers are $N_{\lambda<0}$ and the
background of the points within the coexistence region shows an interpolating
color scale of $N_{\lambda<0}$ for the metallic solution. Lower panel: Same
phase diagram as the upper panel, but showing the distinct
$\bm{\Gamma}_{c}^{\infty}$-lines approximated from the data of the phase
points (see Appendix B). Here, $n_{HM}$ indicates the number of crossed
$\bm{\Gamma}_{c}^{\infty}$-lines of the Hubbard model, coming from
${U\\!=\\!0}$. Dashed red and orange lines (in both panels) mark the
$\bm{\Gamma}_{c}^{\infty}$-lines of the Hubbard atom (HA) according to [9] as
reference, where $n_{HA}$ is the number of crossed lines coming from
${U\\!=\\!0}$.

Our results for the PM solutions of the MIT coexistence region are shown in
the upper panel of Fig. 3. Here the coexistence region is indicated as a blue
framed and shaded area in the $T$-$U$ phase diagram. The two-particle
calculations performed by scanning the MIT starting from the PM side are
marked, respectively, with green diamonds or red squares, depending on whether
a PM or a PI solution is found (the stabilization of the first PI solutions
evidently corresponds to the crossing of the $U_{c2}(T)$ line). The calculated
number $N_{\lambda<0}$ of negative eigenvalues of
$\bm{\chi}_{c}^{\nu\nu^{\prime}}$ is shown right below the corresponding
markers. In the background of the data points, a color scale, interpolating
between the numerical results, indicates the change of $N_{\lambda<0}$ within
the coexistence region. We emphasize here that, in order to make the
comparison between the results obtained in the different parameter regimes
more easily accessible, the same color scale has been used for all intensity-
plots shown in the paper.

In the lower panel of Fig. 3 instead, we show the corresponding vertex
divergence ($\bm{\Gamma}_{c}^{\infty}$) lines, starting with the
$\bm{\Gamma}_{c}^{\infty}$-line associated to
${N_{\lambda<0}\\!=\\!n_{HM}\\!=\\!19}$, and we adopt the same (red and
orange) color-coding [34] as introduced in Sec. II C. As reference for the
Mott insulating phase, the $\bm{\Gamma}_{c}^{\infty}$-lines of the Hubbard
atom (HA) [9, 14] are shown as dashed red and orange lines on the right (PI)
side of the coexistence region, starting with
${N_{\lambda<0}\\!=\\!n_{HA}\\!=\\!26}$.

In both panels of Fig. 3 a clear bending of the
$\bm{\Gamma}_{c}^{\infty}$-lines towards the critical point of the Mott MIT at
${T\\!=\\!0}$ ($U_{c2}^{T=0}$) is observed in the whole PM coexistence region,
corresponding to an increase of $N_{\lambda<0}$ along an increase of $U$ and
$T$. One also notices that the $\bm{\Gamma}_{c}^{\infty}$-lines occur
particularly dense for increasing $U$ at low temperatures. More specifically,
the shape of the bending reminds on the first $\bm{\Gamma}_{c}^{\infty}$-lines
encountered in the correlated metallic regime of the DMFT solution of the HM
[9] (central panel in Fig. 1), as well as the $\bm{\Gamma}_{c}^{\infty}$-lines
found in low-$T$ region of the phase-diagram of the AIM [13]. On a more
physical perspective, the $\bm{\Gamma}_{c}^{\infty}$-lines in the PM
coexistence region display a qualitative similarity with the Mott transition
line $U_{c}$ (solid blue line in Fig. 3), especially for low temperatures,
while no significant match with the effective [23] Kondo temperature
associated to the auxiliary AIM of the DMFT solution (estimated through the
so-called “fingerprints criterion” [21]) can be noticed. These observations
suggest the existence of a direct connection between the vertex divergences
and the Mott MIT itself.

###  III.2 Insulating coexistence region

Figure 4: Phase diagrams of the MIT with PI solution in the coexistence region
(blue-shaded area) for the Hubbard model (HM) on the Bethe lattice. The blue
line $U_{c}$, taken from Ref. [43], denotes the thermodynamic transition.
Upper panel: The data points, red squares for an insulating phase and green
diamonds for a metallic phase, mark the performed DMFT calculations. The
number next to the points shows the corresponding $N_{\lambda<0}$ and the
background of the data points in the coexistence region shows an interpolating
color scale of $N_{\lambda<0}$. Lower panel: Same phase diagram as the upper
panel with approximated $\bm{\Gamma}_{c}^{\infty}$-lines of the HM, $n_{HM}$
is the number of crossed lines coming from ${U\\!=\\!0}$. Dashed red and
orange lines (in both panels) mark the $\bm{\Gamma}_{c}^{\infty}$-lines of the
Hubbard atom (HA) according to [9] as reference, where $n_{HA}$ is the number
of $\bm{\Gamma}_{c}^{\infty}$-lines of the HA.

Our results for the PI solutions of the MIT coexistence region are shown in
the upper panel of Fig. 4. The two-particle calculations performed by scanning
the MIT starting from the PI side are marked, respectively, with red squares
or green diamonds, depending on whether a PI or a PM solution is found (the
stabilization of the first PM solutions evidently corresponds to the crossing
of the $U_{c1}(T)$ line).

Figure 5: $T$-$U$ diagram of the MIT coexistence regions for the Hubbard model
on the Bethe lattice showing the number of negative eigenvalues of the
generalized charge susceptibility $N_{\lambda<0}$ as interpolating color scale
for the metallic solution (left), the insulating solution (middle) and
comparing $N_{\lambda<0}$ at the thermodynamic stable (physical) transition
$U_{c}$ of Ref. [43] (right panel).

As before, the color scale shows the interpolation of $N_{\lambda<0}$ in the
coexistence region. In the lower panel of Fig. 4, the distinct
$\bm{\Gamma}_{c}^{\infty}$-lines (orange/red) are shown. The
$\bm{\Gamma}_{c}^{\infty}$-lines are now straight lines, similar to the dashed
(orange/red) lines of the HA (i.e. the extreme case of an HM with
${t\\!=\\!0}$), shown as reference on the right side of the coexistence
region. Due to the different shape of the divergences lines w.r.t. the PM
solution analyzed before, $N_{\lambda<0}$ still increases with increasing $U$,
but, differently from the PM case of Fig. 3, it increases for decreasing $T$.

At the same time, the qualitative similarity between the
$\bm{\Gamma}_{c}^{\infty}$-line shape in the HA and the PI phase of the HM
solved in DMFT cannot represent, evidently, an identity [44], reflecting the
intrinsic difference between the “perfect” localization of the HA and the
actual one of the Mott PI phase, where finite (though quite small) double-
occupancy is found even in the ground-state. In fact, by comparing the results
for the HA and the PI phase in DMFT, on a quantitative level, (e.g. by
considering the 28-th line of the HA ${n_{HA}\\!=\\!28}$, and the
corresponding line of the HM ${N_{\lambda<0}\\!=\\!n_{HM}\\!=\\!28}$ in the
lower panel of Fig. 4) a systematic shift of
HM-$\bm{\Gamma}_{c}^{\infty}$-lines w.r.t. the
HA-$\bm{\Gamma}_{c}^{\infty}$-lines can be noted. To effectively account for
this shift, one might rescale the interaction of the HA [45], by comparing the
number of $n_{HA}$ to $n_{HM}$ for all phase points within the coexistence
region, which yields an average factor of
${\eta\\!=\\!\frac{n_{HA}}{n_{HM}}\\!=\\!1.19\pm 0.03}$ (see Appendix C for
further details). Hence, at least in terms of vertex divergences, it might be
tempting to “approximate” the PI phase of the HM in a similar spirit as in
Ref. [44], as HA with “reduced” effective interaction $U_{eff}$:

$\mathcal{H}_{PI}=U_{eff}\sum_{i}n_{i\uparrow}n_{i\downarrow}\text{\ \ with\ \
}U_{eff}=\frac{U}{\eta}.$ (7)

Obviously, this analogy cannot be pushed too far: In the limit of
${U\\!\rightarrow\\!\infty}$, the HM-$\bm{\Gamma}_{c}^{\infty}$-lines will
gradually approach the HA-$\bm{\Gamma}_{c}^{\infty}$-lines asymptotically,
i.e., ${U_{eff}\\!\rightarrow\\!U}$ and ${\eta\\!\rightarrow\\!1}$ for
$U\\!\rightarrow\\!\infty$. In practice, ${\eta\approx\\!1.19}$ can be used as
a good approximation for the vertex divergences of the PI phase of the HM, in
the (physically relevant) regime of the MIT [46].

###  III.3  Vertex divergence lines at the phase transition

After separately analyzing how vertex divergences occur in the two (PM and PI)
sets of DMFT solutions in the coexistence region, the question of their
behavior across the thermodynamic transition can be now readily addressed. We
briefly recall that the actual thermodynamic phase transition takes place at
$U\\!=\\!U_{c}(T)$ (cf. Ref. [43]), when the free energy of the PM and the PI
converged DMFT-solutions become equal. The corresponding transition line [43]
marked as white line in all three panels of Fig. 5, thus separates the
thermodynamically stable PI and PM solutions (labeled as physical in Fig. 5)
on the two sides of the coexistence region and is associated, except at its
second-order (quantum) critical endpoints, to an abrupt (first-order) jump of
the physical properties from a PM to a PI behavior, except at its second-order
(quantum) critical endpoints.

After these premises, we summarize in Fig. 5 our results for the vertex
divergences in the proximity of the Mott MIT. In the first two panels, the
respective values $N_{\lambda<0}$ for the PM (left panel) and the PI (middle
panel) solutions are shown as color-intensity plots, making the qualitatively
different shape of the $\bm{\Gamma}_{c}^{\infty}$-lines between the PM and the
PI quite evident. From the direct comparison of the first two panels, the
quantitative difference in $N_{\lambda<0}$ becomes also clearly visible, as
the values of $N_{\lambda<0}$ appear to be systematically higher in the PI
than in the PM realization of the coexistence region.

Figure 6: Left panel: Phase diagram of the inverse temperature
${\beta\\!=\\!(k_{\mathrm{B}}T)^{-1}}$ and interaction $U$ showing the MIT and
the coexistence region with an interpolating color scale of the number of
negative eigenvalues of the generalized charge susceptibility $N_{\lambda<0}$
for the metallic phase, where the thermodynamic transition line $U_{c}$ is
taken from Ref. [43]. Right panel: $N_{\lambda<0}(\beta)$; black markers are
calculated values of $N_{\lambda<0}$ for $U_{c1}$ (triangles), $U_{c}$
(circles) and $U_{c2}$ (squares) according to the left panel. The solid lines
correspond to fits for these data points, yielding
${N_{\lambda<0}^{U_{c1}}(\beta)\\!\cong\\!a_{1}\\!+\\!b_{1}\beta^{-1}}$ [47]
for $\beta\rightarrow\infty$ with the red dashed line as asymptotic value (the
red dotted lines are a guide to the eye),
${N_{\lambda<0}^{U_{c}}(\beta)\\!\propto\\!\ln(\beta)}$ and
${N_{\lambda<0}^{U_{c2}}(\beta)\\!\propto\\!\beta}$. Black empty markers are
estimated values from the fits for ${\beta\\!=\\!300}$ and colored crosses are
the results of corresponding DMFT test calculations. The dashed dotted lines
are intermediate lines ($U^{int}_{i}$) interpolating between $U_{c1}$,
$U_{c}$, and $U_{c2}$.

Hence, when directly considering the vertex divergence behavior for the
thermodynamic transition, as we do in the rightmost panel of Fig. 5,
$N_{\lambda<0}$ displays an evident jump $\Delta N_{\lambda<0}$ between the PM
and PI, at the first-order transition $U_{c}(T)$. On a more quantitative
level, we note that $\Delta N_{\lambda<0}$ between the PM and the PI along
$U_{c}$ increases with decreasing temperature. Only at the second order
critical point at ${T\\!=\\!T_{c}\\!\approx\\!\frac{1}{39}}$ we observe a
continuous transition of $N_{\lambda<0}$ with ${\Delta
N_{\lambda<0}\\!=\\!0}$, as expected. In this regard, $\Delta N_{\lambda<0}$
appears to reflect well the behavior of the order parameter characterizing the
transition. As discussed in the next subsection, however, the same analogy
will not directly apply to the other extreme ($T=0$) of the Mott MIT line, as
this point will be characterized by a divergence of $N_{\lambda<0}$.

###  III.4 Accumulation of vertex divergence lines

Our calculations for $N_{\lambda<0}$ in the coexistence region at finite
temperatures naturally raise the question of how $N_{\lambda<0}$ will behave
at temperatures lower than those accessible by our numerical calculations and,
in particular, for ${T\\!\rightarrow\\!0}$. To this aim, we performed
numerical extrapolations for $N_{\lambda<0}(T)$ from our data and tested the
validity of those by comparing them to additional, numerically heavier
calculations, performed at a “testbed” temperature value (e.g., $\beta=300$
for the PM-phase) lower than the temperature-interval used for the
extrapolation.

Figure 7: Left panel: Phase diagram of the MIT and its coexistence region with
an interpolating color scale of the number of negative eigenvalues
$N_{\lambda<0}$ of the generalized charge susceptibility for the insulating
results. The thermodynamic transition line $U_{c}$ is taken from [43]. Right
panel: $N_{\lambda<0}(\beta)$; black markers are determined values of
$N_{\lambda<0}$ for $U_{c1}$ (triangles), $U_{c}$ (circles) and $U_{c2}$
(squares) according to the left panel. The solid lines mark fits of these
data, yielding linear behavior for the three transition lines. Black empty
markers are estimated values from these fits for ${\beta\\!=\\!200}$ and
colored crosses are the results of corresponding test calculations (a close-up
is shown in the inset on the right). The dashed lines are fits including the
test calculations and the dashed dotted ones present the number of crossed
$\bm{\Gamma}_{c}^{\infty}$-lines of the Hubbard atom along $U_{c1}$, $U_{c}$
and $U_{c2}$

Our results for the PM solution are shown in Fig. 6. Specifically, in the left
panel, we show, as a guidance, $\beta$-$U$ paths in the coexistence region,
along which our extrapolation are performed, superimposed to the corresponding
color-intensity plot for (the numerically interpolated) $N_{\lambda<0}$. In
the right panel, then, we report our numerically calculated data
$N_{\lambda<0}$ as function of the inverse temperature $\beta$ (black markers)
together with the corresponding fits along $U_{c1}$, $U_{c}$, $U_{c2}$ (solid
lines) as well as along intermediate (dash dotted) paths in parameter space
(see Appendix D for further details). The range of the $N_{\lambda<0}(\beta)$
values in the coexistence region is indicated as blue-shaded area. As
mentioned before, the reliability of our extrapolation has been tested by
comparing additional calculations at ${\beta\\!=\\!300}$ (crosses) to the
extrapolated values (empty markers), which yielded a good agreement.

By a closer inspection of the results, we note the following differences in
the ${T\\!\rightarrow\\!0}$ behavior of $N_{\lambda<0}$ along the distinct
paths we selected. For instance, starting with the transition line
$U_{c1}(\beta)$ we find
${N_{\lambda<0}\\!\cong\\!a_{1}\\!+\\!b_{1}\beta^{-1}}$ [47] for
${\beta\\!\rightarrow\\!\infty}$ (${a_{1}\\!\approx\\!6.25}$), which would be
compatible with ${N_{\lambda<0}\\!=\\!6}$ at $U_{c1}(T\\!=\\!0)$. Along the
intermediate lines at ${U_{c1}\\!<\\!U\\!<\\!U_{c}}$ we find
${N_{\lambda<0}\\!\cong\\!a_{i}\\!+\\!b_{i}\beta^{-1}}$ with finite
${a_{i}\\!>\\!a_{1}}$. At the transition lines $U_{c}(\beta)$ and
$U_{c2}(\beta)$, which both end at the second-order critical point at
$U_{c2}(T=0)$, $N_{\lambda<0}(\beta)$ is observed to grow logarithmically and
linearly in $\beta$, respectively. Consistently, the intermediate lines in the
interval ${U_{c}\\!<\\!U\\!<\\!U_{c2}}$ feature a divergent behavior of
${N_{\lambda<0}(\beta)\\!\propto\\!\beta^{x}}$ with ${0\\!<\\!x\\!<\\!1}$.
This indicates that an infinite number of $\bm{\Gamma}_{c}^{\infty}$-lines
need to be crossed before reaching $U_{c2}(T\\!=\\!0)$, i.e., the location of
the Mott-Hubbard MIT in the ground-state.

Our extrapolations for $N_{\lambda<0}(\beta)$ of the PI solution are displayed
Fig. 7. Again, in the left panel of the figure the respective location of the
paths considered in the $\beta$-$U$ coexistence region are shown, superimposed
to the corresponding values of $N_{\lambda<0}(\beta)$ for the PI solutions.
Due to the numerical hurdles of performing two-particle DMFT calculations for
low $T$ in the PI-phase (see Appendix A.1), the extrapolations (solid lines)
and the test calculations (empty markers) have been performed at higher
temperatures than for the PM phase (PI test calculations at
${\beta\\!=\\!200}$), which likely explains the overall less satisfactory
agreement between extrapolation and test calculations w.r.t. the PM case. The
inset in the right panel presents a close-up of the corresponding deviations.
Fits including the $\beta=200$ calculations, which should be regarded as the
most precise results which we could obtain in the PI phase, are displayed as
dashed lines. As comparison, dashed dotted lines present
$N_{\lambda<0}(\beta)$ of the HA along $U_{c1}$, $U_{c}$, and $U_{c2}$, all
showing higher values than the PI solution of the HM in the coexistence
region.

The extrapolations for $N_{\lambda<0}(\beta)$ in the PI (with and without test
calculations) display a linear behavior in $\beta$ for all three transition
lines $U_{c1}$, $U_{c}$, $U_{c2}$, corresponding to an infinite number of
$\bm{\Gamma}_{c}^{\infty}$-lines to be crossed before reaching ${T\\!=\\!0}$.
On a more quantitative level, we note that a comparison between the slopes
$\alpha_{i}$ of the extrapolations of $N_{\lambda<0}(\beta)$ and the
corresponding results of the HA, can be also used to obtain
${\eta\\!=\\!\frac{\alpha_{HA}}{\alpha_{HM}}}$ for an effective HA like
description of the PI solution of the HM in Eq. (7) of Sec. III B, e.g.
${\eta\\!=\\!1.22219}$ for $U_{c1}(T)$ (for details see Appendix C).

From our extrapolations of $N_{\lambda<0}(\beta)$ we can now try to estimate
the overall behavior of $N_{\lambda<0}$ at ${T\\!=\\!0}$ across the MIT. By
crossing the MIT from ${\text{PM}\\!\rightarrow\\!\text{PI}}$ at
${T\\!=\\!0}$, the number of $\bm{\Gamma}_{c}^{\infty}$-lines increases to
infinity at the critical endpoint $U_{c2}(T\\!=\\!0)$. In this respect, at
${T\\!>\\!0}$ the physical transition line $U_{c}(T)$ can bee seen as the
first of the intermediate line-paths in $\beta$,$U$ of Fig. 6 along which
$N_{\lambda<0}(\beta)$ diverges for ${\beta\\!\rightarrow\\!\infty}$, namely
logarithmically. In the metastable metallic phase for
${U_{c}\\!<\\!U\\!<\\!U_{c2}}$ – along the intermediate lines –
$N_{\lambda<0}(\beta)$ can diverge faster than $\log{\beta}$, and finally
${N_{\lambda<0}(\beta)\\!\propto\\!\beta}$ at $U_{c2}$. The same
$\beta$-dependence, namely ${N_{\lambda<0}(\beta)\\!\propto\\!\beta}$, occurs
in the PI coexistence region as well as in the HA.

On the basis of this numerical evidence, we can then conclude that at the
Mott-Hubbard MIT at ${T\\!=\\!0}$ an accumulation point of an infinite number
of vertex divergence lines occurs, a clear link between the Mott MIT at
${T\\!=\\!0}$ and the vertex divergence lines. Clearly, the different paths in
the $\\{U,\beta\\}$ parameter plane, through which the $T\\!=\\!0$ transition
can be approached, e.g. $U_{c}(\beta)$ or $U_{c2}(\beta)$, reflect in
correspondingly different ways in which the divergence of
$N_{\lambda<0}(\beta\rightarrow\infty)$ is eventually achieved 111Evidently,
this also means that the difference of $N_{\lambda<0}(\beta)$ between
different paths approaching the $T\\!=\\!0$ MIT, e.g. $U_{c}(\beta)$ or
$U_{c2}(\beta)$, does not necessarily vanish (and it might even diverge) for
$\beta\rightarrow\infty$, depending on the pair of paths chosen.. However, as
the asymptotic value of $N_{\lambda<0}(\beta\rightarrow\infty)$ always
diverges when approaching $U_{c}(T\\!=\\!0)\\!=\\!U_{c2}(T\\!=\\!0)$
independently of the path, our extrapolated results for
$N_{\lambda<0}(\beta\rightarrow\infty)$ appear consistent 222In particular,
the $T\\!=\\!0$-extrapolated behavior of $N_{\lambda<0}$ could be consistent
with the existence of a unique function describing the increase of
$N_{\lambda<0}$ for increasing $U$ along the $T\\!=\\!0$-axis up to its
divergence at ${U=U_{c2}(T\\!=\\!0)}$, and, hence, does not rule out the
possibility of a smooth evolution of the generalized susceptibility
$\bm{\chi}^{\nu\nu^{\prime}}_{c}$ and of the corresponding physical response
function. with the smooth evolution of several physical quantities (such as
the quasiparticle spectral weight, the double-occupancy, etc.) numerically
reported [50, 51] at the $T=0$-Mott MIT.

##  IV Theoretical and Algorithmic implications

Figure 8:  Phase diagrams of the Hubbard atom (HA) (reconstructed from Ref.
[9]), Hubbard model (HM) and Anderson impurity model (AIM) (reconstructed from
Ref. [13]) with the corresponding $\bm{\Gamma}_{c}^{\infty}$-lines and their
accumulation point (blue dot) at ${T\\!=\\!0}$ (${U_{MIT}\\!\\!0}$ for the HA,
${U_{MIT}\\!=\\!U_{c2}^{T=0}}$ for the HM and
${U_{MIT}\\!\rightarrow\\!\infty}$ for the AIM). The black arrow below
sketches the shift of the accumulation point of the
$\bm{\Gamma}_{c}^{\infty}$-lines between the different models. The
accumulation point at the MIT of the DMFT solution of the Hubbard model at
$U_{c2}(T\\!=\\!0)$ can be seen as intermediate case w.r.t. the “extreme”
cases of the purely insulating Hubbard atom and the purely metallic Anderson
impurity model.

Our analysis of the irreducible vertex divergences in the close proximity of
the Mott MIT has several important implications both of algorithmic as well as
of conceptual nature. We demonstrated the existence of a large number of
divergence lines found in the coexistence region of the MIT, increasing up to
infinity along all calculated parameter paths towards the ${T\\!=\\!0}$ MIT at
${U_{c2}(T\\!=\\!0)}$ evidently poses a huge challenge to the applicability of
all algorithmic approaches using and/or explicitly manipulating irreducible
vertices of DMFT as input, such as the full-fledged, parquet-equation based
version of the dynamical vertex approximation [52, 53, 54, 55, 56] and the
QUADRILEX scheme [57]. In fact, while one could realistically cope with a
situation of a few well-separated divergence lines, it becomes hard to
stabilize the numerical manipulation of irreducible vertices in parameter
regimes, where their frequency structure will display large oscillations, due
to the proximity of several and closely spaced divergence lines.

At the same time, we recall that the occurrence of vertex divergences has been
recently proven to be pivotal [21, 27] to ensure the correct transfer of
physical information among the spin/magnetic and the charge/density or the
pairing channels in the nonperturbative regimes, e.g., to feature the proper
freezing of the local charge (or pairing) fluctuations associated to the pre-
formation of local magnetic moments. Hence, it becomes clear why self-
consistent diagrammatic approaches, whose irreducible vertex functions cannot
diverge per construction [58], such as the truncated functional
renormalization group (fRG) [59] or the parquet approximation [60], will not
be able to yield a consistent picture of the Mott MIT and its related
phenomena.

While it is beyond the scope of this work to address the issue of possible
strategies how to include this relevant piece of physical information in
diagrammatic treatments of the many-electron problems in nonperturbative
regimes, e.g, by means of the merger of DMFT and fRG [61, 62, 63, 64, 65, 66]
or by rewriting ladder diagrammatic resummations or parquet equations in terms
of the full two-particle scattering amplitude of the local problem [67, 31,
54, 68], here we want to discuss some more fundamental theoretical implication
of our results.

Indeed, the identification of the Mott-Hubbard MIT at zero temperature with
the accumulation point of an infinite number of
$\bm{\Gamma}_{c}^{\infty}$-lines does not only clarify the precise relation
linking the occurrence of vertex divergences in the charge channel and the
MIT, highlighting the profoundly nonperturabtive nature of the MIT itself, but
also allows to draw interesting, more general considerations about how the
vertex divergences affect different many-electron problems. In particular, our
results for the HM in DMFT (summarized in the central panel of Fig. 8) can be
put into a broader perspective by directly comparing them to the corresponding
ones for the Hubbard atom [9, 14] (HA, left panel) and the (metallic) Anderson
impurity model of Ref. [13] (AIM, right panel) in Fig. 8, where we used the
same color-coding for denoting the different vertex divergence lines.

Thereby, we note that for the HA, i.e. for an Hubbard model with
${t\\!\equiv\\!0}$, which features a perfect local-moment/insulating behavior
for all ${U\\!>\\!U_{MIT}\\!=\\!0}$, an accumulation point of vertex
divergence lines occurs exactly at the origin of the phase diagram, namely at
${U\\!=\\!T\\!=\\!0}$. On the other hand, by looking at the phase-diagram of
the AIM, where no transition to an insulating ground state occurs for any
value of $U$, we do not observe an accumulation point of vertex divergence
lines in the available data [13] and, in general, we do not expect one at
${T\\!=\\!0}$. We can note, nonetheless, that their low-$T$ distribution
becomes denser by increasing interaction values. One can then assume the
asymptotic value ${U_{MIT}\\!=\\!+\\!\infty}$ as an hypothetical value for the
transition to an insulating state in the AIM. In this regime, the physics
would be eventually dominated by $U$, similarly to the insulating phases of
the other two cases, and the progressively denser distribution of vertex
divergence lines by increasing $U$ could then lead asymptotically to an
accumulation point.

Within this framework [69], the results we obtained for the HM in DMFT can be
naturally interpreted as an “intermediate” case between the two extreme
situations of the HA (with accumulation point at ${U_{MIT}\\!=\\!0}$) and the
AIM (at ${U_{MIT}\\!=\\!+\\!\infty}$), where $U_{MIT}$ yields the finite value
of $U_{c2}(T\\!=\\!0)$, as schematically illustrated by the black arrow
sketched below the three panels of Fig. 8. Heuristically, one could imagine to
start from the AIM case (where ${U_{MIT}\\!=\\!+\\!\infty}$). Then, by
progressively moving the value of $U_{MIT}$ first towards lower finite values
and, subsequently, down to $0$, one could qualitatively “reconstruct” the
other two cases by squeezing the corresponding
$\bm{\Gamma}_{c}^{\infty}$-lines against the corresponding accumulation
points.

On a more formal level, we note that the presence or the absence of an
accumulation point of vertex divergence lines, as those discussed here, might
play an important role in the way vertex divergences may affect the
(${T\\!=\\!0}$)-physical properties of strongly correlated electron systems,
such as, e.g. the validity (or the violation) of the Luttinger theorem [70].
Further, the presence of a increasing number of vertex divergences when
approaching the Mott MIT at $T\\!=\\!0$ from the correlated PM-side appears to
provide a key to reconcile the analytical derivation of [71] with the
numerical evidence [50] of the occurrence of a smooth transition in DMFT
calculations. Indeed, in [71] the impossibility of observing such a smooth
Mott MIT at zero temperature was demonstrated by assuming the convergence of
the (self-consistent) skeleton expansion for all $U\\!<\\!U_{c}(T\\!=\\!0)$,
an assumption which our study (see also the Appendix of [9]) demonstrated to
break down for $U$ much smaller than $U_{c}(T\\!=\\!0)$.

##  V Conclusion and Outlook

In this work, we have investigated the relation between a characteristic
aspect of the breakdown of the self-consistent perturbation expansion [10],
namely the occurrence of divergences of irreducible vertex functions in the
charge channel [1, 9] and the Mott-Hubbard metal-to-insulator transition. To
this aim, by performing DMFT calculations on the two-particle level for the
half-filled Hubbard model on a Bethe lattice, we have systematically studied
the occurrence of irreducible vertex divergences in the coexistence region
adjacent to the MIT. Our results demonstrate how the shape and the number of
irreducible vertex divergence lines ($\bm{\Gamma}_{c}^{\infty}$-lines) in the
coexistence region is significantly different in the PM and the PI solutions,
with an abrupt jump across the thermodynamic first-order transition line,
reflecting the different degree of suppression of on-site charge fluctuations
[8, 10, 21, 27] in the two phases. Further, we could show that, in spite the
evident backward bending displayed by the $\bm{\Gamma}_{c}^{\infty}$-lines in
the correlated PM phase, the number of divergences crossed by approaching the
${T\\!=\\!0}$ MIT at ${U\\!=\\!U_{c2}(T\\!=\\!0)}$ is diverging, making the
location of the MIT itself an accumulation point for
$\bm{\Gamma}_{c}^{\infty}$-lines. This finding establishes a clear connection
between the irreducible-vertex divergences and the occurrence of the MIT,
clarifying the difference with cases where vertex divergences appear in
systems where no MIT takes place, and substantiating the interpretation of the
Mott transition as highly-nonperturbative phenomenon in a more precise
context.

Beyond the possible algorithmic implications of our results, especially
relevant for Feynman diagrammatic approaches to treat intermediate-to-strong
coupling regimes [54, 72], it would be interesting, in the future, to extend
our study by including the effects on non-local correlations, e.g. by means of
cluster extensions [73] of DMFT. In particular, one could investigate, whether
the shift towards lower $U$ values of the Mott MIT, found [73, 74] by
including spatial correlations of progressively larger size in the two-
dimensional Hubbard model on a square lattice as well as the change of slope
in the associated transition line $U_{c}(T)$, would be accompanied by a
corresponding shift of the accumulation point of the
$\bm{\Gamma}_{c}^{\infty}$-lines as well as by a corresponding change of their
bending. This piece of information would be of particular importance for
precisely characterizing the nonperturbative nature of electronic correlations
in two-dimensional systems, for which the location of the MIT, in absence of
geometrical frustration, may be shifted [75] down to ${U\\!=\\!0^{+}}$ in the
low-temperature limit.

###### Acknowledgements.

We thank S. Andergassen, S. Ciuchi, P. Chalupa-Gantner, H. Eßl, D. Fus, G.
Rohringer, G. Sangiovanni and T. Schäfer for very insightful discussions and
P. Chalupa-Gantner also for carefully reading the manuscript. We acknowledge
financial support from the Austrian Science Fund (FWF), within the project
I-5487 (M.R. and A.T.) as well as project I 5868 (Project P01, part of the FOR
5249 [QUAST] of the German Science Foundation, DFG) (S.A.). Matthias Reitner
acknowledges support as a recipient of a DOC fellowship of the Austrian
Academy of Sciences. Calculations have been performed on the Vienna Scientific
Cluster (VSC-4).

##  Appendix A DMFT-calculations of the two-particle Green’s function
$G^{(2)}$

As mentioned in Sec. II D, we use a continuous-time quantum Monte Carlo solver
(CT-QMC) in the hybridization expansion (CT-HYB) [40] of the w2dynamics
package [41, 42] for the DMFT calculations of the two-particle Green’s
function $G^{(2)}$ of the auxiliary impurity model and performed our
calculations on the Vienna scientific cluster (VSC-4). A data set containing
all numerical data and plot scripts used to generate the figures of this
publication is publicly available on the TU Wien Research Data repository
[76].

The detailed calculation process for the two-particle Green’s function for a
specific phase point in $U$, $T$ is the following:

  1. 1.

Calculations of the one-particle Green’s function $G^{(1)}$ with a comparable
low number of measurements (e.g. ${\sim\\!10^{5}}$) at each DMFT iteration, to
converge the DMFT self-consistence cycle. The convergence is verified by
tracking the behavior of $G^{(1)}$ for the first Matsubara frequency and the
double occupancy.

  2. 2.

The first computations are followed by a calculation with a high number of
measurements (e.g. ${\sim\\!10^{7}}$) and a small number of DMFT iterations to
obtain $G^{(1)}$ as function of Matsubara frequencies with a satisfactory
error to noise resolution as starting point for a two-particle calculation.

  3. 3.

One DMFT iteration including a calculation of $G^{(2)}$ with a high number of
measurements (e.g. ${\sim\\!10^{7}}$) using a frequency box size of
${100\\!\times\\!100}$ fermionic frequencies for high temperatures and a box
size up to ${200\\!\times\\!200}$ frequencies for our lowest temperatures.
Thereby, we ensure that the number of negative eigenvalues, originating from
lower frequency structure of $G^{(2)}$ (see e.g. [77, 31, 78, 79]), is not
effected by a further increase of the frequency box size. Each individual
computation required an average of $\sim$10000 CPU hours.

###  A.1 Sampling methods

Within the CT-HYB calculations of the w2dynamics package, different sampling
methods for the computation of one- and two-particle Green’s function can be
used. For most of the parameter regimes we have used partition function
sampling with the superstate-sampling method [42] (segment sampling [80]
yielded no computational advantage).

Figure 9: Left column: real part of the generalized charge susceptibility
$\bm{\chi}^{\nu\nu^{\prime}}_{c}$ for ${\beta\\!=\\!60}$ and ${U\\!=\\!2.46}$
for a PI solution. Right column: corresponding eigenvalues
$\lambda_{\chi_{c}}$ around zero. Upper row is calculated with state sampling,
lower row with worm sampling.

However, for the insulating solution in the coexistence region of our model
and especially at low temperatures, partition function sampling runs into
problems. This manifests in numerical artifacts, which occur as diagonal
stripes in the $\bm{\chi}_{c}$-matrix causing complex pairs of the eigenvalues
(see Fig. 9), which are however prohibited by the particle-hole and
$\operatorname{SU}(2)$-symmetry of the model considered ($\bm{\chi}_{c}$ must
be a real bisymmetric matrix with only real eigenvalues [17]), and therefore
result in unusable data. This is presumably caused by the suppressed
hybridization function $\Delta(\tau)$ in the insulating phase, resulting in a
Monte Carlo estimator with high variance, stemming from functional derivatives
of the partition function w.r.t. the hybridization function [81]. To overcome
these numerical difficulties we used the worm-sampling method of w2dynamics
[41], which uses Monte Carlo sampling in both the partition function and the
Green’s function space. Since, the worm-sampling method is numerically more
expensive than superstate-sampling, we used worm-sampling only for lower
temperatures in the insulating phase (${\beta\\!>\\!60}$).

##  Appendix B Approximation of vertex divergence lines

Our numerical two-particle calculations for the determination of
$N_{\lambda<0}$ have been performed for a finite set of $T$, $U$ parameters.
To extract the points in the phase space, where the vertex divergence occurs,
i.e. ${\lambda\\!=\\!0}$, from our data, we compared different approaches
which are detailed in the following subsections. For our results in Sec. III
we used the polynomial fits described in Appendix B.2. For both methods the
symmetry of the eigenvectors was not taken into account. Hence, a possible
crossing of divergence lines [13] was not investigated. Instead, the line
ordering found in Refs. [9, 17] is assumed.

###  B.1 Approximation via eigenvalues of $\boldmath{\chi}_{c}$

Figure 10: Sketch of the extrapolation scheme to obtain the interaction value
${U\\!=\\!U_{\bm{\Gamma}_{c}^{\infty}}}$ of the zero-crossing of the
eigenvalue $\lambda_{l}$ of $\chi_{c}^{\nu\nu^{\prime}}$ from two-particle
DMFT calculations at different $U_{m}$ for a fixed temperature $T$.

For fixed temperature $T$ scans with two-particle calculations close to the
divergence, one can use a straight forward approach to extract the approximate
interaction value $U_{\bm{\Gamma}_{c}^{\infty}}$, where ${\lambda\\!=\\!0}$
and a vertex divergence occurs: Linearly extrapolating from the first two data
points, where ${\lambda\\!<\\!0}$ is already negative.
$U_{\bm{\Gamma}_{c}^{\infty}}$ is obtained, where the extrapolation equals
zero. This is schematically presented in Fig. 10. The corresponding
$\bm{\Gamma}_{c}^{\infty}$-line of ${\lambda_{l}\\!=\\!0}$ for $U,T$ is
approximated by connecting the obtained results for different $T$. The results
for the $\bm{\Gamma}_{c}^{\infty}$-lines are presented in Fig. 12, where the
left panel corresponds to the extrapolation approach.

###  B.2 Approximation via number $N_{\lambda<0}$ of negative eigenvalues of
$\boldmath{\chi}_{c}$

Figure 11: Polynomial fits of $N_{\lambda<0}(U)$ for several temperatures.
Black markers indicate $U_{c1}(T)$ (triangles), $U_{c}(T)$ from [43]
(circles), and $U_{c2}(T)$ (squares) for the temperatures of the fits. The
inset shows a zoom for low $\beta$ and $U$. Figure 12: Phase diagrams with
approximated $\bm{\Gamma}_{c}^{\infty}$-lines for the PM solution in the
coexistence region, according to the approximations of B.1 and B.2 for several
fixed $T$: (left) extrapolation of the eigenvalues of
$\chi_{c}^{\nu\nu^{\prime}}$ and (right) polynomial fits of $N_{\lambda<0}$
along $U$. While the precise position of the $\bm{\Gamma}_{c}^{\infty}$-lines
in $T,U$ differs between the different approximation schemes, the qualitative
behavior of $N_{\lambda<0}(U,T)$ does not dependent on them.

Alternatively to the previously described extrapolation scheme to approximate
the $\bm{\Gamma}_{c}^{\infty}$-lines, $N_{\lambda<0}$ may be fitted by
polynomial functions $p(U)$ for fixed $T$ by using the least square method.
This approach, using the total number of negative eigenvalues, will in general
yield less accurate positions of the divergence lines as compared to the
approach of Appendix B.1. However, this approach takes the non-linear behavior
of $N_{\lambda<0}$ along the temperature scans into account. The fits of the
metallic coexistence region are shown in Fig. 11. For those fits we used
$p(U)$ with different polynomial degrees. The polynomial degree thereby
increases with decreasing $T$ (increasing $\beta$), due to the rapid non-
linear increase of $N_{\lambda<0}$: ${p(U)\\!\propto\\!U}$ for
${\beta\\!\in\\!\\{40,44,50,60,75\\}}$, ${p(U)\\!\propto\\!U^{2}}$ for
${\beta\\!\in\\!\\{100,133\\}}$ and ${p(U)\\!\propto\\!U^{3}}$ for
${\beta=200}$. The $\bm{\Gamma}_{c}^{\infty}$-lines are approximated by
connecting the points in $T$, $U$ with the same value of $N_{\lambda<0}$. The
results for these lines are shown in the right panel of Fig. 12.

Let us stress that, although, the precise position of the
$\bm{\Gamma}_{c}^{\infty}$-lines in $T,U$ differs between the different
approximation schemes, the overall behavior of $N_{\lambda<0}(U,T)$ and its
asymptotic for ${T\\!\to\\!0}$ does not dependent on them. Eventually, we note
that, to generate $N_{\lambda<0}$ as function of $U$, $T$ for the parameter
space, which is shown as color scale plots in Sec. III, an additional linear
interpolation between different temperatures $T$ has been used.

##  Appendix C Effective Hubbard atom like description of the insulating
Hubbard model

To approximately describe the behavior the divergence lines of the Mott
insulating phase of the Hubbard model with Eq. (7) and to account for the
correct limit of ${U_{eff}\\!\to\\!\infty}$, we can construct a simple
function

$\eta(U)=\frac{U_{c1}^{T=0}}{U}(\eta_{0}-1)+1\ \text{,}$ (8)

where ${\eta_{0}\\!=\\!1.19\\!\pm\\!0.03\\!\approx\\!U^{T=0}_{c1}/2}$ is the
mean shift between $\bm{\Gamma}_{c}^{\infty}$-lines of the HM and the HA
within the coexistence region. In Eq. (8) we interpolate between $\eta_{0}$ at
the beginning of the Mott insulating phase and ${\eta\\!\rightarrow\\!1}$ at
${U\\!\rightarrow\\!\infty}$ to account for the asymptotic behavior of the
$\bm{\Gamma}_{c}^{\infty}$-lines. Hence, for $U_{eff}$ in $\mathcal{H}_{PI}$
of Eq. (7) we get

$\displaystyle U_{eff}$
$\displaystyle=\frac{U}{\eta(U)}=\frac{U^{2}}{U+U_{c1}^{T=0}(\eta_{0}-1)}=$
$\displaystyle\overset{2\eta_{0}\sim
U_{c1}^{T=0}}{=}\frac{U^{2}}{U+{U_{c1}^{T=0}}(\frac{1}{2}U_{c1}^{T=0}-1)}$ (9)

The results for the temperature behavior of $N_{\lambda<0}$ of the PI solution
in Sec. III.4 suggest to introduce another way to evaluate $\eta_{0}$, by
comparing the slopes $\alpha_{i}$ of $N_{\lambda<0}(\beta)$ in Fig. 7 of Sec.
III.2 between the HA and the HM at the same $U$. We can estimate the parameter
$\eta_{0}$ by ${\eta_{0}^{\text{slope}}\\!=\\!\alpha_{HA}/\alpha_{HM}}$. The
resulting values are listed in Tab. 1.

|  $U_{c1}(T)$ |  $U_{c}(T)$ |  $U_{c2}(T)$
---|---|---|---
$\eta_{0}^{\text{slopes}}$ | 1.22 | 1.16 | 1.15
Table 1: Parameter $\eta_{0}$ calculated from the quotient of the slopes of
$N_{\lambda<0}$ along the transition lines $U_{c1}(T)$, $U_{c}(T)$ and
$U_{c2}(T)$ in Fig. 7 for the effective HA description of the HM in the
insulating phase.

In Tab. 2 we compare ${n_{HA}\\!=\\!N_{\lambda<0}}$, calculated from
$U_{eff}(\eta_{0})$ for two different $\eta_{0}$, with the exact
${n_{HM}\\!=\\!N_{\lambda<0}}$ of the HM for two $U$ values at
${\beta\\!=\\!40}$: ${\eta_{0}^{\text{slope}}\\!=\\!1.22}$ at $U_{c1}(T)$ (the
first $U$ value where a Mott insulating phase is possible) and
${\eta_{0}\\!=\\!1.19}$, from our considerations in Sec. III.2. We see that
$n_{HA}$ for $\eta_{0}^{slope}$ provides us a particularly good approximation
for the Mott insulating phase of the Hubbard model.

$\beta=40$ | $n_{HM}$ |  $n_{HA}$ at $U_{eff}$
---|---|---
$\eta_{0}$ | $\eta_{0}^{\text{slope}}$
$U=3.1$ | 31 |  30 | 30
$U=3.5$ | 34 | 36 | 34
Table 2: Test calculation results for the number of crossed vertex divergence
lines and corresponding values according to the Hubbard atom with effective
interaction $U_{eff}$ for ${\eta_{0}\\!=\\!1.19}$ and
${\eta_{0}^{\text{slope}}\\!=\\!1.22}$.

##  Appendix D Calculation of intermediate lines

The dashed dotted lines in the left panel of Fig. 6 are intermediate lines
interpolating between the $U_{c1}(T)$, $U_{c}(T)$, and $U_{c2}(T)$ lines
according to

$\displaystyle U^{int}_{1}(T)=\frac{1}{2}\,\bigl{[}U_{c1}(T)+U_{c}(T)\bigr{]}$
(10) $\displaystyle
U^{int}_{2}(T)=\frac{1}{6}\,\bigl{[}U_{c1}(T)+5\,U_{c}(T)\bigr{]}$ (11)
$\displaystyle
U^{int}_{3}(T)=\frac{1}{4}\,\bigl{[}3\,U_{c}(T)+U_{c2}(T)\bigr{]}$ (12)
$\displaystyle U^{int}_{4}(T)=\frac{1}{2}\,\bigl{[}U_{c}(T)+U_{c2}(T)\bigr{]}$
(13) $\displaystyle
U^{int}_{5}(T)=\frac{1}{4}\,\bigl{[}U_{c}(T)+3\,U_{c2}(T)\bigr{]}.$ (14)

We used the polynomial fits (Fig. 11) to evaluate the corresponding
$N_{\lambda<0}$ along those lines at inverse temperatures
${\beta\\!=\\!\\{40,44,50,60,75,100,133,200\\}}$ and applied the same fitting
routine as along $U_{c1}(T)$, $U_{c}(T)$, and $U_{c2}(T)$ to generate the
corresponding dashed dotted lines in the right panel.

## References

  * Schäfer  _et al._ [2013] T. Schäfer, G. Rohringer, O. Gunnarsson, S. Ciuchi, G. Sangiovanni, and A. Toschi, Divergent precursors of the mott-hubbard transition at the two-particle level, Phys. Rev. Lett. 110, 246405 (2013).
  * Janiš and Pokorný [2014] V. Janiš and V. Pokorný, Critical metal-insulator transition and divergence in a two-particle irreducible vertex in disordered and interacting electron systems, Phys. Rev. B 90, 045143 (2014).
  * Kozik  _et al._ [2015] E. Kozik, M. Ferrero, and A. Georges, Nonexistence of the luttinger-ward functional and misleading convergence of skeleton diagrammatic series for hubbard-like models, Phys. Rev. Lett. 114, 156402 (2015).
  * Stan  _et al._ [2015] A. Stan, P. Romaniello, S. Rigamonti, L. Reining, and J. A. Berger, Unphysical and physical solutions in many-body theories: from weak to strong correlation, New Journal of Physics 17, 093045 (2015).
  * Rossi and Werner [2015] R. Rossi and F. Werner, Skeleton series and multivaluedness of the self-energy functional in zero space-time dimensions, Journal of Physics A: Mathematical and Theoretical 48, 485202 (2015).
  * Ribic  _et al._ [2016] T. Ribic, G. Rohringer, and K. Held, Nonlocal correlations and spectral properties of the falicov-kimball model, Phys. Rev. B 93, 195105 (2016).
  * Rossi  _et al._ [2016] R. Rossi, F. Werner, N. Prokof’ev, and B. Svistunov, Shifted-action expansion and applicability of dressed diagrammatic schemes, Phys. Rev. B 93, 161102(R) (2016).
  * Gunnarsson  _et al._ [2016] O. Gunnarsson, T. Schäfer, J. P. F. LeBlanc, J. Merino, G. Sangiovanni, G. Rohringer, and A. Toschi, Parquet decomposition calculations of the electronic self-energy, Phys. Rev. B 93, 245102 (2016).
  * Schäfer  _et al._ [2016] T. Schäfer, S. Ciuchi, M. Wallerberger, P. Thunström, O. Gunnarsson, G. Sangiovanni, G. Rohringer, and A. Toschi, Nonperturbative landscape of the mott-hubbard transition: Multiple divergence lines around the critical endpoint, Phys. Rev. B 94, 235108 (2016).
  * Gunnarsson  _et al._ [2017] O. Gunnarsson, G. Rohringer, T. Schäfer, G. Sangiovanni, and A. Toschi, Breakdown of traditional many-body theories for correlated electrons, Phys. Rev. Lett. 119, 056402 (2017).
  * Tarantino  _et al._ [2018] W. Tarantino, B. S. Mendoza, P. Romaniello, J. A. Berger, and L. Reining, Many-body perturbation theory and non-perturbative approaches: screened interaction as the key ingredient, Journal of Physics: Condensed Matter 30, 135602 (2018).
  * Vučičević  _et al._ [2018] J. Vučičević, N. Wentzell, M. Ferrero, and O. Parcollet, Practical consequences of the luttinger-ward functional multivaluedness for cluster dmft methods, Phys. Rev. B 97, 125141 (2018).
  * Chalupa  _et al._ [2018] P. Chalupa, P. Gunacker, T. Schäfer, K. Held, and A. Toschi, Divergences of the irreducible vertex functions in correlated metallic systems: Insights from the anderson impurity model, Phys. Rev. B 97, 245136 (2018).
  * Thunström  _et al._ [2018] P. Thunström, O. Gunnarsson, S. Ciuchi, and G. Rohringer, Analytical investigation of singularities in two-particle irreducible vertex functions of the hubbard atom, Phys. Rev. B 98, 235107 (2018).
  * Nourafkan  _et al._ [2019] R. Nourafkan, M. Côté, and A.-M. S. Tremblay, Charge fluctuations in lightly hole-doped cuprates: Effect of vertex corrections, Phys. Rev. B 99, 035161 (2019).
  * Melnick and Kotliar [2020] C. Melnick and G. Kotliar, Fermi-liquid theory and divergences of the two-particle irreducible vertex in the periodic anderson lattice, Phys. Rev. B 101, 165105 (2020).
  * Springer  _et al._ [2020] D. Springer, P. Chalupa, S. Ciuchi, G. Sangiovanni, and A. Toschi, Interplay between local response and vertex divergences in many-fermion systems with on-site attraction, Phys. Rev. B 101, 155148 (2020).
  * Kim and Sacksteder [2020] A. J. Kim and V. Sacksteder, Multivaluedness of the luttinger-ward functional in the fermionic and bosonic system with replicas, Phys. Rev. B 101, 115146 (2020).
  * van Loon  _et al._ [2020] E. G. C. P. van Loon, F. Krien, and A. A. Katanin, Bethe-salpeter equation at the critical end point of the mott transition, Phys. Rev. Lett. 125, 136402 (2020).
  * Reitner  _et al._ [2020] M. Reitner, P. Chalupa, L. Del Re, D. Springer, S. Ciuchi, G. Sangiovanni, and A. Toschi, Attractive effect of a strong electronic repulsion: The physics of vertex divergences, Phys. Rev. Lett. 125, 196403 (2020).
  * Chalupa  _et al._ [2021] P. Chalupa, T. Schäfer, M. Reitner, D. Springer, S. Andergassen, and A. Toschi, Fingerprints of the local moment formation and its kondo screening in the generalized susceptibilities of many-electron problems, Phys. Rev. Lett. 126, 056403 (2021).
  * [22] K. Van Houcke, E. Kozik, R. Rossi, Y. Deng, and F. Werner, Physical and unphysical regimes of self-consistent many-body perturbation theory, arXiv:2102.04508 [cond-mat.str-el] .
  * Mazitov and Katanin [2022a] T. B. Mazitov and A. A. Katanin, Local magnetic moment formation and kondo screening in the half-filled single-band hubbard model, Phys. Rev. B 105, L081111 (2022a).
  * Stepanov  _et al._ [2022] E. A. Stepanov, S. Brener, V. Harkov, M. I. Katsnelson, and A. I. Lichtenstein, Spin dynamics of itinerant electrons: Local magnetic moment formation and berry phase, Phys. Rev. B 105, 155151 (2022).
  * Mazitov and Katanin [2022b] T. B. Mazitov and A. A. Katanin, Effect of local magnetic moments on spectral properties and resistivity near interaction- and doping-induced mott transitions, Phys. Rev. B 106, 205148 (2022b).
  * van Loon [2022] E. G. C. P. van Loon, Two-particle correlations and the metal-insulator transition: Iterated perturbation theory revisited, Phys. Rev. B 105, 245104 (2022).
  * [27] S. Adler, F. Krien, P. Chalupa-Gantner, G. Sangiovanni, and A. Toschi, Non-perturbative intertwining between spin and charge correlations: A ”smoking gun” single-boson-exchange result, arXiv:2212.09693 [cond-mat.str-el] .
  * [28] A. J. Kim and E. Kozik, Misleading convergence of the skeleton diagrammatic technique: when the correct solution can be found, arXiv:2212.14768 [cond-mat.str-el] .
  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Rohringer  _et al._ [2012] G. Rohringer, A. Valli, and A. Toschi, Local electronic correlation at the two-particle level, Phys. Rev. B 86, 125114 (2012).
  * [32] J. Hubbard, Electron correlations in narrow energy bands, Proc. R. Soc. Lond. A. Mathematical and Physical Sciences 276, 238.
  * [33] The specific borders ($U_{c1,c2}(T)$) of the coexistence region have been also explicitly recalculated for this work, finding an essential agreement with the existing results [43], while for $U_{c}(T)$, which is determined by comparing the free energies of the PM and PI phases, literature results [43] have been used.
  * [34] We note here, that in the particle-hole symmetric case considered here, the color of the line indicates not only whether the vertex-divergence occurs only in $\bm{\Gamma}_{c}$ (red) or both $\bm{\Gamma}_{c}$ and $\bm{\Gamma}_{pp}$ (orange), but also whether the eigenvector $V_{i}(\nu)$ associated to the vanishing eigenvalue $\lambda_{i}$ of $\bm{\chi}_{c}$ is anti-symmetric ${V_{i}(\nu)\\!=\\!-V_{i}(-\nu)}$ (red) or symmetric ${V_{i}(\nu)\\!=\\!V_{i}(-\nu)}$ (orange) function of frequency [1, 9, 13]. The symmetry of the eigenvectors was not taken into account in this work, for further information see Appendix B.
  * Nakanishi [1969] N. Nakanishi, A General Survey of the Theory of the Bethe-Salpeter Equation, Progress of Theoretical Physics Supplement 43, 1 (1969).
  * Wilcox [1968] R. M. Wilcox, Bounds for the isothermal, adiabatic, and isolated static susceptibility tensors, Phys. Rev. 174, 624 (1968).
  * Watzenböck  _et al._ [2022] C. Watzenböck, M. Fellinger, K. Held, and A. Toschi, Long-term memory magnetic correlations in the Hubbard model: A dynamical mean-field theory analysis, SciPost Phys. 12, 184 (2022).
  * [38] As for many other properties, at the level of a DMFT analysis, the specific choice of the DOS, e.g., the one of a square lattice or the semicircular one of the Bethe lattice, does not have a qualitative impact on the final results. Furthemore, quantitative discrepancies can be considerably minimized, by choosing, as we do here, DOSes with equal standard deviations.
  * [39] In the particle-particle channel a divergence of the respective ${\bm{\Gamma}_{pp,\uparrow\downarrow}^{\nu\nu^{\prime}\Omega}}$ occurs in a BSE of similar form to the $ph$-case, where ${\bm{\chi}_{pp,\uparrow\downarrow}^{\nu\nu^{\prime}\Omega}\\!=\\!\bm{\chi}_{ph,\uparrow\downarrow}^{\nu\nu^{\prime}(\Omega-\nu-\nu^{\prime})}}$:
${\bm{\chi}_{pp,\uparrow\downarrow}^{\nu\nu^{\prime}\Omega}\\!=\\!-\frac{1}{\beta^{2}}\sum_{\nu_{1}\nu_{2}}\left(\bm{\chi}_{pp,0}^{\nu\nu_{1}\Omega}\\!-\\!\bm{\chi}_{pp,\uparrow\downarrow}^{\nu(\Omega-\nu_{1})\Omega}\right)\,\bm{\Gamma}_{pp,\uparrow\downarrow}^{\nu_{1}\nu_{2}\Omega}\,\bm{\chi}_{pp,0}^{\nu_{2}\nu^{\prime}\Omega}}$,
with
${\bm{\chi}_{pp,0}^{\nu\nu^{\prime}\Omega}\\!=\\!-\beta\,G(\nu)\,G(\Omega-\nu)\,\delta_{\nu\nu^{\prime}}}$.
For further details see Ref. [31, 9].

  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time monte carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011).
  * Wallerberger  _et al._ [2019] M. Wallerberger, A. Hausoel, P. Gunacker, A. Kowalski, N. Parragh, F. Goth, K. Held, and G. Sangiovanni, w2dynamics: Local one- and two-particle quantities from dynamical mean field theory, Computer Physics Communications 235, 388 (2019).
  * Kowalski  _et al._ [2019] A. Kowalski, A. Hausoel, M. Wallerberger, P. Gunacker, and G. Sangiovanni, State and superstate sampling in hybridization-expansion continuous-time quantum monte carlo, Phys. Rev. B 99, 155112 (2019).
  * Blümer [2002] N. Blümer,  _Mott-Hubbard Metal-Insulator Transition and Optical Conductivity in High Dimensions_ , Ph.D. thesis, Universität Mainz (2002).
  * Del Re and Rohringer [2021] L. Del Re and G. Rohringer, Fluctuations analysis of spin susceptibility: Néel ordering revisited in dynamical mean field theory, Phys. Rev. B 104, 235128 (2021).
  * [45] That a rescaling of ${U\\!\to\\!U/\eta}$ accounts for a finite shift of $N_{\lambda<0}$ in the HA, can be readily seen by taking a look at the number of negative eigenvalues with anti-symmetric eigenvectors ${n_{\text{anti-symm.}}\\!=\\!\lfloor\frac{\sqrt{3}\beta U}{4\pi}\rfloor}$ of Ref. [9, 14].
  * [46] Note that in the regime of the MIT ${\eta\\!=\\!1.19\\!\pm\\!0.03}$ is roughly equal to half of the interaction value of $U_{c1}$ at ${T\\!=\\!0}$ (${\frac{1}{2}U_{c1}^{T=0}\\!\approx\\!1.19}$), which is the starting point of a metastable Mott insulating phase. For more details about the effective HA like description of the PI phase, see Appendix C.
  * [47] The applied fitting function used to fit the high $\beta$ range (${\beta\\!\geq\\!75}$) was ${N_{\lambda<0}^{U_{c1}}(\beta)\\!=\\!a_{1}\\!+\\!b_{1}(\beta\\!-\\!\beta_{0})^{-1}}$, which behaves as ${N_{\lambda<0}^{U_{c1}}(\beta)\\!\cong\\!a_{1}\\!+\\!b_{1}\beta^{-1}}$ for ${\beta\\!\rightarrow\\!\infty}$.
  * Note [1] Evidently, this also means that the difference of $N_{\lambda<0}(\beta)$ between different paths approaching the $T\tmspace-{.1667em}=\tmspace-{.1667em}0$ MIT, e.g. $U_{c}(\beta)$ or $U_{c2}(\beta)$, does not necessarily vanish (and it might even diverge) for $\beta\rightarrow\infty$, depending on the pair of paths chosen.
  * Note [2] In particular, the $T\tmspace-{.1667em}=\tmspace-{.1667em}0$-extrapolated behavior of $N_{\lambda<0}$ could be consistent with the existence of a unique function describing the increase of $N_{\lambda<0}$ for increasing $U$ along the $T\tmspace-{.1667em}=\tmspace-{.1667em}0$-axis up to its divergence at ${U=U_{c2}(T\tmspace-{.1667em}=\tmspace-{.1667em}0)}$, and, hence, does not rule out the possibility of a smooth evolution of the generalized susceptibility $\bm{\chi}^{\nu\nu^{\prime}}_{c}$ and of the corresponding physical response function.
  * Karski  _et al._ [2005] M. Karski, C. Raas, and G. S. Uhrig, Electron spectra close to a metal-to-insulator transition, Phys. Rev. B 72, 113110 (2005).
  * Raas and Uhrig [2009] C. Raas and G. S. Uhrig, Generic susceptibilities of the half-filled hubbard model in infinite dimensions, Phys. Rev. B 79, 115136 (2009).
  * Toschi  _et al._ [2007] A. Toschi, A. A. Katanin, and K. Held, Dynamical vertex approximation: A step beyond dynamical mean-field theory, Phys. Rev. B 75, 045118 (2007).
  * Valli  _et al._ [2015] A. Valli, T. Schäfer, P. Thunström, G. Rohringer, S. Andergassen, G. Sangiovanni, K. Held, and A. Toschi, Dynamical vertex approximation in its parquet implementation: Application to hubbard nanorings, Phys. Rev. B 91, 115115 (2015).
  * Rohringer  _et al._ [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory, Rev. Mod. Phys. 90, 025003 (2018).
  * Kauch  _et al._ [2020] A. Kauch, P. Pudleiner, K. Astleithner, P. Thunström, T. Ribic, and K. Held, Generic optical excitations of correlated systems: $\pi$-tons, Phys. Rev. Lett. 124, 047401 (2020).
  * Kaufmann  _et al._ [2021] J. Kaufmann, C. Eckhardt, M. Pickem, M. Kitatani, A. Kauch, and K. Held, Self-consistent ladder dynamical vertex approximation, Phys. Rev. B 103, 035120 (2021).
  * Ayral and Parcollet [2016] T. Ayral and O. Parcollet, Mott physics and collective modes: An atomic approximation of the four-particle irreducible functional, Phys. Rev. B 94, 075159 (2016).
  * [58] We do not consider here, the obvious case of second-order phase transitions, characterized by simultaneous divergences of reducible and irreducible vertices.
  * Metzner  _et al._ [2012] W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Functional renormalization group approach to correlated fermion systems, Rev. Mod. Phys. 84, 299 (2012).
  * Bickers [2004] N. Bickers,  _Theoretical Methods for Strongly Correlated Electrons_ , edited by D. Sénéchal, A.-M. Tremblay, and C. Bourbonnais (Springer-Verlag New York Berlin Heidelberg, 2004) pp. 237–296.
  * Taranto  _et al._ [2014] C. Taranto, S. Andergassen, J. Bauer, K. Held, A. Katanin, W. Metzner, G. Rohringer, and A. Toschi, From infinite to two dimensions through the functional renormalization group, Phys. Rev. Lett. 112, 196402 (2014).
  * Wentzell  _et al._ [2015] N. Wentzell, C. Taranto, A. Katanin, A. Toschi, and S. Andergassen, Correlated starting points for the functional renormalization group, Phys. Rev. B 91, 045120 (2015).
  * Vilardi  _et al._ [2019] D. Vilardi, C. Taranto, and W. Metzner, Antiferromagnetic and $d$-wave pairing correlations in the strongly interacting two-dimensional hubbard model from the functional renormalization group, Phys. Rev. B 99, 104501 (2019).
  * Bonetti [2020] P. M. Bonetti, Accessing the ordered phase of correlated fermi systems: Vertex bosonization and mean-field theory within the functional renormalization group, Phys. Rev. B 102, 235160 (2020).
  * Bonetti  _et al._ [2022] P. M. Bonetti, A. Toschi, C. Hille, S. Andergassen, and D. Vilardi, Single-boson exchange representation of the functional renormalization group for strongly interacting many-electron systems, Phys. Rev. Research 4, 013034 (2022).
  * Chalupa-Gantner  _et al._ [2022] P. Chalupa-Gantner, F. B. Kugler, C. Hille, J. von Delft, S. Andergassen, and A. Toschi, Fulfillment of sum rules and ward identities in the multiloop functional renormalization group solution of the anderson impurity model, Phys. Rev. Research 4, 023050 (2022).
  * Hafermann  _et al._ [2009] H. Hafermann, G. Li, A. N. Rubtsov, M. I. Katsnelson, A. I. Lichtenstein, and H. Monien, Efficient perturbation theory for quantum lattice models, Phys. Rev. Lett. 102, 206401 (2009).
  * Krien  _et al._ [2020] F. Krien, A. Valli, P. Chalupa, M. Capone, A. I. Lichtenstein, and A. Toschi, Boson-exchange parquet solver for dual fermions, Phys. Rev. B 102, 195131 (2020).
  * [69] Here, it is important to note that our study specifically addresses the situation of many-electron models with a purely local and SU(2) invariant-interaction. Hence, our considerations do not necessarily apply to other classes of problems, such as non-interacting electrons in the presence of disorder (e.g. the binary mixture) or in the presence of not-SU(2)-symmetric interactions (e.g. the Falicov Kimball model), where irreducible vertex divergences have been also reported [2, 9] (e.g. in the spin $\uparrow\uparrow$ sector). On the other hand, our conclusions remain fully applicable, mutatis mutandis, to the vertex divergences appearing [17] in fermionic models with a local and SU(2)-invariant attractive interaction (i.e., with ${U\\!<\\!0}$) and, in particular, to their corresponding MITs [82].
  * Fabrizio [2022] M. Fabrizio, Emergent quasiparticles at luttinger surfaces, Nat. Commun. 13, 1561 (2022).
  * Kehrein [1998] S. Kehrein, Density of states near the mott-hubbard transition in the limit of large dimensions, Phys. Rev. Lett. 81, 3912 (1998).
  * Del Re  _et al._ [2019] L. Del Re, M. Capone, and A. Toschi, Dynamical vertex approximation for the attractive hubbard model, Phys. Rev. B 99, 045137 (2019).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Quantum cluster theories, Rev. Mod. Phys. 77, 1027 (2005).
  * Park  _et al._ [2011] H. Park, K. Haule, and G. Kotliar, Magnetic excitation spectra in ${\mathrm{bafe}}_{2}{\mathrm{as}}_{2}$: A two-particle approach within a combination of the density functional theory and the dynamical mean-field theory method, Phys. Rev. Lett. 107, 137007 (2011).
  * Schäfer  _et al._ [2015] T. Schäfer, F. Geles, D. Rost, G. Rohringer, E. Arrigoni, K. Held, N. Blümer, M. Aichhorn, and A. Toschi, Fate of the false mott-hubbard transition in two dimensions, Phys. Rev. B 91, 125109 (2015).
  * Pelz  _et al._ [2023] M. Pelz, S. Adler, M. Reitner, and A. Toschi, Numerical results for ”highly nonperturbative nature of the mott metal-insulator transition: Two-particle vertex divergences in the coexistence region”, (1.0.0) [Data set] TU Wien (2023).
  * Kuneš [2011] J. Kuneš, Efficient treatment of two-particle vertices in dynamical mean-field theory, Phys. Rev. B 83, 085102 (2011).
  * Tagliavini  _et al._ [2018] A. Tagliavini, S. Hummel, N. Wentzell, S. Andergassen, A. Toschi, and G. Rohringer, Efficient bethe-salpeter equation treatment in dynamical mean-field theory, Phys. Rev. B 97, 235140 (2018).
  * Wentzell  _et al._ [2020] N. Wentzell, G. Li, A. Tagliavini, C. Taranto, G. Rohringer, K. Held, A. Toschi, and S. Andergassen, High-frequency asymptotics of the vertex function: Diagrammatic parametrization and algorithmic implementation, Phys. Rev. B 102, 085106 (2020).
  * Werner  _et al._ [2006] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Continuous-time solver for quantum impurity models, Phys. Rev. Lett. 97, 076405 (2006).
  * Gunacker  _et al._ [2015] P. Gunacker, M. Wallerberger, E. Gull, A. Hausoel, G. Sangiovanni, and K. Held, Continuous-time quantum monte carlo using worm sampling, Phys. Rev. B 92, 155102 (2015).
  * Toschi  _et al._ [2005] A. Toschi, P. Barone, M. Capone, and C. Castellani, Pairing and superconductivity from weak to strong coupling in the attractive hubbard model, New Journal of Physics 7, 7 (2005).