# Why the Universal Threshold for Primordial Black Hole Formation is Universal

Alex Kehagias  Physics Division, National Technical University of Athens,
Athens, 15780, Greece  Davide Perrone  0000-0003-4430-4914 Department of
Theoretical Physics and Gravitational Wave Science Center,
24 quai E. Ansermet, CH-1211 Geneva 4, Switzerland  Antonio Riotto
0000-0001-6948-0856 Department of Theoretical Physics and Gravitational Wave
Science Center,
24 quai E. Ansermet, CH-1211 Geneva 4, Switzerland

###### Abstract

We show why the threshold for primordial black hole formation is universal
(independent from the shape of the perturbation) when expressed in terms of
the volume averaged compaction function. The proof is rooted in the self-
similarity of the gravitational collapse phenomenon at criticality.

## 1\. Introduction

The topic of Primordial Black Holes (PBHs) has become much debated in the last
years (see Ref. Bagui  _et al._ (2023) for a recent review) since they might
provide an explanation for some of the signals from binary BH mergers measured
by gravitational wave detectors and an important component of the dark matter
in the universe. Furthermore, the next generation of gravitational wave
experiments will provide an armoury of smoking guns to distinguish PBHs from
their astrophysical counterparts at high redshifts Riotto and Silk (2024).

The most common scenario for the formation of the PBHs is the one where they
are originated from the collapse of large fluctuations formed during inflation
upon horizon re-entry Sasaki  _et al._ (2018). The formation probability, and
therefore their current abundance, is extremely sensitive to the critical
threshold of collapse, that is on the minimum value ${\cal C}_{c}$ that the
amplitude of the compaction function (or the energy overdensity) must have in
order for a fluctuation to form a BH Yoo (2022).

Numerical simulations have shown that the critical threshold depends
sensitively on the curvature at the peak of the compaction function Musco
(2019); Germani and Musco (2019); Escrivà  _et al._ (2020); Musco  _et al._
(2021), that is on the shape of the perturbation. The threshold of the
compaction function varies from 2/5 for broad profiles to 2/3 for very peaked
profiles. This is certainly not a pleasant feature since the theory delivers
only stochastic quantities and one can only calculate the averaged profile of
a perturbation. Furthermore, in realistic cases the critical threshold for
formation is determined by the broadest possible compaction function, and not
by the mean profile. Ianniccari  _et al._ (2024a).

This ambiguity can be eliminated by the interesting observation made in Ref.
Escrivà  _et al._ (2020). Based on a numerical approach, it has been shown
that the threshold to form PBHs from an initial spherically symmetric
perturbation is universal when the compaction function is averaged over a
sphere of radius equal to the location of the maximum of the compaction
function. Its critical value in radiation turns out to be 2/5 independently
from the shape of the perturbation and within the numerical errors Escrivà
_et al._ (2020).

The goal of this short paper is to provide a simple proof of why such
universal threshold for the formation of PBH is indeed universal and
independent of the shape of the collapsing compaction function. As we shall
see, our proof is rooted in the fundamental property of the phenomenon of
gravitational collapse: the existence of attractors within the BH threshold,
that is, attractors of codimension one in phase space, and which are self-
similar. More specifically, the persistence of self-similarity at criticality
will play a key role.

## 2\. Some basic concepts

PBH formation is routinely studied for spherical and time-dependent
configurations. Therefore, we will be concerned with time-dependent
spherically symmetric spacetimes. The latter have a residual symmetry which is
more manifest if the metric is written locally in double-null form as (see,
for instance, Ref. Hayward (1996))

${\rm d}s^{2}=-2e^{-f}{\rm d}\xi^{+}{\rm d}\xi^{-}+R^{2}{\rm d}\Omega^{2},$
(1)

where ${\rm d}\Omega^{2}$ refers to the unit sphere and $R$ and $f$ are
functions of the null coordinates $\xi^{\pm}$. In such double-null form the
residual coordinate freedom consists of all the diffeomorphisms

$\xi^{\pm}\rightarrow\widetilde{\xi}^{\pm}(\xi^{\pm}).$ (2)

The key point is that the areal radius $R$ is a geometrical invariant quantity
(even if $f$ is not) under the residual symmetry, that is it transforms like a
scalar (see, for instance, Ref. Faraoni  _et al._ (2017)). One could also use
one spatial and one temporal direction (which we will adopt in the following)
to show it, but there is no unique choice of such directions and this makes it
more difficult to check coordinate invariance.

The Misner-Sharp spherically symmetric gravitational energy (in units
$G_{N}=1$ and $a=\pm$) Misner and Sharp (1964)

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{1}{2}R\left(1-\nabla_{a}R\,\nabla^{a}R\right)$ (3)
$\displaystyle=$
$\displaystyle\frac{1}{2}R+e^{f}\,R\,\partial_{+}R\,\partial_{-}R$

is consequently an invariant. Finally, we define the compaction function,
which we will use to set the criterion for the PBH formation, to be

$C=\frac{2M}{R}.$ (4)

Being constructed from the Misner-Sharp mass and the areal radius, it is as
well a geometric quantity independent from any spherically symmetric
foliation.

The rescaling freedom (2) can be used to fix the null coordinates

$\sqrt{2}e^{-f/2}{\rm d}\xi^{\pm}=A{\rm d}t\pm B{\rm d}r,$ (5)

such that the metric (1) becomes

${\rm d}s^{2}=-A^{2}{\rm d}t^{2}+B^{2}{\rm d}r^{2}+R^{2}{\rm d}\Omega^{2},$
(6)

