text:
\section{...}

\section*{Acknowledgements}

We thank Jan-Lukas Wynen for helpful discussions on the Hubbard model and software issues. We also thank Michael Kajan for proof reading and for providing a lot of detailed comments.
This work was funded, in part, through financial support from the Deutsche Forschungsgemeinschaft (Sino-German CRC 110 and SFB TRR-55).
During this work E.B. was supported by the U.S. Department of Energy under Contract No. DE-FG02-93ER-40762.
The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA~\cite{jureca} at Forschungszentrum J\"ulich.
We also gratefully acknowledge time on DEEP~\cite{DEEP}, an experimental modular supercomputer at the J\"ulich Supercomputing Centre.
The analysis was mostly done in \texttt{R}~\cite{R_language} using \texttt{hadron}~\cite{hadron}.
We are indebted to Bartosz Kostrzewa for helping us detect a compiler bug and for lending us his expertise on HPC hardware.

\section{Analysis}\label{sec:analysis}

All our results have been extrapolated to the thermodynamic and continuous time limits in a two-dimensional fit simultaneously. In particular the single particle gap and the staggered magnetisation defined by
\begin{equation}
	m_s^2 = \frac{\sum_i S^{ii}_{--}}{V(\kappa\beta)^2}= \frac{2S^{11}_{--}+S^{33}_{--}}{V(\kappa\beta)^2},
	\label{eq:definition_m_s}
\end{equation}
have been fitted using the functional dependences
\begin{align}
	\Delta^2(L,N_t^{}) &=
	\Delta_0^2 + c_0^{\Delta} N_t^{-2}
	+ c_1^{\Delta} L^{-3}\,,\label{eqn_gap_artefacts}\\
	m_s &= m_{s,0} + c^{m_s}_0\,\delta + c^{m_s}_1 \,L^{-2}\,.\label{eqn:extrapolation_m_s}
\end{align}
An example of such a fit can be found in figure~\ref{fig:2d_gap_fit}.

\begin{figure}[t]
	\centering
	\raisebox{0.016\height}{\includegraphics[width=.45\textwidth]{{b_8_u_3.5_m_0_gap_l}.pdf}}
	\includegraphics[width=.45\textwidth]{{b_8_u_3.5_m_0_gap_t}.pdf}
	\caption{Simultaneous two-dimensional fit of $\Delta(N_t,L)$ (in units of $\kappa$) using Eq.~\eqref{eqn_gap_artefacts}, for $\kappa\beta=8$ and $U/\kappa=\num{3.5}$. Data points for $L<9$ have been omitted from the fit,
		but not from the plot. Very small lattices lead to large values of $\Delta$, which are not visible on the scale of the plot.
		This fit has $\chi^2/\text{d.o.f.} \simeq \num{1.1}$, corresponding to a p-value of $\simeq \num{0.34}$.
		\label{fig:2d_gap_fit}}
\end{figure}

\subsection{Results \label{sec:gap_results}}

Our results for $\Delta_0$ and $m_s$ are shown
in \Figref{both_order_parameters} for all values of $U/\kappa$ and $\kappa\beta$, along with an extrapolation (with error band) to zero temperature ($\beta\to\infty$).
For details as to the zero-temperature gap, see Ref.~\cite{semimetalmott}.

\begin{figure}[ht]
	\raisebox{0.06\height}{
		\resizebox{0.5\textwidth}{!}{{\Large
				\input{data/gap}}}}
	\resizebox{0.505\textwidth}{!}{
			\input{data/m_s_intensive}}}
	\caption{All quantities in units of $\kappa$ and after the thermodynamic and continuum limit extrapolations. $\beta$ is the inverse temperature.
		The single-particle gap $\Delta_0(U,\beta)$ (left) and the AFMI order parameter (staggered magnetization) $m_s$ (right).
		We also show $\Delta_0(U,\beta = \infty)$ and $m_s(U,\beta = \infty)$ as solid black lines with error band. The legend from the left plot applies to both.}
	\label{fig:both_order_parameters}
\end{figure}

Finally, we perform a simultaneous data collapse fit to
\begin{align}
	\beta\Delta_0 &= F(\beta^\mu(U-U_c))\,,\label{eq:collapse_delta}\\
	\beta^{\mu\upbeta} m_s &= G(\beta^\mu(U-U_c))\label{eq:collapse_m_s}
\end{align}
with some universal functions $F$ and $G$. Here $\mu=1/(z\nu)$ with the dynamical exponent $z=1$. Thus we obtain the critical parameters $U_c=\ensuremath{\num{3.835(14)}}$, $\nu=\critNu$ and $\upbeta=\ensuremath{\num{0.898(37)}}$. The collapse plots are shown in figure~\ref{figure:structure factors}. The critical coupling $U_c$ and the critical exponent $\mu$ appear in both contributions~\eqref{eq:collapse_delta} and~\eqref{eq:collapse_m_s} to the collapse fit proving that gap and AFM order do indeed emerge simultaneously.

\begin{figure*}[ht]
	\includeepslatex{data/gap_collapse0}
	\includeepslatex{data/m_s_collapse}
	\caption{Data collapse plots with the optimal parameters of $U_c$, $\nu$, and $\upbeta$ obtained from a simultaneous collapse fit to the gap $\Delta$ (left panel) and the order parameter $m_s$ (right panel).
	Note that the ``outliers'' are due to particularly small $\beta$ and are excluded from the analysis.
	 All quantities are expressed in units of $\kappa$.}
	\label{figure:structure factors}
\end{figure*}

\section{Conclusions}\label{sec:conclusion}

Our work represents the first instance where the grand canonical Brower-Rebbi-Schaich (BRS) algorithm has been applied to the hexagonal Hubbard model
(beyond mere proofs of principle), and we have found highly promising results.

We calculated the single particle gap $\Delta$ as well as all operators that contribute to the antiferromagnetic (AFM), ferromagnetic (FM), and charge density wave (CDW) order parameters of the Hubbard Model on a honeycomb lattice (all the details in Refs.~\cite{semimetalmott,more_observables}).
Furthermore we provide a comprehensive analysis of the temporal continuum, thermodynamic and zero-temperature limits for all these quantities. The favorable scaling of the HMC
enabled us to simulate lattices with $L > 100$ and to
perform a highly systematic treatment of all three limits.

The semimetal-antiferromagnetic Mott insulator (SM-AFMI) transition falls into the Gross-Neveu (GN)-Heisenberg
$SU(2)$ universality class~\cite{PhysRevLett.97.146401}.
The GN-Heisenberg critical exponents have been studied by means of
multiple methods.
In Table~\ref{tab:crit_values}, we give an up-to-date comparison of these calculations with our results. Our value for $U_c/\kappa$ is in overall agreement with previous Monte Carlo (MC) simulations.
For the critical exponents $\nu$ and $\upbeta$, the situation is less clear. Our results for $\nu$ (assuming $z = 1$ due to Lorentz invariance~\cite{PhysRevLett.97.146401}) agree
best with the MC calculation (in the Blankenbecler-Sugar-Scalapino (BSS) formulation) of Ref.~\cite{Buividovich:2018yar}, followed by the FRG and large~$N$ calculations.
On the other hand, our critical exponent $\nu$ is systematically larger than most projection Monte Carlo (PMC) calculations and first-order $4-\epsilon$ expansion results. The agreement appears to be significantly
improved when the $4-\epsilon$ expansion is taken to higher orders, although the discrepancy between expansions for $\nu$ and $1/\nu$ persists.
Finally, our critical exponent $\upbeta$ does not agree with any results previously derived in the literature. They have been clustering in two regions. The PMC methods and first order $4-\epsilon$ expansion yielding values between $0.7$ and $0.8$, the other methods predicting values larger than 1. Our result lies within this gap at approximately $0.9$ and our uncertainties do not overlap with any of the other results.

\input{section/other-results}

Thus, though we are confident to have pinned down the nature of the phase transition and to have performed a thorough analysis of the critical parameters, the values of $\nu$ and $\upbeta$ remain ambiguous. This is mostly due to a large spread of incompatible results existing in the literature prior to this work. With our values derived by an independent method we add a valuable confirmation for the critical coupling and some estimations of $\nu$ as well as a new but plausible result for $\upbeta$.

Our progress sets the stage for future high-precision calculations of additional observables of
the Hubbard model and its extensions, as well as other Hamiltonian theories of strongly correlated electrons.
We anticipate the continued advancement of calculations with ever increasing system sizes, through the leveraging of additional state-of-the-art techniques from lattice QCD,
such as multigrid solvers on GPU-accelerated architectures. We are actively pursuing research along these lines.

\section{Analysis \label{sec:analysis}}

We now analyze the single-particle gap $\Delta_0$ as a function of coupling $U-U_c$ and inverse temperature $\beta$, in order to determine the critical coupling
$U_c/\kappa$ of the quantum (AFMI) phase transition, along with some of the associated critical exponents. Our MC results for $\Delta_0$ as a function of $U$ are shown in \Figref{delta_against_U}.
We make use of the standard finite-size scaling (FSS) method~\cite{newmanb99, dutta_aeppli, Shao_2016}, whereby a given observable $\mathfrak{O}$ is described by
\begin{equation}
\mathfrak{O} = (U-U_c)^q \mathcal{F}_{\mathfrak{O}}^{}(L/\xi, \beta/\xi_t^{}),
\label{eqn:scaling_1}
\end{equation}
in the vicinity of $U_c$, where $q$ is the relevant critical exponent. The scaling function $\mathcal{F}$ accounts for the effects of finite spatial system size $L$ and inverse temperature $\beta$.
The spatial and temporal correlation lengths $\xi$ and $\xi_t$ are
\begin{equation}
\xi \sim (U-U_c)^{-\nu},
\quad
\xi_t^{} \sim (U-U_c)^{-\nu_t^{}} = (U-U_c)^{-z\nu} \sim \xi^z,
\label{eqn:corr}
\end{equation}
such that in the thermodynamic limit we can express the FSS relation as
\begin{equation}
\mathfrak{O} = \beta^{-q/(z\nu)} F_{\mathfrak{O}}^{}(\beta^{1/(z\nu)}(U-U_c)),
\label{eqn:scaling_5}
\end{equation}
which is analogous to \Refs{Assaad:2013xua, Otsuka:2015iba}, apart from the scaling argument here being
$\beta$ instead of $L$, and the correlation length exponent $z\nu$ picked up a factor of the dynamical exponent $z$.

We assume the standard scaling behavior~\cite{Herbut:2009qb} for the single-particle gap
\begin{equation}
\Delta_0^{} \sim \xi_t^{-1} \sim (U-U_c)^{z\nu},
\end{equation}
as extracted from the asymptotic behavior of the correlator $C(t)$. Therefore, Eq.~\eqref{eqn:scaling_5} gives
\begin{equation}
\Delta_0^{} = \beta^{-1}
F(\beta^{1/(z\nu)}(U-U_c)),
\label{eqn:scaling_c2}
\end{equation}
as the FSS relation for $\Delta_0$. Our treatment of Eq.~\eqref{eqn:scaling_c2}
is similar to the data-collapsing procedure of Refs.~\cite{beach2005data,Campostrini:2014lta}.
We define a universal (not explicitly $\beta$-dependent), smooth function $F$ such that
\begin{align}
u & \coloneqq \beta^\mu (U-U_c^{}), \quad \mu:= 1/(z\nu), \\
f & \coloneqq \beta \Delta_0^{},\\
f & \overset{!}{=} F(u),
\label{eq:scaling_for_collapse}
\end{align}
and by adjusting the parameters $\mu$ and $U_c$, we seek to minimize the dependence of the observable $\Delta_0$ on the system size~\cite{19881,Campostrini:2014lta}, in our
case the inverse temperature $\beta$. In practice, these parameters are determined such that all points of the (appropriately scaled) gap $\Delta_0$ lie on a single line in a $u$-$f$ plot.

In this work, we have not computed the AFMI order parameter (staggered magnetization) $m_s$. However, we can make use of the findings of \Ref{Assaad:2013xua} to
provide a first estimate of the associated critical exponent $\beta$ (which we denote by $\tilde\beta$, to avoid confusion with the inverse temperature $\beta$). Specifically,
\Ref{Assaad:2013xua} found that describing $\Delta_0/U$ according to the FSS relation
\begin{equation}
\Delta_0^{}/U = \beta^{-\tilde\beta/(z\nu)}
G(\beta^{1/(z\nu)}(U-U_c)),
\label{eqn:scaling_gg}
\end{equation}
produced a scaling function $G$ which was indistinguishable from the true scaling function for $m_s$, in spite of the \textit{ansatz} $m_s \sim \Delta/U$ being a
mean-field result. In our notation, this corresponds to
\begin{align}
g \coloneqq \beta^\zeta \Delta_0^{}/U & \overset{!}{=} G(u), \quad \zeta:= \tilde\beta/(z\nu),
\label{eq:scaling_for_collapse_magn}
\end{align}
similarly to Eq.~\eqref{eq:scaling_for_collapse}.

