# New inflationary exact solution from Lie symmetries

Andronikos Paliathanasis<EMAIL_ADDRESS>Institute of Systems Science,
Durban University of Technology, Durban 4000, South Africa  Instituto de
Ciencias Físicas y Matemáticas, Universidad Austral de Chile, Valdivia
5090000, Chile

###### Abstract

For the inflaton field we determine a new exact solution by using the Lie
symmetry analysis. Specifically, we construct a second-order differential
master equation for arbitrary scalar field potential by assuming that the
spectral index for the density perturbations $n_{s}$ and the scalar to tensor
ratio $r$ are related as $n_{s}-1=h\left(r\right)$. Function $h\left(r\right)$
is classified according to the admitted Lie symmetries for the master
equation. The possible admitted Lie symmetries form the $A_{2}$, $A_{3,2}$,
$A_{3,3}$ and $sl\left(3,R\right)$ Lie algebras. The new inflationary solution
is recovered by the Lie symmetries of the $A_{3,3}$ algebra. Scalar field
potential is derived explicitly, while we compare the resulting spectral
indices with the observations.

Inflation; scalar field; exact solutions; Lie symmetries

##  I Introduction

Inflation is the main mechanism to explain the homogeneity and isotropy of the
observable universe at present time. During the inflationary period, the
universe was dominated by the inflaton. The inflaton drove the dynamics such
that to provide the expansion guth . Nevertheless, the inflationary models are
mainly defined on homogeneous spacetimes, or on background spaces with small
inhomogeneities st1 ; st2 .

In w1  it was found that the presence of a positive cosmological constant in
Bianchi cosmologies leads to expanding Bianchi spacetimes, evolving toward an
rapidly expanding de Sitter universe. That was the first result to support the
cosmic “no-hair” conjecture nh1 ; nh2 . The conjecture states that all
expanding universes with a positive cosmological constant admit as asymptotic
expanding solution. The necessity of the rapid expansion is that it provides a
rapid expansion for the size of the universe such that the latter effectively
loses its memory on the initial conditions, which means that the rapid
expansion solves the “flatness”, “horizon” and monopole problem f1 ; f2 .

Because inflation occurred in a finite time period during the cosmological
evolution, instead of the cosmological constant, the inflaton is assumed to be
described by a scalar field, which provides a dynamical behavior for the
evolution of the field equations. In the single scalar field inflationary
models the expansion appears when the scalar field potential dominates. In
these models, the inflationary period is determined explicitly by the nature
of the scalar field potential. In the literature, there is a plethora of
scalar field potentials which have been proposed the last decades, see for
instance ref1a ; ref1 ; ref2 ; ref3 ; newinf ; ref4 ; ref5a ; ref5 ; ref6a ;
ref6 ; ref7 ; ref8 ; ref9 ; ref10 ; ref11 .

Moreover, single scalar field gravitational models attribute the additional
degrees of freedom provided by modified theories of gravity, such in the
quadratic gravity sa1 . The quadratic gravity is an geometric approach for the
construction of the scalar field potential, through a conformal transformation
with the use of a Lagrange multiplier sot1 . The quadratic gravity, is also
known as Starobinsky model for inflation and it is the model which is mainly
supported by the recent cosmological observations pin1 . There is a plethora
of inflationary models which have been proposed the last years, multifield
inflationary models or inflationary models which follow from modifications of
the Einstein-Hilbert Action Integral, see for instance Aref2 ; Aref3 ; Aref3b
; Aref5 ; Aref6 ; Aref6b ; Aref7 ; ar1 ; ar2 ; ar3 ; ar4 . A catalogue with
the viable inflationary models published in ec1 , while the updated version
published later in ec2  presents the viable inflationary models after the
release of Planck 2013 data.

The field equations for the inflaton model are that of a minimally coupled
scalar field in the context of General Relativity with a spatially flat
Friedmann–Lemaître–Robertson–Walker (FLRW) as a background space. The field
equations are nonlinear differential equations of second-order and form a
singular Lagrangian system with dynamical variables the scalar factor of the
metric tensor. The existence of exact solutions for the field equations
depends on the functional form of the scalar field potential ell1 ; ell2
which drives the scalar field dynamics. In sym1 ; sym2 , the scalar field
potential has been constraint such that the field equations to admit
conservation laws. The existence of the conservation laws indicates the
Liouville integrability for the field equations such that exact and analytic
solutions to be determined. However, an innovative approach for the
determination of analytic solutions in scalar field cosmology applied in ns1 .
Because the field equations form a singular dynamical system, there exist
infinity number of nonlocal conservation laws, which can be used to reduce the
order for the field equations. Indeed, with the use of the nonlocal
conservation laws, the generic algebraic solution has been found, for
arbitrary potential function. The results of ns1  applied for the derivation
of new inflationary solutions in ref10 , while this approach applied for the
reconstruction of the inflationary potential in ref11 .

Specifically, in ref11  it was assumed that spectral index for the density
perturbations, $n_{s}$, and the tensor-to-scalar ratio, $r$, are related by a
function such that $n_{s}-1=h\left(r\right)$. where $h\left(r\right)$ is
either constant, linear or quadratic function in $r$. By applying the results
of ns1  in each case a master differential equation has been defined which
gives the resulting inflationary potential. The resulting inflationary
solutions fits the cosmological observations, while the closed-form
expressions for the scalar field potential were found. In this piece of work,
we focus on the classification of function $h\left(r\right)$ according to the
admitted Lie symmetries for the master equations. This geometric selection
rule is inspired by the results of the applications of the symmetry analysis
in sym1 . The classification scheme that we follow, it was established by
Ovsiannikov on the classification of the unknown function for a nonlinear
Schrödinger equation ovs1 . A similar reconstruction approach for the slow-
roll inflation can be found in the series of studies ff1 ; ff2  where the
scalar-to-tensor ratio has been assumed to be a specific function with
independent variable the number of e-folds. The corresponding scalar-tensor
theory and the resulting $f\left(R\right)$-theory determined. The plan of the
paper is as follows.

In Section II we define the gravitational model for a single scalar field
inflationary theory and we present Einstein’s field equations. Moreover, we
present the generic algebraic solution for arbitrary potential found in ns1 ,
and we define the master equation of our analysis. The Lie symmetry analysis
for the master equation of our consideration is performed in Section III. In
Section IV we present the new inflationary exact solution as also we calculate
the spectral indices $n_{s}$ and $r$. Finally, in Section V we draw our
conclusions.

##  II The inflaton field

Consider the Action Integral of Einstein-Hilbert Action with a scalar field
minimally coupled to gravity, that is,

$S=\int
dx^{4}\sqrt{-g}\left[R-\frac{1}{2}g^{\mu\nu}\phi_{;\mu}\phi_{;\nu}+V(\phi)\right]$
(1)

in which $R$ is the Ricci-scalar for the background space with metric tensor
$g_{\mu\nu}\left(x^{\kappa}\right)$, $\phi\left(x^{\kappa}\right)$ is the
scalar field and $V\left(\phi\left(x^{\kappa}\right)\right)$ is the scalar
field potential. Variation with respect to the metric tensor provides the
Einstein-field equations

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=T_{\mu\nu}^{\left(\phi\right)}$ (2)

where $R_{\mu\nu}$ is the Ricci tensor, and $T_{\mu\nu}^{(\phi)}$ is the
energy-momentum tensor which corresponds to the scalar field.
$T_{\mu\nu}^{\left(\phi\right)}$ is defined as

$T_{\mu\nu}^{\left(\phi\right)}=\phi_{;\mu}\phi_{;\nu}-g_{\mu\nu}\left(\frac{1}{2}\phi^{;\kappa}\phi_{;\kappa}-V\left(\phi\right)\right).$
(3)

Furthermore, variation with respect to the scalar field in the Action
Integral, provides the equation of motion for the scalar field, that is, the
Klein-Gordon equation
$T_{~{}~{}~{}~{}~{}~{}~{}~{}~{};\nu}^{\left(\phi\right)\mu\nu}=0,$ that is,

$g^{\mu\nu}\phi_{;\mu\nu}-V\left(\phi\right)=0.$ (4)

For the background space we consider a spatially flat FLRW universe with line
element

$ds^{2}=-dt^{2}+a^{2}\left(t\right)\left(dx^{2}+dy^{2}+dz^{2}\right).$ (5)

Function $a\left(t\right)$ is the scale factor while the Hubble function is
defined as $H=\frac{\dot{a}\,}{a}$, an overdot denote differentials with
respect to comoving proper time, $t$.

Hence, the Einstein’s field equations (2) are

$3H^{2}=\frac{1}{2}\dot{\phi}^{2}+V(\phi),$ (6)
$2\dot{H}+3H^{2}=-\frac{1}{2}\dot{\phi}^{2}+V(\phi),$ (7)

while the Klein-Gordon equation reads

$\ddot{\phi}+3H\dot{\phi}+V_{,\phi}=0.$ (8)

We have assumed that field $\phi$ inherits the symmetries of the background
space, that is, $\phi=\phi\left(t\right)$. Inflation occurs, when the scalar
field potential dominates, $3H^{2}\simeq V\left(\phi\right)$ and the scalar
field approaches a stationary point
$\dot{\phi}\simeq-\frac{V_{,\phi}}{3H}~{}$l01 .

Thus in order to understand the existence of an inflationary era provided by a
potential function, the potential slow-roll parameters (PSR)

$\varepsilon_{V}=\left(\frac{V_{,\phi}}{2V}\right)^{2}~{}\,,~{}\eta_{V}=\frac{V_{,\phi\phi}}{2V},$
(9)

have been introduced l01 .

The condition an inflationary era to exists is expressed as,
$\varepsilon_{V}<<1$, while in order for the inflationary phase to last long
enough it is required require the second PSR parameter also to be small, i.e.
$\eta_{V}<<1$.

Some inflationary potentials which have been proposed in the literature are
presented. The quadratic $R+\left(\frac{R}{6M}\right)^{2}$ inflationary model
with potential function in the Einstein-frame sa1

$V_{A}\left(\phi\right)=V_{0}\left(1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M}}\right)^{2},$
(10)

the intermediate inflationary potential oll0

$V_{B}\left(\phi\right)=\frac{48A^{2}}{\left(\Delta+4\right)^{2}}\left(\frac{\phi-\phi_{0}}{\sqrt{2A\Delta}}\right)^{-\Delta}-\frac{4A\Delta}{\left(\Delta+4\right)^{2}}\left(\frac{\phi-\phi_{0}}{\sqrt{2A\Delta}}\right)^{-\Delta-2}~{},$
(11)

the natural inflation oll1

$V_{C}\left(\phi\right)=V_{0}\left(1+\cos\left(\frac{\phi}{f}\right)\right),$
(12)

the hyperbolic inflationary potential ref7

$V_{D}\left(\phi\right)=V_{0}\sinh^{q}\left(p\right)~{},$ (13)

and many others. For an extended list of the proposed inflationary models we
refer the reader in the reviews ec1 ; ec2 .