where now $R$, $A$ and $B$ are functions of $(r,t)$. Denoting the derivatives
with respect to $r$ with a prime and the derivative with respect to $t$ with a
dot, one can write the equations of motion for the Misner-Sharp mass given a
perfect fluid with energy density $\rho$ and pressure density $P$ Misner and
Sharp (1964)

$\displaystyle M^{\prime}$ $\displaystyle=$ $\displaystyle 4\pi R^{2}\rho
R^{\prime},$ $\displaystyle\dot{M}$ $\displaystyle=$ $\displaystyle-4\pi
R^{2}P\dot{R}.$ (7)

Near a regular center $R=0$, one has $R^{\prime}={\cal O}(R^{0})$ and, if
$\partial/\partial t$ is tangent to the centre, $\dot{R}={\cal O}(R)$. As a
consequence, if $\rho$ is ${\cal O}(R^{0})$, then one can conclude in full
generality that Harada and Maeda (2001)

$M(r\sim 0,t)=\frac{4\pi}{3}\rho(r\sim 0,t)R^{3}(r\sim 0,t)+{\cal O}(R^{4}),$
(8)

for any given time. It is important to notice that this behavior is completely
general and gauge-independent. It does not require the solution to be self-
similar, but it only depends on the regularity of the solution at the origin.
We will use this fundamental result in the following.

## 3\. The self-similarity at criticality

The first description of critical phenomena in gravitational collapse goes
back to the seminal references Choptuik (1993); Evans and Coleman (1994). The
BH threshold in the space of initial data for general relativity shows a
surprising simple structure characterized by universality and power-law
scaling of the BH mass (for a review, see Ref. Gundlach (2003)). This is
explained by the existence of exact solutions which are self-similar
attractors near the center of collapse.

The dynamics is the following. A given overdensity of the radiation fluid has
an initial comoving scale larger than the comoving Hubble radius. At some
point in time the latter catches the size of the overdensity and it is at this
moment that numerical simulations provide the threshold value for PBH
formation. It is also the moment when the metric and the fluid quickly
approach a self-similar behaviour Musco and Miller (2013) which depends only
on the variable

$z=\frac{r}{(-t)},\,\,\,\,\,\,t<0$ (9)

and independent from the time variable

$\tau=-\ln(-t).$ (10)

At later times, self-similarity is broken, leading eventually to the formation
of a BH if the evolution is super-critical, that is if the compaction function
at its maximum is larger than a critical value (for a review, see Ref. Bagui
_et al._ (2023)).

The key point for us is that near, but not precisely critical, spacetimes
develop a self-similar region for $r$ much smaller than the initial size of
the configuration, but then exit self-similarity on even smaller, but finite
scales, in order to form a BH.

On the contrary, the critical configuration for BH formation is the one that
approaches the self-similar solution upon horizon re-entry and retains the
self-similar dependence even when $r\rightarrow 0$ and $t\rightarrow 0$. In
other words, by definition, critical configurations out of which one deduces
the threshold for BH formation, are defined to be those that stay always self-
similar.

It is this property of never-ending self-similarity for the critical
configurations which we will use in the next section to define the most
suitable threshold for PBH formation.

## 4\. Why the universal threshold for PBH formation is universal

Let us consider a PBH formation during the radiation phase of the early
universe. We define the volume averaged compaction function at a given time
$t_{0}$ as

$\displaystyle\overline{C}(R)=\frac{3}{R^{3}}\int_{0}^{R}{\rm
d}x\,x^{2}\,C(x),$ (11)

where $C(R(r,t_{0}))=C(r,t_{0})$. The time $t_{0}$ has to be thought of as the
time a given perturbation enters the horizon and the self-similar behaviour
takes over. The quantity (11) is manifestly gauge-invariant being constructed
from gauge-invariant quantities.

We now rescale the spatial and time coordinates performing a dilatation
transformation

$r\to\lambda r\,\,\,\,{\rm and}\,\,\,\,t\to\lambda t,$ (12)

such that the variable $z$ is invariant

$z\to z.$ (13)

Since the spacetime is self-similar, it admits a conformal Killing vector,
which for the metric (6) is given by

$\displaystyle\xi=t\partial_{t}+r\partial_{r},$ (14)

and satisfies

$\nabla_{\mu}\xi_{\nu}+\nabla_{\nu}\xi_{\mu}=2g_{\mu\nu}.$ (15)

It is easy to see that the compaction function defined in Eq. (4) is invariant
under the dilatation (12)

$\delta C=t\partial_{t}C+r\partial_{r}C=0,$ (16)

and depends only on the zooming coordinate $z$. Similarly, the averaged
compaction function is also invariant since

$\delta\overline{C}(R)=3C(R)-3\overline{C}(R)+\frac{3}{R^{3}}\int_{0}^{R}{\rm
d}x\,x^{2}\,t\partial_{t}C(x),$ (17)

and thus, after using Eq. (16), one finds

$\delta\overline{C}(R)=0.$ (18)

This implies that one can calculate $\overline{C}(R)$ at any $r$, and in
particular around $r=0$, by appropriate rescalings of both $r$ and $t$. Since
the areal radius $R$ scales like $r$ (it can be always be expressed as $r$
times a function of $z$ Harada and Maeda (2001)), that is

$R\to\lambda R,$ (19)

then the volume averaged compaction function scales as

$\displaystyle\overline{C}(R)$ $\displaystyle=$
$\displaystyle\frac{3}{(\lambda R)^{3}}\int_{0}^{\lambda R}{\rm d}(\lambda
x)\,({\lambda x})^{2}\,C(\lambda x)$ (20) $\displaystyle=$
$\displaystyle\frac{3}{\overline{R}^{3}}\int_{0}^{\overline{R}}{\rm
d}\overline{x}\,\overline{x}^{2}\,C(\overline{x}),$

