in every time-step. In our simulations, we adopt the gradient-based probing technique combined with the kernel polynomial method (KPM) [@barros13; @wang18; @chern18]. This allows for an efficient estimation of the electron density matrix $\rho^{\rm qp}$ at a computational cost that scales linearly with system size. All simulations discussed below were performed on triangular lattice of unprecedented size of $N = 150^2$ sites.

In this work we focus on the intrinsic Mott transition in the absence of magnetic ordering. We consider only non-magnetic Gutzwiller solutions, $p_{i, \uparrow} = p_{i, \downarrow} = p_i$, and $n_{i, \uparrow} = n_{i, \downarrow} = n_i$, in all simulations discussed below. The following model parameters are used. $t^{(0)}_{\rm nn1} = -1$, $g = 0.55$, $K_0 = K_1 = 0.35$, $m = 2$, $\gamma = 0.2$, and the time step $\Delta t = 0.02$. The number of the Chebyshev polynomials $M = 2000$. In order to assist transitions between the local minima of the local order parameter $\Phi_i$, a small noise of the order of $\delta\Phi_i = 0.1$ is added to the slave-boson amplitudes obtained from the previous time-step at the beginning of the new Gutzwiller iteration.

Thermal Mott transition {#sec:thermal-mott}
=======================

![(Color online) \[fig:phase-diagram\] (a) Schematic phase 