Similarly, with the PSR parameters we can define the corresponding Hubble
slow-roll parameters (HSR) as l01

$\varepsilon_{H}=-\frac{d\ln H}{d\ln a}=\left(\frac{H_{,\phi}}{H}\right)^{2},$
(14) $\eta_{H}=-\frac{d\ln H_{,\phi}}{d\ln a}=\frac{H_{,\phi\phi}}{H}.$ (15)

The two different set of parameters, the HSR and PSR parameters are related by
the expressions

$\varepsilon_{V}\simeq\varepsilon_{H}~{}\text{and~{}}\eta_{V}\simeq\varepsilon_{H}+\eta_{H}.$
(16)

Therefore, inflation occurs when $\varepsilon_{H}<<1$. The limit
$\varepsilon_{H}=1$, is called the end of inflation where we have the exit
from the acceleration era. In the following we work with the Hubble slow-roll
parameters.

As far as the observable values for the spectral indices $n_{s}$ and $r$, are
concerned, from the Planck 2018 collaboration pin2018  follow that the
spectral index for the density perturbations is constraint as $n_{s}=0.9649\pm
0.0042,$ while the tensor to scalar ratio, $r$ is constraint as $r<0.10$.

These indices are related with the HSR parameters in the first-order
approximation as

$r=10\varepsilon_{H}~{},$ (17) $n_{s}=1-4\varepsilon_{H}+2\eta_{H}.$

In the second-order approximation we shall consider the additional slow-roll
parameter $\xi_{H}\equiv\frac{H_{\phi}H_{\phi\phi\phi}}{H^{2}}$, such that
$n_{s}$ and $n_{s}^{\prime}$ are expressed as

$n_{s}\equiv
1-4\varepsilon_{H}+2\varepsilon_{H}-8\left(\varepsilon_{H}\right)^{2}\left(1+2C\right)+\varepsilon_{H}\eta_{H}\left(10C+6\right)-2C\xi_{H},$
(18) $n_{s}^{\prime}\equiv 2\varepsilon_{H}\eta_{H}-2\xi_{H}.$ (19)

where $C=\gamma_{E}+\ln 2-2\simeq-0.7296$ and the range of the scalar spectral
index is $n_{s}^{\prime}=-0.005\pm 0.013$ pin2018 .

In the following we consider the spectral indices in the first-order
approximation.

###  II.1 Algebraic solution

In ns1 , the generic algebraic solution for the field equations (6)-(8) is
presented for arbitrary potential function. Indeed, with the use of new
variables the unknown potential function can be introduced inside the metric
tensor, while the dynamical variables for the field equations are expressed in
terms of an arbitrary function. The field equations are reduced in algebraic
equation. Such a solution is called algebraic solution.

We define the new variable
$dt=\exp\left(\frac{F\left(\omega\right)}{2}\right)d\omega$ with $\omega=6\ln
a$. The FLRW line element becomes

$ds^{2}=-e^{F\left(\omega\right)}d\omega^{2}+e^{\omega/3}(dx^{2}+dy^{2}+dz^{2}).$
(20)

with Hubble function $H\left(\omega\right)=\frac{1}{6}e^{-\frac{F}{2}}$.

Hence, from the field equations (6)-(8) we find that then the scalar field is
expressed as ns1

$\phi(\omega)=\pm\frac{\sqrt{6}}{6}\int\\!\\!\sqrt{F^{\prime}(\omega)}d\omega$
(21)

while the scalar field potential reads

$V(\omega)=\frac{1}{12}e^{-F(\omega)}\left(1-F^{\prime}(\omega)\right).$ (22)

in which
$F^{\prime}\left(\omega\right)=\frac{dF\left(\omega\right)}{d\omega}$.

Function $F\left(\omega\right)$ can be constraint by various approaches. In
ref10  closed-form scalar field solutions were found by assuming specific
functional forms for the equation of state parameter
$P_{\phi}=P_{\phi}\left(\rho_{\phi}\right)$ of the scalar field,
$P_{\phi}=\frac{1}{2}\dot{\phi}^{2}-V\left(\phi\right)$,
$\rho_{\phi}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)$.

In the new set of variables, the HSR parameters are expressed as follows ref10

$\varepsilon_{H}=3F^{\prime}~{},~{}\eta_{H}=3\frac{\left(F^{\prime}\right)^{2}-F^{\prime\prime}}{F^{\prime}},$
(23)

while the number of e-folds is defined as

$N_{e}=\int_{t_{i}}^{t_{f}}H\left(t\right)dt=\ln\frac{a_{f}}{a_{i}}=\frac{1}{6}\left(\omega_{f}-\omega_{i}\right).$
(24)

In ref10  it was considered that the spectral indices are constraint as

$n_{s}-1=h\left(r\right),$ (25)

where $h\left(r\right)$ has been assumed to be
$h\left(r\right)=h_{1}+h_{2}r+h_{3}r^{2}.~{}$For this latter assumption, new
exact and analytic inflationary solutions were found.

In the first-order approximation, from (25) and (23) the master equation (25)
reads

$F^{\prime\prime}+G\left(F^{\prime}\right)=0,~{}G\left(F^{\prime}\right)=-\frac{1}{6}h\left(F^{\prime}\right)F^{\prime}+\left(F^{\prime}\right)^{2}.$
(26)

In the following we constrain function $G\left(F^{\prime}\right)$, i.e.
$h\left(F^{\prime}\right)$, by the requirement the master equation (26) to
admit Lie point symmetries. This selection rule has geometric characteristics.
Symmetries are geometric objects in the space where the differential equation
is defined. Such a geometric selection rule is in agreement with the geometric
nature of gravitational physics, for a discussion we refer to sym1 .

We proceed with the presentation of the symmetry analysis for the master
equation (26).

##  III Lie symmetry analysis

The theory of symmetries of differential equations established at the end of
the 19th century by Sophus Lie kumei . The novelty of Lie’s approach is that
the transformations group which leave invariant a differential equation can be
used to simplify the given equation. In the case of ordinary differential
equations the Lie point symmetries are used to reduce the order of the
differential equation.

Let us demonstrate this on our master equation (26). By assuming the new
variables $f=F^{\prime}$, the master equation (26) can be written in the
equivalent form $f^{\prime}=G\left(f\right)$, that is
$\int\frac{df}{G}=\omega-\omega_{0}$. This is the most common reduction
process which holds for autonomous dynamical systems where the autonomous
point symmetry vector $\partial_{\omega}$ exists.

The derivation of Lie (point) symmetries for a given dynamical system is based
on a simple algorithm. We briefly discuss the algorithm.

Consider the infinitesimal point transformation

$\displaystyle\bar{\omega}$
$\displaystyle=\omega+\varepsilon\xi\left(\omega,F\right)~{},~{}$ (27)
$\displaystyle\bar{F}$
$\displaystyle=F+\varepsilon\eta\left(\omega,F\right)~{},$ (28)

with generator the vector field
$X=\xi\left(\omega,F\right)\partial_{\omega}+\eta\left(\omega,F\right)\partial_{F}$.
Then, the master equation (26) is invariant under the Action of the later
point transformation if and only if

$X^{\left[2\right]}\left(F^{\prime\prime}-G\left(F^{\prime}\right)\right)=0$
(29)

in which
$X^{\left[2\right]}=X+\eta^{\left[1\right]}\partial_{F^{\prime}}+\eta^{\left[2\right]}\partial_{F^{\prime\prime}}$
is the second extension of $X$ in the tangent space, where
$\eta^{\left[1\right]}$ and $\eta^{\left[2\right]}$ are defined as

$\displaystyle\eta^{\left[1\right]}$
$\displaystyle=\eta^{\prime}-F^{\prime}\xi^{\prime},$ (30)
$\displaystyle\eta^{\left[2\right]}$
$\displaystyle=\eta^{\left[1\right]^{\prime}}-F^{\prime\prime}\xi^{\prime}.$
(31)

When for a given point transformation the symmetry condition (29) is valid,
the vector field $X$ will be called a Lie symmetry.

The symmetry condition (29) for the master equation (26) is expressed as
follow

$\displaystyle 0$
$\displaystyle=\eta_{,\omega\omega}-F^{\prime}\xi_{,\omega\omega}+\left(F^{\prime}\right)^{2}\left(\eta_{,FF}-2\xi_{,\omega
F}\right)-\left(F^{\prime}\right)^{3}\xi_{,FF}+2F^{\prime}\eta_{,\omega
F}+G_{,F^{\prime}}\eta_{,\omega}+$
$\displaystyle+\left(G_{,F^{\prime}}F^{\prime}-G\right)\eta_{,F}+\left(3GF^{\prime}-G_{,F^{\prime}}\left(F^{\prime}\right)^{2}\right)\xi_{,F}+\left(2G-G_{,F^{\prime}}F^{\prime}\right)\xi_{,\omega}.$
(32)

From the later equation we can define systems of partial differential
equations which shall constraint the unknown functions
$\xi\left(\omega,F\right)$ and $\eta\left(\omega,F\right)$. However the
constraint system it depends on the functional form of
$G\left(F^{\prime}\right)$. Indeed, for arbitrary function
$G\left(F^{\prime}\right)$ we find the generic solution
$\xi\left(\omega,F\right)=\alpha_{1}$ and
$\eta\left(\omega,F\right)=\alpha_{2}$, which means that the generic symmetry
vector is $X=\alpha_{1}\partial_{\omega}+\alpha_{2}\partial_{F}$. Hence, the
two independent symmetries, $X_{1}=\partial_{\omega}$ and
$X_{2}=\partial_{F}$, follow, with commutator $\left[X_{1},X_{2}\right]=0.$

However, there are special forms of $G\left(F^{\prime}\right)$ in which the
master equation (26) admits additional Lie symmetries. This classification
problem is that we attempt to solve in this study. It is important to mention
that the Lie point symmetries for second-order differential equations have
been widely studied in the literature see for instance cs1 ; cs2 .

It is well known that a second-order differential equation can admits $0,$
$1,~{}2,~{}3$ or $8$ Lie symmetries. When eight Lie point symmetries exist,
the the differential equation is called maximal symmetric and there exist a
similarity transformation such that the equation to be written in the form of
the free particle.

For arbitrary function $G\left(F^{\prime}\right)$ the master equation (26)
admits two Lie symmetries, thus we investigate for specific forms of
$G\left(F^{\prime}\right)$ in which the master equation admits $3$ or $8$ Lie
symmetries.

Thus, the classification scheme provides for equation (26) gives the three
functional forms
$G_{A}\left(F^{\prime}\right)=G_{0}\left(F^{\prime}\right)^{\nu+1},~{}\nu\neq-1,0,1,2$,
$G_{B}\left(F^{\prime}\right)=G_{0}\exp\left(\mu F^{\prime}\right),~{}\mu\neq
0$ and
$G_{C}\left(F^{\prime}\right)=G_{3}\left(F^{\prime}\right)^{3}+G_{2}\left(F^{\prime}\right)^{2}+G_{1}F^{\prime}+G_{0}$,
where $G_{0-3}$ and $\mu~{}$are constants parameters