where $\overline{R}=\lambda R$.

The key point is that we can now take $\lambda$ as small as we wish in order
to bring $\lambda R$ as close as desired to the central origin. For what
explained in the previous section, for the critical configuration the self-
similar solution persists in this limit independently from how close we are at
$R\sim r\sim 0$.

We can use Eqs. (8) and (4) to express Eq. (20) as

$\displaystyle\overline{C}(R)$ $\displaystyle=$
$\displaystyle\frac{3}{\overline{R}^{3}}\int_{0}^{\overline{R}}{\rm
d}\overline{x}\,\frac{8\pi}{3}\,\overline{x}^{4}\,\rho(0,\lambda t_{0})\,$
(21) $\displaystyle=$
$\displaystyle\frac{3}{\overline{R}^{3}}\int_{0}^{\overline{R}}{\rm
d}\overline{x}\frac{8\pi}{3}\lambda^{-2}\overline{x}^{4}\rho(0,t_{0}),$

where we have used the scaling $\rho(0,\lambda
t_{0})=\lambda^{-2}\rho(0,t_{0})$, which comes from the fact that in full
generality $t^{2}\rho(0,t)=$ constant Harada and Maeda (2001). From Eq. (21)
we find that

$\displaystyle\overline{C}(R)$ $\displaystyle=$
$\displaystyle\frac{3}{5}\left(\frac{8\pi}{3}\rho(0,t_{0})R^{2}(0,t_{0})\right)$
(22) $\displaystyle=$ $\displaystyle\frac{3}{5}\left(\lim_{r\to
0}\frac{2M(r,t_{0})}{R(r,t_{0})}\right)$ $\displaystyle=$
$\displaystyle\frac{3}{5},$

where in the last passage we have used the fact that for critical solution the
formation of the BH imposes

$\displaystyle\lim_{r\to 0}\frac{2M(t_{0},r)}{R(t_{0},r)}=1$ (23)

since the critical BH has zero mass and therefore zero horizon radius. In
general for a fluid of equation of state $\omega$ one has

$\displaystyle\overline{\mathcal{C}}=f(\omega)\overline{C},\,\,\,f(\omega)=\frac{3(1+\omega)}{5+3\omega},$
(24)

such that for radiation ($\omega=1/3$) we finally obtain

$\displaystyle\overline{\mathcal{C}}_{c}=\frac{2}{5}.$ (25)

This is the value obtained by numerical simulations. We have thus proven that
the volume averaged compaction function provides a universal threshold for the
PBH formation because the critical solution exhibits persistent self-
similarity near the center which is rooted in the critical behavior of the
gravitational collapse phenomenon. It will be interesting to understand if our
result has an interepretation in view of the correspondence between the
formation of PBHs and the stability of circular null geodesics around the
collapsing perturbation Ianniccari  _et al._ (2024b), where self-similarity
plays also a crucial role.

###### Acknowledgements.

We thank V. De Luca for useful comments. A.R. acknowledges support from the
Swiss National Science Foundation (project number CRSII5_213497) and from the
Boninchi Foundation for the project “PBHs in the Era of GW Astronomy”.

## References

  * Bagui  _et al._ (2023) E. Bagui _et al._ (LISA Cosmology Working Group), (2023), arXiv:2310.19857 [astro-ph.CO] .
  * Riotto and Silk (2024) A. Riotto and J. Silk, (2024), arXiv:2403.02907 [astro-ph.CO] .
  * Sasaki  _et al._ (2018) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Class. Quant. Grav. 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Yoo (2022) C.-M. Yoo, Galaxies 10, 112 (2022), arXiv:2211.13512 [astro-ph.CO] .
  * Musco (2019) I. Musco, Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Germani and Musco (2019) C. Germani and I. Musco, Phys. Rev. Lett. 122, 141302 (2019), arXiv:1805.04087 [astro-ph.CO] .
  * Escrivà  _et al._ (2020) A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D 101, 044022 (2020), arXiv:1907.13311 [gr-qc] .
  * Musco  _et al._ (2021) I. Musco, V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. D 103, 063538 (2021), arXiv:2011.03014 [astro-ph.CO] .
  * Ianniccari  _et al._ (2024a) A. Ianniccari, A. J. Iovino, A. Kehagias, D. Perrone, and A. Riotto, (2024a), arXiv:2402.11033 [astro-ph.CO] .
  * Hayward (1996) S. A. Hayward, Phys. Rev. D 53, 1938 (1996), arXiv:gr-qc/9408002 .
  * Faraoni  _et al._ (2017) V. Faraoni, G. F. R. Ellis, J. T. Firouzjaee, A. Helou, and I. Musco, Phys. Rev. D 95, 024008 (2017), arXiv:1610.05822 [gr-qc] .
  * Misner and Sharp (1964) C. W. Misner and D. H. Sharp, Phys. Rev. 136, B571 (1964).
  * Harada and Maeda (2001) T. Harada and H. Maeda, Phys. Rev. D 63, 084022 (2001), arXiv:gr-qc/0101064 .
  * Choptuik (1993) M. W. Choptuik, Phys. Rev. Lett. 70, 9 (1993).
  * Evans and Coleman (1994) C. R. Evans and J. S. Coleman, Phys. Rev. Lett. 72, 1782 (1994).
  * Gundlach (2003) C. Gundlach, Phys. Rept. 376, 339 (2003), arXiv:gr-qc/0210101 .
  * Musco and Miller (2013) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013), arXiv:1201.2379 [gr-qc] .
  * Ianniccari  _et al._ (2024b) A. Ianniccari, A. J. Iovino, A. Kehagias, D. Perrone, and A. Riotto, (2024b), arXiv:2404.02801 [astro-ph.CO] .