It should be noted that Refs.~\cite{Toldin:2014sxa, Otsuka:2015iba} also considered sub-leading corrections to the FSS relation for $m_s$.
Here, data points outside of the scaling region have been omitted instead. When taken together with additional fit parameters,
they do not increase the significance of the fit. On the contrary, they reduce its stability. As a cutoff, we take $\beta U < 8$ for Eq.~\eqref{eq:scaling_for_collapse} and $\beta U < 10$ for
Eq.~\eqref{eq:scaling_for_collapse_magn}.
The reason for this limited scaling region can be understood as follows. The gap vanishes at $U=0$ regardless of $\beta$, thus $\Delta_0(u=-\beta^\mu U_c)=0$, and therefore data points with small
$\beta$ and $U$ no longer collapse onto $F(u)$. Similarly, as we show in Appendix~\ref{sec:thermal_gap}, for $U \ll U_c$ we have $\Delta_0\sim U\beta^{-2}$, which implies $g\sim\beta^{\zeta-2}\approx\beta^{-1}$, and represents a strong deviation from the expected scaling. With decreasing $\beta$, the effect increases and materialises at larger $U$.
We shall revisit the issue of sub-leading corrections in a future MC study of $m_s$.

\begin{figure}[ht]
	\input{data/gap_collapse0}
	\input{data/gap_collapse1}
	\caption{$f = \beta\Delta_0$ (left panel) and $g = \beta^\zeta\Delta_0/U$ (right panel) as a function of $u = \beta^\mu(U-U_c)$ according to Eqs.~\eqref{eq:scaling_for_collapse} and~\eqref{eq:scaling_for_collapse_magn}, respectively. All quantities are plotted with the parameters for optimal data collapse, and expressed in units of $\kappa$.}
	\label{fig:fully_collapsed_data}
\end{figure}

We perform the data-collapse analysis by first interpolating the data using a smoothed cubic spline for each
value of $\beta$. Next we calculate the squared differences between the interpolations for each pair of two different temperatures and integrate over these squared differences. Last the sum over all the integrals is minimised. The resulting optimal data collapse plots are shown
in \Figref{fully_collapsed_data}. The optimization of Eq.~\eqref{eq:scaling_for_collapse} is performed first (see \Figref{fully_collapsed_data}, left panel). This
yields $\mu=\ensuremath{\num{0.847(31)}}$, therefore $z\nu=\critNu$, and $U_c/\kappa=\ensuremath{\num{3.835(14)}}$, where the errors have been determined using parametric bootstrap.
Second, the data collapse is performed for Eq.~\eqref{eq:scaling_for_collapse_magn} (see \Figref{fully_collapsed_data}, right panel)
which gives $\zeta=\ensuremath{\num{0.761(5)}}$, such that $\tilde\beta=\ensuremath{\num{0.898(37)}}$. While the $u$-$f$ collapse is satisfactory, the $u$-$g$ collapse does not
materialize for $u \ll 0$. This is not surprising, because our scaling \textit{ansatz} does not account for the thermal gap, which we discuss in \Appref{thermal_gap}.
Moreover, these findings are consistent with the notion that the AFMI state predicted by mean-field theory is destroyed by quantum fluctuations at weak coupling, at which point
the relation $m_s \sim \Delta/U$ ceases to be valid.
Also, note that in \Ref{Assaad:2013xua} the plot corresponding to our \Figref{fully_collapsed_data} starts at $u=-2$, a region where our $u$-$g$ collapse works
out as well (for sufficiently large $\beta$).

We briefly describe our method for determining the errors of the fitted quantities and their correlation matrix. We use a mixed error propagation scheme, consisting of a parametric bootstrap part
as before, and an additional influence due to a direct propagation of the bootstrap samples of $\mu$ and $U_c$. For every $(\mu,U_c)$-sample we generate a new data set, following the
parametric sampling system. Then, the variation of the composite samples should mirror the total uncertainty of the results.
The correlation matrix\footnote{We give enough digits to yield percent-level matching to our full numerical results when inverting the correlation matrix.}
\begin{align}
	\mathrm{corr}(U_c,\, z\nu,\, \tilde\beta) &= \matr{S[round-mode=places,round-precision = 4]S[round-mode=places,round-precision = 4]S[round-mode=places,round-precision = 4]}{
	 1.00000000 & -0.01987849 & -0.2330682\\
	-0.01987849 &  1.00000000 &  0.8649537\\
	-0.23306820 &  0.86495372 &  1.0000000
		}
\end{align}
clearly shows a strong correlation between $z\nu$ and $\tilde\beta$ (due to the influence of $\mu$ on both of them), whereas $U_c$ is weakly correlated.

From the definitions of the scaling functions $F$ and $G$, one finds that they become independent of $\beta$ at $U = U_c$. Therefore, all lines in a
$U$-$f$ and $U$-$g$ plot should cross at $U=U_c$. The average of all the crossing points (between pairs of lines with different $\beta$) yields an estimate of $U_c$ as well.
We obtain $U_c^{f}/\kappa=\ensuremath{\num{3.877(46)}}$ from the crossings in the $U$-$f$ plot (see \Figref{crossing}, left panel) and $U_c^{g}/\kappa=\ensuremath{\num{3.876(79)}}$ from the $U$-$g$ plot
(see \Figref{crossing}, right panel). The errors are estimated from the standard deviation of all the crossings $\sigma_c$ and the number of crossings $n_c$ as $\sigma_c/\sqrt{n_c}$.
We find that all three results for $U_c$ are compatible within errors. However, the result from the data collapse analysis is much more precise than the crossing analysis because the data collapse analysis performs a global fit over all the MC data points, not just the points in the immediate vicinity of $U_c$.

\begin{figure}[ht]
	\input{data/gap_crossing}
	\input{data/gap_crossing1}
	\caption{Single-particle gap $\Delta_0$ scaled by $\beta$ (left panel) and by $\beta^\zeta/U$ (right panel), as a function of $U$. Note that $\Delta_0$ has been extrapolated to infinite volume
	and $\delta \to 0$ (continuum limit). Each line represents a given inverse temperature $\beta$. Here, $\zeta=\ensuremath{\num{0.761(5)}}$ as obtained from the data collapse fit.
	The averages of all the crossing points $U_c^f$ and $U_c^g$ have been marked by vertical lines. For $U_c^g$, data with $\beta\le 4$ have been omitted due to the thermal gap discussed
	in \Appref{thermal_gap}. All quantities are given in units of $\kappa$.
	\label{fig:crossing}}
\end{figure}

\subsection{Zero temperature extrapolation}\label{sec:zero temperature extrapolation}

Let us first consider the scaling properties as a function of $\beta$ for $U = U_c$. As the scaling function reduces to a constant, we have
\begin{align}
\Delta_0^{} \propto \beta^{-1}\,.
\label{eq:scaling_beta}
\end{align}
For $U > U_c$ and $\beta \to \infty$, we have $\mathfrak{O} \sim (U-U_c)^q$ by construction
and therefore
\begin{align}
\Delta_0^{} & = \begin{cases}
0 & U \le U_c \\
c_U^{} (U-U_c)^{z\nu} & U > U_c
\end{cases}\label{eq:scaling_U}
\end{align}
where $c_U^{}$
is a constant of proportionality.

We are now in a position to perform a simple extrapolation of $\Delta_0$ to $\beta \to \infty$, in order to visualize the quantum phase transition.
Let us assume that the scaling in Eq.~\eqref{eq:scaling_beta} holds approximately for $U$ slightly above $U_c$, up to a constant shift $\Delta_0^\infty$,
the zero-temperature gap. Thus we fit
\begin{align}
\Delta_0^{} & = \Delta_0^\infty + c_\beta^{}\beta^{-1},
\end{align}
at a chosen value of $U/\kappa = 4$. The resulting value of $\Delta_0^\infty/\kappa=\ensuremath{{\color{red}\Delta\Delta\Delta}}$ then allows us to determine the
coefficient $c_U = \ensuremath{{\color{red}uuu}}$ (in units of $\kappa$) of Eq.~\eqref{eq:scaling_U}. This allows us to plot the
solid black line with error bands in \Figref{delta_against_U}. Such an extrapolation should be regarded as valid only in the immediate
vicinity of the phase transition. For $U \gg U_c$ the data seem to approach the mean-field result $\tilde\beta = 1/2$~\cite{Sorella_1992}. Furthermore,
we note that an inflection point has been observed in~\Ref{Assaad:2013xua} at $U/\kappa\approx 4.1$, though this effect is thought to be an artifact of the extrapolation of the
MC data in the Trotter error $\delta$. The error estimation for the extrapolation $\beta \to \infty$ follows the same scheme as the one for the $u$-$g$ data collapse.
We obtain the error band in \Figref{delta_against_U} as the area enclosed by the two lines corresponding to the lower bound of $U_c$ and the upper bound
of $z\nu$ (on the left) and vice versa (on the right). This method conservatively captures any correlation between the different parameters.

\section{Formalism
\label{sec:formalism}}

The Hubbard model is a theory of interacting fermions, that can hop between nearest-neighbor sites. We focus on the two-dimensional honeycomb lattice,
which is bipartite in terms of $A$~sites and $B$~sites. The Hamiltonian is given by
\begin{equation}
H := - \sum_{xy}
(a\ensuremath{{}^{\dagger}}_x h_{xy}^{} a_y^{} + b\ensuremath{{}^{\dagger}}_x h_{xy}^{} b_y^{})
+ \frac{1}{2} \sum_{xy} \rho_x^{} V_{xy}^{} \rho_y^{},
\label{Ham_def}
\end{equation}
where we have applied a particle-hole transformation to a theory of spin \ensuremath{\uparrow}\ and spin \ensuremath{\downarrow}\ electrons. Here,
$a^\dag$ and $a$ are creation and annihilation operators for particles (spin-up electrons), and $b^\dag$ and $b$ are similarly
for holes (spin-down electrons). As usual, the signs of the $b$ operators have been switched for $B$~sites.
The matrix $h_{xy} := \kappa \delta_{\langle x,y \rangle}$ describes nearest-neighbor hopping, while $V_{xy}$ is the potential between particles on different sites,
and $\rho_x := b\ensuremath{{}^{\dagger}}_x b_x - a\ensuremath{{}^{\dagger}}_x a_x$ is the charge operator. In this work, we study the Hubbard model with on-site interactions only, such that $V_{xy} = U\delta_{xy}$;
the ratio $U/\kappa$ determines whether we are in a strongly or weakly coupled regime.

Hamiltonian theories such as~(\ref{Ham_def}) have for a long time been studied with lattice MC methods~\cite{Paiva2005,Beyl:2017kwp}, as this allows for
a fully \textit{ab initio} stochastic evaluation of the thermal trace, or Grassmann path integral. There is a large freedom of choice in the construction of lattice MC algorithms,
including the discretization of the theory, the choice of Hubbard-Stratonovich (or auxiliary field) transformation, and the algorithm used to update the auxiliary field variables.
This freedom can be exploited to optimize the algorithm with respect to a particular computational feature. These pertain to the scaling of the computational effort with system (lattice)
size $L$, inverse temperature $\beta$, number of time slices $N_t=\beta/\delta$, interaction strength $U/\kappa$, and electron number density (for simulations away from half filling).
Hamiltonian theories are often simulated with an exponential (or compact) form of both the kinetic and
potential energy contributions to the partition function (or Euclidean time projection amplitude), and with random Metropolis updates of the auxiliary fields (which may be either
discrete or continuous). In condensed matter and atomic physics, such methods are referred to as the Blankenbecler-Sugar-Scalapino (BSS)
algorithm~\cite{Blankenbecler:1981jt}.

In Lattice Quantum Chromodynamics (QCD), the high dimensionality of the theory and the need to precisely approach the continuum limit have led to the development of specialized algorithms
which optimize the computational scaling with $L$. These efforts have culminated in the HMC algorithm, which combines elements of the
Langevin, Molecular Dynamics (MD), and Metropolis algorithms~\cite{Duane:1987de}. The application of HMC to the Hubbard model~(\ref{Ham_def}) has proven to be
surprisingly difficult, due to problems related to ergodicity, symmetries of the Hamiltonian, and the correct approach to the (temporal) continuum limit.
For a thorough treatment of these from the point of view of HMC, see \Refs{Luu:2015gpl,Wynen:2018ryx,Ostmeyer:thesis}.
In order to realize the expected $\sim V^{5/4}$ computational scaling (where $V = 2L^2$), a suitable conjugate gradient (CG)
method has to be found for the numerical integration of the MD equations of motion.
The Hasenbusch preconditioner~\cite{hasenbusch} from Lattice QCD has recently been found to work for the Hubbard model as well~\cite{2018arXiv180407195K}.
The resulting combination of HMC with the Hubbard model is referred to as the Brower-Rebbi-Schaich (BRS) algorithm~\cite{Brower:2011av,Brower:2012zd}, which is closely related to the
BSS algorithm. The main differences are the linearized kinetic energy (or nearest-neighbor hopping) term, and the purely imaginary auxiliary field, which is updated using HMC moves in the
BRS algorithm.