For $G_{A}\left(F^{\prime}\right)$, equation (26) admits the three Lie
symmetries
$\left\\{X_{1},X_{2},X_{3}=\nu\omega\partial_{\omega}+\left(\nu-1\right)F\partial_{F}\right\\}$
with commutators $\left[X_{1},X_{2}\right]=0,~{}\left[X_{1},X_{2}\right]=\nu
X_{1}$ and $\left[X_{2},X_{3}\right]=\left(\nu-1\right)X_{2}$.

For the exponential function $G_{B}\left(F^{\prime}\right)$, the Lie
symmetries of equation (26) are
$\left\\{X_{1},X_{2},\bar{X}_{3}=\mu\omega\partial_{x}+\left(F\mu-\omega\right)\partial_{F}\right\\}$
with commutators $\left[X_{1},\bar{X}_{3}\right]=\mu X_{1}-X_{2}$ and
$\left[X_{2},\bar{X}_{3}\right]=\mu X_{2}$.

Finally, for $G_{C}\left(F^{\prime}\right)$ equation (26) is maximal symmetric
and admits eight Lie symmetries. The representation of the admitted Lie
algebra depends on the values of the free parameters thus we omit it. The
later case has been widely studied before in ref11 . Moreover, the similarity
transformations which connect the different inflationary models which belong
to the family of $G_{C}\left(F^{\prime}\right)$ are presented in sc1 .

For the inflationary models provided by $G_{A}\left(F^{\prime}\right)$ and
$G_{B}\left(F^{\prime}\right)$ we derive the corresponding functions
$h\left(r\right)$ are
$h_{A}\left(r\right)=-6\left(h_{0}r^{\nu}-\frac{1}{10}r\right)$ and
$h_{B}\left(r\right)=-6\left(h_{0}r^{-1}e^{\mu r}-\frac{1}{10}r\right)$.

From the cosmological observations we know that $n_{s}-1\simeq 0$, while
$r<0.11$. Thus, from these two models, $h_{A}\left(r\right)$ and
$h_{B}\left(r\right)$, only model $h_{A}\left(r\right)$ provides a behaviour
$\lim_{r\rightarrow 0}\left(h_{A}\right)\simeq 0$ for $\nu>0.~{}$On the other
hand, $h_{A}\left(r\right)<<1$ when $r<<1$, and for $\nu<-1$ if and only if
$h_{0}\simeq\frac{1}{r^{\nu+1}}$, that is, for very large values of the free
parameter $h_{0}$.

In the following section we focus our analysis on model $h_{A}\left(r\right)$,
where we investigate the closed-form solution, we discuss the physical
properties of the model and we investigate the evolution for the inflationary
parameters.

##  IV New inflationary exact solution

For $G=G_{A}\left(F^{\prime}\right)$, the master equation reads

$F^{\prime\prime}+G_{0}\left(F^{\prime}\right)^{\nu+1}=0,$ (33)

with closed-form solution

$F\left(\omega\right)=\frac{\left(G_{0}\nu\omega\right)^{1-\frac{1}{\nu}}}{G_{0}\left(\nu-1\right)}\text{.}$
(34)

Thus, the physical parameters for the inflaton field, such as the energy
density $\rho_{\phi}$, the pressure $P_{\phi}$ and the equation of state
parameter $w_{\phi}=\frac{P_{\phi}}{\rho_{\phi}}$ can be constructed
analytical

$\rho_{\phi}\left(\omega\right)=\frac{1}{12}\exp\left(\frac{\left(G_{0}\nu\omega\right)^{1-\frac{1}{\nu}}}{G_{0}\left(\nu-1\right)}\right)~{},$
(35)
$P_{\phi}\left(\omega\right)=\frac{1}{12}\exp\left(\frac{\left(G_{0}\nu\omega\right)^{1-\frac{1}{\nu}}}{G_{0}\left(\nu-1\right)}\right)\left(-1+2\left(G_{0}\nu\omega\right)^{-\frac{1}{\nu}}\right)~{},$
(36)

and

$w_{\phi}\left(\omega\right)=\left(-1+2\left(G_{0}\nu\omega\right)^{-\frac{1}{\nu}}\right).$
(37)

We remark that for $\nu>0$ and for large values of
$\omega,~{}w_{\phi}\left(\omega\right)\rightarrow-1$, otherwise when $\nu<0$,
$w_{\phi}\left(\omega\right)\rightarrow-1$ for small values of $\omega$.

As far as the scalar field and the scalar field potential are concerned, we
derive the following expressions

$\phi\left(\omega\right)=\sqrt{\frac{2}{3}}\frac{\left(G_{0}\nu\omega\right)^{1-\frac{1}{2\nu}}}{G\left(2\nu-1\right)}~{},~{}\nu\neq\frac{1}{2},$
(38)
$\phi\left(\omega\right)=\sqrt{\frac{2}{3}}\frac{\ln\omega}{G_{0}}~{},~{}\nu=\frac{1}{2}~{},$
(39)

and

$V\left(\omega\right)=\frac{1}{12}\left(1-\left(G_{0}\nu\omega\right)^{-\frac{1}{\nu}}\right)\exp\left(\frac{\left(G_{0}\nu\omega\right)^{1-\frac{1}{\nu}}}{G_{0}\left(\nu-1\right)}\right)\text{.}$
(40)

In terms of the scalar filed $\phi$, the potential function
$V\left(\phi\right)$ is given by the following functional form

$V\left(\phi\right)=\frac{1}{12}\left(1-\left(\frac{3}{2}\right)^{\frac{1}{1-2\nu}}\left(\left(2\nu-1\right)G_{0}\phi\right)^{\frac{2}{1-2\nu}}\right)\exp\left(\frac{\left(\frac{2}{3}\right)^{\frac{1-\nu}{2\nu-1}}}{G_{0}\left(\nu-1\right)}\left(\left(2\nu-1\right)G_{0}\phi\right)^{1+\frac{1}{2\nu-1}}\right)~{},~{}\nu\neq\frac{1}{2}$
(41)

or

$V\left(\phi\right)=\frac{1}{12}\left(1-\frac{4}{G_{0}}\exp\left(-\sqrt{6}G_{0}\phi\right)\right)\exp\left(-\frac{4}{G_{0}}\exp\left(-\sqrt{\frac{3}{2}}G_{0}\phi\right)\right)~{},~{}\nu=\frac{1}{2}\text{.}$
(42)

In Fig. 1 we present the qualitative evolution for the scalar field potential
(41) for $G_{0}=1$ and for different values of the free parameter $\nu$.

Figure 1: Qualitative evolution for the scalar field potential (41) for
$G_{0}=1$. Solid line is for $\nu=\frac{3}{2}$, dotted line is for
$\nu=\frac{5}{2}$ and dashed line is for $\nu=4$.

According to our knowledge, this inflationary potential has not been presented
before in the literature.

###  IV.1 Spectral indices

For the closed-form solution of $F\left(\omega\right)$, the definition of the
slow-roll parameters from expressions (23) we can define the spectral indices
$n_{s}$ and $r$ in terms of the number of e-fold $N_{e}$. Recall, that the end
of inflation occurs when $\varepsilon\left(\omega_{f}\right)=1$, that is,
$\omega_{f}=\frac{3^{\nu}}{G_{0}\nu}$.

Consequently, the spectral indices by using the HSR parameters are written as
follows

$n_{s}=1-6\left(\left(G_{0}\nu\right)\left(6N_{e}+\frac{3^{\nu}}{G_{0}\nu}\right)\right)^{-\frac{1}{\nu}}+\frac{6G_{0}}{3^{\nu}+6G_{0}N_{e}\nu}~{},$
(43)
$r=30\left(\left(G_{0}\nu\right)\left(6N_{e}+\frac{3^{\nu}}{G_{0}\nu}\right)\right)^{-\frac{1}{\nu}}.$
(44)

We observe that the indices $n_{s}$ and $r$ depend on the number of e-fold
$N_{e}$, and on the free parameters $\left(\nu,G_{0}\right).~{}$From the
observations the number of e-fold it is considered to be in the range
$N_{e}=\left(50,60\right)$.

In Figs. 2, 3 and in 4 we present the qualitative evolution of the spectral
indices according to the different values of the free parameters. In Fig. 2 we
present the contour plot for the spectral indices $n_{s}$ and $r~{}$ in the
first-order approximation, in the space of the free parameter
$\left(\nu,G_{0}\right)$ and for $N_{e}=50$, $N_{e}=55$ and $N_{e}=60$. we
observe that as power $\nu$ increases then the model fits the observations as
$G_{0}$ increases exponentially. In Figs. 3 and 4 we present the parametric
plot in the space $\left(n_{s},r\right)$ for varying parameter the number of
e-fold $N_{e}$, and for specific values of the free parameters of the model
$\left(\nu,G_{0}\right)$. In Fig. 4 we assume the that
$G_{0}=\bar{G}_{0}\nu^{-1}$ and plots are for specific values of $\nu$.

From the results of Figs. 2, 3 and in 4 it is clear that the new inflationary
model provided by the symmetry analysis provide values for the spectral
indices according to the cosmological observations, see Section II, for the
values of the Planck 2018 collaboration, Table 3 in pin1 .

Figure 2: Qualitative evolution for the spectral indices
$\left(n_{s},r\right)$ in the first approximation in the two-dimensional space
for the free variables $\left(\nu,G_{0}\right)$, for $N_{e}=50$ (left figs.)
$N_{e}=55~{}$(center figs.) and $N_{e}=60~{}$(right figs.) Figure 3:
Parametric plots for the spectral indices $\left(n_{s},r\right)$ in the first
approximation varying number of e-folds $N_{e}=\left(50,55,60\right)$ and
different values of $G_{0}$ and $\nu=\nu_{0}+\varepsilon$. Left fig. is for
$\left(G_{0},\nu_{0}\right)=\left(10,1.5\right)$, center fig. is for
$\left(G_{0},\nu_{0}\right)=\left(10^{2},2\right)$ and right fig. is for
$\left(G_{0},\nu_{0}\right)=\left(10^{3},3\right)$. Solid lines are for
$\varepsilon=0.01,~{}$dashed lines are for $\varepsilon=0.02$ and dotted lines
are for $\varepsilon=0.03$.  Figure 4: Parametric plots for the spectral
indices $\left(n_{s},r\right)$ in the first approximation varying number $\nu$
and different values of $G_{0}$ and $N_{e}$. Left fig. is for
$\left(G_{0},N_{e}\right)=\left(1/\nu,\left\\{50,55,60\right\\}\right)$,
center fig. is for
$\left(G_{0},N_{e}\right)=\left(10^{3}/\nu,\left\\{50,55,60\right\\}\right)$
and right fig. is for
$\left(G_{0},N_{e}\right)=\left(10^{4}/\nu,\left\\{50,55,60\right\\}\right)$.
Solid lines are for $N_{e}=50,~{}$dashed lines are for $N_{e}=55$ and dotted
lines are for $N_{e}=60$. Left fig. is for $\nu=\left(1.01,1.2\right)$, center
fig. is for $\nu=\left(2.01,2.5\right)$, right fig. is for
$\nu=\left(2.5,3\right)$.