UTHEP-760, UTCCS-P-139

# Metal-insulator transition in (2+1)-dimensional Hubbard model with tensor
renormalization group

Shinichiro Akiyama<EMAIL_ADDRESS>Graduate School of Pure and
Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan
Yoshinobu Kuramashi<EMAIL_ADDRESS>Center for Computational
Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan  Takumi
Yamashita111Present address: 28-20 Kojogaoka, Otsu, Shiga 520-0821 Japan
<EMAIL_ADDRESS>Center for Artificial Intelligence
Research, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan

###### Abstract

We investigate the doping-driven metal-insulator transition of the
(2+1)-dimensional Hubbard model in the path-integral formalism with the tensor
renormalization group method. We calculate the electron density $\langle
n\rangle$ as a function of the chemical potential $\mu$ choosing three values
of the Coulomb potential with $U=80$, $8$, and $2$ as representative cases of
the strong, intermediate, and weak couplings. We have determined the critical
chemical potential at each $U$, where the Hubbard model undergoes the metal-
insulator transition from the half-filling plateau with $\langle n\rangle=1$
to the metallic state with $\langle n\rangle>1$. Our results indicate that the
model exhibits the metal-insulator transition over the vast region of the
finite coupling $U$.

I44

##  1 Introduction

The Hubbard model, which is a simple theoretical model to describe electron
systems with repulsive Coulomb interactions, is expected to have rich phase
structures so that it has been attracting the interest of not only condensed
matter physicists but also particle physicists. It has been widely known that
the Hubbard model has a similar path-integral form to the Nambu–Jona-Lasinio
(NJL) model Nambu:1961tp ; Nambu:1961fr , which is a low energy effective
theory in Quantum Chromodynamics (QCD): Both consist of a hopping term and a
four-fermi interaction term. Their similarity, unfortunately, leads to sharing
the so-called sign problem, which is a notorious difficulty in the numerical
analyses based on the Monte Carlo approach.

Recently the authors have successfully applied the tensor renormalization
group (TRG) method222In this paper the TRG method or the TRG approach refers
to not only the original numerical algorithm proposed by Levin and Nave
Levin:2006jai  but also its extensions PhysRevB.86.045139 ; Shimizu:2014uva ;
Sakai:2017jwp ; Adachi:2019paf ; Kadoh:2019kqk ; Akiyama:2020soe ;
adachi2020bondweighted ; Kadoh:2021fri . to investigate the phase transition
of the four-dimensional ($4d$) NJL model at high density and very low
temperature Akiyama:2020soe . This work was followed by the application of the
TRG method to analyze the metal-insulator transition of the $(1+1)d$ Hubbard
model by calculating the electron density $\langle n\rangle$ as a function of
the chemical potential $\mu$ Akiyama:2021xxr . Our results for the critical
chemical potential $\mu_{\rm c}$ and the critical exponent $\nu$ are
consistent with an exact solution based on the Bethe ansatz
PhysRevLett.20.1445 ; LIEB20031 .

In this paper, we apply the TRG method to investigate the doping-driven metal-
insulator transition in the $(2+1)d$ Hubbard model. 333The model has also been
investigated by the tensor network method based on the Hamiltonian formalism,
like a fermionic PEPS, which is also free from the sign problem. For a recent
study, see Ref. Schneider:2021hqj , for example. The TRG method, which was
originally proposed to study two-dimensional (2$d$) classical spin systems
Levin:2006jai , has been developed to study wide varieties of fermionic models
in particle physics Shimizu:2014uva ; Shimizu:2014fsa ; Shimizu:2017onf ;
Takeda:2014vwa ; Sakai:2017jwp ; Yoshimura:2017jpk ; Kadoh:2018hqq ;
Kadoh:2019ube ; Kuramashi:2019cgs ; Akiyama:2020ntf ; Akiyama:2020soe ;
PhysRevD.101.094509 . It is also confirmed that the TRG method does not suffer
from the sign problem by studying various quantum field theories
Shimizu:2014uva ; Shimizu:2014fsa ; Shimizu:2017onf ; Takeda:2014vwa ;
Kawauchi:2017dnj ; Kadoh:2018hqq ; Kadoh:2019ube ; Kuramashi:2019cgs ;
Akiyama:2020ntf ; Akiyama:2020soe ; Akiyama:2021xxr ; Bloch:2021mjw ;
Nakayama:2021iyp . We calculate the electron density $\langle n\rangle$ as a
function of $\mu$ with three choices of $U=80$, $8$ and $2$. The $\mu$
dependence of $\langle n\rangle$ allows us to determine the critical chemical
potential $\mu_{\rm c}$ at the doping-driven metal-insulator transition from
the half-filling plateau with $\langle n\rangle=1$ to the metallic state with
$\langle n\rangle>1$. Our results at $U=80$, 8 and 2 show that $|\mu_{\rm
c}-U/2|$ monotonically diminishes as $U$ decreases and seems to converge on
$|\mu_{\rm c}-U/2|=0$ at $U=0$. This indicates the possibility that the model
exhibits the metal-insulator transition over the wide region of the finite
coupling.

This paper is organized as follows. In Sec. 2 we define the Hubbard model in
the path-integral formalism and give a brief description of the numerical
algorithm. In Sec. 3 we present the $\mu$ dependence of the electron density
and determine the critical chemical potential $\mu_{\rm c}$ at the doping-
driven metal-insulator transition. Section 4 is devoted to summary and
outlook.

##  2 Formulation and numerical algorithm

###  2.1 (2+1)-dimensional Hubbard model in the path-integral formalism