The BSS algorithm has preferentially been used within the canonical ensemble~\cite{Meng2010,Assaad:2013xua,Otsuka:2015iba}. This entails a
projection Monte Carlo (PMC) calculation, where particle number is conserved and one is restricted to specific many-body Hilbert spaces.
PMC is highly efficient at accessing zero-temperature (or ground state) properties, especially when the number of particles is constant, for
instance the $A$ nucleons in an atomic nucleus. For the Hubbard model on the honeycomb lattice at half-filling, the fully anti-symmetric trial wave function
encodes a basis of $2L^2$ electrons to be propagated in Euclidean time. Due to this scaling of the number of trial wave functions, PMC and grand canonical
versions of the BSS algorithm both exhibit $\sim V^3$ scaling (with random, local Metropolis updates). In contrast to PMC, the grand canonical
formalism resides in the full Fock space, and no trial wave function is used. Instead, Boltzmann-weighted thermal expectation values of observables are extracted.
At low temperatures and large Euclidean times, spectral observables are measured relative to the ground state of the Fock space, which is the half-filling state
(an explicit example is given in \Secref{results}).
With HMC updates, such an algorithm has been found to scale as $\sim V^{5/4}$~\cite{Creutz:1988wv,2018arXiv180407195K}.
A drawback of the grand canonical ensemble is the explicit inverse temperature $\beta$. Thus, the limit $\beta \to \infty$ is taken by extrapolation, or
more specifically by finite-size scaling. This $\beta$-dependence may be considerable, though observable-dependent. While PMC simulations are
not free of similar effects (due to contamination from excited state contributions),
they are typically less severe due to the absence of backwards-propagating states in Euclidean time.

For a number of reasons, HMC updates have proven difficult for the BSS algorithm. First, the exponential form of the fermion operator $M$ causes $\det M$
to factorize into regions of positive and negative sign. Though this does not imply a sign problem at half filling (the action $S \sim |\det M|^2$),
it does introduce boundaries in the energy landscape of the theory, which HMC trajectories in general cannot cross without special and computationally very expensive
methods~\cite{Fodor:2003bh,Fodor:2004wx,Cundy:2005pi}. For $\beta \to \infty$, this fragmentation effect increases dramatically, and causes an ergodicity problem with HMC.
Second, while this problem can be circumvented by a complex-valued auxiliary field, the resulting ``complexified'' HMC algorithm shows poor (roughly cubic) scaling with $V$~\cite{Beyl:2017kwp}.
It is interesting to note how the BRS algorithm avoids this ergodicity problem. Due to the linearization of the hopping term in the fermion operator (with imaginary auxiliary field),
the boundaries impassable to HMC are reduced in dimension and can be avoided~\cite{Wynen:2018ryx}. Naturally, the BSS and BRS formulations become equivalent in the temporal
continuum limit. Then, the ergodicity problem would eventually be recovered when $\delta \to 0$ (where MC simulations are in any case not practical).
A drawback specific to BRS is that spin symmetry is explicitly broken for $\delta \neq 0$, due to the linearization of the hopping term~\cite{Buividovich:2016tgo}.
Hence, the choice of BSS versus BRS represents a tradeoff between the retention of more symmetries at finite $\delta$, and faster convergence to the
continuum limit (BSS), or improved ergodicity and computational scaling with $V$ (BRS).