##  V Conclusions

In this work we investigated the construction of a new exact inflationary
models in the inflaton theory by using Lie symmetry analysis for the master
equation (26). The master equation of our consideration is constructed by the
assumption the spectral index $n_{s}$ and the tensor to scalar ratio $r$, to
be related by a function $G$. From this assumption the second-order
differential equation (26) follows.

We performed a classification for the function which relates $n_{s}$ and $r$
by assuming the master equation to admit Lie symmetries. By applying the
symmetry condition we found that the master equation admits two Lie symmetries
for arbitrary function $G$, while for $G=G_{A}$ and $G=G_{B}$ admits three Lie
symmetries and for $G=G_{C}$ the admitted Lie symmetries form the
$sl\left(3,R\right)$ Lie algebra and the master equation is maximally
symmetric. The case $G_{C}$ was investigated in details in a previous study
ref11 , hence in this work we focused on $G_{A}$ and $G_{B}$.

For $G_{A}$ and $G_{B}$, the master equation (26) admits three Lie symmetries,
which form the Lie algebras $A_{3,3}$ and $A_{3,2}$ in the Patera et al.
classification pat1 , respectively. However function $G_{B}$ does not provide
a relation for the indices $n_{s}$,$~{}r$ as provided by the observations,
thus, we focused our analysis on the function $G_{A}$ and the Lie algebra
$A_{3,3}$.

We solved the master equation and we were able to derive the closed-form
solution for the scalar field and for all the physical variables of the
inflaton model. Furthermore, we wrote the closed-form expressions for the
spectral indices $n_{s}$, $r$ , in the first-order approximation by using the
HSR parameters. We found that the spectral indices depend on the number of
e-fold and on two free parameters. By presenting the qualitative evolution of
the spectral indices we remark that they can fit in the cosmological
observations.

In Fig. 5, we present the qualitative evolution for the spectral indices
$n_{s}$ and $n_{s}^{\prime}$ in the second approximation for different values
of the free variables $\left(\nu,G_{0}\right)$. We observe that
$n_{s}^{\prime}$ is always positive valued, and takes small values when
$n_{s}-1\simeq 0$. According the cosmological observations, the running index
is constraint $\eta_{s}^{\prime}=-0.005\pm 0.013$, therefore, small positive
values are inside the $1\sigma$ region.

Figure 5: Qualitative evolution for the spectral indices
$\left(n_{s},n_{s}^{\prime}\right)$ in the second approximation in the two-
dimensional space for the free variables $\left(\nu,G_{0}\right)$, for
$N_{e}=50$ (left figs.) $N_{e}=55~{}$(center figs.) and $N_{e}=60~{}$(right
figs.)

We conclude that the Lie symmetry approach is a powerful method for the
derivation of exact solutions and in this study the symmetry method applied
for the derivation of a new inflationary exact solution. A natural question
which raised is if this approach can be applied in $f\left(R\right)$-theory or
in scalar-tensor gravity ff3 . These three theories are related under
conformal transformations. The analytical solution determined in ns1  for the
scalar field can be easily used to write the corresponding analytic solutions
for the conformal equivalent theories, and extend the present symmetry
analysis in the conformal frame. However, such an analysis overpass the scopus
of the present work and will be discussed elsewhere.

In a future study we plan to investigate further this scalar field solution
model and we want to use this model as dark energy candidate and investigate
if it solves the $H_{0}$-tension h01 .

###### Acknowledgements.

This work is based on the research supported in part by the National Research
Foundation of South Africa (Grant Number 131604).

## References

  * (1) A. Guth, Phys. Rev. D 23, 347 (1981)
  * (2) V. Muller, H.-J. Schmidt and A.A. Starobinsky, Phys. Lett. B 202, 2, 198 (1988)
  * (3) L.A Kofman, A.D. Linde and A.A. Starobinsky, Phys. Lett. B 157, 5-6, 361 (1985)
  * (4) R. Wald, Phys. Rev. D 28, 2118 (1983)
  * (5) G.W. Gibbons and S.W Hawking, Phys. Rev. D 15, 2738 (1977)
  * (6) S.W. Hawking and J.G. Moss. Phys. Lett. B 110, 35 (1982)
  * (7) K. Sato, MNRAS 195, 467 (1981)
  * (8) J.D Barrow and A. Ottewill, J. Phys. A 16, 2757 (1983)
  * (9) A.D. Linde, Phys. Lett. B 108, 389 (1982)
  * (10) A.D. Linde, Phys. Lett. B 129, 177 (1983)
  * (11) A.R. Liddle, Phys. Lett. B 220, 502 (1989)
  * (12) T. Charters, J.P. Mimoso and A. Nunes, Phys. Lett. B 472, 21 (2000)
  * (13) J.D. Barrow, Phys. Rev. D 48, 1585 (1993)
  * (14) J.D. Barrow and P. Saich, Class. Quantum Grav. 10, 279 (1993)
  * (15) J.D. Barrow and P. Parsons, Phys. Rev. D 52, 5576 (1995)
  * (16) N. Deruelle, C. Gundlach and D. Langlois, Phys. Rev. D 46, 5337 (1992)
  * (17) H. Sheikhahmadi, A. Mohammadi, A. Aghamohammadi, T. Harko, R. Herrera, C. Corda, A. Abebe and K. Saaidi, EPJC 79, 1038 (2019)
  * (18) R. Kallosh and A. Linde, JCAP 13, 027 (2013)
  * (19) S. Basilakos and J.D. Barrow, Phys. Rev. D 91, 103517 (2015)
  * (20) J. de Haro, J. Amorós and S. Pan, Phys. Rev. D 93, 084018 (2016)
  * (21) L.A. Urena-Lopez, J. Phys. Conf. 761, 012076 (2016)
  * (22) J.D. Barrow and A. Paliathanasis, Phys. Rev. D 94, 083518 (2016)
  * (23) J.D. Barrow and A. Paliathanasis, Gen. Rel. Gravit. 50, 82 (2018)
  * (24) A.A. Starobisnky, Phys. Lett. B 91, 99 (1980)
  * (25) T.P. Sotiriou, Gravity and Scalar Fields. In: Papantonopoulos E. (eds) Modifications of Einstein’s Theory of Gravity at Large Distances. Lecture Notes in Physics, vol 892, Springer (2015)
  * (26) Planck Collaboration: Y. Akrami et al. A&A 641, A10 (2020)
  * (27) J. Gacia-Bellido, Nucl. Phys. B 423, 221 (1994)
  * (28) Y. Gong, Phys. Rev. D 59, 083507 (1999)
  * (29) A. Beesham, S.V. Chervon and S.D. Maharaj, Class. Quantum Grav. 26, 075017 (2009)
  * (30) M. Yamaguchi, Class. Quantum Gravit. 28, 103001 (2011)
  * (31) J. Ellis and N.E. Mavromatos, Phys. Rev. D 88, 085029 (2013)
  * (32) A. Kehagias, A.M. Dizgah and A. Riotto, Phys. Rev. D 89, 043527 (2014)
  * (33) S. Hellerman, J. Kehayias and T.T. Yanagida, Phys. Lett. B 742, 390 (2015)
  * (34) D. Wands, Lect. Notes Phys. 738, 275 (200)
  * (35) S.D. Odintsov, V.K. Oikonomou and F.P. Fronimos, Nucl. Phys. B 963, 115299 (2021)
  * (36) S.D. Odintsov and V.K. Oikonomou, Phys. Rev. D 98, 044039 (2018)
  * (37) D. Benisty, E.I. Guendelman and E.N. Saridakis, EPJC 80, 480 (202)
  * (38) J. Martin, C. Ringeval and V. Vennin, Phys. Dark Univ. 5-6, 75 (2014)
  * (39) J. Martin, C. Ringeval, R. Trotta and V. Vennin, JCAP 03, 039 (2014)
  * (40) G.F.R. Ellis and M.S. Madsen, Class. Quantum Grav. 8, 667 (1991)
  * (41) J.G. Russo, Phys. Lett. B 600, 185 (2004)
  * (42) S. Basilakos, M. Tsamparlis and A. Paliathanasis, Phys. Rev. D 83, 103512 (2011)
  * (43) A. Paliathanasis, M. Tsamparlis, S. Basilakos and J.D. Barrow, Phys. Rev. D 91, 123535 (2015)
  * (44) N. Dimakis, A. Karagiorgos, A. Zampeli, A. Paliathanasis, T. Christodoulakis and P.A. Terzis, Phys. Rev. D 93, 123518 (2016)
  * (45) L. V. Ovsiannikov, Group analysis of differential equations, Academic Press, New York, (1982)
  * (46) S.D. Odintsov and V.K. Oikonomou, Nucl. Phys. B 929, 79 (2018)
  * (47) S.D. Odintsov and V.K. Oikonomou, Annals of Physics 388, 267 (2018)
  * (48) A.R. Liddle, P. Parson and J.D. Barrow, Formalising the Slow-Roll Approximation in Inflation, Phys. Rev. D 50, 7222 (1994)
  * (49) J.D. Barrow, Phys. Lett. B 235, 40 (1990)
  * (50) K. Freese, J.A. Frieman and A.V. Olinto, Phys. Rev. Lett. 65, 3233 (1990)
  * (51) Planck Collaboration: Y. Akrami et al. Planck 2018 results. X. Constraints on inflation, A&A 641, A10 (2020)
  * (52) G.W. Bluman and S. Kumei, Symmetries and Differential Equations, Springer-Verlag, New York, (1989)
  * (53) F.M. Mahomed, Symmetry group classification of ordinary differential equations: Survey of some results, Math. Meth. Appl. Sci. 30, 1995 (2007)
  * (54) V. Dorodnitsyn, R. Kozlov and P. Winternitz, Lie group classification of second-order ordinary difference equations, J. Math. Phys. 41, 480 (2000)
  * (55) A. Paliathanasis, Symmetry Analysis In Inflationary Cosmology, AIP Conference Proceedings 2153, 020015 (2019)
  * (56) J. Patera, R.T. Sharp and P. Winternitz, Invariants of real low dimension Lie algebras, J. Math. Phys. 17, 986 (1976)
  * (57) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D.F. Mota, A.G. Riess and J. Silk, In the Realm of the Hubble tension - a Review of Solutions, Class. Quantum Grav. 38, 153001 (2021)
  * (58) S. Nojiri, S.D. Odintsov and V.K. Oikonomou, Phys. Rept. 692, 1 (2017)


# Itinerant magnetism and magnetic polarons in the triangular lattice Hubbard
model

Ivan Morera  Institute for Theoretical Physics, ETH Zurich, Wolfgang-Pauli-
Str. 27, 8093 Zurich, Switzerland  Eugene Demler  Institute for Theoretical
Physics, ETH Zurich, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland

###### Abstract