We consider the partition function of the Hubbard model in the path-integral
formalism on an anisotropic lattice with the physical volume $V=L_{x}\times
L_{y}\times\beta$, whose spatial extension is defined as
$L_{\sigma}=aN_{\sigma}~{}(\sigma=x,y)$ with $a$ the spatial lattice spacing.
$\beta$ denotes the inverse temperature, which is divided as
$\beta=1/T=\epsilon N_{\tau}$. Following Ref. Akiyama:2021xxr , the path-
integral expression of the partition function is given by

$\displaystyle
Z=\int\left(\prod_{n\in\Lambda_{2+1}}\prod_{s=\uparrow,\downarrow}{\rm
d}\bar{\psi}_{s}(n){\rm d}\psi_{s}(n)\right){\rm e}^{-S},$ (1)

where $n=(n_{x},n_{y},n_{\tau})\in\Lambda_{2+1}(\subset\mathbb{Z}^{3})$
specifies a position on the lattice $|\Lambda_{2+1}|=N_{x}\times N_{y}\times
N_{\tau}$. Since the Hubbard model describes the spin-1/2 fermions, they are
labeled by $s=\uparrow,\downarrow$, corresponding to the spin-up and spin-
down, respectively. Introducing the notation,

$\displaystyle\psi(n)=\left(\begin{array}[]{c}\psi_{\uparrow}(n)\\\
\psi_{\downarrow}(n)\end{array}\right),~{}\bar{\psi}(n)=\left(\bar{\psi}_{\uparrow}(n),\bar{\psi}_{\downarrow}(n)\right),$
(4)

the action $S$ is defined as

$\displaystyle S$ $\displaystyle=\sum_{n_{\tau},n_{x},n_{y}}\epsilon
a^{2}\left\\{\bar{\psi}(n)\left(\frac{\psi(n+{\hat{\tau}})-\psi(n)}{\epsilon}\right)\right.$
$\displaystyle\left.-t\sum_{\sigma=x,y}\left(\bar{\psi}(n+{\hat{\sigma}})\psi(n)+\bar{\psi}(n)\psi(n+{\hat{\sigma}})\right)+\frac{U}{2}\left(\bar{\psi}(n)\psi(n)\right)^{2}-\mu\bar{\psi}(n)\psi(n)\right\\}.$
(5)

The kinetic terms in the spatial directions contain the hopping parameter $t$.
The four-fermi interaction term represents the Coulomb repulsion of electrons
at the same lattice site. The chemical potential is denoted by the parameter
$\mu$. Note that the half-filling is realized at $\mu=U/2$ in the current
definition. We assume the periodic boundary condition in the spatial
direction, $\psi(N_{x}+1,n_{y},n_{\tau})=\psi(1,n_{y},n_{\tau})$ and
$\psi(n_{x},N_{y}+1,n_{\tau})=\psi(n_{x},1,n_{\tau})$, while the anti-periodic
one in the temporal direction,
$\psi(n_{x},n_{y},N_{\tau}+1)=-\psi(n_{x},n_{y},1)$. In the following
discussion, we always set $a=1$.

###  2.2 Numerical algorithm

Based on Ref. Akiyama:2020sfo , the tensor network representation of Eq. (1)
is immediately obtained. Set $d=2$ in the Appendix of Ref. Akiyama:2021xxr
and one can find out the Grassmann tensor which generates the Grassmann tensor
network of Eq. (1). The resulting Grassmann tensor
$\mathcal{T}_{\Psi_{x}\Psi_{y}\Psi_{\tau}\bar{\Psi}_{\tau}\bar{\Psi}_{y}\bar{\Psi}_{x}}$
is of rank 6 and we evaluate $\mathrm{gTr}[\prod_{n}\mathcal{T}]$ with the
anisotropic TRG (ATRG) algorithm Adachi:2019paf  whose extension to the
Grassmann integrals, referred as the Grassmann ATRG (GATRG), is given in Ref.
Akiyama:2020soe . We also follow the coarse-graining procedure employed in the
previous study of $(1+1)d$ Hubbard model Akiyama:2021xxr . Firstly, we carry
out $m_{\tau}$ times of renormalization along with the temporal direction,
which can be seen as the imaginary time evolution of the local Grassmann
tensor. Secondly, $3d$ ATRG procedure is applied as the spacetime coarse-
graining. As in the case of $(1+1)d$ Hubbard model, we have found that the
optimal $m_{\tau}$ satisfies the condition $\epsilon 2^{m_{\tau}}\sim
O(10^{-1})$ in the sense of preserved tensor norm.

##  3 Numerical results

###  3.1 Algorithmic-parameter dependence

The partition function of Eq. (1) is evaluated using the numerical algorithm
explained above on lattices with the physical volume $V=L_{x}\times
L_{y}\times\beta=N_{x}\times N_{y}\times\epsilon N_{\tau}$ with
$(N_{x},N_{y},N_{\tau})=(2^{m},2^{m},2^{m_{\tau}+m})~{}(m,m_{\tau}\in\mathbb{N})$.
We employ $U=80$, 8 and 2 for the four-fermi coupling with $t=1$ for the
hopping parameter. In Fig. 1 we plot the $\mu$ dependence of the thermodynamic
potential $\ln Z/V$ at $U=8$ on $V=4096^{2}\times 1677.7216$ with the bond
dimension $D=80$ in the GATRG algorithm choosing $\epsilon=2^{12}\times
10^{-4},2^{8}\times 10^{-4},2^{4}\times 10^{-4},10^{-4}$. For each value of
$\epsilon$, $m_{\tau}$ is chosen via the condition $\epsilon
2^{m_{\tau}}=2^{12}\times 10^{-4}=O(10^{-1})$ following Ref. Akiyama:2021xxr .
We find clear discretization effects for the coarsest case of
$\epsilon=2^{12}\times 10^{-4}$. On the other hand, the results with
$\epsilon=2^{4}\times 10^{-4}$ and $10^{-4}$ show good consistency. This means
that the discretization effects with $\epsilon=10^{-4}$ are negligible. We
employ $\epsilon=10^{-4}$ in the following calculations.