Here, we apply the BRS algorithm with HMC updates to the Hubbard model~(\ref{Ham_def}). This entails the stochastic evaluation of a path integral over a Hubbard-Stratonovich
field $\phi$.
The exact form of the fermion operator $M$ depends on the choice of discretization
for time derivatives. We adopt the conventions of \Ref{Luu:2015gpl}, which used a ``mixed differencing'' operator, with forward differencing in time for $A$ sites, and backward
for $B$ sites. For reference, we note that this scheme does not introduce a fermion doubling problem. For the non-interacting theory, mixed differencing gives $\ordnung{\delta^2}$ discretization
errors (per time step). Numerically, mixed differencing has been shown to approach the limit $\delta \to 0$ faster than pure forward or backward differencing when $U > 0$.
The explicit form of $M$ is
\begin{equation}
\begin{split}
M^{AA}_{(x,t)(y,t')} & = \delta_{xy}^{} \left(\delta_{t+1,t^\prime}^{} - \delta_{t,t^\prime}^{} \exp(-i\ensuremath{\tilde{\phi}}_{x,t}^{}) \right), \\
M^{BB}_{(x,t)(y,t')} & = \delta_{xy}^{} \left(\delta_{t,t^\prime} - \delta_{t-1,t'} \exp(-i\ensuremath{\tilde{\phi}}_{x,t}^{}) \right), \\
M^{AB}_{(x,t)(y,t')} = M^{BA}_{(x,t)(y,t^\prime)} & =- \ensuremath{\tilde{\kappa}}\,\delta_{\erwartung{x,y}} \delta_{t,t^\prime},
\end{split}
\label{eqn_ferm_op}
\end{equation}
where the dependence on $A$ and $B$ sites has been written out. All quantities multiplied by $\delta$ have been denoted by a tilde. While the hopping term
in Eq.~(\ref{eqn_ferm_op}) has been linearized, the auxiliary field $\tilde\phi$ enters through exponential ``gauge links'' familiar from Lattice QCD. As first discussed
in Refs.~\cite{Brower:2011av,Brower:2012zd}, such gauge links contain a ``seagull term'', which needs to be correctly handled in order to recover the physical $\delta \to 0$ limit.
This condition is satisfied by Eq.~(\ref{eqn_ferm_op}), and further details are given in \Appref{mixed_bias}.

\section{Introduction}

The Fermi-Hubbard model---a prototypical model of electrons hopping between lattice sites---has a rich phenomenology of strongly-correlated electrons, requiring nonperturbative
treatment.
On a honeycomb lattice, where it provides a basis for studying electronic properties of carbon nanosystems like graphene and nanotubes, the Hubbard model is expected to exhibit a second-order quantum phase transition between a (weakly coupled) semi-metallic (SM) state and an antiferromagnetic Mott insulating (AFMI) state as a function of the electron-electron coupling.

It is noteworthy that a fully-controlled \emph{ab initio} characterization of the SM-AFMI transition has not yet appeared; Monte Carlo (MC) calculations have not provided unique,
generally accepted values for the critical exponents~\cite{Otsuka:2015iba,semimetalmott}. Such discrepancies are primarily attributed to the adverse scaling of MC algorithms with spatial system
size $L$ and inverse temperature $\beta$.
The resulting systematic error is magnified by an incomplete understanding of the extrapolation of operator expectation values to the thermodynamic and temporal continuum limits.

In this work we summarise the recent progress made, using Lattice Monte Carlo (LMC) techniques, in removing he systematic uncertainties which affect determinations of the critical exponents of
the SM-AFMI transition in the honeycomb Hubbard model~\cite{semimetalmott,more_observables}.
We present a unified, comprehensive, and systematically controlled treatment of both the single particle gap and the staggered magnetisation, the former being the electric and the latter the anti-ferromagnetic (AFM) order parameters of the phase transition.
We thus confirm the AFM nature of the transition from first principles. The simultaneous opening of the gap and emergence of AFM order is illustrated in figure~\ref{fig:opening_gap}.
We find that the critical coupling $U_c/\kappa=\ensuremath{\num{3.835(14)}}$ and
the critical exponents---expected to be the exponents of the SU(2) Gross-Neveu, or chiral Heisenberg, universality class~\cite{heisenberg_gross_neveu}---to be $\nu=\critNu$ and $\upbeta=\ensuremath{\num{0.898(37)}}$.

\begin{figure}[ht]
	\input{tikz/opening_gap}
	\caption{Left: The two energy bands (in multiples of the hopping $\kappa$) of the non-interacting Hubbard model as a function of the momentum $k$ normalised by the lattice spacing $a$. Center: Inset showing the Dirac cones. A band gap $\Delta$ separating the bands opens in the phase transition, once a critical coupling $U_c$ is surpassed. The bottom figure is only a qualitative visualisation, not the exact result. Right: The sublattice symmetry is broken at the same critical coupling and the disordered state (a superposition of all possibilities) transitions to an antiferromagnetic order. We show in the following (and in~\cite{semimetalmott,more_observables}) that the transitions happen simultaneously.}
	\label{fig:opening_gap}
\end{figure}

We formulate the grand canonical Hubbard model at half filling in the particle-hole basis and without a bare staggered mass.
Its Hamiltonian reads
\begin{equation}
	H=-\kappa \sum_{\erwartung{x,y}}\left(p^\dagger_{x}p^{\phantom{\dagger}}_{y}+h^\dagger_{x}h^{\phantom{\dagger}}_{y}\right)+\frac{U}{2}\sum_{x}\rho_x\rho_x\,,\label{eqn_particle_hole_hamiltonian}
\end{equation}
where $p$ and $h$ are fermionic particle and hole annihilation operators, $\kappa$ is the hopping amplitude, $U$ the on-site interaction, and
\begin{equation}
	\rho_x = p^\dagger_x p_x -h^\dagger_x h_x,
	\label{eq:charge_operator}
\end{equation}
is the charge operator.
Using Hasenbusch-accelerated HMC with the BRS formulation~\cite{2018arXiv180407195K,Brower:2012zd} and a
mixed time differencing~\cite{semimetalmott} which has favorable computational scaling properties,
we generate ensembles of auxiliary field configurations for different linear spatial extents $L$ with a maximum of $L$=$102$ corresponding to 20,808 lattice sites, interaction strengths $U$ (with fixed hopping $\kappa$), inverse temperatures $\beta = 1/T$, and $N_t$ Trotter steps;\footnote{Throughout this work, we use an upright $\upbeta$ for the critical exponent and a slanted $\beta$ for the inverse temperature.}
See \Ref{semimetalmott} for full details.
\section{Subtleties of the mixed differencing scheme
\label{sec:mixed_bias}}

As in Ref.~\cite{Luu:2015gpl}, we have used a mixed-differencing scheme for the $A$ and $B$ sublattices, which was first suggested by
Brower \textit{et al.} in Ref.~\cite{Brower:2012zd}. While mixed differencing does not cancel the linear Trotter error completely, it does diminish it significantly,
as shown in Ref.~\cite{Luu:2015gpl}. We shall now discuss some fine points related to the correct continuum limit when mixed differencing is used.

\subsection{Seagull term}

Let us consider the forward differencing used in Ref.~\cite{Luu:2015gpl} for the $AA$ contribution to fermion operator,
\begin{equation}
M^{AA}_{(x,t)(y,t^\prime)} = \delta_{xy}^{} \left\{
-\delta_{t,t^\prime}^{} + \left[\exp(i\ensuremath{\tilde{\phi}}_{x,t}^{})-\ensuremath{\tilde{m}_s}\right] \delta_{t+1,t^\prime}^{} \right\},
\label{ferm_op}
\end{equation}
with the gauge links and an explicit staggered mass $\ensuremath{\tilde{m}_s}$ (not to be confused with the AFMI order parameter) on the time-off-diagonal. We recall
that a tilde means that the corresponding quantity is multiplied by $\delta=\beta/N_t$. An expansion in $\delta$ gives
\begin{align}
M^{AA}_{(x,t)(y,t^\prime)}
& = \delta_{xy}^{} \left\{
-\delta_{t,t^\prime} + \left[1 + i\ensuremath{\tilde{\phi}}_{x,t}^{} - \frac{\ensuremath{\tilde{\phi}}_{x,t}^2}{2} - \ensuremath{\tilde{m}_s}\right] \delta_{t+1,t^\prime}^{} \right\} + \ordnung{\delta^3}, \\
& = \delta_{xy}^{} \left\{
\delta_{t+1,t^\prime} - \delta_{t,t^\prime} + \left[ i\ensuremath{\tilde{\phi}}_{x,t}^{} - \frac{\ensuremath{\tilde{\phi}}_{x,t}^2}{2} - \ensuremath{\tilde{m}_s}\right] \delta_{t+1,t^\prime}^{} \right\} + \ordnung{\delta^3}, \\
& = \delta_{xy}^{} \left\{
\delta\partial_t^{} +\frac{\delta^2\partial^2_t}{2} + \left[ i\ensuremath{\tilde{\phi}}_{x,t}^{} - \frac{\ensuremath{\tilde{\phi}}_{x,t}^2}{2} - \ensuremath{\tilde{m}_s}\right]
(\delta_{t,t^\prime}^{} + \delta\partial_t^{})\right\} + \ordnung{\delta^3}, \\
& = \delta_{xy}^{} \left\{
\delta\partial_t^{} + ( i\ensuremath{\tilde{\phi}}_{x,t}^{} - \ensuremath{\tilde{m}_s}^\prime ) \delta_{t,t^\prime} - \delta\ensuremath{\tilde{m}_s}\partial_t^{} \right\} + \ordnung{\delta^3},
\end{align}
where terms proportional to $\delta$ remain the continuum limit. In the last step we defined
\begin{equation}
\ensuremath{\tilde{m}_s}^\prime := \ensuremath{\tilde{m}_s} + \frac{\ensuremath{\tilde{\phi}}_{x,t}^2}{2} - \frac{\delta^2\partial^2_t}{2} - i\delta\ensuremath{\tilde{\phi}}_{x,t}^{}\partial_t^{},
\label{eqn_effective_staggered_mass}
\end{equation}
as the ``effective'' staggered mass at finite $\delta$.
The same calculation can be performed for the $BB$ contribution of Ref.~\cite{Luu:2015gpl}. This gives
\begin{align}
M^{BB}_{(x,t)(y,t^\prime)}
& = \delta_{xy}^{} \left\{
\delta\partial_t^{} + ( i\ensuremath{\tilde{\phi}}_{x,t}^{} + \ensuremath{\tilde{m}_s}^\prime ) \delta_{t,t^\prime} - \delta\ensuremath{\tilde{m}_s}\partial_t^{} \right\} + \ordnung{\delta^3},
\end{align}
where the effective staggered mass has the opposite sign, as expected.

Let us discuss the behavior of $\ensuremath{\tilde{m}_s}^\prime$ when $m_s \to 0$ and $\delta \to 0$. First,
$\partial_t^2$ is negative semi-definite, $\ensuremath{\tilde{\phi}}_{x,t}^2$ is positive semi-definite, and $i\ensuremath{\tilde{\phi}}_{x,t}\partial_t$ is indefinite,
but one typically finds that $\ensuremath{\tilde{m}_s}^\prime \ge \ensuremath{\tilde{m}_s}$. For a vanishing bare staggered mass $\ensuremath{\tilde{m}_s} \to 0$, this creates a non-vanishing
bias between the sublattices at $\delta \neq 0$, which is due to the mixed differencing scheme. Numerically, we find that this effect prefers $\erwartung{m_A-m_B}>0$.
Second, the ``seagull term'' $\ensuremath{\tilde{\phi}}_{x,t}^2$ is not suppressed in the continuum limit, as first noted in Ref.~\cite{Brower:2012zd}. This happens because,
in the vicinity of continuum limit, the Gaussian part of the action becomes narrow, and
$\ensuremath{\tilde{\phi}}$ is approximately distributed as $\ensuremath{\tilde{\phi}} \sim \mathcal{N}(0,\sqrt{\delta U})$. Because
$\ensuremath{\tilde{\phi}}$ scales as $\sqrt{\delta}$, the seagull term is not in $\ordnung{\delta^2}$ but effectively linear in $\delta$, denoted as $\efford{\delta}$.
The seagull term contains important physics, and should be correctly generated by the gauge links.

\subsection{Field redefinition
\label{sec_field_redef}}

Following Brower \textit{et al.} in \Ref{Brower:2012zd}, the seagull term can be absorbed by means of a redefinition of the Hubbard-Stratonovich field,
at the price of generating the so-called ``normal-ordering term'' of Ref.~\cite{Brower:2012zd}, which is of
physical significance. Let us briefly consider how this works in our case. The field redefinition of \Ref{Brower:2012zd} is
\begin{equation}
\phi_x^{} := \varphi_x^{} - \frac{\delta}{2} V_{xy}^{} \varphi_y^{} \psi^*_y \psi_y^{},
\label{phi_shift}
\end{equation}
in terms of the field $\psi$ on which the fermion matrix $M$ acts. For backward differencing of $M$, the full Hamiltonian
includes the terms
\begin{align}
& \frac{1}{2} \sum_t \tilde\phi_{x,t}^{} \tilde V_{xy}^{-1} \tilde \phi_{y,t}^{} +
\sum_{t,t^\prime} \psi^*_{x,t} \bigg[ \delta_{t,t^\prime}^{}
-\exp(-i\tilde\phi_{x,t}^{}) \delta_{t-1,t^\prime}^{} + \tilde m \delta_{t-1,t^\prime}^{} \bigg] \psi_{x,t^\prime}^{}
\nonumber \\
& \qquad =
\ { \frac{1}{2} \int dt \, \phi_x^{} V_{xy}^{-1} \phi_y^{}}
\ { + \int dt \, \psi^*_x \partial_t^{} \psi_x^{}}
\ { + \, i \int dt \, \psi^*_x \phi_x^{} \psi_x^{}}
\ { + \, m \int dt \, \psi^*_x \psi_x^{}}
\ { + \,\frac{\delta}{2} \int dt \, \psi^*_x \phi_x^2 \psi_x^{}}
+ \efford{\delta},
\label{lattice_B}
\end{align}
where the Gaussian term generated by the Hubbard-Stratonovich transformation is included, and $\tilde m$ is left unspecified for the moment.
If we apply the redefinition~\eqref{phi_shift} to the Gaussian term, we find
\begin{align}
\frac{1}{2} \int dt \, \phi_x^{} V_{xy}^{-1} \phi_y^{} & = \frac{1}{2} \int dt \, \varphi_x^{} V_{xy}^{-1} \varphi_y^{}
- \frac{\delta}{2} \int dt \, \varphi_x^{} V_{xz}^{-1} V_{zy}^{} \varphi_y^{} \psi^*_y \psi_y^{} + \efford{\delta},
\nonumber \\
& = \frac{1}{2} \int dt \, \varphi_x^{} V_{xy}^{-1} \varphi_y^{}
- \frac{\delta}{2} \int dt \, \psi^*_x \varphi_x^2 \psi_x^{} + \efford{\delta},
\end{align}
and along the lines of \Ref{Brower:2012zd}, we note that the $\varphi_x^2$ seagull term cancels the $\phi_x^2$
seagull term of Eq.~(\ref{lattice_B}) to leading order in $\delta$. As we are performing a path integral over the Hubbard-Stratonovich field,
we need to account for the Jacobian of the field redefinition, which is
\begin{align}
\int \mathcal{D}\tilde\phi =
\int \mathcal{D}\tilde \varphi \, \det \left[\frac{\partial\tilde\phi_{x,t}^{}}{\partial\tilde\varphi_{y,t^\prime}^{}}\right]
& = \int \mathcal{D}\tilde\varphi \, \exp \left[\mathrm{Tr} \log\left(\delta_{xy}^{}\delta_{t,t^\prime}^{} - \frac{\delta}{2}
V_{xy}^{} \psi^*_{y,t} \psi_{y,t}^{} \delta_{t,t^\prime}^{}\right) \right]
\nonumber \\
& \simeq \int \mathcal{D}\varphi \,
\exp\left(- \frac{1}{2} \int dt \,  V_{xx}^{} \psi^*_x \psi_x^{}\right),
\end{align}
where in the last step we used $\log(1+\delta z) = \delta z+\ordnung{\delta^2}$ before taking the continuum limit.
We conclude that the seagull term in the expansion of the gauge links has the correspondence
\begin{equation}
\frac{\delta}{2} \int dt \, \psi^*_x \phi_x^2 \psi_x^{}
\longleftrightarrow \frac{1}{2} \int dt \, V_{xx}^{} \psi^*_x \psi_x^{},
\label{seagull_norm}
\end{equation}
which is exactly the normal-ordering term proportional to $V_{xx}^{}/2$ of \Ref{Smith:2014tha}.
Hence, as argued in \Ref{Brower:2012zd}, the normal-ordering term should be omitted when gauge links are used, as an equivalent term is dynamically
generated by the gauge links.
This statement is valid when backward differencing is used for both sublattices. In \Appref{alternative derivative}, we discuss
how this argument carries over to the case of forward and mixed differencing.

\subsection{Alternative forward difference
\label{sec:alternative derivative}}

In case the backward differencing of Eq.~(\ref{lattice_B}) is used for both sublattices as in \Ref{Smith:2014tha}, then simply taking
the usual staggered mass term
\begin{equation}
\tilde m = \ensuremath{\tilde{m}_s} \quad (x \in A),
\qquad
\tilde m = -\ensuremath{\tilde{m}_s} \quad (x \in B),
\end{equation}
suffices to get the correct Hubbard Hamiltonian, as both sublattices receive a dynamically generated normal-ordering term with coefficient $V_{xx}^{}/2$.
However, the mixed-difference lattice action in \Ref{Luu:2015gpl} produces a ``staggered'' normal-ordering term, with $-V_{xx}^{}/2$ for sublattice $A$ and
$V_{xx}^{}/2$ for sublattice $B$. Hence, with the mixed-difference operator of \Ref{Luu:2015gpl} (forward for sublattice $A$, backward for sublattice $B$),
we should instead take
\begin{equation}
\tilde m = \tilde V_{00}^{} + \ensuremath{\tilde{m}_s} \quad (x \in A),
\qquad
\tilde m = -\ensuremath{\tilde{m}_s} \quad (x \in B),
\end{equation}
in order to again obtain the physical Hubbard Hamiltonian with normal-ordering and staggered mass terms.
Therefore, in our current work we adopt the alternative forward differencing
\begin{align}
M^{AA}_{(x,t)(y,t^\prime)} & = \delta_{xy}^{} \left\{ \delta_{t+1,t^\prime}^{}
- \left[\exp(-i\ensuremath{\tilde{\phi}}_{x,t}^{}) + \ensuremath{\tilde{m}_s}\right] \delta_{t,t^\prime}^{} \right\},
\label{MAA_modified}
\end{align}
instead of Eq.~\eqref{ferm_op}, which again yields a normal-ordering term $V_{xx}^{}/2$ for sublattice $A$. As in \Ref{Smith:2014tha}, we thus
retain the desirable feature of a completely dynamically generated normal-ordering term. In our actual numerical simulations, we set the bare
staggered mass $\ensuremath{\tilde{m}_s} = 0$. In our CG solver with Hasenbusch preconditioning, we work with finite $\ensuremath{\tilde{m}_s}$~\cite{hasenbusch}.
The spectrum of the operator~\eqref{MAA_modified} lacks conjugate reciprocity, which causes an ergodicity problem~\cite{Wynen:2018ryx}.

\section{Observables}\label{sec:observables}

\subsection{The single particle gap}

The single particle gap measures the distance between the two energy bands which is closest at the Dirac momenta $K$ and $K'$. It can readily be obtained from the single particle correlator
\begin{align}
	C(t) \coloneqq \erwartung{p_{K,t}^{\phantom{\dagger}} p\ensuremath{{}^{\dagger}}_{K,0}} \equiv \erwartung{p_{K',t}^{\phantom{\dagger}} p\ensuremath{{}^{\dagger}}_{K',0}}
\end{align}
at said momenta in imaginary time by a fit of the effective mass $\ensuremath{m_\text{eff}}\xspace$ to the functional form
\begin{align}
	C(\tau) = a \cosh\left(\ensuremath{m_\text{eff}}\xspace\,\delta\left(\tau-N_t/2 \right)\right),
	\label{eq:cosh fit}
\end{align}
with $a$ and $\ensuremath{m_\text{eff}}\xspace$ as fit parameters in a reasonable scaling region. Here $\tau=t/\delta$ has integer values and we symmetrised the correlator about $t=\beta/2$. For an example of such a fit see figure~\ref{fig:plateau_fit} and for more details see Ref.~\cite{semimetalmott}.
The single particle gap is simply given by
\begin{equation}
	\Delta = 2 \ensuremath{m_\text{eff}}\xspace,
\end{equation}
for a specific value of $L$ and $\delta$.

\begin{figure}[t]
	\centering
	\includegraphics[width=.45\textwidth]{{b_8_l_15_t_64_u_3.5_plateau}.pdf}
	\includegraphics[width=.45\textwidth]{{b_12_l_6_t_72_u_3.85_plateau}.pdf}
	\caption{Examples of effective mass determinations from the {single-particle correlator}s, extracted from ensembles of auxiliary field configurations.
		The blue line with error band gives $\ensuremath{m_\text{eff}}\xspace$ with statistical error, obtained from a hyperbolic cosine effective mass \eqref{eq:cosh fit}.
		The length of the blue band indicates the fitting region. For comparison, a constant fit to the effective local mass is
		shown by the dashed orange line. The dot-dashed red line shows the estimation of the systematic error~\cite{semimetalmott}.
		Note that the red and orange lines have been extended outside of the fitting region, for clearer visibility.
		Left panel: $\kappa\beta=8$, $L=15$, $N_t=64$, $U/\kappa=\num{3.5}$. Right panel: $\kappa\beta=12$, $L=6$, $N_t=72$, $U/\kappa=\num{3.85}$.
		The timeslice index $\tau$ is integer-valued. The effective mass \ensuremath{m_\text{eff}}\xspace is given in units of $\kappa$.
		\label{fig:plateau_fit}}
\end{figure}

\subsection{Spin structure factors}

At half-filling, we compute expectation values of one-point and two-point functions of bilinear local operators,
the spins
\begin{equation}
	S^i_{x} = \half (p^\dagger_x,\, (-1)^x h_x) \sigma^i (p_x,\, (-1)^x h^\dagger_x)\ensuremath{{}^{\top}}
\end{equation}
and the charge~\eqref{eq:charge_operator}
where the $\sigma^i$ are Pauli matrices and $(-1)^x$ provides a
minus sign depending on the triangular sublattice of the honeycomb to which the site $x$ belongs; the sign originates in the particle-hole transformation.
In this work we focus on the staggered sum of spins
\begin{equation}
	S^i_{-} = \sum_x (-1)^x S^i_{x},
\end{equation}
which computes the difference between the sublattices, though the uniform sum of spins and the charge operators are readily available too and have been investigated in Ref.~\cite{more_observables}. We omit them here because they are no order parameters of the phase transition of interest.

Furthermore we define the correlator
\begin{equation}
	S^{ii}_{--} = \average{S^i_- S^i_-}.
\end{equation}
Because the dynamical exponent $z=1$~\cite{PhysRevLett.97.146401} we can obtain intensive order parameters at any finite temperature by taking the thermodynamic limit of these extensive quantities and dividing them by $\beta^2$, rather than the spatial volume $V=L^2$.

\section{Finding a plateau
\label{sec:plateau_estimation}}

Here we present an automatized, deterministic method that reliably finds the optimal plateau in a given data set (such as the effective mass $m(\tau)$). Specifically, our method finds
the region of least slope and fluctuations, and checks whether this region is a genuine plateau without significant drift. If a given time series does not exhibit an acceptable plateau,
our method returns an explicit error message.

Apart from the time series $m(\tau)$ expected to exhibit a plateau, the algorithm requires two parameters to be chosen in advance.
The first is the minimal length $\lambda$ a plateau should have. The second is an ``analysis window'' of width $\mu\le\lambda$.
This controls how many data points are considered in the analysis of local fluctuations. We find that
\begin{align}
\lambda & = \frac{N_t^{}}{6}, \\
\mu & = \log_2^{}(N_t^{}),
\end{align}
are in most cases good choices.

Algorithm~\ref{alg_plateau} describes the procedure in detail. The idea is to find a balance between least statistical and systematic fluctuations.
Statistical fluctuations decrease with increasing plateau length. This is why we seek to choose the plateau as long as possible, without running into a region
with large systematic deviations. This property can also be used to our advantage. If we calculate the mean from a given time $\tau_2$ to all the
previous times, the influence of another point compatible with the mean will decrease with the distance from $\tau_2$. Thus the local fluctuation of the
running mean decreases, until it reaches a point with significant systematic deviation. This local fluctuation minimum marks the optimal $\tau_1$. The plateau
then ranges from $\tau_1$ to $\tau_2$. We check that it does not exhibit significant drift, by fitting a linear function and checking if the first order term deviates
from zero within twice its error. By repeating the analysis for all possible values of $\tau_2$, the globally best plateau can be found, as determined by
least local fluctuations of the running mean.

	\begin{algorithm}
		\SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
		\Input{$N_t$, $m[0,\dots,\,N_t-1]$, $\lambda$, $\mu$}
		\Output{$\tau_1$, $\tau_2$}
		\BlankLine
		\For{$\tau'=\mu-1,\dots,\,N_t-1$}{
			\For{$\tau=0,\dots,\, \tau'$}{
				$\overline{m}[\tau;\tau'] = \mathrm{mean}\!\left(m[\tau,\dots,\tau']\right)$\;
			}
			\For{$\tau=0,\dots,\, \tau'-\mu+1$}{
				$\sigma[\tau;\tau'] = \mathrm{sd}\!\left(\overline{m}[\tau,\dots,\tau+\mu-1;\tau']\right)$\;
			}
			$\tau_1^*[\tau'] = \underset{\tau\in \left\{0,\dots\,, \tau'-\mu\right\}}{\mathrm{argmin}} \left(\sigma[\tau;\tau']\right)$\;
		}
		$\Lambda_{\phantom{0}} = \left\{\left(\tau,\tau'\right)\,|\:\tau=\tau_1^*[\tau'],\,\tau'-\tau\ge\lambda\right\}$\;
		$\Lambda_0 = \left\{\left(\tau,\tau'\right)\in \Lambda\,|\:m[\tau,\dots,\tau'] \text{ has no significant drift}\right\}$\;
		\uIf{$\Lambda_0 \neq \varnothing$}{
			$\left(\tau_1,\tau_2\right) = \underset{(\tau,\tau')\in \Lambda_0}{\mathrm{argmin}} \left(\sigma[\tau;\tau']\right)$\;
		}\Else{
			No acceptable plateau of requested length found.
		}
		\caption{Finding a fit range for a plateau in a time series. \label{alg_plateau}}
	\end{algorithm}

As every range in the set $\Lambda_0$ from Algorithm~\ref{alg_plateau} is a valid plateau, it allows us to estimate the systematic error due to the choice of plateau.
We simply repeat the calculation of the relevant observable for all ranges in $\Lambda_0$, and interpret the standard deviation of the resulting set of values as a systematic uncertainty.

\section{The gap
\label{sec:results}}

\subsection{The single-particle correlator}

We shall now describe the procedure of obtaining the single-particle gap $\Delta$ as a function of $U/\kappa$ and $\beta$, from a given ensemble of auxiliary-field
configurations. We recall that $a$ and $a\ensuremath{{}^{\dagger}}$ are annihilation and creation operators for quasiparticles, and similarly $b$ and $b\ensuremath{{}^{\dagger}}$ for (quasi-)holes.
For instance, by creating and destroying quasiparticles at different locations and times, we obtain the correlator
\begin{equation}
C_{xy}^{}(t)
:= \left\langle a_{x,t}^{\phantom{\dagger}} a\ensuremath{{}^{\dagger}}_{y,0} \right\rangle
= \frac{1}{\ensuremath{\mathcal{Z}}} \int \ensuremath{\mathcal{D} }\phi\; M[\phi] \inverse_{(x,t),(y,0)} \exp(-S[\phi])
= \left\langle M[\phi] \inverse_{(x,t),(y,0)} \right\rangle,
\end{equation}
as an ensemble average,
where $S[\phi]$ is the Euclidean action and $\ensuremath{\mathcal{Z}}$ is the partition function---the integral without $M$.
We have used Wick contraction to replace the operators with the fermion propagator.

We now move to the Heisenberg picture and express the correlator
as a thermal trace
\begin{equation}
C_{xy}(t) = \left\langle a_{x,t}^{\phantom{\dagger}} a\ensuremath{{}^{\dagger}}_{y,0} \right\rangle
= \frac{1}{\ensuremath{\mathcal{Z}}} \tr\bigg\{a_{x,t}^{\phantom{\dagger}} a\ensuremath{{}^{\dagger}}_{y,0} \exp(-\beta H) \bigg\}
= \frac{1}{\ensuremath{\mathcal{Z}}} \tr\bigg\{\exp(-H(\beta-t)) \, a_x^{\phantom{\dagger}} \exp(-Ht) a\ensuremath{{}^{\dagger}}_y \bigg\},
\end{equation}
and by inserting the identity (resolved in the interacting basis), we find the spectral decomposition
\begin{align}
\label{eq:spectral decomposition}
C_{xy}^{}(t) = \left\langle a_{x,t}^{\phantom{\dagger}} a\ensuremath{{}^{\dagger}}_{y,0} \right\rangle & =
\frac{1}{\sum_i \exp(-\beta E_i)}
\sum_{m,n} \exp(-\beta E_m) \exp(-(E_n-E_m)t)
\, z_{mxn}^{\vphantom{*}} z\ensuremath{{}^{*}}_{myn}, \\
\label{eq:overlap factor}
z_{mxn}^{} & :=
\left\langle m \middle| a_x^{} \middle| n \right\rangle,
\end{align}
where the summation indices $m$ and $n$ denote interacting eigenstates with energies $E_m^{}$ and $E_n^{}$, respectively, and the $z_{ijk}^{}$ are referred to as ``overlap factors''.
At large $\beta$ and in the limit of large Euclidean time, the correlator decays exponentially,
\begin{equation}
\lim_{t\ensuremath{\rightarrow}\infty} \lim_{\beta\ensuremath{\rightarrow}\infty}
C_{xy}^{}(t) \simeq \exp(-(E_1-E_0)t) \, z_{0x1}^{\vphantom{*}} z\ensuremath{{}^{*}}_{0y1},
\end{equation}
where the energy $E_1^{}$ is measured relative to the energy $E_0$ of ground state of the Fock space, assuming that the associated
overlap factors are non-zero.

To continue it is convenient to Fourier transform the fermion propagator, which is a function of the lattice coordinates $x$ and $y$, to a function of the (lattice) momenta $k$ and $p$\footnote{We note that under the ensemble average the correlators are diagonal in $k$ and $p$ due to translational invariance.},
\begin{equation}
C_{kp}^{}(t)=\frac{1}{L^4}\sum_{x,y}\exp(-ikx-ipy) \, C_{xy}^{}(t)\ .
\end{equation}
Repetition of this procedure yields an ensemble of ``measurements'', the average of which is our MC estimate of $C_{kp}(t)$.
As there are two Dirac points $K$ and $K^\prime$, by symmetry we have\footnote{Because of the underlying sublattices $A$ and $B$,
there are in principle two independent correlators for each momentum $k$~\cite{Luu:2015gpl}.  However, these two correlators
are degenerate at each Dirac point $K$ and $K^\prime$; we average them to construct $C_{KK}$ and $C_{K^\prime K^\prime}$, respectively.}
\begin{equation}
C_{KK}^{}(t) = C_{K^\prime K^\prime}^{}(t),
\end{equation}
which holds for the expectation values, but not on a configuration-by-configuration basis. By choosing to Fourier transform $x$ and $y$ to the Dirac momenta $K$ or $K^\prime$, we adjust the overlap factors so that we can take
$E_1^{}$ to refer to the energy at the Dirac point. We define the ``effective mass''
\begin{equation}
\ensuremath{m_\text{eff}}\xspace := E_1^{}-E_0^{},
\end{equation}
which can be extracted from the correlator according to
\begin{equation}
\label{eq:usual effective mass}
\ensuremath{m_\text{eff}}\xspace(t) = -\lim_{\beta\ensuremath{\rightarrow}\infty} \partial_t \ln C_{KK}^{}(t),
\qquad
\ensuremath{m_\text{eff}}\xspace = \lim_{t\ensuremath{\rightarrow}\infty} \ensuremath{m_\text{eff}}\xspace(t),
\end{equation}
such that the single-particle gap is given by
\begin{equation}
\Delta = 2 \ensuremath{m_\text{eff}}\xspace,
\end{equation}
for a specific value of $L$ and $\delta$. We discuss the extrapolations in these variables in Section~\ref{sec:gap_results}.

By symmetry, the two Dirac points are indistinguishable, and the correlator expectation values are symmetric in time around
$\beta/2$ (or $\ensuremath{N_t}\xspace/2$ in units of $\tau=t/\delta$). We therefore average and fold the correlator according to
\begin{equation}
\label{eq:symmetrized}
C(t) := \frac{1}{4} \big(
C_{KK}^{}(t) + C_{K^\prime K^\prime}(t) +
C_{KK}^{}(\beta-t) + C_{K^\prime K^\prime}(\beta-t)\big),
\end{equation}
on a configuration-by-configuration basis,
which increases our numerical precision without the need for generating additional MC configurations.
In most cases, thermal effects due to backwards-propagating states cannot be completely neglected, and the isolation of an unambiguous exponential decay is difficult.
Instead of obtaining $\Delta$ from
\begin{align}
\label{eq:const fit}
\ensuremath{m_\text{eff}}\xspace(\tau)\,\delta = \ln\frac{C(\tau)}{C(\tau+1)},
\end{align}
we use the symmetry of the correlator at the Dirac point about $\tau=\beta/2$ and calculate
\begin{equation}
\label{eq:cosh fit}
\ensuremath{m_\text{eff}}\xspace(\tau)\, \delta = \cosh^{-1}\left(\frac{C(\tau+1)+C(\tau-1)}{2C(\tau)}\right)\ ,
\end{equation}
where $C(\tau)$ is the folded and symmetrized correlator \eqref{eq:symmetrized}.
This hyperbolic cosine form \eqref{eq:cosh fit} is found to be more accurate than the exponential form \eqref{eq:const fit},
as the quasiparticle masses are rather small, and the backward-propagating contributions non-negligible (due to the finite extent in $\beta$).
Once the optimal region from which to extract the effective mass has been found (for details, see \Appref{plateau_estimation}),
we fit the correlator in this region to the form
\begin{align}
C(\tau) = a \cosh\left(\ensuremath{m_\text{eff}}\xspace\,\delta\left(\tau-N_t/2 \right)\right),
\end{align}
with $a$ and $\ensuremath{m_\text{eff}}\xspace$ as fit parameters. In \Figref{plateau_fit}, we show the difference between the hyperbolic cosine fit and a direct fit of $\ensuremath{m_\text{eff}}\xspace(\tau)$. The former is shown
by a blue line and band (statistical error), and the latter by an orange dashed line. While the results are in general agreement, the direct fit of $\ensuremath{m_\text{eff}}\xspace(\tau)$ tends to overestimate the
gap for small {effective mass}es and high statistical noise levels. The estimated systematic error due to the choice of the fit range (as explained in \Appref{plateau_estimation}) is
shown by the red dot-dashed line. The statistical errors are obtained by a bootstrap procedure, and the total error has been estimated by adding the statistical and
systematic uncertainties in quadrature. The data analysis described in this and the following sections has mostly been performed in the \texttt{R} language~\cite{R_language},
in particular using the \texttt{hadron} package~\cite{hadron}.

\begin{figure}[t]
\centering
\includegraphics[width=.45\textwidth]{{b_8_l_15_t_64_u_3.5_plateau}.pdf}
\includegraphics[width=.45\textwidth]{{b_12_l_6_t_72_u_3.85_plateau}.pdf}
\caption{Examples of effective mass determinations from the {single-particle correlator}s, extracted from ensembles of auxiliary field configurations.
The blue line with error band gives $\ensuremath{m_\text{eff}}\xspace$ with statistical error, obtained from a hyperbolic cosine effective mass \eqref{eq:cosh fit}.
The length of the blue band indicates the fitting region. For comparison, a constant fit to the effective mass \eqref{eq:usual effective mass} is
shown by the dashed orange line. The dot-dashed red line shows the estimation of the systematic error, as explained in \Appref{plateau_estimation}.
Note that the red and orange lines have been extended outside of the fitting region, for clearer visibility.
Left panel: $\kappa\beta=8$, $L=15$, $N_t=64$, $U/\kappa=\num{3.5}$. Right panel: $\kappa\beta=12$, $L=6$, $N_t=72$, $U/\kappa=\num{3.85}$.
The timeslice index $\tau$ is integer-valued. The effective mass \ensuremath{m_\text{eff}}\xspace is given in units of $\kappa$.
\label{fig:plateau_fit}}
\end{figure}

\subsection{Lattice artifacts \label{sec:lattice_artefacts}}

Once we have determined the single-particle gap $\Delta$, we still need to consider the limits $\delta \to 0$, $L \to \infty$, and $\beta \to \infty$. We shall discuss the
first two limits here, and return to the issue of finite-size scaling in $\beta$ later.

In general, grand canonical MC simulations of the Hubbard model are expected to show very mild finite-size effects which vanish exponentially with $L$, as found by
\Ref{Wang:2017ohn}. This situation is more favourable than in canonical simulations, where observables typically scale as a power-law in $L^{-1}$~\cite{Wang:2017ohn}.
Let us briefly consider the findings of other recent MC studies.
For the extrapolation in $L$, \Refs{Assaad:2013xua,Meng2010}
(for $\Delta$ and the squared staggered magnetic moment $m_s^2$), \Ref{Otsuka:2015iba} (for $m_s^2$) and \Ref{Buividovich:2018yar}
(for the square of the total spin per sublattice), find a power-law dependence of the form $a+bL\inverse+cL^{-2}$.
On the other hand, \Ref{Ulybyshev_2017} found little or no dependence on $L$ for the
conductivity.
A side effect of a polynomial dependence on $L^{-1}$ is that (manifestly positive) extrapolated quantities may become
negative in the limit $L \to \infty$.

The continuous time limit $\delta \to 0$ was taken very carefully in \Ref{Otsuka:2015iba},
by simulations at successively smaller $\delta$ until the numerical results stabilized. With the exponential (or compact)
kinetic energy term used in \Ref{Otsuka:2015iba}, the Trotter error of observables should scale as $\ordnung{\delta^2}$. As shown, for instance in
\Ref{Wynen:2018ryx}, discretization errors of observables for our linearised kinetic energy operator are in general of $\ordnung\delta$. Even with an
exponential kinetic energy operator, some extrapolation in $\delta$ is usually required~\cite{Wynen:2018ryx}.
For further details concerning the limit $\delta \to 0$, see \Appref{mixed_bias}.

In \Appref{thermal_gap}, we argue that the residual modification of $\Delta$ due to the finite lattice extent $L$ should be proportional to $L^{-3}$.
This is not expected for all observables, but only for those that satisfy two conditions. First, the observable should not (for single MC configurations)
have errors proportional to $L^{-2}$, which are not cancelled by an ensemble average. This condition is satisfied by the correlation functions at the Dirac points,
but not in general for (squared) magnetic or other locally defined quantities. For example, $m_s^2$ exhibits $\sim L^{-2}$ fluctuations which are
positive for every MC configuration. As the average of these positive contributions does not vanish, the error is effectively proportional to $L^{-2}$.
Second, the correlation length $\xi$ has to fulfil $\xi \ll L$, such that the error contribution $\sim\exp(-L/\xi)$ as in \Ref{Wang:2017ohn} remains suppressed.
While we expect $L^{-3}$ scaling to hold far from phase transitions, this may break down
close to criticality, where $\xi \to \infty$. We find that deviations from inverse cubic scaling are small when
\begin{align}
L \gg \frac{\kappa\beta}{\pi},
\end{align}
where $\pi/\beta$ is the minimum Matsubara frequency dominating the correlation length ($\kappa\beta/\pi\gtrsim\xi$). We employed $\kappa \beta\le 12$ implying the requirement $L\gg 4$. Numerically, we find that for $L \ge 9$ our observed dependence on $L$ is entirely governed by inverse cubic scaling (see Figures~\ref{fig:2d_gap_fit} and~\ref{fig:finite_size_always_cubic}).

\subsection{Extrapolation method}

In this study, we have chosen to perform a simultaneous extrapolation in $\delta$ and $L$, by a two-dimensional chi-square minimization.
We use the extrapolation formula
\begin{align}
\mathfrak{O}(L,N_t^{}) =
\mathfrak{O} + a_1^{} L^{-3} + a_2^{} N_t^{-2},
\label{eqn_white_2d}
\end{align}
where $\mathfrak{O}$ is an observable with expected Trotter error of $\ordnung{\delta^2}$, and $a_1, a_2$ are fit parameters.
This is similar to the procedure of \Ref{White:1989zz} for expectation values of the Hamiltonian.
Before we describe our fitting procedure in detail, let us note some advantages of Eq.~(\ref{eqn_white_2d}) relative to a method where observables
are extracted at fixed $(U,\beta)$ by first taking the temporal continuum limit $\delta \to 0$,
\begin{align}
\mathfrak{O}(L,N_t^{}) =
\mathfrak{O}(L) + a_2^{}(L) N_t^{-2},
\label{eqn_white_1d_step1}
\end{align}
where each value of $\mathfrak{O}(L)$ and $a_2^{}(L)$ is obtained from a separate chi-square fit, where $L$ is
held fixed. This is followed by
\begin{align}
\mathfrak{O}(L) =
\mathfrak{O} + a_1^{} L^{-3},
\label{eqn_white_1d_step2}
\end{align}
as the final step. Clearly, a two-dimensional chi-square fit using Eq.~(\ref{eqn_white_2d}) involves a much larger
number of degrees of freedom relative to the number of adjustable parameters. This feature makes it easier to judge the quality
of the fit and to identify outliers, which is especially significant when the fit is to be used for extrapolation. While we have presented arguments for the
expected scaling of $\Delta$ as a function of $L$ and $\delta$, in general the true functional dependence on these variables is not \textit{a priori} known. Thus, the only unbiased check on the
extrapolation is the quality of each individual fit. This criterion is less stringent, if the data for each value of $L$ is individually extrapolated
to $\delta \to 0$.

The uncertainties of the fitted parameters have been calculated via
parametric bootstrap, which means that the bootstrap samples have been generated by drawing from independent normal distributions, defined by every single
value of $\Delta(L,N_t^{})$ and its individual error.

\subsection{Results \label{sec:gap_results}}

\begin{figure}[t]
\centering
\raisebox{0.016\height}{\includegraphics[width=.45\textwidth]{{b_8_u_3.5_m_0_gap_l}.pdf}}
\includegraphics[width=.45\textwidth]{{b_8_u_3.5_m_0_gap_t}.pdf}
\caption{Simultaneous two-dimensional fit of $\Delta(N_t,L)$ (in units of $\kappa$) using Eq.~\eqref{eqn_gap_artefacts}, for $\kappa\beta=8$ and $U/\kappa=\num{3.5}$.
Note that Eq.~\eqref{eqn_gap_artefacts} only incorporates effects of $\ordnung{L^{-3}}$ and $\ordnung{N_t^{-2}}$. Data points for $L<9$ have been omitted from the fit,
but not from the plot. Very small lattices lead to large values of $\Delta$, which are not visible on the scale of the plot.
This fit has $\chi^2/\text{d.o.f.} \simeq \num{1.1}$, corresponding to a p-value of $\simeq \num{0.34}$.
\label{fig:2d_gap_fit}}
\end{figure}

Given the fermion operator described in \Secref{formalism}, we have used Hasenbusch-accelerated HMC~\cite{hasenbusch,2018arXiv180407195K}
to generate a large number of ensembles at different values of $L$, $U/\kappa$, $\kappa\beta$, and $\kappa\delta=\kappa\beta/ N_t$.
We have used six inverse temperatures $\kappa\beta\in\{3, 4, 6, 8, 10, 12\}$\footnote{
These correspond to a highest temperature of $T\approx\SI{1.04e4}{\kelvin}$ and a lowest temperature of $T\approx\SI{2.6e3}{\kelvin}$. According to \Ref{temperature_phase_trans}
this range is well below the critical temperature $T_c\approx\SI{1.3e4}{\kelvin}$, above which the sublattice symmetry is restored.

} and a number of couplings in the range $U/\kappa \in [1.0, 5.5]$, which is expected to bracket the critical coupling $U_c/\kappa$ of the AFMI transition. For each temperature and
coupling, we scanned a large range in $L$ and $\kappa\delta$, to provide reliable extrapolations to the physical limits. These are
$L \in [3, 102]$ \footnote{
	We generated several very large lattices for a few sets of parameters to check the convergence.
	As we observed a nearly flat dependence of $\Delta$ on $L\gtrsim12$, we did most of the analysis with medium-sized lattices ($L\le 21$) in order to conserve
	computation time.
} and $\kappa\delta$ from~$1/4$ down to~$1/40$. Our values of $L$ were chosen to be 0 (mod 3) so that the
non-interacting dispersion has momentum modes exactly at the Dirac points. In other words, we select the lattice geometry to ensure that $E_1-E_0 = 0$ in the absence of
an interaction-induced (Mott) gap.

Let us briefly summarize our procedure for computing $\Delta$. We have calculated the single-particle correlator\ $C(t)$
as an expectation value of the inverse fermion matrix, at both independent Dirac-points $K, K^\prime$, and averaged over them to increase our statistics.
For each set of parameters simulated, those correlators have been fitted with an effective mass in order to extract $\Delta$.
We are then in a position to take the limits $\delta \to 0$ and $L \to \infty$, which we accomplish by a simultaneous fit to the functional form
\begin{align}
\Delta^2(L,N_t^{}) =
\Delta_0^2 + c_0^{} N_t^{-2}
+ c_1^{} L^{-3}\,,
\label{eqn_gap_artefacts}
\end{align}
where the fitted quantities are the (infinite-volume, temporal continuum) gap $\Delta_0$, and the leading corrections proportional to $c_0^{}$ and $c_1^{}$.
Around $\delta = 0$, we find the leading $L$-dependence
\begin{align}
\Delta(L) = \Delta_0^{} + \frac{c_1^{}}{2\Delta_0^{}}L^{-3} + \ordnung{L^{-6}},
\end{align}
where $c_1^{}$ is numerically found to be very small. Also, around $L = \infty$, we have
\begin{align}
\Delta(N_t^{}) = \Delta_0^{} + \frac{c_0^{}}{2\Delta_0^{}} N_t^{-2} + \ordnung{N_t^{-4}}\,.
\end{align}
It should be noted that inclusion of a term of order $N_t^{-1}$ did not improve the quality of the fit. This observation
can be justified by the expected suppression of the linear term due to our mixed differencing scheme.
Though the term of $\ordnung\delta$ is not removed analytically, it appears small enough to be numerically unresolvable.

An example of such a fit is shown in \Figref{2d_gap_fit}. We find that Eq.~\eqref{eqn_gap_artefacts} describes our data well with only a minimal set of parameters for large $L$ and
small $\delta$. Lattices with $L\le6$ have been omitted from the fit, as such data does not always lie in the scaling region where the data shows
$\ordnung{L^{-3}}$ convergence, as explained in \Secref{lattice_artefacts}.
Our results for $\Delta_0$ are shown
in \Figref{delta_against_U} for all values of $U/\kappa$ and $\kappa\beta$, along with an extrapolation (with error band) to zero temperature ($\beta\to\infty$).
For details as to the zero-temperature gap, see \Secref{analysis}.

\begin{figure}[t]
\input{data/gap}
\caption{The single-particle gap $\Delta_0(U,\beta)$, with all quantities in units of $\kappa$, after the thermodynamic and continuum limit extrapolations.
We also show $\Delta_0(U,\beta \to \infty)$ as a solid black line with error band (see \Secref{zero temperature extrapolation}). For
$U< U_c \simeq \ensuremath{\num{3.835(14)}}$ the zero-temperature gap vanishes.
\label{fig:delta_against_U}}
\end{figure}

\section{Thermal gap}\label{sec:thermal_gap}

It is useful to consider the influence of the inverse temperature $\beta$ on the single-particle gap, in order to provide a better understanding
of the scaling of $\Delta$ with $\beta$. Naturally, we are not able to solve the entire problem analytically, so we shall consider small
perturbations in the coupling $U$, and assume that the dispersion relation of graphene is not significantly perturbed by the interaction (which is
expected to be the case when $U$ is small). Let us now compute the expectation value of the number of electrons excited from the ground state.
As we consider exclusively the conduction band, we assume that the particle density follows Fermi-Dirac statistics.  We take the positive-energy part of
\begin{equation}
\omega_k^{} := \kappa\sqrt{\tilde\omega_k^2},
\qquad
\tilde\omega_k^2 =
3 + 4 \cos(3a k_x^{}/2) \cos(\sqrt{3}a k_y^{}/2) + 2\cos(\sqrt{3} a k_y^{}),
\end{equation}
where $a \simeq 1.42$~{\AA} is the nearest-neighbor lattice spacing and we assume that every excited electron contributes an energy $E(U)$ to a ``thermal gap'' $\Delta(\beta)$.
These considerations yield the gap equation
\begin{equation}
\Delta(\beta) = E(U) a^2 f_{\text{BZ}}^{} \int_{k\in\text{BZ}}
\frac{d^2k}{(2\pi)^2} \frac{1}{1+\exp(\beta\omega_k^{})},
\qquad
f_{\text{BZ}}^{} := \frac{3\sqrt{3}}{2},
\label{eqn_therm_gap_inf_L}
\end{equation}
where the factor $f_{\text{BZ}}^{}$ is due to the
hexagonal geometry of the first Brillouin zone (BZ). It should be emphasized that the thermal gap is not the interaction-driven Mott gap we are studying here,
even though it is not numerically distinguishable from the latter. A thermal gap can occur even if the conduction and valence bands touch or overlap.
The physical interpretation of the thermal gap (as explained above) is a measure of the degree of
excitation above the ground state, based on the number of excited states that are already occupied in thermal equilibrium.

\subsection{Finite temperature}

\begin{figure}[t]
	\input{data/gap_norm}
	\caption{Illustration of the thermal gap in the weakly coupled regime, as given by Eq.~\eqref{eqn:thermal_gap}. Our MC data for the single-particle gap $\Delta$ from
	\Figref{delta_against_U} is shown multiplied by $\beta^2$. All quantities are expressed in appropriate units of $\kappa$.}
	\label{fig:beta_sq_delta_against_U}
\end{figure}

Let us evaluate Eq.~(\ref{eqn_therm_gap_inf_L}) under the assumption that $\beta$ is large.
Then, the integrand only contributes in the region where $\omega_k^{} \approx 0$, in other words
near the Dirac points $K$ and $K^\prime$ located at momenta $k_D^{}$.
In the vicinity of a Dirac point, the dispersion relation reduced to the well-known Dirac cone
\begin{equation}
\omega_k^{} \simeq v_F^{} | k-k_D^{} |,
\qquad
v_F^{} := 3\kappa a/2,
\end{equation}
with Fermi velocity $v_F$. Within this approximation, we may
sum over the two Dirac points and perform the angular integral, which gives
\begin{align}
\Delta(\beta) & \approx
2 E(U) a^2 f_{\text{BZ}}^{} \int_{k\in\mathbb{R}^2}
\frac{d^2k}{(2\pi)^2}
\frac{1}{1+\exp(\beta v_F^{} | k | )}
= E(U) a^2 f_{\text{BZ}}^{}
\int_{0}^\infty \frac{dk}{\pi}
\frac{k}{1+\exp(\beta v_F^{} k)}, \\
\label{eqn:thermal_gap}
& = E(U) a^2 f_{\text{BZ}}^{} \times \frac{\pi}{12} (\beta v_F^{})^{-2}
= \frac{\sqrt{3}\pi}{18} \, E(U) (\beta\kappa)^{-2}
\approx \num{0.3023} \, E(U) (\beta\kappa)^{-2},
\end{align}
where the Fermi-Dirac integral has been evaluated in terms of the polylogarithm function. We expect the
error of this approximation to be exponentially suppressed in $\beta$.

In \Figref{beta_sq_delta_against_U}, we validate Eq.~\eqref{eqn:thermal_gap} using our data for $\Delta$ shown in \Figref{delta_against_U}.
We find that the prediction of quadratic scaling in $\beta$ and the linear approximation in $U$ are quite accurate.
A fit of Eq.~\eqref{eqn:thermal_gap} to our MC data for $\Delta$ in the weak-coupling regime $U/\kappa \le 2$ gives
\begin{align}
E(U) = \num{5.95\pm.15}\,U,
\end{align}
under the (perturbative) assumption $E(U)\propto U$. This fit is shown in \Figref{beta_sq_delta_against_U}, where we plot $\beta^2\Delta$ as a function
of $U$ (in proper units of $\kappa$). In the weakly coupled regime, the MC data for $\beta^2\Delta$ coincide and fall on a straight line. Once the critical coupling $U_c$ is
approached, the points for various $\beta$ separate. As expected, the linear dependence on $U$ persists longest for small $\beta$, as temperature effects
are dominant over interaction effects. The quadratic scaling with $\beta$ is most accurate for large $\beta$, which is in line with the expected
exponential convergence stated above.

\subsection{Finite lattice size}

Let us also consider the leading correction to the thermal gap due to finite lattice size $L$.
The discretized form of Eq.~(\ref{eqn_therm_gap_inf_L}) is
\begin{equation}
\Delta(L,\beta)=\frac{E(U)}{L^2}\sum_{k\in \text{BZ}}\frac{1}{1+\exp(\beta\omega_k^{})}
\label{eqn_therm_gap_discrete}
\end{equation}
which in general has a very complicated convergence behavior.
Because of periodic boundary conditions, Eq.~(\ref{eqn_therm_gap_discrete}) is
an effective trapezoidal approximation to Eq.~(\ref{eqn_therm_gap_inf_L}), thus the
convergence is \textit{a priori} expected to scale as $\ordnung{L^{-2}}$.

We shall now obtain a precise leading-order error estimation. Let us discretize the
first BZ is discretized into a regular triangular lattice, with lattice spacing $h\propto L^{-1}$.
We integrate a function $f(x,y)$ over a single triangle, spanned by the coordinates $(\pm h/2, 0)$
and $(0, \sqrt{3}h/2)$,
\begin{equation}
I := \int_{-h/2}^{h/2} dx \int_0^{b(x)} dy \, f(x,y),
\qquad
b(x) := \sqrt{3}h/2-\sqrt{3} | x |,
\end{equation}
and subtract the average of $f(x,y)$ over the corner points multiplied by the area of the triangle,
\begin{equation}
\hat{I} := \frac{\sqrt{3}h^2}{4} \times \frac{1}{3}
\left[ f(-h/2,0)+f(h/2,0)+f(0,\sqrt{3}h/2) \right],
\end{equation}
which gives the (local) error
\begin{equation}
\delta I :=
I-\hat{I}=-\frac{\sqrt{3}}{64} \left[\delsqu{f}{x}(0)+\delsqu{f}{y}(0)\right] h^4 + \ordnung{h^5},
\end{equation}
due to discretization. The global error is obtained by summing over the complete BZ,
\begin{align}
\sum_{k\in \text{BZ}} \delta I(k) & =
-\frac{\sqrt{3}}{64}\sum_{k\in \text{BZ}} \left[\delsqu{f}{x}(k)+\delsqu{f}{y}(k)\right] h^4 + \ordnung{L^2h^5}, \\
&\propto \frac{1}{L^4} \sum_{k\in \text{BZ}} \left[\delsqu{f}{x}(k)+\delsqu{f}{y}(k)\right] + \ordnung{L^{-3}}, \\
&\propto \frac{1}{L^2} \int_{k\in \text{BZ}} d^2k \left[\delsqu{f}{x}(k)+\delsqu{f}{y}(k)\right] + \ordnung{L^{-3}},
\end{align}
which equals
\begin{align}
\sum_{k\in \text{BZ}} \delta I(k) & \propto
\frac{1}{L^2} \oint_{k\in \partial\text{BZ}} \nabla f(k) \cdot d\vec k + \ordnung{L^{-3}},
\label{eqn_integral_1st_bz} \\
& \propto \ordnung{L^{-3}},
\end{align}
where Gauss's theorem has been applied in Eq.~\eqref{eqn_integral_1st_bz}. Hence,
the projection of $\nabla f$ onto the normal of the BZ is integrated over the boundary of the BZ.
As every momentum-periodic function takes the same values on the opposite edges of the BZ, the result sums up to zero.
Surprisingly, one then finds that the second order error term in $L$ vanishes. For the special case of
$f(k)\propto 1/(1+\exp(\beta\omega_k^{}))$, the gradient in BZ-normal direction vanishes everywhere on the boundary, and the integral in Eq.~\eqref{eqn_integral_1st_bz}
is trivially zero.

\begin{figure}[t]
\centering
\includegraphics[width=.45\textwidth]{{b_10_u_3.9_m_0_gap_l}.pdf}
\includegraphics[width=.45\textwidth]{{b_4_u_5_m_0_gap_l}.pdf}
\caption{Simultaneous two-dimensional fit of $\Delta(N_t,L)$ (in units of $\kappa$) using Eq.~\eqref{eqn_gap_artefacts},
for $\kappa\beta=10$ and $U/\kappa=\num{3.9}$ (left panel) and $\kappa\beta=4$ and $U/\kappa=\num{5.0}$ (right panel). Only the extrapolations in $L$ are shown.
Data points for $L<9$ have been omitted from the fits, but not from the plots.
These fits have $\chi^2/\text{d.o.f.} \simeq \num{0.83}$ and p-value of $\simeq \num{0.62}$ (left panel), and
$\chi^2/\text{d.o.f.} \simeq \num{1.1}$ and p-value of $\simeq \num{0.36}$ (right panel).
\label{fig:finite_size_always_cubic}}
\end{figure}

Higher orders in $U$ influencing the thermal gap are not as easy to calculate, but can in principle be dealt with using
diagrammatic techniques in a finite-temperature Matsubara formalism~\cite{bruus_qft}. As an example, we
know from \Ref{Giuliani2009} that $v_F^{}$ is influenced (at weak coupling) only at $\ordnung{U^2}$.
Let us finally provide some further numerical evidence for the expected cubic finite-size effects in $L$.
In \Figref{2d_gap_fit}, we have already shown that cubic finite-size effects are a good approximation for $U < U_c$, as expected for
states with small correlation lengths. In \Figref{finite_size_always_cubic}, we show that the cubic behavior in $L$ still holds for $U \simeq U_c$
and $U > U_c$.

redpajama_set_name: RedPajamaArXiv


I'm absolutely delighted to announce that my debut collection of short stories, This Taste for Silence, will be published in mid-2019 by University of Queensland Press (UQP). It will be great to see all my current stories - new, and previously published - in the one book. UQP has a brilliant reputation and I'm proud to be joining them.


Details of the menu will be released closer to the time. Expect nothing less than a Christmas menu with all the most delicious trimmings and equally fabulous vegetarian alternatives, using only the very best local produce.
Somerset based rock and pop band, The IOUs, will be back performing live on the night, playing a wide and varied range of popular music to keep you entertained throughout the evening, so don't forget those dancing shoes!
Call us for more details and to book on 01460 249758.
We look forward to welcoming you to a very festive Frogmary Green Farm this Christmas.
Join us for lots of festive cheer, delicious food and live music! For more information, call us on 01460 249758.


Be warned. This is one of those rambling, little-bit-of-everything posts.
First the boats… which were a theme of the League of Vermont Writers summer picnic at the Lake Champlain Maritime Museum last weekend. The League invited me to be one of two speakers for the event, and we had a great time talking about research for both historical fiction and for more contemporary pieces. I was so glad I stayed to hear poet Daniel Lusk speak after lunch, too.
He was charming and read from his latest work-in-progress, a collection of poems inspired by what's under the surface of Lake Champlain, from fish to shipwrecks, to our elusive lake monster, Champ. And for those of you planning to attend the Burlington Book Festival this fall, take note… Daniel will be reading from his work at Burlington's Waterfront Theatre on September 27th. He wasn't sure what time yet, but if you check the festival website, I'm sure they'll have a schedule posted before long.
In other news this week, my husband and I have spent a fair amount of time rescuing our dock pieces from waves as Lake Champlain's water level continues to rise with all our summer rain.
I think naked docks look funny, don't you?
After we finished lugging all the pieces around, I plopped down on the sun porch and devoured my ARC of this book in a single day.
Don't miss this one when it comes out in October. L.K. Madigan's FLASH BURNOUT has one of the best teen boy voices in any book I've ever read. I can't begin to tell you how many times the dialogue made me laugh out loud. So, so perfect and funny. But not fluffy either – this book has a way of being hysterically funny and deeply profound all at once, and it tackles some tough issues. I loved it and was sad when I finished. Up next on my reading list (and saving me from too much pouting) is THE EVOLUTION OF CALPURNIA TATE, which I've heard great things about, too.
I haven't had much reading time today yet because I spent the afternoon at the public library in Alburg, Vermont for a presentation about my second Lake Champlain historical novel, CHAMPLAIN AND THE SILENT ONE. No matter how many presentations I give, I'm always moved by a welcome like this.
Those balloons have had me smiling all day. And I loved the group of people who came out for this event, a mix of all ages, including a high school freshman who just had her own poem about the Holocaust published. It was fun to talk about writing with her and her sister!
My other big highlight of the week was a research appointment to help me write a new scene in SUGAR ON SNOW, my middle grade ice skating novel that's coming out with Walker Books in 2010. The main character's friend takes care of the bees in her family's apple orchard, and I've decided to set the new scene there. (If you read the book later on, you will see why – it just fits the mood of the conversation perfectly.) So yesterday, in 87-degree weather, I met up with a delightful area beekeeper (thank you, Dave!), zipped myself into a borrowed beekeeper's outfit, and got to work watching and listening to 100,000 honeybees and learning how to care for them.
It's probably good that you can't see in the photo how sweaty I was. Despite the heat, it was a fantastic afternoon, and I learned some amazing things about honeybees, some of which are perfect for the new scene and some of which are fascinating but don't belong in my book, and so I'll just quietly enjoy those. For those of you who write & research, it's a challenge, isn't it? Knowing what has to be left out, even when it's information you love?
The Alburg Library was my last author event until Labor Day weekend, when we're throwing a launch party for THE BRILLIANT FALL OF GIANNA Z. More on that soon, but for now, know that you're invited to Flying Pig Bookstore in Shelburne, VT on Saturday, September 5th at 11:00. Nonna's famous funeral cookies from the book will make their debut, too!
For now…it's back to the revision, with some camping, family visits, and general summer frolicking sprinkled in, too. I hope you've had a terrific July!


Congratulations, you've just accepted an offer for the job you've been chasing for months. You're excited and relieved...but how is your boss is going to react when you resign...?
Once upon a time, you handed in your notice and looked forward to your leaving party; unfortunately it's not quite that simple anymore. Today when you hand in your notice, there is a very realistic chance that you will receive a counter offer; but beware, once the initial flattery and extra cash has lost its appeal, you may still be left feeling exactly the same as when you started looking for a new job.
What inconvenient timing, I was about to go on annual leave for 3 weeks!
How am I going to look in front of the board/senior management for losing this person?
How can I get him/her to stay?
'Loving' you liked they've never 'loved' you before – be suspicious of this.
Do your reasons for wanting to leave still exist? You may have a number of reasons; salary was too low, no promotion in sight, don't like your boss. You may be offered more money to stay, which can be tempting, but if you still have other issues outstanding, you'll probably end up leaving anyway.
Before accepting a counter offer, ask yourself why your employer has made the offer. There is a strong possibility that the cons will outweigh the pros and you will realise that your decision to resign was right after all.


text:
Pretty vacant sex pistols hqcl in Tennessee
Dirty sex acts funny names in Sunnyvale
Most liberal states for sex offenders in Maitland
Charles darrell toliver sex offender in Fresno
Hung tin ming wife sexual dysfunction in McKinney
Bryan anthony mouton sex offender in Baltimore
Civil commitment of sex offenders new york in Exeter
Proponents argue that such provisions offer an important community protection safeguard by incapacitating a high risk subgroup of sex offenders. What happens if someone, who is under strict and intensive supervision and treatment SISTviolates the conditions of his or her parole supervision?
After the case is referred to the Office of the Attorney General by the case review team, the Attorney General reviews the case, taking into consideration any supervision that the offender will be under after he or she is released from prison. The Office of Mental Health and the Attorney General's staff are able to civil commitment of sex offenders new york in Exeter records concerning the sex offender, but these records are confidential and cannot be made available to the public.
A program should establish internal and external multi-disciplinary bodies to provide oversight of program, staff, and client issues.
ATSA recommends that states provide a contemporary, properly designed prison-based treatment program to enable incarcerated offenders to receive treatment in order to promote risk-reduction and successful community reintegration.
Most of the sex offenders are represented by Mental Hygiene Legal Service. Risk assessment is a process for determining the likelihood of reoffense, utilizing specific factors, with specific scoring rules, that places a specific sex offender within a certain risk group. However, the sex offender may make a request that the case be moved to the court where the sex offender was convicted of the sex offense.
Under the supervision of the Division of Parole, the individual must submit to sex offender treatment and other conditions. What happens after the civil commitment of sex offenders new york in Exeter is filed?
Ждать civil commitment of sex offenders new york in Exeter этим столкнулся
If he follows the rules for the next 10 years, he can apply to be removed from the registry. Fifteen minutes later, Gamer Gurl replied that she was The vulgarities and snide tone seemed too adult. But one donor stands out. By contrast, some other states, such as Minnesota, are facing legal challenges to their civil commitment laws because they lack regular evaluations of offenders.
He was in shock, naked. Or a sly adult?
This website uses cookies to give you the best experience when you visit our website. By continuing to use this website, you consent to our use of these cookies.
In New York, some persons who have committed sex offenses are considered too dangerous to live the community. These individuals may be held indefinitely in secure institutions under laws allowing for civil commitment.
The final decision of where the case should be heard is up to the judge. First, the legal mechanism by which the offenders are detained when civilly committed depends on clinical criteria primarily created or defined by legislative bodies rather than by the scientific or mental health communities.
Typically, these laws provide a legal mechanism for the confinement of a limited number of adult sexual offenders in a secure treatment facility after incarceration when a court determines they are likely to engage in future acts of sexual violence.
The court holds a hearing after the petition is filed, without a jury, to determine if there is probable cause to believe that the sex offender is a sex offender requiring civil management. The sex offender cannot appeal the probable cause finding.
jackson county missouri sex offenders list in North Yorkshire
Aug 26, · Offenders were spending more time in prison, and the number of offenders on the sex registry increased. In Florida, sentences are often shorter. Peter Aiken, a Florida defense lawyer, has. Apr 03, · In New York, civil management of sex offenders is codified in Article 10 of the Mental Hygiene Law. It took effect in as part of the Sex Offender Management and Treatment Act ("SOMTA"). A "dangerous sex offender requiring confinement" is a person who committed one or more sex offenses, who is suffering from a mental abnormality.
are same sex marriages legal in oregon in Grand Rapids
A "dangerous sex offender requiring confinement" is a person who committed one or more sex offenses, who is suffering from a mental. At an annual cost of $, per civilly-committed sex offender, New York's civil commitment program is the second most expensive in the.
same sex marriage cakes in Charters Towers
Last year, the Supreme Court of the United States upheld a Kansas law that permits the government to detain for indefinite periods of time persons who have. Frequently Asked Questions regarding New York State's Civil Management Laws. State law permits civil management of certain sex offenders after they have completed their What is a dangerous sex offender requiring confinement?
ask sue sex education in Bunbury
have been enacted in response to crimes, committed in the late s and early s, with violation of sex offenders' constitutional rights (American Civil Liberties Union of The New York State sex offender registration act was created in educational facilities, and health facilities among others (Witten, Exeter, &. Typically, these laws provide a legal mechanism for the confinement of a limited number of adult sexual offenders in a secure treatment facility after incarceration.
love sex and candy lyrics in Cincinnati
According to the New York Times, it costs about $, per resident west, Ohio, does not have a civil commitment statute for sex offenders. nd Street, 15thFloor, New York, NY The Civil Commitment of Sexual Predators: A Policy Review. Andrew J. Harris.
single sex classrooms facts in Jersey City 50231 | 50232 | 50233 | 50234 | 50235 robin thicke sex therapy album free download in Newcastle
I am a sex addict film online in Québec
Benefits of sex offender registry in Busselton
Michael boone sex offender in Lewisville
Darsteller sex and the city der film sissi in Vaughn
Michelle sex and the city in Oshawa
Female sex offender registry ny in Pittsburgh
Public registry of sex offenders in Gladstone-Tannum Sands
Taslima nasrin sex scandal news in Pickering
Single sex schools test scores in Escondido
Opposite sex friendship boundaries in Leonora
Dramaalert runescape sex offender returns for sale in Saguenay
Loving v. virginia same sex marriage in Escondido
Same sex parenting usa in Milton
Single sex education pros and con in Stamford
Sam consiglio sex predator in Geelong
Same sex domestic violence cases in washington in Oklahoma City
Do gingers have higher sex drives in Esperance
Have better sex my husband in Bunbury
Berger sex and the city imdb in Savannah
Registered sex offenders in lansing michigan in South Carolina
Wa state sentencing guidelines for juvenile sex offenders in Darlington
Masters of sex recap involuntary head in Blainville
Federally registered sex offenders in Whyalla
Thomas franklin mullins sex offender in Irving
Christian sex therapy london in Mississippi
Nd attorney general sex offender registration in Corona
Shawn thornton sex registry in New Westminster

redpajama_set_name: RedPajamaCommonCrawl


text:
Valli Hilaire — April 9, 2008
Ok, so I'm really excited about this. I was given a copy of the new 2-disc DVD set Daytona 500: 50 Years of the "Great American Race."
I watched it last night and it gives you some really cool insight into the history, legacy and champions of the Daytona 500, NASCAR's biggest and most significant race as we all know. The first disc is an interactive experience that allows you to watch this year's Daytona 500 from the driver's perspective. It's pretty trippy to feel like you're in the car with them. The featured drivers include Kasey Kahne, Tony Stewart, Denny Hamlin, Dale Earnhardt Jr. and of course the race winner Ryan Newman. You get to hear the team radio communications and you also have the option of switching to the broadcast coverage as well.
The second disc gives you a round table discussion (or square table as Michael Waltrip put it) with past Daytona 500 winners like David Pearson, Junior Johnson and Mario Andretti to name a few. The conversation is great and full of little interesting anecdotes, plus great footage from the archives of the race's history. Oh and I just want to say that it's pretty uncanny how much IndyCar driver Marco Andretti looks like his grandfather Mario when he was in his twenties. But I digress. The second disc also features some interviews with drivers from today including Jeff Gordon, Jimmie Johnson and Kevin Harvick.
I highly recommend this DVD set to anyone that is a serious NASCAR fan. It's a cool way to experience the race all over again. This DVD set is available for purchase online and you can also win a copy of the set from me! I was given an extra copy to give away to one of my awesome readers! Click here to enter the contest. The contest is over! Thanks!
The 50th running of the Daytona 500 is one of NASCAR®'s most unforgettable and highly-anticipated milestone events. Already the most prestigious race in motorsports, this landmark running is sure to be remembered and cherished by millions of race fans. On April 15, join A&E Home Video and NASCAR® Collections to celebrate the golden anniversary of "The Great American Race™" with DAYTONA 500: THE 50TH ANNIVERSARY, a highly-collectible 2-DVD set that puts fans in the driver's seat like never before! Available just weeks after the 50th Daytona® 500 (February 17), this is THE official NASCAR® collectible of this historic race.
As part of DAYTONA® 500: THE 50TH ANNIVERSARY, the 2008 Daytona 500 will be presented in the ground-breaking "My Race" Interactive DVD format. This innovative user experience allows fans to choose the view of the race from six different driver perspectives plus the main broadcast feed, along with driver and announcer audio feeds — putting them right in the middle of the action. The second disc features the great stories and heroes that have defined the Daytona 500 through two exclusive NASCAR® programs: "A Sports Century" features a dynamic roundtable discussion of the moments and personalities that defined Daytona with a host and panel of NASCAR® experts including Richard Petty and Jeff Gordon, while "The Top 50 Daytona® 500 Moments" is the definitive telling of the story of NASCAR®'s Greatest single greatest event. Additional DVD extras include official NASCAR® PC wallpapers and time-lapse screensavers.
With a history of legend-defining and myth-making spectaculars, the Daytona® 500 winners circle has embraced NASCAR® heroes including: Richard Petty, Mario Andretti, Jeff Gordon, and Dale Earnhardt Senior and Junior. And DAYTONA® 500: THE 50TH ANNIVERSARY, an unprecedented and revolutionary DVD release, only adds to the legend of NASCAR®'s biggest race of the year.
Tags: ContestsDale Earnhardt Jr.David PearsonDenny HamlinJeff GordonJimmie JohnsonJunior JohnsonKasey KahneKevin HarvickMarco AndrettiMario AndrettiMichael WaltripRyan NewmanTony Stewart
Next post when tony met kid rock
Previous post give a little bit: autism awareness month