We use density matrix renormalization group to investigate the phase diagram
of the Fermi Hubbard model on a triangular lattice with densities above half-
filling, $1\leq n<2$. We discuss the important role of kinetic magnetism and
magnetic polarons. For strong interactions and low doublon dopings, attractive
interaction between polarons results in phase separation between the fully
polarized state at finite doping and the commensurate spin-density wave state
at half-filling. For intermediate interaction strength and small doping,
competition between antiferromagnetic superexchange and kinetic magnetism
gives rise to the incommensurate spin density wave (I-SDW) phase. Fully
polarized ferromagnetic (FPF) phase for weak interactions is limited to
dopings close to the van Hove singularity in the density of states. With
increasing interactions the FPF phase expands to lower dopings. For strong
interactions it reaches the low doping regime and is better understood as
arising from proliferation of Nagaoka-type ferromagnetic polarons. Other
phases that we find include a partially polarized phase, another type of I-SDW
at high densities, and Müller-Hartmann ferromagnetism close to the band
insulating regime.

Introduction. Recent experimental studies of electrons in van der Waals moiré
materials revealed unusual character of magnetism in triangular lattices.
Experiments with heterobilayer transition metal dichalcogenide (TMDc)
compounds indicated that magnetic correlations are absent in the Mott
insulating phase at half-filling ($n=1$) but appear upon doping [1, 2, 3].
This can be contrasted to square lattices, in which strongest magnetic order
is observed in the Mott state at $n=1$, while doping tends to suppress it.
These findings provide motivation for theoretical analysis of magnetic phases
in the Hubbard model on a triangular lattice with finite doping (for
theoretical studies of magnetism at filling $n=1$, see Refs. [4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]). Prior
work suggested several unusual phenomena upon doping the Mott insulator state
[27, 28, 29, 30, 31, 32, 33, 34, 35]. Delocalization of holes strengthens the
120∘ Neel order [28, 36], whereas doublons tend to favor ferromagnetic
correlations in accordance with Nagaoka-Thouless theorem [37, 38]. Other types
of itinerant ferromagnetism are also expected to be relevant in the triangular
lattice such as Mielke’s flat-band ferromagnetism [39, 40, 41] or Müller-
Hartmann ferromagnetism [42] and the presence of a van Hove singularity
suggests the appearance of a Stoner instability at weak coupling [43, 5, 44,
45].

Magnetic polarons play a crucial role for understanding the emergence of
magnetic correlations in a wide range of materials [46]. In seminal works by
Nagaev [47, 48, 49], it was proposed that conduction electrons doped onto
antiferromagnetic semiconductors would generate significant ferromagnetic
regions around them, the so-called “ferrons”, wherein the electrons would
become self-trapped. Similar phenomena were subsequently discussed for dilute
doped magnetic semiconductors, such as Ga1-xMnxAs and Ge1-xMnx. In these
materials random potential constrains the motion of doped holes to regions of
finite size and ferromagnetic bubbles form around each dopant. With increasing
dopant density, such bubbles begin to overlap, and ferromagnetic transition
can be understood as a percolation transition of localized magnetic polarons
[50]. Magnetic polarons have also been investigated within the Hubbard model
[51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 65, 66, 67, 68,
69, 70], and the emergence of ferromagnetic polarons, commonly referred to as
Nagaoka’s polarons, has also been examined in square lattices [62]. In this
case Nagaoka polarons emerge only in the very-strong coupling regime
$U/t\gtrsim 200$.

Figure 1: Schematic magnetic phase diagram of the Hubbard model in the
triangular lattice in (a). The top (bottom) panel represents results obtained
with the $t-J$ (Hubbard) model. We encounter a commensurate spin density wave
(C-SDW) corresponding to 120∘ order (green regions), a phase separated (PS)
regime (blue-to-white regions) where the system spatially phase separates in a
ferromagnetic doublon-rich region and an undoped region with C-SDW order, a
fully polarized (FP) state (full blue regions), a partially polarized (PP)
state (light blue region), and a flat-band ferromagnet (FBF) state (red
regions). We encounter three critical points (CP1, CP2 and CP3) where
different critical lines merge. In (b) we show the spin structure factors
corresponding to the different phases.

Recent experimental progress of quantum simulators made it possible to observe
magnetic polarons directly. In quantum dot systems, Nagaoka-Thouless
ferromagnetism has been observed in a small $2\times 2$ plaquette [71].
Experiments with cold atoms in triangular optical lattices have demonstrated
the tendency of the system towards ferromagnetism (antiferromagnetism) upon
particle (hole) doping of the Mott insulating state [72]. Additionally, two
different cold atom experiments have reported direct imaging of
antiferromagnetic (Haerter-Shastry’s) and ferromagnetic (Nagaoka’s) polarons
in the triangular lattice [73, 74]. These recent experimental advances provide
strong motivation to explore theoretically the phase diagram of the doped
Fermi Hubbard model. In this paper, we focus on fermion densities higher than
1, for which we find interesting competition between kinetic ferromagentism
driven by charge carriers (doublons) and antiferromagnetic exchange
interactions. While previous theoretical studies have explored the emergence
of ferromagnetism in the triangular lattice with variational wavefunctions
[43, 75, 5, 76, 77], a full picture of the phase diagram is still unclear.

Itinerant magnetism in the triangular lattice Hubbard model. In this work we
study the Hubbard model in the triangular lattice

$\displaystyle\hat{H}$ $\displaystyle=-t\sum_{\langle
i,j\rangle,\sigma}\left(\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{j,\sigma}+\textrm{h.c.}\right)+U\sum_{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow},$
(1)

where $\hat{c}_{i,\sigma}^{\dagger}$ ($\hat{c}_{i,\sigma}$) creates (destroys)
a fermion with spin $\sigma$ at site $i$,
$\hat{n}_{i,\sigma}=\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$
corresponds to the particle number operator, $\langle i,j\rangle$ denotes
nearest-neighbor sites, $t>0$ is the hopping strength and $U$ is the on-site
interaction strength.

Figure 2: In (a) we represent the total spin $S$ as a function of the
superexchange interaction $J/t$ for a $15\times 4$ system with two doublons.
Continuous black line corresponds to Eq. (3). In the inset panel we show the
doublon density profile for $J/t=0.06$ and $N_{D}=2$. Two doublons bunch
together in space denoting an effective attraction. In (b) we show the energy
dependence with the total magnetization for two doublons and $J/t=1/20,1/50$.
The dashed line denotes the total spin of the ground state for $S^{z}=0$. In
(c) we show the polarization energy $\Delta$ as a function of the doping
$\delta$ for two different interaction strengths. The continuous black line
denotes the linear dependence $\Delta/(tN_{s})\sim 4t/U-1.15\delta$.

To obtain the ground state of Eq. (1) we perform density matrix
renormalization group (DMRG) simulations implementing the U(1)$\times$U(1)
symmetry in YC4 and YC6 cylinders [21] with a maximum number of sites
$N_{s}=60$ and a maximum bond dimension $\chi\sim 6000$. Moreover, to study
the strongly interacting regime $U/t\gtrsim 40$ we perform simulations of the
effective $t-J$ model with a superexchange interaction $J=4t^{2}/U$ and
neglect three-site terms.

In Fig. 1 we present the magnetic phase diagram of the Hubbard model in the
triangular lattice as a function of doping $\delta=n-1$ and on-site
interaction $U/t$. For small hole doping ($\delta<0$) and strong interactions,
we observed enhancement of the 120∘ commensurate spin density wave order
(C-SDW). The mechanism of this enhancement has been discussed elsewhere [34].
For small doublon doping ($\delta>0$) and for $U/t>40$, each doped doublon
develops a Nagaoka-type ferromagnetic bubble around itself. In the ground
state we do not find a homogeneous Nagaoka’s polaron gas state, but phase
separation between the doublon rich fully polarized state and a half-filled
C-SDW state. In the phase separated regime, the system as a whole exhibits
partial polarization and the total spin increases with the doublon density.
For weaker interactions ($U/t<40$), we observe that a small doublon doping
results in the formation of an incommensurate spin density wave (I-SDW). In
this phase, peaks of the spin structure factor $\chi(q)$, which is defined in
Eq. (2), are shifted from the $K$ points towards zero momentum by an amount
that increases monotonically with increasing doublon density. We will refer to
momenta corresponding to the peaks of the spin structure factor as $q^{*}$. We
point out two similarities between the I-SDW phase and the stripe phase of
cuprates: 1) when there are incommensurate peaks, we no longer see
commensurate peaks at $K$; 2) for small doping we have linear scaling
$|q^{*}-K|$ with doping $\delta$ (similar to “Yamada plot” in cuprates [78]).
However, we do not observe sharp peaks in the density structure factor as
expected for a stripe phase. Therefore, we conclude that the I-SDW state is a
spiral phase with no clear density structure. In Fig. 1 this phase corresponds
to point (II).

At a critical doublon doping, we obtain a transition towards a fully polarized
(FP) state for all values of $U/t$ considered, see (III) in Fig. 1.
Additionally, the critical doublon doping decreases with increasing on-site
repulsion and it vanishes in the limit $U/t\rightarrow\infty$. A robust FP
state appears in this limit for a wide range of doublon dopings
$0<\delta\lesssim 0.7$ in our widest cylinder. Furthermore, we observe that
ferromagnetism is centered around doping $\delta\sim 1/2$ in the intermediate
interaction regime $U/t\sim 20$ signaling the relevance of the van Hove
singularity and the possible connection to a Stoner instability. The
transition towards the FP state can be preceded by an I-SDW, a PS state or a
partially polarized (PP) phase depending on the strength of the on-site
interaction. The PP phase precedes the FP one by lowering the on-site
interaction. It is characterized by a peak in the spin structure factor
centered at zero momentum and a finite non-saturated magnetization.

We identify three critical points (CPs) in the phase diagram in which two
critical lines merge, see CP1, CP2 and CP3 in Fig. 1. At infinitesimal doublon
doping we obtain a CP (CP1) characterizing the transition from Nagaoka’s
polaron to Brinkman-Rice-Shraiman-Siggia type polaron [52, 55, 62]. At higher
doublon dopings we encounter two different critical points (CP2 and CP3) which
signal that the PP and PS states cannot be smoothly connected. Note that in
each of the critical points CP2 and CP3 three different phases meet. We cannot
rule out the possibility that with increasing width of the cylinders, critical
points CP2 and CP3 start approaching each other and in the two-dimensional
system we will have a single critical point in which four different phases
meet.

Figure 3: Normalized spin structure factor $S(\mathbf{q})$ at null momentum in
(a), at momentum $\mathbf{q}^{*}$ in (b) and magnitude of the momentum
$|\mathbf{q}^{*}|$ in (c) as a function of the doping factor $\delta$ for
different interaction strengths computed in a YC4 cylinder.