Figure 1: Thermodynamic potential at $U=8$ on $V=4096^{2}\times 1677.7216$
lattice. $\beta$ is divided with $\epsilon=2^{12}\times 10^{-4}=0.4096$,
$2^{8}\times 10^{-4}=0.0256$, $2^{4}\times 10^{-4}=0.0016$ and
$10^{-4}=0.0001$. The bond dimension is set to be $D=80$.

We investigate the convergence behavior of the thermodynamic potential by
defining the quantity

$\displaystyle\delta=\left|\frac{\ln Z(D)-\ln Z(D=80)}{\ln Z(D=80)}\right|$
(6)

on $V=4096^{2}\times 1677.7216$ lattice with $\epsilon=10^{-4}$. In Fig. 2, we
plot the $D$ dependence of $\delta$ at $(U,t)=(8,1)$ with the choices of
$\mu=6.0,7.5$ and 8.5. As we will see below, $\mu=6.0$ corresponds to $\langle
n\rangle\approx 1.0$ and $\mu=8.5$ does to $\langle n\rangle\approx 1.5$. We
observe that $\delta$’s at these values of $\mu$ decrease as a function of
$D$, though some of them are fluctuating.

Figure 2: Convergence behavior of thermodynamic potential with $\delta$ of Eq.
(6) at $U=8$ with $\mu=6.0,7.5$ and 8.5 as a function of $D$ on
$V=4096^{2}\times 1677.7216$ lattice.

###  3.2 Strong coupling limit

We first consider the atomic limit at $(U,t)=(8,0)$. This case is analytically
solvable. The electron density $\langle n\rangle$ is obtained by the numerical
derivative of the thermodynamic potential in terms of $\mu$:

$\displaystyle\langle n\rangle=\frac{1}{V}\frac{\partial\ln
Z(\mu)}{\partial\mu}\approx\frac{1}{V}\frac{\ln Z(\mu+\Delta\mu)-\ln
Z(\mu-\Delta\mu)}{2\Delta\mu}.$ (7)

In Fig. 3 we compare the numerical and exact results for the $\mu$ dependence
of $\langle n\rangle$ at $(U,t)=(8,0)$. Note that we set $m_{\tau}=24$ because
this case is equivalent to the model defined on $V=1\times\beta$ lattice.
Thanks to the vanishing hopping structure in the spatial direction, we can
check the validity of the imaginary time evolution explained in Sec. 2.2. The
agreement of our numerical result with the exact solution shows that the
imaginary time evolution carried out by the GATRG works precisely. The
deviation from the exact value defined by

$\displaystyle\delta_{\rm exact}(D)=\left|\frac{\ln Z(D)-\ln Z_{\rm
exact}}{\ln Z_{\rm exact}}\right|$ (8)

is at most $O(10^{-4})$ in the range of $0\leq\mu\leq 16$. For $\mu<0$, the
exact thermodynamic potential is equal to zero and the results obtained by the
TRG are also equal to zero within a level of double precision.

Since the phase diagram of the metal-insulator transition in the (2+1)$d$
Hubbard model is not well known so far, we investigate the $\mu$ dependence of
$\langle n\rangle$ choosing $U=80$ as a representative case in the strong
coupling region. In Fig. 4 we plot the electron density $\langle n\rangle$ as
a function of $\mu$ in the vicinity of $\mu\sim U$ with $D=80$. We have
checked that the convergence behavior of $\delta$ at $U=80$ is better than
that at $U=8$. We observe that the electron density starts to increase from
$\langle n\rangle=1$ at $\mu=77.0(2)$ or $\mu/U=0.9625(25)$ and reaches
$\langle n\rangle=2$ with $\mu\gtrsim 83.0$ or $\mu/U\gtrsim 1.04$. The $\mu$
dependence of $\langle n\rangle$ is smooth and continuous so that there is no
signal of the first-order phase transition. We expect that the critical
chemical potential $\mu_{\rm c}$ at the doping-driven metal-insulator
transition approaches to $\mu_{\rm c}/U=1$ toward the atomic limit and the
transition from $\langle n\rangle=1$ to $\langle n\rangle=2$ becomes a step-
function as a function of $\mu/U$.

Figure 3: Electron density $\langle n\rangle$ in the $(U,t)=(8,0)$ case at
$\beta=1677.7216$ with $\epsilon=10^{-4}$ as a function of $\mu$. The solid
line shows the exact solution and the blue circles are the results obtained by
the TRG. Figure 4: Electron density $\langle n\rangle$ as a function of $\mu$
at $(U,t)=(80,1)$ on $V=4096^{2}\times 1677.7216$ with $\epsilon=10^{-4}$. The
bond dimension is set to be $D=80$.

###  3.3 Critical chemical potential at $U=8$ and $2$