Finally, in the large doublon doping regime ($0.8>\delta>0.7$), we encounter
an I-SDW’ phase which presents dips at the $K$ points in the spin structure
factor distinguishing it from the I-SDW phase at low doublon doping, see (IV)
in Fig. 1. Moreover, we observe that a ferromagnetic state develops in the
limit of doping with two holes the band insulator at $\delta=1$, see (V) in
Fig. 1. We name this phase flat-band ferromagnetic (FBF) phase and we
interpret it as Müller-Hartmann ferromagnetism [42]. It can also be understood
as Mielke’s flat band ferromagnetism [39, 40, 41] in the limit of
infinitesimal doping. We point out however, that we only find this regime for
two holes regardless of the system size. Hence in the thermodynamic limit it
corresponds to the vanishing hole density, see [79] for a discussion of
ferromagnetism stability at finite density.

Having discussed the main features of the phase diagram, we now present the
observables that we use to identify the different phases. One of the key
characteristics that we study is the spin structure factor

$\displaystyle\chi(\mathbf{q})=\sum_{ij}\left(\langle\mathbf{S}_{i}\mathbf{S}_{j}\rangle-\langle\mathbf{S}_{i}\rangle\langle\mathbf{S}_{j}\rangle\right)e^{i\mathbf{q}\cdot\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)}.$
(2)

We note that the total spin of the system can be related to the same
correlation function using the relation $\chi(\mathbf{0})=S\left(S+1\right)$.
We use analysis of the total spin of the system to identify PP and FP phases.
To confirm the nature of the PP phases we also compute ground state energy of
the system for different values of the total magnetization $S^{z}$. Such
magnetization-energy relation $E(S^{z})$ has a simple physical interpretation.
When the system has a tendency to develop ferromagentic order, fixing $S^{z}$
to a value that is smaller than the preferred magnetization, means that there
will be additional magnetization in the perpendicular direction. Hence for a
PP state we expect energy to be independent on magnetization until
magnetization reaches the preferred total spin of the system. When $S^{z}$
increases beyond the preferred magnetization, energy will start increasing
with increasing magnetization. Another experimentally relevant quantity that
we analyze is the energy difference between the fully polarized state,
$(S^{z}=S_{M})$, and the unpolarized state with $(S^{z}=0)$, which denotes the
energy cost to fully polarize the system. Finally, to detect the appearance of
phase separation we compute the doublon density
$\langle\hat{n}_{i}^{D}\rangle\equiv\langle\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}\rangle$
and inspect where the density is localized in a finite region of the system.

Magnetic polarons and effective interaction between them. To understand the
low-doping regime we first examine the two-doublon problem. We focus on the
regime $U/t\gg 1$, when the $t-J$ model accurately describes the Hubbard
model. Superexchange interaction $J$ favors a 1200 C-SDW configuration of
spins which, based on our simulations, leads to the ground state energy
$E_{J}/N_{s}\sim-1.25J$, where $N_{s}$ is the number of sites. However, a
single doublon minimizes its constrained kinetic energy by inducing a
ferromagnetic background around it as a result of the constructive pattern of
interference of different propagation paths. The competition between these two
terms leads to a Nagaoka’s polaron of finite radius $R$. The doublon’s
propagation is confined within the ferromagnetic region, and therefore, its
kinetic energy can be approximated as $E_{t}/t\sim(2.4)^{2}/R^{2}$ in the
continuum limit, as shown in [49, 62]. Considering both the kinetic and
magnetic energy of the Nagaoka’s polaron, we obtain the total spin $S=\pi
R^{2}/2$,

$\displaystyle S\sim 1.95N_{D}\sqrt{t/J},$ (3)

where we assumed that each Nagaoka’s polaron contributes to the total spin of
the system in the same way. The size of the Nagaoka’s polarons increases as
the superexchange interaction is reduced and therefore, the total spin of the
system increases. When polarons cover the entire system a long-range
ferromagnetic state is formed in agreement with Nagaoka-Thouless’ theorem.

To verify the formation of Nagaoka’s polarons we numerically obtain the total
spin as aforementioned. As superexchange interaction, $J/t$, is decreased, the
total spin $S$ starts to increase until it reaches its saturation value
$S_{M}=(N_{s}-N_{D})/2$ for $J/t\rightarrow 0$, see Fig. 2 (a). Moreover, two
polarons experience an effective interaction induced by the magnetic
environment. If two doublons share their ferromagnetic polaron bubbles, they
can reduce the superexchange energy gain by reducing the area of the shared
ferromagnetic bubble, but in doing so, they incur a higher kinetic energy due
to the Pauli exclusion principle. In the regime $J/t<0.1$, two polarons have a
tendency to come close to each other and share their ferromagnetic bubbles.
This signals an attractive polaron interaction, see Fig. 2 (a). However, for
$J/t>0.1$, the total spin vanishes and two Nagaoka’s polarons effectively
repeal each other.

To confirm the partially polarized nature of the PS regime we compute its
magnetization-energy relation $E(S^{z})$, see Fig. 2 panel (b). The ground
state exhibits no energy dependence with the magnetization until a critical
magnetization is reached at which the energy starts to increase. Moreover, the
critical magnetization coincides with the total spin of the ground state at
zero magnetization confirming that each Nagaoka’s polaron creates a finite
region of space that costs no energy to polarize in any direction.

Density-induced ferromagnetism and phase separation. We now consider the
regime of low doping, $\delta<0.1$, and large on-site interaction, $J/t<0.1$,
see Fig. 1. In this case we find phase separation between a doublon-rich
ferromagnetic region and an undoped region with 120o C-SDW order. In this
regime, the total spin increases linearly with the number of doublons denoting
that the ferromagnetic region increases with doublon doping. Consequently the
region with C-SDW order diminishes leading to a state which exhibits partial
polarization. By increasing the doublon density the ferromagnetic puddle
starts to increase, eventually occupying the entire system. At this point,
doublons delocalize throughout the entire system and the system transitions to
a long-range ferromagnetic state.

In the intermediate interacting regime ($40>U/t>20$) we perform simulations of
the Hubbard model Eq. (1) in YC4 cylinders. We do not observe signatures of a
PS state at $U/t\leq 40$ for small doublon dopings. Instead, we note that
doublons establish an I-SDW at finite momenta, see Fig. 1 (II) and the
associated momenta decrease with increasing doublon doping, see Fig. 3. We
observe that the momenta suddenly jump from a finite to a vanishing value at a
critical doublon doping, denoting a first order phase transition between the
I-SDW and the FP phases. In Fig. 2 (c) we present the polarization energy for
the Hubbard model at different interaction strengths. The polarization energy
linearly decreases (increases) upon particle (hole) doping, signaling that
dopants favor ferromagnetism or antiferromagnetism in the particle or hole
doping regimes, respectively. Furthermore, the linear relation suggests that
each dopant contributes to the magnetic order, reinforcing the idea that
doublons trigger the formation of ferromagnetic regions. When polarization
energy vanishes, a transition towards a long-range ferromagnet occurs.

Stability of long-range ferromagnetism at large doublon dopings. In the doping
regime $\delta\sim 0.4-0.6$ we observe a strong ferromagnetic phase in both
YC4 and YC6 cylinders, that persists up to the intermediate interacting regime
$U/t\sim 20$, see Fig. 3. The ferromagnetic phase in the intermediate
interacting regime is centered around $\delta\sim 1/2$ signaling the
importance of the van Hove singularity and suggesting that the mechanism
leading to a ferromagnetic state is likely to be a Stoner instability. We note
that there is no global distinction between Stoner-type ferromagnetism and the
Nagaoka-type kinetic ferromagnetism aforementioned, since they correspond to
the same phase. Notably, our simulations in YC6 cylinders reveal that these
two mechanisms can be continuously connected by increasing $U/t$.

Long-range ferromagnetism becomes unstable around $1>\delta\gtrsim 0.7$ in
both geometries (YC4 and YC6) even in the limit $U/t\rightarrow\infty$. In
this regime, we note that the spin structure factor exhibits minima at the
$K$, $K^{\prime}$ and $\Gamma$ points and an I-SDW’ is formed, see Fig. 1
(IV). Moreover, we observe that the polarization energy per site becomes very
small in the I-SDW’ phase, see Fig. 2 (c). Furthermore, we observe that the
spin structure factor becomes almost independent on the interaction strength
$U/t$ in this high doping regime, see Fig.3.

Finally, let us discuss the vanishing doping situation. Specifically, we
explore the case of two holes introduced on top of the band insulator
$\delta\rightarrow 1^{-}$. By performing a low-energy expansion around the $K$
and $K^{\prime}$ points we obtain a ground state which is a spin triplet and a
valley (orbital) singlet, see [42, 80] for a similar discussion. In this
limit, the ground state exhibits maximum spin and the spin structure factor
exhibits features at the $K$, $K^{\prime}$ and $\Gamma$ points, see Fig. 1
(V), reflecting that the holes populate the $K$ and $K^{\prime}$ valleys in an
antisymmetric way. This has to be contrasted with the aforementioned long-
range ferromagnetism, which is characterized by a single peak in the spin
structure factor at the $\Gamma$ point. Moreover, flat-band ferromagnetism
survives for any interaction strength $U/t$. Therefore, our results suggest
that Müller-Hartmann ferromagentism appears in the vanishing doping limit
($\delta\rightarrow 1^{-}$), but it becomes unstable at finite doping and it
is not connected with the ferromagentic phase at intermediate dopings
($0.7>\delta>0$) even in the limit $U/t\rightarrow\infty$. Although we cannot
entirely rule out the possibility of finite density Müller-Hartmann
ferromagnetism being present in the full 2D limit, our simulations conducted
in our widest cylinder do not support its existence.

Conclusions and outlook. We have investigated the ground state magnetic
properties of the Hubbard model in the triangular lattice at finite doping. In
the low doublon doped regime and for $U/t>40$ we find that individual doublons
form Nagaoka’s ferromagentic polarons. Attractive interaction between such
polarons results in a phase separation between a half-filled state and a
doublon rich phase with full polarization. However, for $U/t<40$ we observe an
incommensurate spin density wave. By increasing particle doping we find a
transition towards a long-range ferromagnetism for all interaction strengths
considered. Therefore, our results show that long-range ferromagnetism is
stable in the particle doped triangular lattice Hubbard model for a wide range
of parameters.

In the intermediate interacting regime $U/t\sim 10$ one also expects the
emergence of a chiral spin liquid [21, 25]. However, the fate of the chiral
spin liquid upon particle doping is still unclear and we are currently limited
by system size to speculate on the chirality of the ground state.
Nevertheless, recent works have pointed out that a chiral metallic phase
arises when the spin liquid is hole doped [36]. Consequently, a comprehensive
exploration of the doublon-doped system in this intermediate interacting
regime will provide deeper insights on the chiral nature of the ground state.
Furthermore, hole spectral functions have been shown to be a useful observable
that reveals the chiral spin nature of the ground state and the appearance of
dynamical Nagaoka’s polarons [81, 82]. In combination with the recent
photoexcitation probes discussed in [83] the nature of the Nagaoka’s polarons
could be explored in current cold atoms experiments.