Now let us investigate the metal-insulator transition in the intermediate
coupling region at $(U,t)=(8,1)$. There are a lot of previous work to
investigate a possible superconducting phase expected in this coupling region.
Since we are interested in the thermodynamic and zero-temperature limit, we
first check the volume dependence of the electron density $\langle n\rangle$.
In Fig. 5 we plot the $\mu$ dependence of $\langle n\rangle$ at $U=8$ changing
the lattice sizes with $\epsilon=10^{-4}$, $m_{\tau}=12$ and $D=80$. We
observe that the size of $(N_{x},N_{y},N_{\tau})=(2^{12},2^{12},2^{24})$,
which corresponds to $V=4096^{2}\times 1677.7216$, is sufficiently large to be
identified as the thermodynamic and zero-temperature limit. We observe the
$\langle n\rangle=0$ plateau for $\mu\lesssim-4$ and the $\langle n\rangle=2$
one for $12\lesssim\mu$. The half-filling state is characterized by the
plateau of $\langle n\rangle=1$ in the range of $2\lesssim\mu\lesssim 6$.
These plateaus yield the vanishing compressibility $\kappa=\partial\langle
n\rangle/\partial\mu$ indicating the insulating states.

Figure 5: Electron density $\langle n\rangle$ at $U=8$ on two lattice sizes,
$V=256^{2}\times 104.8576$ and $4096^{2}\times 1677.7216$, as a function of
$\mu$. The bond dimension is set to be $D=80$. Figure 6: Electron density
$\langle n\rangle$ at $U=8$ on $V=4096^{2}\times 1677.7216$ with
$\epsilon=10^{-4}$ as a function of $\mu$. The bond dimensions are $D=80$,
$72$, $64$ and $56$. Fit results are drawn by the solid lines for each bond
dimension.

In order to determine the critical chemical potential $\mu_{\rm c}$ in the
limit of $D\rightarrow\infty$ at $U=8$ on $V=4096^{2}\times 1677.7216$ lattice
we make a global fit of $\langle n\rangle$ with $D=80$, 72, 64 and 56 in the
metallic phase near the transition point. In Fig. 6 we plot the results of
$\langle n\rangle$ at $D=80$, $72$, $64$ and $56$ with a much finer resolution
of $\Delta\mu$ than Fig. 5 focusing on the range of $6.0\leq\mu\leq 8.2$,
which covers the region of $1\leq\langle n\rangle\leq 1.5$. This figure
provides us a closer look at the $\mu$ dependence of $\langle n\rangle$ around
$\mu_{\rm c}$. The results at $D=80$, $72$, $64$ and $56$ are almost
degenerate indicating the small $D$ dependence. For the global fit we employ
the following quadratic fitting function:

$\displaystyle\langle n\rangle=1+\alpha\left(\mu-\mu_{\rm
c}(D)\right)+\beta\left(\mu-\mu_{\rm c}(D)\right)^{2}$ (9)

with $\mu_{\rm c}(D)=\mu_{\rm c}(D=\infty)+\gamma/D$, where $\alpha$, $\beta$,
$\gamma$ and $\mu_{\rm c}(D=\infty)$ are the fit parameters. The solid curves
in Fig. 6 represent the fit results over the range of $6.3\leq\mu\leq 8.0$. We
obtain $\mu_{\rm c}(D=\infty)=6.43(4)$, which is presented in Table 1 together
with other fitting results. It may be instructive to compare Figs. 5, 6, and
the estimated location of $\mu_{\rm c}(D=\infty)$ with numerical data in Refs.
PhysRevB.67.085103 ; PhysRevB.87.035110 , though their calculations are
carried out on very small lattice sizes and at low but finite temperatures.

Table 1: Critical chemical potential $\mu_{\rm c}(D=\infty)$ at $U=8$ and $2$ determined by the global fit of the data at $D=80$, $72$, $64$ and $56$. $U$ | 8 | 2
---|---|---
fit range | [6.3, 8.0] | [1.2, 3.4]
$\mu_{\rm c}(D=\infty)$ | 6.43(4) | 1.30(6)
$\alpha$ | 0.372(9) | 0.39(1)
$\beta$ |  $-$0.051(6) |  $-$0.054(5)
$\gamma$ |  $-$7(2) |  $-$13(4)

We repeat the same analysis for the weak coupling case at $U=2$. We apply the
fit function of Eq. (9) to four data sets with the bond dimensions of $D=80$,
72, 64 and 56. Fit results are depicted in Fig. 7 and their numerical values
are summarized in Table 1. Our results show that the deviation of $|\mu_{\rm
c}(D)-U/2|$ diminishes as the Coulomb potential $U$ decreases. It is likely
that $|\mu_{\rm c}(D)-U/2|$ vanishes only at $U=0$. This means that the model
exhibits the metal-insulator transition over the wide regime of the finite
coupling, including the weak coupling region. This conclusion may provide us a
different scenario of the phase diagram from that predicted by the dynamical
mean-field theory (DMFT) RevModPhys.68.13 ; there exists some $U_{\rm c}$ such
that no metal-insulator transition occurs with $U<U_{\rm c}$.

Figure 7: Electron density $\langle n\rangle$ at $U=2$ on $V=4096^{2}\times
1677.7216$ with $\epsilon=10^{-4}$ as a function of $\mu$. The bond dimensions
are $D=80$, $72$, $64$ and $56$. Fit results are drawn by the solid lines for
each bond dimension.

##  4 Summary and outlook

We have investigated the doping-driven metal-insulator transition of the
(2+1)$d$ Hubbard model in the path-integral formalism employing the TRG
method. The electron density $\langle n\rangle$ is calculated in the wide
range of $\mu$ corresponding to $0\leq\langle n\rangle\leq 2$. We have also
determined the critical chemical potential $\mu_{\rm c}$ at three values of
$U$. Our results indicate that the deviation $|\mu_{\rm c}-U/2|$ vanishes only
at $U=0$. This means that the model exhibits the metal-insulator transition
over the vast regime of the finite coupling $U$. As a next step, it would be
interesting to investigate the metal-insulator transition of the (3+1)$d$
Hubbard model.