Finally, our results reveal the occurrence of phase separation within the low
doublon-doped regime, signaling the presence of an effective attraction
between Nagaoka’s polarons. However, phase separation can be suppressed upon
inclusion of non-local Coulomb interaction, which can be relevant for
experiments with TMDCs. We envision that Coulomb interaction could stabilize a
bipolaron without leading to phase separation which could pave the way to
explain the emergence of superconductivity in doped triangular Mott
insulators. Furthermore, thermal fluctuations can also overcome attraction
between Nagaoka’s polarons and lead to a uniform polaronic phase.

Note added. When this paper was being finalized, we learned about paper [84]
that addressed the phase diagram for the particular filling factor
$\delta=1/2$.

Acknowledgments. We acknowledge useful discussions with Waseem S. Bakr, Utso
Bhattacharya, Annabelle Bohrdt, Anant Kale, Youqi Gang, Markus Greiner, Fabian
Grusdt, Lev Haldar Kendrick, Wen-Wei Ho, Atac Imamoglu, Martin Lebrat, Max L.
Prichard, Benjamin M. Spar, Muqing Xu, Zoe Z. Yan. We acknowledge support from
the ARO grant number W911NF-20-1-0163, the SNSF project 200021_212899, and the
Swiss State Secretariat for Education, Research and Innovation (SERI) under
contract number UeM019-1. Tensor network computations have been performed
using TeNPy [85].

## References

  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Simulation of hubbard model physics in wse2/ws2 moiré superlattices, Nature 579, 353 (2020).
  * Ciorciaro  _et al._ [2023] L. Ciorciaro, T. Smoleński, I. Morera, N. Kiper, S. Hiestand, M. Kroner, Y. Zhang, K. Watanabe, T. Taniguchi, E. Demler, and A. İmamoğlu, Kinetic magnetism in triangular moiré materials, Nature 623, 509 (2023).
  * Tao  _et al._ [2023] Z. Tao, W. Zhao, B. Shen, P. Knüppel, K. Watanabe, T. Taniguchi, J. Shan, and K. F. Mak, Observation of spin polarons in a frustrated moiré hubbard system, arXiv:2307.12205 (2023).
  * Morita  _et al._ [2002] H. Morita, S. Watanabe, and M. Imada, Nonmagnetic insulating states near the mott transitions on lattices with geometrical frustration and implications for $\kappa$-(et)2cu2(cn)3, Journal of the Physical Society of Japan 71, 2109 (2002), https://doi.org/10.1143/JPSJ.71.2109 .
  * Honerkamp [2003] C. Honerkamp, Instabilities of interacting electrons on the triangular lattice, Phys. Rev. B 68, 104510 (2003).
  * Shimizu  _et al._ [2003] Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, and G. Saito, Spin liquid state in an organic mott insulator with a triangular lattice, Phys. Rev. Lett. 91, 107001 (2003).
  * Koretsune  _et al._ [2007] T. Koretsune, Y. Motome, and A. Furusaki, Exact diagonalization study of mott transition in the hubbard model on an anisotropic triangular lattice, Journal of the Physical Society of Japan 76, 074719 (2007).
  * Sahebsara and Sénéchal [2008] P. Sahebsara and D. Sénéchal, Hubbard model on the triangular lattice: Spiral order and spin liquid, Phys. Rev. Lett. 100, 136402 (2008).
  * Clay  _et al._ [2008] R. T. Clay, H. Li, and S. Mazumdar, Absence of superconductivity in the half-filled band hubbard model on the anisotropic triangular lattice, Phys. Rev. Lett. 101, 166403 (2008).
  * Watanabe  _et al._ [2008] T. Watanabe, H. Yokoyama, Y. Tanaka, and J. Inoue, Predominant magnetic states in the hubbard model on anisotropic triangular lattices, Phys. Rev. B 77, 214505 (2008).
  * Galanakis  _et al._ [2009] D. Galanakis, T. D. Stanescu, and P. Phillips, Mott transition on a triangular lattice, Phys. Rev. B 79, 115116 (2009).
  * Yoshioka  _et al._ [2009] T. Yoshioka, A. Koga, and N. Kawakami, Quantum phase transitions in the hubbard model on a triangular lattice, Phys. Rev. Lett. 103, 036401 (2009).
  * Yang  _et al._ [2010] H.-Y. Yang, A. M. Läuchli, F. Mila, and K. P. Schmidt, Effective spin model for the spin-liquid phase of the hubbard model on the triangular lattice, Phys. Rev. Lett. 105, 267204 (2010).
  * Kokalj and McKenzie [2013] J. Kokalj and R. H. McKenzie, Thermodynamics of a bad metal–mott insulator transition in the presence of frustration, Phys. Rev. Lett. 110, 206402 (2013).
  * Yamada [2014] A. Yamada, Magnetic properties and mott transition in the hubbard model on the anisotropic triangular lattice, Phys. Rev. B 89, 195108 (2014).
  * Li  _et al._ [2014] G. Li, A. E. Antipov, A. N. Rubtsov, S. Kirchner, and W. Hanke, Competing phases of the hubbard model on a triangular lattice: Insights from the entropy, Phys. Rev. B 89, 161118 (2014).
  * Tocchio  _et al._ [2014] L. F. Tocchio, C. Gros, R. Valentí, and F. Becca, One-dimensional spin liquid, collinear, and spiral phases from uncoupled chains to the triangular lattice, Phys. Rev. B 89, 235107 (2014).
  * Laubach  _et al._ [2015] M. Laubach, R. Thomale, C. Platt, W. Hanke, and G. Li, Phase diagram of the hubbard model on the anisotropic triangular lattice, Phys. Rev. B 91, 245125 (2015).
  * Mishmash  _et al._ [2015] R. V. Mishmash, I. González, R. G. Melko, O. I. Motrunich, and M. P. A. Fisher, Continuous mott transition between a metal and a quantum spin liquid, Phys. Rev. B 91, 235140 (2015).
  * Shirakawa  _et al._ [2017] T. Shirakawa, T. Tohyama, J. Kokalj, S. Sota, and S. Yunoki, Ground-state phase diagram of the triangular lattice hubbard model by the density-matrix renormalization group method, Phys. Rev. B 96, 205130 (2017).
  * Szasz  _et al._ [2020] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Chiral spin liquid phase of the triangular lattice hubbard model: A density matrix renormalization group study, Phys. Rev. X 10, 021042 (2020).
  * Szasz and Motruk [2021] A. Szasz and J. Motruk, Phase diagram of the anisotropic triangular lattice hubbard model, Phys. Rev. B 103, 235132 (2021).
  * Cookmeyer  _et al._ [2021] T. Cookmeyer, J. Motruk, and J. E. Moore, Four-spin terms and the origin of the chiral spin liquid in mott insulators on the triangular lattice, Phys. Rev. Lett. 127, 087201 (2021).
  * Wietek  _et al._ [2021] A. Wietek, R. Rossi, F. Šimkovic, M. Klett, P. Hansmann, M. Ferrero, E. M. Stoudenmire, T. Schäfer, and A. Georges, Mott insulating states with competing orders in the triangular lattice hubbard model, Phys. Rev. X 11, 041013 (2021).
  * Chen  _et al._ [2022] B.-B. Chen, Z. Chen, S.-S. Gong, D. N. Sheng, W. Li, and A. Weichselbaum, Quantum spin liquid with emergent chiral order in the triangular-lattice hubbard model, Phys. Rev. B 106, 094420 (2022).
  * Downey  _et al._ [2023] P.-O. Downey, O. Gingras, J. Fournier, C.-D. Hébert, M. Charlebois, and A.-M. S. Tremblay, Mott transition, widom line, and pseudogap in the half-filled triangular lattice hubbard model, Phys. Rev. B 107, 125159 (2023).
  * Haerter and Shastry [2005] J. O. Haerter and B. S. Shastry, Kinetic antiferromagnetism in the triangular lattice, Phys. Rev. Lett. 95, 087202 (2005).
  * Haerter  _et al._ [2006a] J. O. Haerter, M. R. Peterson, and B. S. Shastry, Finite-temperature properties of the triangular lattice $t\text{$-$}j$ model and applications to ${\mathrm{na}}_{x}\mathrm{Co}{\mathrm{o}}_{2}$, Phys. Rev. B 74, 245118 (2006a).
  * Haerter  _et al._ [2006b] J. O. Haerter, M. R. Peterson, and B. S. Shastry, Strong correlations produce the curie-weiss phase of ${\mathrm{na}}_{x}{\mathrm{coo}}_{2}$, Phys. Rev. Lett. 97, 226402 (2006b).
  * Sposetti  _et al._ [2014] C. N. Sposetti, B. Bravo, A. E. Trumper, C. J. Gazza, and L. O. Manuel, Classical antiferromagnetism in kinetically frustrated electronic models, Phys. Rev. Lett. 112, 187204 (2014).
  * Zhang  _et al._ [2018] S.-S. Zhang, W. Zhu, and C. D. Batista, Pairing from strong repulsion in triangular lattice hubbard model, Phys. Rev. B 97, 140507 (2018).
  * Morera  _et al._ [2021] I. Morera, A. Bohrdt, W. W. Ho, and E. Demler, Attraction from frustration in ladder systems 10.48550/arxiv.2106.09600 (2021).
  * Lee  _et al._ [2023] K. Lee, P. Sharma, O. Vafek, and H. J. Changlani, Triangular lattice hubbard model physics at intermediate temperatures, Phys. Rev. B 107, 235105 (2023).
  * Morera  _et al._ [2023a] I. Morera, M. Kanász-Nagy, T. Smolenski, L. Ciorciaro, A. m. c. Imamoğlu, and E. Demler, High-temperature kinetic magnetism in triangular lattices, Phys. Rev. Res. 5, L022048 (2023a).
  * Davydova  _et al._ [2023] M. Davydova, Y. Zhang, and L. Fu, Itinerant spin polaron and metallic ferromagnetism in semiconductor moiré superlattices, Phys. Rev. B 107, 224420 (2023).
  * Zhu  _et al._ [2022] Z. Zhu, D. N. Sheng, and A. Vishwanath, Doped mott insulators in the triangular-lattice hubbard model, Phys. Rev. B 105, 205110 (2022).
  * Nagaoka [1966] Y. Nagaoka, Ferromagnetism in a narrow, almost half-filled $s$ band, Phys. Rev. 147, 392 (1966).
  * Thouless [1965] D. J. Thouless, Exchange in solid 3He and the Heisenberg Hamiltonian, Proceedings of the Physical Society 86, 893 (1965).
  * Mielke [1991a] A. Mielke, Ferromagnetic ground states for the hubbard model on line graphs, Journal of Physics A: Mathematical and General 24, L73 (1991a).
  * Mielke [1991b] A. Mielke, Ferromagnetism in the hubbard model on line graphs and further considerations, Journal of Physics A: Mathematical and General 24, 3311 (1991b).
  * Mielke [1992] A. Mielke, Exact ground states for the hubbard model on the kagome lattice, Journal of Physics A: Mathematical and General 25, 4335 (1992).
  * Müller-Hartmann [1995] E. Müller-Hartmann, Ferromagnetism in hubbard models: Low density route, Journal of Low Temperature Physics 99, 349 (1995).
  * Hanisch  _et al._ [1995] T. Hanisch, B. Kleine, A. Ritzl, and E. Müller-Hartmann, Ferromagnetism in the hubbard model: instability of the nagaoka state on the triangular, honeycomb and kagome lattices, Annalen der Physik 507, 303 (1995).
  * Gao  _et al._ [2007] M. Gao, S. Zhou, and Z. Wang, Itinerant and localized magnetism on the triangular lattice: Sodium-rich phases of ${\mathrm{na}}_{x}{\mathrm{coo}}_{2}$, Phys. Rev. B 76, 180402 (2007).
  * Martin and Batista [2008] I. Martin and C. D. Batista, Itinerant electron-driven chiral magnetic ordering and spontaneous quantum hall effect in triangular lattice models, Phys. Rev. Lett. 101, 156402 (2008).
  * Alexandrov and Devreese [2010] A. S. Alexandrov and J. T. Devreese,  _Advances in polaron physics_ (Springer Heidelberg, Heidelberg, 2010).
  * lordanskil and Smirnov [1967] S. V. lordanskil and A. V. Smirnov, Hole concentration dependence of the curie temperature of (ga,mn)sb in a field-effect structure, Zh. Eksp. Teor. Fiz. 6, 484 (1967).
  * lordanskil and Smirnov [1968] S. V. lordanskil and A. V. Smirnov, Hole concentration dependence of the curie temperature of (ga,mn)sb in a field-effect structure, Zh. Eksp. Teor. Fiz. 54, 228 (1968).
  * Nagaev [1992] E. Nagaev, Self-trapped states of charge carriers in magnetic semiconductors, Journal of Magnetism and Magnetic Materials 110, 39 (1992).
  * Kaminski and Das Sarma [2002] A. Kaminski and S. Das Sarma, Polaron percolation in diluted magnetic semiconductors, Phys. Rev. Lett. 88, 247202 (2002).
  * L. N. Bulaevskii and Khomskii [1968] E. L. N. L. N. Bulaevskii and D. I. Khomskii, A new type of auto-localized state of a conduction electron in an antiferromagnetic semiconductor, Zh. Eksp. Teor. Fiz. 54, 1562 (1968).
  * Brinkman and Rice [1970] W. F. Brinkman and T. M. Rice, Single-particle excitations in magnetic insulators, Phys. Rev. B 2, 1324 (1970).
  * Trugman [1988] S. A. Trugman, Interaction of holes in a hubbard antiferromagnet and high-temperature superconductivity, Phys. Rev. B 37, 1597 (1988).
  * Schmitt-Rink  _et al._ [1988] S. Schmitt-Rink, C. M. Varma, and A. E. Ruckenstein, Spectral function of holes in a quantum antiferromagnet, Phys. Rev. Lett. 60, 2793 (1988).
  * Shraiman and Siggia [1988] B. I. Shraiman and E. D. Siggia, Mobile vacancies in a quantum heisenberg antiferromagnet, Phys. Rev. Lett. 61, 467 (1988).
  * Sachdev [1989] S. Sachdev, Hole motion in a quantum néel state, Phys. Rev. B 39, 12232 (1989).
  * Kane  _et al._ [1989] C. L. Kane, P. A. Lee, and N. Read, Motion of a single hole in a quantum antiferromagnet, Phys. Rev. B 39, 6880 (1989).
  * Dagotto  _et al._ [1989] E. Dagotto, A. Moreo, and T. Barnes, Hubbard model with one hole: Ground-state properties, Phys. Rev. B 40, 6721 (1989).
  * Boninsegni and Manousakis [1991] M. Boninsegni and E. Manousakis, Quasihole excitation in a quantum antiferromagnet: Variational monte carlo calculation, Phys. Rev. B 43, 10353 (1991).
  * Boninsegni and Manousakis [1992] M. Boninsegni and E. Manousakis, Variational description of a quasihole excitation in a quantum antiferromagnet, Phys. Rev. B 45, 4877 (1992).
  * White and Scalapino [1997] S. R. White and D. J. Scalapino, Hole and pair structures in the t-j model, Phys. Rev. B 55, 6504 (1997).
  * White and Affleck [2001] S. R. White and I. Affleck, Density matrix renormalization group analysis of the nagaoka polaron in the two-dimensional $t-j$ model, Phys. Rev. B 64, 024411 (2001).
  * Kung  _et al._ [2015] Y. F. Kung, E. A. Nowadnick, C. J. Jia, S. Johnston, B. Moritz, R. T. Scalettar, and T. P. Devereaux, Doping evolution of spin and charge excitations in the hubbard model, Phys. Rev. B 92, 195108 (2015).
  * Wang  _et al._ [2015] Y. Wang, K. Wohlfeld, B. Moritz, C. J. Jia, M. van Veenendaal, K. Wu, C.-C. Chen, and T. P. Devereaux, Origin of strong dispersion in hubbard insulators, Phys. Rev. B 92, 075119 (2015).
  * Grusdt  _et al._ [2018] F. Grusdt, M. Kánasz-Nagy, A. Bohrdt, C. S. Chiu, G. Ji, M. Greiner, D. Greif, and E. Demler, Parton theory of magnetic polarons: Mesonic resonances and signatures in dynamics, Phys. Rev. X 8, 011046 (2018).
  * Koepsell  _et al._ [2019] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Imaging magnetic polarons in the doped fermi–hubbard model, Nature 572, 358 (2019).
  * Soriano and Katsnelson [2020] D. Soriano and M. I. Katsnelson, Magnetic polaron and antiferromagnetic-ferromagnetic transition in doped bilayer ${\mathrm{cri}}_{3}$, Phys. Rev. B 101, 041402 (2020).
  * Jiang  _et al._ [2020] Y.-F. Jiang, J. Zaanen, T. P. Devereaux, and H.-C. Jiang, Ground state phase diagram of the doped hubbard model on the four-leg cylinder, Phys. Rev. Res. 2, 033073 (2020).
  * Ji  _et al._ [2021] G. Ji, M. Xu, L. H. Kendrick, C. S. Chiu, J. C. Brüggenjürgen, D. Greif, A. Bohrdt, F. Grusdt, E. Demler, M. Lebrat, and M. Greiner, Coupling a mobile hole to an antiferromagnetic spin background: Transient dynamics of a magnetic polaron, Phys. Rev. X 11, 021022 (2021).
  * Koepsell  _et al._ [2021] J. Koepsell, D. Bourgund, P. Sompet, S. Hirthe, A. Bohrdt, Y. Wang, F. Grusdt, E. Demler, G. Salomon, C. Gross, and I. Bloch, Microscopic evolution of doped mott insulators from polaronic metal to fermi liquid, Science 374, 82 (2021), https://www.science.org/doi/pdf/10.1126/science.abe7165 .
  * Dehollain  _et al._ [2020] J. P. Dehollain, U. Mukhopadhyay, V. P. Michal, Y. Wang, B. Wunsch, C. Reichl, W. Wegscheider, M. S. Rudner, E. Demler, and L. M. K. Vandersypen, Nagaoka ferromagnetism observed in a quantum dot plaquette, Nature 579, 528 (2020).
  * Xu  _et al._ [2023] M. Xu, L. H. Kendrick, A. Kale, Y. Gang, G. Ji, R. T. Scalettar, M. Lebrat, and M. Greiner, Frustration- and doping-induced magnetism in a fermi–hubbard simulator, Nature 620, 971 (2023).
  * Lebrat  _et al._ [2023] M. Lebrat, M. Xu, L. H. Kendrick, A. Kale, Y. Gang, P. Seetharaman, I. Morera, E. Khatami, E. Demler, and M. Greiner, Observation of nagaoka polarons in a fermi-hubbard quantum simulator, arXiv:2308.12269 (2023).
  * Prichard  _et al._ [2023] M. L. Prichard, B. M. Spar, I. Morera, E. Demler, Z. Z. Yan, and W. S. Bakr, Directly imaging spin polarons in a kinetically frustrated hubbard system, arXiv:2308.12951 (2023).
  * Koretsune and Ogata [2002] T. Koretsune and M. Ogata, Resonating-valence-bond states and ferromagnetic correlations in the doped triangular mott insulator, Phys. Rev. Lett. 89, 116401 (2002).
  * Watanabe  _et al._ [2004] T. Watanabe, H. Yokoyama, Y. Tanaka, J.-i. Inoue, and M. Ogata, Variational monte carlo studies of pairing symmetry for the t–j model on a triangular lattice, Journal of the Physical Society of Japan 73, 3404 (2004), https://doi.org/10.1143/JPSJ.73.3404 .
  * Weber  _et al._ [2006] C. Weber, A. Läuchli, F. Mila, and T. Giamarchi, Magnetism and superconductivity of strongly correlated electrons on the triangular lattice, Phys. Rev. B 73, 014519 (2006).
  * Yamada  _et al._ [1998] K. Yamada, C. H. Lee, K. Kurahashi, J. Wada, S. Wakimoto, S. Ueki, H. Kimura, Y. Endoh, S. Hosoya, G. Shirane, R. J. Birgeneau, M. Greven, M. A. Kastner, and Y. J. Kim, Doping dependence of the spatially modulated dynamical spin correlations and the superconducting-transition temperature in ${\mathrm{la}}_{2\mathrm{-}\mathit{x}}{\mathrm{sr}}_{x}{\mathrm{cuo}}_{4}$, Phys. Rev. B 57, 6165 (1998).
  * Pieri  _et al._ [1996] P. Pieri, S. Daul, D. Baeriswyl, M. Dzierzawa, and P. Fazekas, Low-density ferromagnetism in the hubbard model, Phys. Rev. B 54, 9250 (1996).
  * Li  _et al._ [2023] C. Li, M.-G. He, C.-Y. Wang, and H. Zhai, Frustration induced itinerant ferromagnetism of fermions in optical lattice, arXiv:2305.01682 (2023).
  * Kadow  _et al._ [2022] W. Kadow, L. Vanderstraeten, and M. Knap, Hole spectral function of a chiral spin liquid in the triangular lattice hubbard model, Phys. Rev. B 106, 094417 (2022).
  * Kadow  _et al._ [2023] W. Kadow, H.-K. Jin, J. Knolle, and M. Knap, Single-hole spectra of kitaev spin liquids: From dynamical nagaoka ferromagnetism to spin-hole fractionalization, arXiv:2309.15157 (2023).
  * Morera  _et al._ [2023b] I. Morera, C. Weitenberg, K. Sengstock, and E. Demler, Exploring kinetically induced bound states in triangular lattices with ultracold atoms: spectroscopic approach, arXiv:2312.00768 (2023b).
  * He  _et al._ [2023] Y. He, R. Rausch, M. Peschke, C. Karrasch, P. Corboz, N. Bultinck, and S. A. Parameswaran, Itinerant magnetism in the triangular lattice hubbard model at half-doping: Application to twisted transition-metal dichalcogenides, (2023), arXiv:2311.10146 [cond-mat.str-el] .
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy), SciPost Phys. Lect. Notes , 5 (2018).