## Acknowledgment

Numerical calculation for the present work was carried out with the
supercomputer Fugaku provided by RIKEN (Project ID: hp200314) and also with
the Oakforest-PACS (OFP) computer under the Interdisciplinary Computational
Science Program of Center for Computational Sciences, University of Tsukuba.
This work is supported in part by Grants-in-Aid for Scientific Research from
the Ministry of Education, Culture, Sports, Science and Technology (MEXT) (No.
20H00148) and JSPS KAKENHI Grant Number JP21J11226 (S.A.).

## References

  * (1) Yoichiro Nambu and G. Jona-Lasinio, Phys. Rev., 122, 345–358 (1961).
  * (2) Yoichiro Nambu and G. Jona-Lasinio, Phys. Rev., 124, 246–254 (1961).
  * (3) Michael Levin and Cody P. Nave, Phys. Rev. Lett., 99(12), 120601 (2007), arXiv:cond-mat/0611687.
  * (4) Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Phys. Rev. B, 86, 045139 (Jul 2012).
  * (5) Yuya Shimizu and Yoshinobu Kuramashi, Phys. Rev., D90(1), 014508 (2014), arXiv:1403.0642.
  * (6) Ryo Sakai, Shinji Takeda, and Yusuke Yoshimura, PTEP, 2017(6), 063B07 (2017), arXiv:1705.07764.
  * (7) Daiki Adachi, Tsuyoshi Okubo, and Synge Todo, Phys. Rev. B, 102(5), 054432 (2020), arXiv:1906.02007.
  * (8) Daisuke Kadoh and Katsumasa Nakayama (2019), arXiv:1912.02414.
  * (9) Shinichiro Akiyama, Yoshinobu Kuramashi, Takumi Yamashita, and Yusuke Yoshimura, JHEP, 01, 121 (2021), arXiv:2009.11583.
  * (10) Daiki Adachi, Tsuyoshi Okubo, and Synge Todo (11 2020), arXiv:2011.01679.
  * (11) Daisuke Kadoh, Hideaki Oba, and Shinji Takeda (7 2021), arXiv:2107.08769.
  * (12) Shinichiro Akiyama and Yoshinobu Kuramashi, Phys. Rev. D, 104(1), 014504 (2021), arXiv:2105.00372.
  * (13) Elliott H. Lieb and F. Y. Wu, Phys. Rev. Lett., 20, 1445–1448 (Jun 1968).
  * (14) Elliott H. Lieb and F.Y. Wu, Physica A: Statistical Mechanics and its Applications, 321(1), 1–27, Statphys-Taiwan-2002: Lattice Models and Complex Systems (2003).
  * (15) Manuel Schneider, Johann Ostmeyer, Karl Jansen, Thomas Luu, and Carsten Urbach, Phys. Rev. B, 104(15), 155118 (2021), arXiv:2106.13583.
  * (16) Yuya Shimizu and Yoshinobu Kuramashi, Phys. Rev., D90(7), 074503 (2014), arXiv:1408.0897.
  * (17) Yuya Shimizu and Yoshinobu Kuramashi, Phys. Rev., D97(3), 034502 (2018), arXiv:1712.07808.
  * (18) Shinji Takeda and Yusuke Yoshimura, PTEP, 2015(4), 043B01 (2015), arXiv:1412.7855.
  * (19) Yusuke Yoshimura, Yoshinobu Kuramashi, Yoshifumi Nakamura, Shinji Takeda, and Ryo Sakai, Phys. Rev., D97(5), 054511 (2018), arXiv:1711.08121.
  * (20) Daisuke Kadoh, Yoshinobu Kuramashi, Yoshifumi Nakamura, Ryo Sakai, Shinji Takeda, and Yusuke Yoshimura, JHEP, 03, 141 (2018), arXiv:1801.04183.
  * (21) Daisuke Kadoh, Yoshinobu Kuramashi, Yoshifumi Nakamura, Ryo Sakai, Shinji Takeda, and Yusuke Yoshimura, JHEP, 02, 161 (2020), arXiv:1912.13092.
  * (22) Yoshinobu Kuramashi and Yusuke Yoshimura, JHEP, 04, 089 (2020), arXiv:1911.06480.
  * (23) Shinichiro Akiyama, Daisuke Kadoh, Yoshinobu Kuramashi, Takumi Yamashita, and Yusuke Yoshimura, JHEP, 09, 177 (2020), arXiv:2005.04645.
  * (24) Nouman Butt, Simon Catterall, Yannick Meurice, Ryo Sakai, and Judah Unmuth-Yockey, Phys. Rev. D, 101, 094509 (May 2020).
  * (25) Hikaru Kawauchi and Shinji Takeda, EPJ Web Conf., 175, 11015 (2018), arXiv:1710.09804.
  * (26) Jacques Bloch, Raghav G. Jha, Robert Lohmayer, and Maximilian Meister, Phys. Rev. D, 104(9), 094517 (2021), arXiv:2105.08066.
  * (27) Katsumasa Nakayama, Lena Funcke, Karl Jansen, Ying-Jer Kao, and Stefan Kühn (7 2021), arXiv:2107.14220.
  * (28) Shinichiro Akiyama and Daisuke Kadoh, JHEP, 10, 188 (2021), arXiv:2005.07570.
  * (29) J. Bonča and P. Prelovšek, Phys. Rev. B, 67, 085103 (Feb 2003).
  * (30) Kaden R. A. Hazzard, Ana Maria Rey, and Richard T. Scalettar, Phys. Rev. B, 87, 035110 (Jan 2013).
  * (31) Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg, Rev. Mod. Phys., 68, 13–125 (Jan 1996).