\section{Introduction}
\label{sec:introduction}

The use of identically prepared systems, or ensembles, has been
essential to our understanding of equilibrium and far-from-equilibrium
properties of classical and quantum systems. Traditionally, three
types of ensembles are used---(a) the microcanonical ensemble, which
involves systems with fixed energy and particle number, (b) the
canonical ensemble (CE), which involves systems with fixed particle number
in contact with a large reservoir (at temperature $T$) with which they
can exchange energy, and (c) the grand canonical ensemble (GE), which
involves systems in contact with a reservoir with which they can
exchange energy and particles (in equilibrium, the average particle
number is determined by the chemical potential $\mu$). Whereas these
three ensembles pose fundamentally different physical constraints,
it can be shown that they are equivalent in the thermodynamic limit
(provided, of course, that temperatures and chemical potentials are
selected appropriately). Being technically easier to deal with, the
canonical and grand canonical ensembles are the most commonly used
ensembles in the literature. Several texts on statistical mechanics
cover these topics in detail; see e.g.,~Ref.~\cite{huang87}.

In finite systems, differences appear between calculations carried
out using the three ensembles. These differences, dubbed finite-size
effects, have to do with the effect of energy and particle number
fluctuations, and with boundary effects. For example, to describe an
isolated system with mean energy $E$, it is most appropriate to
use the microcanonical ensemble with that energy. However,
one can also use a canonical ensemble at a temperature $T$ for which
the mean energy is $E$. Since the systems used to construct the
canonical ensemble have different energies from the ones used to
construct the microcanonical ensemble, one finds differences in the
predictions of each ensemble.  Remarkably, one can show that energy
fluctuations in the canonical ensemble typically scale as the square
root of the volume of the system, whereas the average energy scales
as the volume of the system.  Hence, the ratio between the energy
fluctuations and the average energy scales as the inverse of the
square root of the volume, and vanishes in the thermodynamic
limit. One then finds that differences between the predictions of
each ensemble decrease polynomially with increasing volume
(at fixed density). The same applies if one considers the grand
canonical ensemble, where particle number fluctuations typically
scale with the square root of the volume of the system. Indeed,
explicit calculations in one dimensional (1D) lattices have shown
that the differences between the predictions of the canonical and
grand-canonical ensembles for various observables decrease with the
inverse of the number of particles (or lattice sites) in the system
\cite{lebowitz61,rigol_05}.

Experiments usually deal with thermodynamically large systems, whereas
numerical analyses of many-body interacting systems can generally be
done for only (relatively) much smaller system sizes.  Hence, when
trying to theoretically predict/reproduce the outcome an experimental
measurement, a question of much relevance is: \emph{which ensemble should
one use to minimize finite size effects and obtain the ``thermodynamic
limit'', or, experimental result?} From the previous discussion about
the differences between ensembles, one might naively conclude that finite
size effects always scale polynomially with system size and that, therefore,
the best one can do theoretically is to optimize exponents and prefactors.

In this article we show that this is not the case. There is a
preferred ensemble (the grand canonical ensemble) and preferred boundary
conditions (periodic boundary conditions, so that the system is translationally
invariant) for which finite-size effects are exponentially small in the
system size. This holds if the system of interest is in an unordered
(i.e., without long or quasi-long range order) phase at finite temperature.
We also consider a different approach to calculating finite-temperature
properties of many-particle systems, namely, numerical linked cluster
expansions (NLCEs) \cite{rigol_bryant_06,rigol_bryant_07,rigol_bryant_07b}.
We show that NLCEs not only exhibit exponential convergence with
increasing system size but generally outperform grand canonical
ensemble calculations in systems with periodic boundary conditions.

The paper is organized as follows. In Sec.~\ref{sec:gener-cons}
we argue, based on a high temperature expansion of the partition
function, that grand canonical ensemble calculations in
translationally invariant systems have exponentially small finite-size
errors. In Sec.~\ref{sec:verification}, we discuss analytically
solvable examples, the 1D and 2D Ising models, that substantiate the
arguments in Sec.~\ref{sec:gener-cons}. In Sec.~\ref{sec:pert-th}, we present
a proof that finite-size errors are indeed exponentially small in the grand
canonical ensemble for translationally invariant noninteracting systems
and that, within perturbation theory, the same scaling applies to interacting
systems. We then study numerically, in Sec.~\ref{sec:numerical-tests},
three examples where we systematically compare results from canonical and
grand canonical ensemble calculations, each for open boundary conditions (OBC)
and periodic boundary conditions (PBC), and NLCEs. We
summarize our results and conclude in Sec.~\ref{sec:conclusions}.

\section{General considerations}
\label{sec:gener-cons}

In this section, we argue that the GE for a translationally invariant
system (we abbreviate the CE [GE] with open and
  periodic boundary conditions as CE-O [GE-O] and CE-P [GE-P]
  respectively) has exponentially small finite size corrections.  For
that, we make use of a $\beta$ expansion of the free energy, where
$\beta=(k_{B}T)^{-1}$ is the inverse temperature and $k_{B}$ is the
Boltzmann constant. This kind of expansion has been used extensively
in the literature to compute partition functions for various
models \cite{domb60a,domb60b,oitmaa_hamer_book_06}.

Consider the Taylor expansion of the grand partition function $Z\equiv
\Tr e^{-\beta \hat{H}}$ (we set $\mu=0$ for brevity; all the arguments
below are valid for nonzero $\mu$, which will be required for bosons
to prevent Bose-Einstein condensation):
\begin{equation}
  Z(\beta) = \Tr(1) - \beta \Tr(\hat{H}) + \frac{\beta^{2}}{2!}\Tr(\hat{H}^{2}) + \ldots
\end{equation}
We are interested in $\ln Z$, from which thermodynamic quantities can
be obtained by taking suitable $\beta$ or $\mu$ derivatives,
\begin{multline}
  \label{eq:9}
  \ln Z(\beta) = \ln\Tr(1) - \beta \frac{\Tr(\hat{H})}{\Tr(1)} + \\
  + \frac{\beta^{2}}{2}
  \left[\frac{\Tr(\hat{H}^{2})}{\Tr(1)}-\frac{\Tr(\hat{H})^{2}}{\Tr(1)^{2}}\right]
  + \cdots
\end{multline}
We note that
\begin{equation}
  \label{eq:8}
  \frac{\Tr(\hat{H}^{n})}{\Tr(1)} = \frac{\Tr(\hat{H}^{n}e^{-0\cdot \hat{H}})}
  {\Tr(e^{-0\cdot \hat{H}})},
\end{equation}
is an infinite temperature expectation value. At infinite temperature
all unconnected parts of the system, however close to each other, are
uncorrelated. Therefore, the expansion in Eq.~\eqref{eq:9} reduces to
a sum over only the connected graphs that can be embedded in the
finite system \cite{domb60a,sykes66,huang87}. In the CE, the particle
number constraint, i.e., that the total particle number is fixed,
implicitly correlates unconnected pieces of a graph, and therefore
this simplification does not occur.

For a system that has no ordered phase or, equivalently, where order
appears only at $T=0$, the above expansion must converge beyond a
certain order (because the correlation length is finite).  Since this
must occur for any temperature $\beta < \infty$, the convergence must
come from the coefficients of the $\beta$-expansion, i.e., from the
traces in Eq.~\eqref{eq:9}.  In other words, for the expansion $(\ln
Z)/N = a_{0}+a_{1}\beta + a_{2}\beta^{2}+\cdots$, the coefficient
$a_{n}$ must fall faster than $e^{-n}$ for the series to converge for
any $\beta$. The convergence cannot come from cancellation of terms of
opposite sign, since any such cancellation can work only at some
fine-tuned value of $\beta$.

For a system that has a phase transition between an unordered high
temperature phase and an ordered low temperature phase at a finite
critical temperature $\beta_{c}$, the coefficients do not exhibit this
behavior -- in the critical phase the correlation length is infinite
and all orders of this expansion are relevant. The convergence of the
series for $\beta<\beta_{c}$ instead comes from the fact that
$\beta/\beta_{c} < 1$.  We verify these arguments in the 1D and 2D
Ising models.

From here on, we assume that we are in a phase where the
$\beta$-expansion converges. We will now show that with PBC (when the
system is translationally invariant), all orders of the expansion
Eq.~\eqref{eq:9} up to the system size (to be properly defined below)
are identical to those in the thermodynamic limit. We will further
show that this is not the case with open boundary conditions.

\subsection{Periodic boundary conditions}

Consider a system with $N$ sites and periodic boundary conditions (a
system that is translationally invariant). The $\beta$ expansion of
$\ln Z$ is shown in Eq.~\eqref{eq:9}, in which each term can be
represented by a graph embedded on the finite system. We will call
these graphs clusters \cite{huang87}.  First note that each cluster
has $N$ equivalent positions on the lattice since the system is
translationally invariant.  That gives a factor of $N$ that we move to
the left hand size in Eq.~\eqref{eq:9} to get $(\ln Z)/N$, i.e., an
intensive quantity. Let us now consider a cluster with $c$
sites. First, since we have a cumulant expansion, as discussed above,
only connected clusters enter (see, e.g.,
Refs.~\cite{huang87,domb60a,oitmaa_hamer_book_06} for details).  If
the extent of our cluster in each direction is less than the system
size in that directions (say $L$), then the cluster has open boundary
conditions. Furthermore, even if the system size is increased in any
direction, this cluster is present. Hence, this cluster is present in
the thermodynamic limit. In general, \emph{every cluster with $c$
  sites that, in the finite lattice with $N$ sites, does not wrap
  around any boundary appears in the infinite system, and vice versa}.
Therefore the contribution of this cluster in a finite system is
exactly the same as its contribution in the thermodynamic limit. On
the other hand, a cluster with $L$ sites in any given direction wraps
around a boundary, i.e., it does not appear in the thermodynamic
limit. As a result, clusters that wrap around boundaries give
contributions that are not present in the thermodynamic limit
\cite{lebowitz61}. Hence, the difference between results in the
thermodynamic limit and in finite-size periodic systems is
$O(\beta^{L-p})$, $p$ being determined by the Hamiltonian.  We note
that $p$ is $O(1)$ for local Hamiltonians, which are the ones of
interest here. Further, based on the earlier argument, we must have
that the coefficient at this order falls faster than $e^{-(L-p)}$ or
that the expansion parameter ($\beta/\beta_{c}$) is smaller than one.
Therefore, finite-size errors in a GE calculation of a translationally
invariant system at any temperature in the unordered phase are
smaller than $O(e^{-L})$, for systems with linear dimension $L$.

\subsection{Open boundary conditions}

For a system with OBC, one immediately realizes that clusters do not
have $N$ equivalent positions on the lattice. As a result, even if a
given cluster in the finite system appears in the thermodynamic limit,
its contribution in the finite system will differ from the
thermodynamic limit. For example, for a lattice model in which the
Hamiltonian is a sum of terms involving only nearest neighbor sites,
the term linear in $\beta$ in Eq.~\eqref{eq:9} for a system with OBC
has a correction $O(A/2N)$ relative to the result for PBC, where $A$
is the number of sites in the boundary.  This correction vanishes as
$1/L$ with increasing system size.  Complicated geometric and
combinatorial factors appear at higher orders, all of which approach
the thermodynamic limit result with increasing system size. Hence,
none of the coefficients of a $\beta$ expansion for a finite system
with OBC match the result in the thermodynamic limit, and
finite size errors in $(\ln Z)/N$ are $O(1/L)$.

\subsection{Numerical linked cluster expansions}

Rather than making calculations of finite systems with periodic or
open boundary conditions, and then extrapolating the results to the
thermodynamic limit, another way to calculate finite-temperature
properties of lattice systems in the thermodynamic limit is to use
NLCEs
\cite{rigol_bryant_06,rigol_bryant_07,rigol_bryant_07b,tang_khatami_13}.
The idea in this case is to directly use the linked cluster expansion
of the infinite (translationally invariant) system, for which any
extensive quantity ${\cal O}$ per site can be computed as the sum
\begin{equation}
  \label{eq:nlce1}
  \frac{\cal O}{N} = \sum_{c}M(c)\times W_{\cal O}(c),
\end{equation}
over all connected clusters $c$ that can be embedded in the infinite
lattice. In Eq.~\eqref{eq:nlce1}, $M(c)$ is the \emph{multiplicity} of
cluster $c$, namely, the number of ways per site in which cluster $c$
can be embedded on the lattice, and $W_{\cal O}(c)$ is the
\emph{weight} of the cluster $c$ for observable ${\cal O}$.
$W_{\cal O}(c)$ is calculated by an inclusion-exclusion principle, one
systematically subtracts contributions from the connected subclusters
of $c$ \cite{sykes66}
\begin{equation}\label{eq:nlce2}
  W_{\cal O}(c) = \mathcal{O}(c) - \sum_{s\subset c}W_{\cal O}(s).
\end{equation}
${\cal O}(c)$ is the value of the observable evaluated on the cluster
$c$. In NLCEs, ${\cal O}(c)$ is obtained using a full exact
diagonalization of the Hamiltonian for cluster $c$.

Due to computational limitations, only a finite number of clusters can
ultimately be calculated in Eq.~\eqref{eq:nlce1}. Nevertheless, as
shown in
Refs.~\cite{rigol_bryant_06,rigol_bryant_07,rigol_bryant_07b}, NLCEs
can converge at lower temperatures than high temperature expansions,
and sometimes all the way to the ground state for systems with
unordered ground states. Also, NLCEs can provide very accurate
results for temperatures at which exact diagonalization results for
systems with periodic boundary conditions suffer from very large
finite-size effects. A pedagogical introduction to implementing NLCEs
can be found in Ref.~\cite{tang_khatami_13}.

In what follows, we compare NLCE results with those obtained in
calculations in finite systems with different boundary conditions.
Our goal is to find how each of them converges to the thermodynamic
limit result and which converges the fastest. For NLCEs, the accuracy
of the results is determined by the size of the largest clusters
considered in the sum in Eq.~\eqref{eq:nlce1} and the model under
consideration.

\section{Verification in Ising models}
\label{sec:verification}

In this section, we verify the arguments given in
Sec.~\ref{sec:gener-cons} in the 1D and 2D Ising models, both of which
can be solved analytically.

\subsection{1D Ising model}
\label{sec:1d-ising}

In the thermodynamic limit, the log of the partition function per
site, $\Omega$, can be obtained using the transfer matrix method
\cite{huang87}, and is given by (we set $J=1$)
\begin{equation}
  \label{eq:1}
  \Omega(\beta) = \ln(e^{\beta} + e^{-\beta}).
\end{equation}
The expansion in powers of $\beta$ of this result is
\begin{equation}
  \label{eq:4}
  \Omega(\beta) = \ln2 + \frac{\beta^{2}}{2}-\frac{\beta^{4}}{12}
  + \frac{\beta^{6}}{45}-\frac{17\beta^{8}}{2520}
  + \frac{31\beta^{10}}{14175}+ \cdots
\end{equation}

The result for finite systems with periodic boundary conditions is
given by
\begin{equation}
  \label{eq:2}
  \begin{split}
    \Omega_L(\beta) &= \frac{1}{L}\ln \left[(e^{\beta} +
      e^{-\beta})^{L}
      + (e^{\beta} - e^{-\beta})^{L}\right] \\
    &= \Omega(\beta) + \frac{\tanh^{L}\beta}{L} + \cdots
  \end{split}
\end{equation}
Since $0\leq \tanh\beta<1$ for $0\leq \beta<\infty$, the finite size
error is indeed faster than exponential. Quantities like the energy,
which are derivatives of the free energy, converge exponentially fast
with $L$.

Expanding Eq.~\eqref{eq:2} in powers of $\beta$ for different values
of $L$, we get
\begin{equation}
  \label{eq:3}
  \begin{split}
    \Omega_2(\beta) &= \ln2 + \beta^{2} + \ldots \\
    \Omega_3(\beta) &= \ln2 + \frac{\beta^{2}}{2} + \frac{\beta^{3}}{3} + \ldots \\
    \Omega_4(\beta) &= \ln2 + \frac{\beta^{2}}{2} + \frac{\beta^{4}}{6} + \ldots \\
    \Omega_5(\beta) &= \ln2 + \frac{\beta^{2}}{2} -
    \frac{\beta^{4}}{12} +
    \frac{\beta^{5}}{5} + \ldots \\
    \Omega_6(\beta) &= \ln2 + \frac{\beta^{2}}{2} -
    \frac{\beta^{4}}{12} + \frac{17\beta^{6}}{90} + \ldots
  \end{split}
\end{equation}
As one can see, the results are exact to $O(\beta^{L-1})$. Naturally,
as one goes to lower temperatures, the correlation length increases
and larger systems are required to capture the relevant powers of
$\beta$.  It is easy to verify that, with open boundary conditions,
the corrections are always $O(1/L)$.

As discussed in Sec.~\ref{sec:gener-cons}, although the coefficients
in Eq.~\eqref{eq:3} are exact up to $O(\beta^{L-1})$, the convergence
for \emph{all} temperatures comes from the fact that they fall off
rapidly (faster than exponential) with increasing expansion
order. Figure \ref{fig:ising-fit} shows a plot of the coefficients of
expansion in Eq.~\eqref{eq:4} along with a fit to $a b^{-L}/L$ that
demonstrates the faster-than-exponential behavior.

To conclude our discussion of the 1D Ising model, we evaluate the
first few orders of the NLCE for this model (instead of numerical
exact diagonalization of the clusters, we obtain these results
analytically).  First, we evaluate the partition function [$\ln
Z_L(\beta)$] on finite clusters with OBC
\begin{equation}
  \label{eq:ising-nlce-prop}
  \begin{split}
    \ln Z_1(\beta) &= \ln2, \\
    \ln Z_2(\beta) &= \ln2 + \ln\left(e^{\beta} + e^{-\beta}\right),\\
    \ln Z_3(\beta) &= \ln2 + \ln\left(e^{2\beta}+e^{-2\beta} +
      2\right),
  \end{split}
\end{equation}
and then carry out the subtractions. The weights are given by [see
Eq.~\eqref{eq:nlce2}]
\begin{equation}
  \label{eq:ising-nlce-weights}
  \begin{split}
    W_{1} &= \ln Z_{1}(\beta)  = \ln2, \\
    W_{2} &= \ln Z_{2}(\beta)-2W_{1} = \ln\left(e^{\beta} + e^{-\beta}\right)- \ln2,\\
    W_{3} &= \ln Z_{3}(\beta)-2W_{2}-3W_{1}= 0.
  \end{split}
\end{equation}
Hence, the result for $\Omega$ obtained in calculations including up
to $n$ sites, $(\Omega)_n$, is given by [see Eq.~\eqref{eq:nlce1}]
\begin{equation}
  \label{eq:ising-nlce-sums}
  \begin{split}
    (\Omega)_1 &= W_{1}=\ln2, \\
    (\Omega)_2 &= W_{1}+W_{2}= \ln\left(e^{\beta} + e^{-\beta}\right), \\
    (\Omega)_3 &= W_{1}+W_{2}+W_{3} = \ln\left(e^{\beta} +
      e^{-\beta}\right).
  \end{split}
\end{equation}
It can be verified that the last result is valid at all higher orders
in the ``NLCE''.  The thermodynamic limit result is therefore obtained
by just considering clusters with one and two sites. This is an
infinite improvement over the use of the grand canonical ensemble with
periodic boundary conditions.  Whereas this infinite gain is specific
to the 1D Ising model --- the model can, after all, be solved using a
two dimensional transfer matrix, we show in what follows that the fact
that NLCEs outperform exact calculations in finite systems appears to
be generic.

\begin{figure}[!tb]
  \centering
  \includegraphics[width=0.95\columnwidth]{ising-new.eps}
  \caption{(Color online) Coefficients (absolute value) of the
    $\beta$-expansion of the free energy vs. the expansion order
    $l$. We show the coefficients in Eq.~\eqref{eq:4} for the 1D Ising
    model, and a fit to the function $a b^{-l}/l$ with $a=2.0000$ and
    $b=1.5708$ (the difference between the exact result and
    $2(\pi/2)^{-l}/l$ vanishes exponentially fast with $l$), and from
    Eq.~\eqref{eq:7} for the 2D Ising model. The coefficients in the
    latter case do not fall off exponentially fast. Because of this,
    the $\beta$-expansion only converges for $\beta < \beta_{c}$. It
    is only in this regime that finite-size errors in grand-canonical
    ensemble calculations of translationally invariant systems are
    exponentially small in system size. \textbf{Inset:}
      Shows the rational part of the coefficients in
      Eq.~\eqref{eq:7}. They decrease at first but increase after
      $O(\eta^{24})$. See text for further discussion.}
  \label{fig:ising-fit}
\end{figure}

\subsection{2D Ising model}
\label{sec:2d-ising}

For the 2D Ising model, $\Omega$ is given by \cite{huang87,Onsager_44}
\begin{multline}
  \label{eq:5}
  \Omega(\beta) = \ln[2\cosh(2\beta)] + \\ +
  \int_{0}^{\pi}\frac{d\phi}{2\pi}
  \ln\left[\frac{1+\sqrt{1-\frac{4\sin^{2}\phi}{\cosh^{2}(2\beta)\coth^{2}(2\beta)}}}{2}
  \right].
\end{multline}
The $\beta$-expansion of this result is given by
\begin{equation}
  \label{eq:6}
  \Omega(\beta) = \ln2 + \beta^{2} + \frac{5\beta^{4}}{6} + \frac{32\beta^{6}}{45}
  + \frac{425\beta^{8}}{252} + \ldots
\end{equation}
One can see that the coefficients become larger than 1 for higher
orders. The correct expansion parameter for models with a
finite-temperature transition is $\beta/\beta_{c}$.  For the classical
2D Ising model on a square lattice, $\beta_{c} = \ln(1+\sqrt{2})/2$.
This gives (with $\eta \equiv \beta/\beta_{c}$),
\begin{equation}
  \label{eq:7}
  \Omega(\beta) = \ln2 + \frac{a^{2}\eta^{2}}{4} + \frac{5a^{4}\eta^{4}}{96} +
  \frac{a^{6}\eta^{6}}{90} + \frac{425a^{8}\eta^{8}}{64512} + \ldots
\end{equation}
where $a\equiv\ln(1+\sqrt{2})<1$. The coefficients of the $\eta$-expansion
are plotted versus the order of the expansion in Fig.~\ref{fig:ising-fit}.
They do not fall off faster than exponential, or exponentially.
We note that the rational part of the coefficients of the
first few orders of the $\eta$-expansion reported in Eq.~\eqref{eq:7} is
deceiving. They decrease with increasing order of
expansion. This, together with the fact that $a<1$, suggests that the
coefficients of the $\eta$-expansion should fall off
faster than exponentially. However, as shown in the inset in
Fig.~\ref{fig:ising-fit}, the aforementioned rational part \emph{increases}
with increasing order of expansion after $O(\eta^{24})$. Because of this,
convergence in the $\eta$ expansion is only expected for $\eta<1$ and
does not come from the coefficients.

We also calculate $\ln Z_{L\times L}$ for small systems with $N=L\times
L$ sites in the GE for both periodic and open boundary
conditions. With PBCs
\begin{eqnarray}
  \label{eq:2d-ising-g-p}
  \ln Z^{\rm 2D-P }_{1\times1} &=& \ln2, \nonumber\\
  \ln Z^{\rm 2D-P}_{2\times 2} &=& \ln2 + \ln\left(e^{4\beta} + e^{-4\beta}+2\right), \\
  \ln Z^{\rm 2D-P}_{3\times 3} &=& \ln2 + \ln\bigl(e^{-18 \beta
  }+9 e^{-10 \beta }+24 e^{-6 \beta }+ \nonumber\\
  && +99 e^{-2 \beta }+72 e^{2
    \beta }+51 e^{6 \beta }\bigr),\nonumber
\end{eqnarray}
whereas, with OBCs,
\begin{eqnarray}
  \label{eq:2d-ising-g-o}
  \ln Z^{\rm 2D-O}_{1\times1} &=& \ln2 \nonumber\\
  \ln Z^{\rm 2D-O}_{2\times 2} &=& \ln2 + \ln\left(e^{4\beta} + e^{-4\beta}+2\right) \\
  \ln Z^{\rm 2D-O}_{3\times 3} &=& \ln2 + \ln\bigl(e^{-12 \beta
  }+4 e^{-8 \beta }+16 e^{-6 \beta }+\nonumber\\
  &&+23 e^{-4 \beta }+48 e^{-2
    \beta }+
  48 e^{2 \beta }+23 e^{4 \beta }+\nonumber\\
  &&+ 16 e^{6 \beta }+4 e^{8
    \beta }+ e^{12 \beta }+72 \bigr)\nonumber.
\end{eqnarray}

For the $\beta$-expansion of the $3\times 3$ systems, up to the first
term that differs from Eq.~\eqref{eq:6}, we obtain
\begin{equation}
  \label{eq:2d-ising-finite-expansion}
  \begin{split}
    \Omega^{\rm 2D-O}_{3\times 3} &= \ln2 + \frac{2\beta^{2}}{3} +
    \ldots \\
    \Omega^{\rm 2D-P}_{3\times 3} &= \ln2 + \beta^{2} +
    \frac{2\beta^{3}}{3} + \ldots
  \end{split}
\end{equation}
We see that whereas for OBC the coefficient of the second order term is
incorrect, for PBC it is correct, i.e., once again GE-P gives results
that are correct to $O(\beta^{L-1})$, with $L=3$ in this case.

For the NLCE calculation with clusters with up to four sites, we
obtain
\begin{eqnarray}
  \label{eq:2d-ising-nlce}
  &&(\Omega)_4 = 20 \ln \left(e^{-\beta }+e^{\beta }\right)+  54
  \ln \left(e^{- \beta } e^{2 \beta }+1\right)\nonumber\\
  && \ -38
  \ln \left(e^{-2 \beta }+e^{2 \beta }+2\right)+\ln \left(e^{-4
      \beta }+e^{4 \beta
    }+6\right).\
\end{eqnarray}
Expanding in powers of $\beta$, and reporting terms up to the first
one that differs from Eq.~\eqref{eq:6}, we get
\begin{equation}\label{eq:2dinlcee}
  (\Omega)_4 = \ln2 + \beta^{2} + \frac{5\beta^{4}}{6} -
  \frac{58\beta^{6}}{45} + \ldots
\end{equation}
The above result is correct to $O(\beta^{5})$. We must stress that
Eq.~\eqref{eq:2dinlcee} was obtained in an expansion in which the
largest cluster has $N=4$, while
Eqs.~\eqref{eq:2d-ising-finite-expansion} are for systems with
$N=9$. The gain is evident.

\section{Finite temperature perturbation theory to all orders}
\label{sec:pert-th}

In the case of bosons or fermions on a lattice that can be treated by
finite-temperature perturbation theory (with a noninteracting theory
as the unperturbed starting point), a proof that finite-size errors are
exponentially small to all orders in perturbation theory can be made based
on the momentum-space representation of the Hamiltonian. The proof is
essentially identical for bosons and fermions, so we focus on the former.

We consider a generic massive scalar field theory in a 1D lattice,
with unit lattice spacing, $L$ sites, and PBC:
\begin{equation}
  \label{H}
  \hat{H}_{s}=\frac12 \sum_{j=1}^L\left[\hat{\pi}_j^2
    +(\hat{\varphi}_{j+1}-\hat{\varphi}_j)^2 + m^2\hat{\varphi}_j^2\right],
\end{equation}
where $[\hat{\varphi}_j,\hat{\pi}_{j'}]=i\delta_{jj'}$ and
$\hat{\varphi}_{j+L}\equiv\hat{\varphi}_j$,
$\hat{\pi}_{j+L}\equiv\hat{\pi}_j$.
This Hamiltonian is diagonalized via
\begin{equation}
  \begin{split}
    \hat{\varphi}_j &= \frac{1}{\sqrt{2L}}\sum_{n=0}^{L-1}\omega_n^{1/2}
    \left[e^{\frac{2\pi i n j}{L}}\hat{a}_n + e^{-\frac{2\pi i n
          j}{L}}\hat{a}^\dagger_n\right],\\
    \hat{\pi}_j &= -\frac{i}{\sqrt{2L}}\sum_{n=0}^{L-1}\omega_n^{-1/2}
    \left[e^{\frac{2\pi i n j}{L}}\hat{a}_n - e^{-\frac{2\pi i n j}{L}}\hat{a}^\dagger_n\right],
  \end{split}
\end{equation}
so that
\begin{equation} \label{Hk} \hat{H}_{s} = \sum_{n=0}^{L-1}\omega_n\hat{a}^\dagger_n
  \hat{a}_n + \text{constant}\hfill,
\end{equation}
where $[\hat{a}^{}_n,\hat{a}^\dagger_{n'}] = \delta_{nn'}$, $\omega_n =\omega(k_n)$, $k_n = 2\pi n/L$,
$n\in[0,L-1]$, and
\begin{equation}
  \label{wk}
  \omega(k) = \sqrt{2(1-\cos k)+m^2} =
  \sqrt{4\sin^2(k/2)+m^2}.
\end{equation}

We now want to compute the grand canonical partition function
\begin{equation}
  Z(\beta,\mu)\equiv \mathop{\rm Tr}e^{-\beta (\hat{H} -\mu \hat{N})}\,,
  \label{Z}
\end{equation}
where $\hat{N}=\sum_n \hat{a}^\dagger_n \hat{a}_n$ is the total number operator.
We take the trace in the Fock basis of eigenstates of each $\hat{a}^\dagger_n a_n$,
\begin{equation}
  Z_L(\beta,\mu) = \prod_{n=0}^{L-1}\sum_{N_n=0}^\infty
  e^{-\beta(\omega_n-\mu)N_n} = \prod_{n=0}^{L-1} \frac{1}{1-e^{-\beta(\omega_n-\mu)}}.
  \label{Z1}
\end{equation}
Equivalently,
\begin{equation}
  \Omega_L(\beta,\mu) = \frac{1}{L}\ln Z_L(\beta,\mu) =
  \frac{1}{L}\sum_{n=0}^{L-1} F(k_n),
  \label{OL}
\end{equation}
where we have defined
\begin{equation}
  \label{F}
  F(k) \equiv -\ln[1-e^{-\beta(\omega(k)-\mu)}].
\end{equation}
In the thermodynamic limit, the sum over $n$ becomes an integral,
\begin{equation}
  \Omega(\beta,\mu)\equiv\lim_{L\to\infty}\Omega_L(\beta,\mu)=
  \frac{1}{2\pi}\int_{0}^{2\pi}dk\,F(k)\,.
  \label{Oi}
\end{equation}

We now wish to show that $|\Omega_L(\beta,\mu)-\Omega(\beta,\mu)|$ is
exponentially small in $L$.

We first note that $F(k)$ is periodic in $k$ with period $2\pi$.
Therefore its Fourier expansion takes the form
\begin{equation}
  F(k)=\sum_{j=-\infty}^{+\infty}\tilde F_j\,e^{i jk}\,,
  \label{fsum}
\end{equation}
where the Fourier coefficients are given by
\begin{equation}
  \tilde F_j = \frac{1}{2\pi} \int_{0}^{2\pi}dk\,e^{-i jk}F(k)\,.
  \label{tildef}
\end{equation}
Using \eqs{Oi} and (\ref{tildef}), we get
\begin{equation}
  \Omega(\beta,\mu) = \tilde F_0\,.
  \label{Oi1}
\end{equation}
Using \eqs{OL} and (\ref{fsum}), we get
\begin{eqnarray}
  \Omega_L(\beta,\mu) &=& \frac{1}{L}\sum_{n=0}^{L-1} F(2\pi n/L)
  = \frac{1}{L}\sum_{n=0}^{L-1}\sum_{j=-\infty}^{+\infty}\tilde
  F_j\,e^{2\pi i jn/L} \nonumber\\
  &=& \sum_{j=-\infty}^{+\infty}\tilde F_j\biggl[\frac{1}{L}
  \sum_{n=0}^{L-1}e^{2\pi i jn/L}\biggr]  \label{OL5}
  \\& =& \sum_{j=-\infty}^{+\infty}\tilde
  F_j\delta_{j\,\rm{mod}\,L,\,0}
  = \sum_{j'=-\infty}^{+\infty}\tilde F_{j'\!L}\,.\nonumber
\end{eqnarray}
Subtracting \eq{Oi1} from \eq{OL5}, we get \begin{equation}
  O_L(\beta,\mu)-\Omega(\beta,\mu) = \sum_{j\ne 0}\tilde F_{jL}\,.
  \label{Odiff}
\end{equation}
Examining \eqs{wk} and (\ref{F}), we see that if $m>0$ and $\mu<m$
(necessary to avoid Bose condensation), then $F(k)$ is continuous and
infinitely differentiable for all real $k$.  It then follows from a
general theorem of Fourier series \cite{tolstov76} that $\tilde F_{jL}$ goes
to zero faster than any power of $|j|L$ as $L\to\infty$.  The sum over
$j$ in \eq{Odiff} will then be dominated by the $j=\pm1$ terms.  We
conclude that $|O_L(\beta,\mu)-\Omega(\beta,\mu)|$ is exponentially
small in $L$ if $m>0$ and $\mu<m$.

We can verify this explicitly. Again examining \eqs{wk} and (\ref{F}),
we see that $F(k)$ is changing most rapidly near $k=0$. For $m>0$ and
small enough $k$, $F(k)$ can be approximated via
\begin{equation}
  F(k)\simeq -\ln[1-e^{-\beta(m-\mu)}e^{-\beta k^2/2m}].
  \label{F0}
\end{equation}
Assuming $\mu<m$, we have
\begin{equation}
  F(k)\simeq \sum_{n=1}^\infty\frac{1}{n}e^{-n\beta(m-\mu)}e^{-n \beta k^2/2m}\,.
  \label{F2}
\end{equation}
In this approximation, we get
\begin{equation}
  \tilde F_{jL} \simeq
  \sqrt{\frac{m}{2\pi\beta}} \sum_{n=1}^\infty \frac{1}{n^{3/2}}
  e^{-n\beta(m-\mu)}e^{-j^2 L^2 m/2 \beta n}\,.
  \label{tF2}
\end{equation}
The sum over $n$ can be approximated by steepest descent.  For
$j^2L^2\gg1/m(m-\mu)$, we find
\begin{equation}
  \tilde F_{jL} \simeq  \frac{1}{|j|L} \exp\Bigl[-[2m(m-\mu)]^{1/2}|j|L\Bigr]\,,
  \label{tF3}
\end{equation}
This is exponentially small in $L$, as expected from the general
theorem.  The sum over $j$ in \eq{Odiff} is then dominated by
$j=\pm1$, and we have
\begin{equation}
  |\Omega_L(\beta,\mu)-\Omega(\beta,\mu)|
  \simeq \frac{2}{L}\exp\Bigl[-[2m(m-\mu)]^{1/2}L\Bigr]\,,
  \label{Odiff2}
\end{equation}
which is the behavior observed for the 1D Ising model in
Fig.~\ref{fig:ising-fit}.

This proof extends straightforwardly to higher dimensions, assuming
periodic boundary conditions in each dimension.

Now consider adding an interaction term to $\hat{H_s}$, such as
$g\sum_j\varphi_j^4$.  We can compute $\Omega_L(\beta,\mu)$ order by order
in finite temperature perturbation theory.  Each term is represented by
a connected Feynman diagram \cite{fetter_bk_03}. The expression for a diagram with $p$
propagators and $v$ vertices takes the form
\begin{equation} \frac{1}{L}\sum_{n_1=0}^{L-1}
  \ldots \frac{1}{L}\sum_{n_{p-v+1}=0}^{L-1} g^v
  F(k_{n_1},\ldots,k_{n_{p-v+1}})\,.
  \label{feyn}
\end{equation}
For $m>0$ and $\mu<m$, $F$ is infinitely differentiable in each
$k_{n_i}$. Hence, one can once again apply the general theorem that says
that the difference between \eq{feyn} and the $L\to\infty$ limit must
be exponentially small in $L$.

Other models of bosons and/or fermions can be analyzed in exactly the same
way, with only the form of $F(k)$ changing. As long as $F(k)$ is infinitely
differentiable, the argument holds.

\section{Numerical tests}
\label{sec:numerical-tests}

In this section, we discuss numerical results that support the
generality of the conclusions reached so far. The results are obtained
using full exact diagonalization calculations for fermionic systems in
one dimension. We study spinless fermions with nearest-neighbor
hoppings that are either noninteracting or that have nearest and
next-nearest-neighbor interactions. They are described by the
following generic Hamiltonian
\begin{multline}
  \label{eq:gen_f_ham}
  \hat{H}_{\rm f} =
  \sum_{i}\Biggl[-t(\hat{c}^{\dag}_{i}\hat{c}^{}_{i+1}+\text{H.c.})
  + V \left(\hat{n}_{i}-\frac12\right)\left(\hat{n}_{i+1}-\frac12\right) \\
  + V'\left(\hat{n}_{i}-\frac12\right)\left(\hat{n}_{i+2}-\frac12\right)\Biggl],
  \qquad\quad
\end{multline}
where $\hat{n}^{}_{i}=\hat{c}^{\dag}_{i}\hat{c}^{}_{i}$ is the number
operator.  We present results for different values of the parameters
$V$ and $V'$, and using both the canonical and grand canonical
ensembles with open and periodic boundary conditions. We also report
NLCE results.

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.98\columnwidth]{FF_v02.eps}
  \caption{(Color online) Comparison of results of canonical and grand canonical
    ensembles with open and periodic boundary conditions, and NLCEs,
    for noninteracting fermions hopping on a one dimensional lattice
    (corresponds to Eq.~\eqref{eq:gen_f_ham} with $t=1$,
    $V=V'=0$). The energy difference $\delta E_l$ (see text) is
    plotted for two different temperatures in (a) and (b). Panel (c)
    shows the occupation of the $k=0$ mode, or, equivalently, the
    correlation $L^{-1}\sum_{ij}\langle c^{\dag}_{i}c_{j}\rangle$.
    $\delta n_{k=0}$ is zero for the grand-canonical ensemble with
    PBCs, see text.}
  \label{fig:free-ferm}
\end{figure}

In all cases, we calculate the energy per site $E_{l}$. For
noninteracting fermions, we also compute the occupation of the
momentum $k=0$ mode, $n^l_{k=0}\equiv \sum_{i,j}\langle
\hat{c}^{\dag}_{i}\hat{c}^{}_{j}\rangle$.  For the interacting models,
we compute the the nearest-neighbor single-particle correlation
function $K_l = \sum_{i}\langle
\hat{c}^{\dag}_{i}\hat{c}^{}_{i+1}+{\rm H.c.}\rangle$. For the results
reported, by $l$ we mean the number of sites of the finite system or
the number of sites of the largest cluster in the NLCE. We plot the
finite size errors $\delta E_{l}$, $\delta n^l_{k=0}$, and $\delta
K_{l}$, with $\delta O_{l}\equiv (O_{l}-O)/O$ where $O$
($E,\,n_{k=0},\,K$) is either the exact analytic result in the
thermodynamic limit, when known, or the highest order result from a
numerical linked cluster expansion.

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.98\columnwidth]{t1V1_v01.eps}
  \caption{(Color online) (a) $\delta E_{l}$ and (b) $\delta K_{l}$ for interacting
    fermions [Eq.~\eqref{eq:gen_f_ham}] with $t=1$ (unit of energy),
    $V=1.0,\,V'=0$, and $T=1.0$.}
  \label{fig:t1v1}
\end{figure}
\begin{figure}[!b]
  \centering
  \includegraphics[width=0.98\columnwidth]{t1V1Vp1_v01.eps}
  \caption{(Color online) (a) $\delta E_{l}$ and (b) $\delta K_{l}$ for interacting
    fermions [Eq.~\eqref{eq:gen_f_ham}] with $t=1$ (unit of energy),
    $V=1.0,\,V'=1.0$, and $T=1.0$.}
  \label{fig:t1v1vp1}
\end{figure}

Figure \ref{fig:free-ferm} reports results for $\delta E_{l}$ and
$\delta n^l_{k=0}$ for noninteracting fermions
[Eq.~\eqref{eq:gen_f_ham} with $t=1$ and $V=V'=0$].  Panels (a) and
(b) show the finite-size errors of the energy at two temperatures, and
(c) shows the finite-size errors of the zero-mode occupation (the sum
of all one-particle correlations). The results for the grand-canonical
energy in finite systems can be obtained analytically
\begin{equation}
  \label{eq:10}
  E_{l} = \sum_{n=0}^{l-1} \frac{\epsilon_{n} e^{-\beta\epsilon_{n}}}
  {1+e^{-\beta\epsilon_{n}}},
\end{equation}
where $\epsilon_{n}$ are the single-particle eigenenergies.  The
momentum distribution function with periodic boundary conditions is
given by the Fermi-Dirac distribution
\begin{equation}
  \label{eq:11}
  n_{k} = \frac{1}{1+e^{-\beta\epsilon_{k}}}.
\end{equation}
We note that $n_{k}$ for finite systems (for the values of $k$ allowed)
is exactly the same as in the thermodynamic limit. Hence, $\delta
n^l_{k=0}=0$ for the grand canonical ensemble with PBC.  Therefore,
no error is reported for the GE-P in Fig.~\ref{fig:free-ferm}(c).  The results in
the canonical ensemble are obtained as described in
Ref.~\cite{rigol_05}.

For the two temperatures and two quantities reported in
Fig.~\ref{fig:free-ferm}, one can see that the errors of the GE-O,
CE-P, and CE-O results decrease as $1/l$. This is expected and is made
apparent in the plots by comparing those results to the $1/l$ plots
(dashed lines) depicted for reference. On the other hand, the GE-P
(for the energy) and NLCE errors can be seen to decrease exponentially
with increasing $l$. For the energy, the NLCE errors are much smaller
than the GE-P ones, showing once again that NLCEs generally ($n_{k}$
for noninteracting fermions being a counterexample) outperform
grand canonical calculations in finite translationally invariant
systems.

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.96\columnwidth]{2dmodels_v01.eps}
  \caption{(Color online) $\delta E_{l}$ in NLCE calculations of the Heisenberg model
    in four different two-dimensional lattices, as indicated. The
    approach to the thermodynamic limit result is again exponentially
    fast as one increases the order of the expansion. The results presented
    here were taken from
    Refs.~\cite{rigol_bryant_07,tang_khatami_13,tang_paiva_12,tang_paiva_13}.
    $J=1$ is taken to be the unit of energy.}
  \label{fig:2dmodels}
\end{figure}

Figures \ref{fig:t1v1} and \ref{fig:t1v1vp1} show the results for: (a)
$\delta E_{l}$ and (b) $\delta K_{l}$ for interacting fermions
[Eq.~\eqref{eq:gen_f_ham}] with $(t,V,V')=(1,1,0)$ and $(1,1,1)$,
respectively.  In all cases, it is apparent that only the GE-P and
NLCE errors decrease exponentially fast with $l$, while the GE-O and
CE errors decrease as $1/l$. Once again, these results (now for
interacting systems) show that the NLCE errors are significantly
smaller than those of the GE-P. Here, we have used the highest order
of the NLCE ($l=18$) as the estimate for the
thermodynamic limit \cite{rigol_14}. All results for the interacting
models were obtained using full exact diagonalization of the
Hamiltonian.

Finally, although we cannot make an analysis in 2D equivalent to that
presented for interacting models in 1D (because of the exponential
scaling of the computational cost combined with the fast increase of $N$),
we can still verify that NLCE calculations
have exponentially small errors. In Fig.~\ref{fig:2dmodels}, we show
results for $\delta E_{l}$ for the 2D Heisenberg model in four
different lattice geometries --- square \cite{tang_khatami_13},
honeycomb \cite{tang_paiva_12,tang_paiva_13}, kagome, and triangular
\cite{rigol_bryant_07}. In all cases, the error is once again seen to
decrease exponentially fast with increasing the number of sites in the
clusters considered.

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

We have shown that grand canonical ensemble calculations in
translationally invariant systems that are in unordered phases at
finite temperatures have exponentially small finite size errors, whereas
canonical ensemble calculations have errors that are power law in system
size. Hence, while full exact diagonalization calculations in the
canonical ensemble are computationally less demanding than grand canonical
ones, if one is interested in accurately computing quantities in an
unordered finite-temperature phase, using the grand canonical ensemble is
preferable. The additional computational cost incurred by
diagonalizing all particle sectors (around a factor two in fermionic
systems at half filling) for the GE is far less than the speed up
gained by having to study smaller systems (exponential). Furthermore, we have shown
that numerical linked cluster expansions generally have even smaller
errors than grand canonical ensemble calculations in translationally
invariant systems. The additional cost incurred in the diagonalization
of many clusters with a given size is far less than the speed up
gained by having to study clusters that are much smaller than finite
systems with periodic boundary conditions. The benefit of using NLCEs
is most striking in two-dimensional lattices.

The specific system sizes (cluster sizes)
required to observe exponential convergence in grand canonical calculations
of systems with periodic boundary conditions (in numerical linked cluster expansions)
depend on details such as the model under consideration and the temperatures
of interest, which set the relevant correlation length. Otherwise, our
conclusions are completely general.

\section*{Acknowledgments}
We thank Anders Sandvik and Cristian Batista for stimulating
discussions.  This work was supported by the Office of Naval Research
(D.I. and M.R.), and by NSF Grant PHY13-16748 (M.S.).


\section{Introduction}

Let $M$ be a real orthogonal matrix of size $n\times n$. We are interested in the smallest and largest modulus among the entries of $M$:
\begin{equation}\label{lmum}
l_M:=\min_{1\le i,j\le n} |m_{i,j}|, \ \textrm{and} \ u_M:=\max_{1\le i,j\le n} |m_{i,j}|
\end{equation}
and, more precisely, in estimating the maximal possible value for $l_M$ and the minimal possible value of $u_M$.
In other words, we want to estimate the following quantities:
\begin{equation}\label{lnun}
l_n:=\max_{M\in O(n)}l_M  \quad\textrm{and} \quad
u_n:=\min_{M\in O(n)}u_M,
\end{equation}
where $O(n)$ denotes the group of orthogonal matrices.

\medskip

In order to determine $l_n$ (respectively, $u_n$) one has to control the lowest (resp. uppermost)
absolute value in an orthogonal matrix $M$. It is also natural to ask whether one can control both
quantities simultaneously. For this purpose we introduce the measure of
``flatness'' of an $n\times n$ orthogonal matrix $M$ as
\begin{equation}\label{fm}
f_M=\min \{\e>0 : \frac{1-\e}{\sqrt{n}}\le |m_{i,j}|\le  \frac{1+\e}{\sqrt{n}}, \ \ \textrm{for all} \ 1\le i,  j\le n \},
\end{equation}
and we are interested in how flat an orthogonal matrix can be:
\begin{equation}\label{fn}
f_n:=\min_{M\in O(n)}f_M.
\end{equation}

It is trivial to see that $l_n\le \frac{1}{\sqrt{n}}$, $u_n\ge \frac{1}{\sqrt{n}}$ for all $n$, and the famous Hadamard conjecture states that if $4|n$ then $l_n=u_n=\frac{1}{\sqrt{n}}$. It is also easy to see that
$f_n\geq \frac{c}{n}$, unless there exists an $n\times n$ Hadamard matrix (Proposition \ref{smallerror} below).
This note was motivated by the natural question of L. Baratchart (personal communication), as to whether $l_n\ge \frac{c}{\sqrt{n}}$ for some absolute constant $c$. The answer to this question is positive and provided by Theorem \ref{thm:block} below. We will also prove $u_n\le \frac{c}{\sqrt{n}}$ with another absolute constant $c$ in Theorem \ref{thm:un}. However, we list here some more restrictive questions on flat orthogonal matrices, all of which we have only partial answers to.

\begin{problem}\label{probln}
Is it true that $l_n=(1+o(1))\frac{1}{\sqrt{n}}$?
\end{problem}

\begin{problem}\label{probun}
Is it true that $u_n=(1+o(1))\frac{1}{\sqrt{n}}$?
\end{problem}
\noindent Can we control $l_n$ and $u_n$ simultaneously?
\begin{problem}\label{probfn}
Is it true that $f_n=o(1)$?
\end{problem}

These questions can be seen as relaxations of the Hadamard problem,
and the matrices leading to such bounds could be called almost Hadamard matrices.
However, a different notion of {\em almost Hadamard matrices} was already introduced and considered in \cite{banica1, banica2} (in those papers the emphasis is on various matrix norms and not on the entries of the matrix). Therefore, to make a clear distinction, we prefer to use the terminology {\it flat orthogonal matrices} here.

Finally, circulant matrices play a special role in applications, therefore we can add this as an additional constraint.

\begin{problem}\label{probcircfn}
Let $f_n^{circ}:=\min_{M\in O(n)\cap Circ(n)}f_M$, where $Circ(n)$ denotes the set of $n\times n$ circulant matrices and $f_M$ was defined in \eqref{fm}. Is it true that $f_n^{circ}=o(1)$?
\end{problem}

We firmly believe that the answer to the first three problems is positive, while we are undecided as to the last one. Recall that a conjecture of Ryser asserts that there are no $n\times n$ circulant Hadamard matrices if $n>4$.

\medskip

The remaining of this paper is split into two sections. The first one is devoted to general constructions that lead to bounds on $u_n$ and $l_n$ valid in arbitrary dimension. We then devote the last section to improved bounds when the size of the matrix have various arithmetic properties.

\section{General constructions}

We begin by a simple result which shows that the flatness parameter $f_n$ cannot be expected to be very small. In other words, requiring $f_n$ to be very small is equivalent to requiring a Hadamard matrix of order $n$ to exist.

\begin{proposition}\label{smallerror}
Let $n\geq 3$ and let $\varepsilon>0$ be such that
$$
\varepsilon <\begin{cases}
                 1/n&\mbox{if $n$ is odd}\\
                 2/n&\mbox{if $n$ is even}
                \end{cases}.
$$
 Assume there exists an orthogonal matrix $M$ such that
 for every $j,k=1,\dots, n$,
 \begin{equation}
  \label{eq:almosthadamard}
  \left(\frac{1-\varepsilon}{n}\right)^{1/2}\leq |m_{j,k}|\leq \left(\frac{1+\varepsilon}{n}\right)^{1/2}.
 \end{equation}
Then $n$ is a multiple of $4$ and there exists a Hadamard matrix of order $n$.
\end{proposition}
\begin{proof}
Let $S$ denote the matrix defined by the sign of the entries of $M$, i.e. $s_{j,k}=\frac{1}{\sqrt{n}}\textrm{sign}\,m_{j,k}$. Consider two rows $\vs_j$ and $\vs_k$ of $S$. Let $r$ denote the number of columns where the entries in $\vs_j$ and $\vs_k$ match, and $n-r$ where they differ. Then, for the corresponding rows $\vm_j, \vm_k$ of $M$ we have
\begin{eqnarray*}
0&=&\langle \vm_j,\vm_k\rangle\\
&&\begin{cases}\dst\leq \frac{r}{n}(1+\varepsilon)- \frac{n-r}{n}(1-\varepsilon)
   =\frac{2r-n+\varepsilon n}{n}\\
  \dst \geq \frac{r}{n}(1-\varepsilon)- \frac{n-r}{n}(1+\varepsilon)
   =\frac{2r-n-\varepsilon n}{n}
  \end{cases}.
\end{eqnarray*}
However, if $n$ is odd the interval  $(2r-n-\varepsilon n, 2r-n+\varepsilon n)$ does not contain zero (or any even number), a contradiction. If $n$ is even, the interval  $(2r-n-\varepsilon n, 2r-n+\varepsilon n)$ contains zero if and only if $r=\frac{n}{2}$, in which case the corresponding rows $\vs_j, \vs_k$ are also orthogonal, and we conclude that $S$ is a Hadamard matrix.
\end{proof}

We continue with a simple block construction which proves that $l_n$ is at least as large as $\frac{1}{2\sqrt{n}}$.

\begin{theorem}\label{thm:block}
For any dimension $n$ we have $l_n\ge \frac{1}{2\sqrt{n}}$.
\end{theorem}
\begin{proof}
Assume $n=2^r+q$ where $q<2^r$, and introduce the notation $s=2^r-q$ for brevity. Let $H$ be a Hadamard matrix of order $2^r$ (of course, such a matrix exists, e.g. the tensorial power $F_2^{\otimes r}$). We use the normalization that $H$ has entries $\pm\frac{1}{\sqrt{2^r}}$ instead of $\pm 1$. Let $\tilde{M}$ be the extension of $H$ by an identity matrix of order $q$ in the lower right. Split $\tilde{M}$ into blocks of size $s$ and $q$ and $q$ as follows (the indexes simply indicating the sizes of the blocks):
$$
\tilde{M}=\left[
\begin{array}{c|c|c}
H_{s,s} & H_{s,q} & 0\\ \hline
H_{q,s} & H_{q,q} & 0\\ \hline
0 & 0 & I_q
\end{array}\right].
$$

Let $U$ denote the following orthogonal block-matrix:
$$
U=\left[
\begin{array}{c|c|c}
I_{s} & 0 & 0\\ \hline
0 & \frac{1}{\sqrt{2}}I_{q} & -\frac{1}{\sqrt{2}}I_{q}\\ \hline
0 & \frac{1}{\sqrt{2}}I_{q} & \frac{1}{\sqrt{2}}I_{q}
\end{array}\right],
$$
and let $M:=U^T \tilde{M} U$. A direct calculation shows that
$$
M=\left[
\begin{array}{c|c|c}
H_{s,s} & \frac{1}{\sqrt{2}}H_{s,q} & -\frac{1}{\sqrt{2}}H_{s,q}\\ \hline
\frac{1}{\sqrt{2}}H_{q,s} & \frac{1}{2}(H_{q,q}+I_q) & \frac{1}{2}(-H_{q,q}+I_q)\\ \hline
-\frac{1}{\sqrt{2}}H_{q,s} & \frac{1}{2}(-H_{q,q}+I_q) & \frac{1}{2}(H_{q,q}+I_q)
\end{array}\right].
$$
The smallest modulus among the entries of $M$ is $\frac{1}{2\sqrt{2^r}}\ge \frac{1}{2\sqrt{n}}$, which proves the theorem. (Note, however, that the largest appearing modulus is $\frac{1}{2}+\frac{1}{\sqrt{2^r}}$, so this construction gives no indication with respect to Problem \ref{probfn}.)
\end{proof}

Proving that $u_n\le \frac{c}{\sqrt{n}}$ is also fairly easy, as one can use a block-diagonal construction.

\begin{theorem}\label{thm:un}
For any dimension $n$ we have $u_n\le \frac{2+o(1)}{\sqrt{n}}$.
\end{theorem}
\begin{proof}
If $A_1\in O(n_1)$ and $A_2\in O(n_2)$ are orthogonal matrices then the block-diagonal matrix $A_1\oplus A_2$ is an orthogonal matrix of order $n_1+n_2$. This implies $u_{n+m}\le \max\{u_n, u_m\}$. Also, Propositions \ref{had4km1} and \ref{n4kp1} below show that $u_p=\frac{1+o(1)}{\sqrt{p}}$ whenever the dimension $p$ is a prime. For general $n$ we must invoke the following weak-Goldbach type result from number theory: for every $\e>0$, every large enough odd number $n$ can be written as a sum of three primes $n=p_1+p_2+p_3$ where each $p_i$ lies in the interval $[(1-\e)n/3, (1+\e)n/3]$, while every large enough even number $n$ can be written as a sum of four primes $n=p_1+p_2+p_3+p_4$ where each $p_i$ lies in the interval $[(1-\e)n/4, (1+\e)n/4]$. This follows from the method of Vinogradov \cite{vin}. Hence, the block-diagonal construction together with this weak-Goldbach type result implies $u_n\le \frac{\sqrt{3}+o(1)}{\sqrt{n}}$ if $n$ is odd, and $u_n\le \frac{2+o(1)}{\sqrt{n}}$ if $n$ is even.
\end{proof}

\begin{remark}\rm
In connection with $l_n$ and $u_n$ it is natural to examine how the entries of a typical random orthogonal matrix behave.
It is well known that for any coordinate $z$ of a random unit vector we have $Pr(|z|>\frac{t}{\sqrt{n}})\le e^{-t^2/2}$
(this can be seen by upper bounding the area of a spherical cap of radius $r$ by that of a sphere of radius $r$).
This implies, via a simple calculation, that $u_M\le \frac{c\sqrt{\log n}}{\sqrt{n}}$ for a random orthogonal matrix $M$,
with high probability. Therefore, a random orthogonal matrix is typically ``not far'' from the bound $u_M\le \frac{c}{\sqrt{n}}$
given in Theorem \ref{thm:un}. On the contrary, the lowest absolute value in a random unit vector is smaller than $\frac{c}{n}$
with high probability (this can be seen by generating a random unit vector as the normalized vector of $n$ independent
Gaussian random variables). The same holds, a fortiori, for a random orthogonal matrix, which shows that the lowest
entry is typically very far from the optimal bound given in Theorem \ref{thm:block}.
\hfill $\square$
\end{remark}

\begin{remark}\rm
If we assume that the Hadamard conjecture holds then $l_n=u_n=\frac{1}{\sqrt{n}}$ for all $n$ divisible by 4. By the simple construction of Proposition \ref{had4km1} we also conclude that for $n\equiv 3 \ (mod \ 4)$ the quantities $l_n$ and $u_n$ are both of the magnitude $(1+o(1))\frac{1}{\sqrt{n}}$. However, we could not prove such a statement for the case $n\equiv 1, \ 2 \ (mod \ 4)$.
\end{remark}

\section{Specific constructions}

In the rest of this note we give some positive partial results with respect to Problems \ref{probln}, \ref{probfn}, \ref{probcircfn}. As noted in \cite[Section 3]{banica1} symmetric balanced incomplete block designs give rise to orthogonal matrices with two entries. If the parameters of the block design are suitable then the entries will be close to $\pm \frac{1}{\sqrt{n}}$. Namely, we have the following special case.

\begin{proposition}\label{had4km1}
(i) If the dimension $n$ is such that a Hadamard matrix $H$ of size $(n+1)\times (n+1)$ exists, then there exists an orthogonal matrix $M$ of size $n\times n$ with all entries having modulus $(1+o(1))\frac{1}{\sqrt{n}}$. In particular, this is the case if $n=p^r$ where $p=4k-1$ is a prime and $r$ is odd.

(ii) If the dimension $n$ is a prime of the form $4k-1$, then $M$ can be chosen to be circulant.
\end{proposition}
\begin{proof}
(i) We can assume that the Hadamard matrix $H$ is in standard form, i.e. the first row and column of $H$ consist of 1's. Delete the first row and column of $H$, and in the remaining matrix $H'$ replace the entries $+1$ and $-1$ by the variables $x$ and $y$, respectively. Orthogonality of any two rows of $H'$ is now equivalent to $\frac{n-3}{4}x^2+\frac{n+1}{4}y^2+ \frac{n+1}{2}xy=0$, and the unit length of the rows of $H'$ is ensured by $\frac{n-1}{2}x^2+\frac{n+1}{2}y^2=1$. This system of equations admits the (non-unique) solution
$$x=\frac{-1}{\sqrt{n+1}-1}$$
$$y=\frac{2}{\sqrt{n+1}}-\frac{1}{\sqrt{n+1}-1},$$
and we define $M$ as the $n\times n$ matrix with these values. Note that both $x$ and $y$ have the order of magnitude $(1+o(1))\frac{1}{\sqrt{n}}$. In fact, the error term $o(1)$ here has the order of magnitude $O(\frac{1}{\sqrt{n}})$. We remark that $M$ corresponds to the {\it Hadamard design} associated to $H$, and this construction is a special case of the one described in \cite[Section 3]{banica1}. If $n=p^r$ where $p=4k-1$ is a prime and $r$ is odd, then a Hadamard matrix of size $(n+1)\times (n+1)$ exists by the Paley construction.

(ii) If $n$ is a prime of the form $4k-1$ then the Paley construction essentially leads to a circulant matrix. Namely, consider the following circulant matrix $M$ corresponding to the quadratic character of $\mathbb{F}_p$: $[M]_{i,j}=x$ if $i-j$ is a quadratic residue, and $[M]_{i,j}=y$ if $i-j$ is a non-residue or zero. The above values of $x,y$ ensure orthogonality of $M$.
\end{proof}

When $n$ is a prime of the form $4k+1$, the construction is much less trivial, as described below.

\begin{proposition}\label{n4kp1}
If the dimension is a prime $p=4k+1$, then there exists a circulant orthogonal matrix $M$ of size $p\times p$ with all entries having modulus $(1+o(1))\frac{1}{\sqrt{p}}$.
\end{proposition}
\begin{proof}
For this proof it will be convenient to identify the cyclic group $\ZZ_p$ with the numbers $(\frac{-p+1}{2}, \dots, \frac{p-1}{2})$, and use coordinates of vectors accordingly. Also, let $Q, NQ\subset \ZZ_p$ denote the set of quadratic residues and non-residues, respectively (0 is not included in either $Q$ or $NQ$). Consider the vector $\vv=(v_\frac{-p+1}{2}, \dots, v_\frac{p+1}{2})$ given by the quadratic character, i.e. $v_0=0$, and for nonzero $j$ we have $v_j=\pm 1$ according to whether $j\in Q$ or $j\in NQ$. Note that $\vv$ is symmetric, $v_j=v_{-j}$, because $p=4k+1$. Note also that $\hat \vv(k)=\sum_{j} v_j e^{2\pi ijk/p}=\pm \sqrt{p}$ for $k\ne 0$, and $\hat \vv(0)=0$. We will need a random modification of $\vv$. The construction is analogous to the one given in \cite[Theorem 9.2]{ruzsmat}.

\medskip

Fix $\e>0$, and let $\rho=\frac{1}{\sqrt{p}}$. As $NQ$ is symmetric, it can be written as a disjoint union of a set $H$ and its negative, $NQ=H\cup -H$, in other words $NQ=\cup_{y\in H} \{y, -y\}$. For each $y\in H$ consider independent random variables $\xi_y$ such that $\p(\xi_y=1/2)=\rho$, $\p(\xi_y=0)=1-\rho$. For $y=0$ set $\p(\xi_y=1/4)=\rho$, $\p(\xi_y=0)=1-\rho$. Consider the random vector $\vw=(w_\frac{-p+1}{2}, \dots, w_\frac{p-1}{2})$ given by $w_y=\xi_y$ if $y\in H$ or $y=0$, $w_y=\xi_{-y}$ if $-y\in H$, and $w_y=0$ if $y\in Q$. Let us evaluate the Fourier transform of the random vector $\vw$.

\begin{equation}
\E(\hat \vw(k))=\frac{1}{2}\rho \left(\frac{1}{2}+\sum_{y\in NQ}e^{2\pi i yk/p}\right)= \begin{cases}
    \sqrt{p}/4 & \text{if } k=0, \\
    \pm 1/4 & \text{if } k\ne 0
  \end{cases}
\end{equation}
and, for all $k$,
\begin{equation}
\D^2(\hat \vw(k))=\frac{1}{4}\rho(1-\rho) \left(\frac{1}{4}+\sum_{y\in H}(e^{2\pi i yk/p}+e^{-2\pi i yk/p})^2\right)\le
    \frac{1}{4}\sqrt{p}.
\end{equation}

We invoke here a large deviation inequality of Chernov as stated in \cite[Theorem 1.8]{taovu06}: let $X_1, \ldots, X_n$ be independent random variables satisfying $| X_i - \E(X_i)| \leq 1 $ for all $i$. Put
  $ X = X_1 + \ldots + X_n$ and let $\sigma^2$ be the variance of $X$. For any $t>0$ we have
   \[ \p( |X - \E(X)| \geq t\sigma) \leq 2 \max \left( e^{-t^2/4}, e^{-t\sigma/2} \right) .\]

\medskip

Using this estimate with $t=p^\e$ we obtain

\begin{equation}
\p\left(|\hat \vw(0)-\sqrt{p}/4|\ge \frac{1}{2}p^{\frac{1}{4}+\e}\right)\le 2e^{-\frac{p^{2\e}}{4}},
\end{equation}

\begin{equation}
\p\left(|\hat \vw(k) \mp \frac{1}{4}|\ge \frac{1}{2}p^{\frac{1}{4}+\e}\right)\le 2e^{-\frac{p^{2\e}}{4}} \ \text{for all } k\ne 0.
\end{equation}

Therefore, with high probability none of the above events occur, and we have $|\hat \vw(0)-\sqrt{p}/4|\le O(p^{\frac{1}{4}+\e})$ and $|\hat \vw(k)|\le O(p^{\frac{1}{4}+\e})$. Fix such a favourable vector $\vw$. We can also assume without loss of generality that $\xi_0=1/4$ and hence $w_0=1/4$ (we are free to change $w_0$ from 0 to $1/4$, if necessary, without altering the order of magnitude of $\hat \vw$).

\medskip

Finally, consider the vector $\vz=\vv+4\vw$, and let $\vu=\frac{1}{p}\hat{\vz}$. The vector $\vz$ is unimodular (some of the $-1$ entries in $\vv$ were changed to $+1$ and the value at 0 was changed to $+1$). Therefore the circulant matrix $M$ with first row $\vu$ is orthogonal. Also, the entries of $\vu$ are all of absolute value $\frac{1}{\sqrt{p}}(1+p^{-\frac{1}{4}+\e})$, by construction.

\medskip

We remark that a similar construction works if the dimension $p$ is of the form $4k-1$ but the result is inferior to Proposition \ref{had4km1} in the sense that the error term is larger.
\end{proof}

By combining the results of the propositions above we can answer Problem \ref{probfn} and \ref{probcircfn} for dimensions $n$ which are composed of large prime factors.

\begin{corollary}\label{cor1}
For any fixed $m$ let $\NN_m$ denote the set of positive integers $n=2^s p_1 \dots p_r$ such that each odd prime factor $p_j\ge n^{1/m}$ (the primes may appear with multiplicity). \\
(i) For every $m$ and every $n\in \NN_m$ there exists an orthogonal matrix $M$ of size $n\times n$ with all entries having modulus $(1+o_m(1))\frac{1}{\sqrt{n}}$. \\
(ii) If all the odd primes $p_j$ appearing in the factorization of $n$ are distinct and $s=0$ or $s=2$ then there exists a circulant orthogonal matrix $M$ of size $n\times n$ with all entries having modulus $(1+o_m(1))\frac{1}{\sqrt{n}}$.
\end{corollary}
\begin{proof}
(i) Let $n=2^sp_1 \dots p_r\in \NN_m $. Then $r\le m$ by the definition of $\NN_m$. Let $H$ denote a Hadamard matrix of size $2^s$ (such matrix exists, of course, e.g. $H=F_2^{\otimes s}$), and let $M_{p_j}$ denote the matrices corresponding to the primes $p_j$, as constructed in Propositions \ref{had4km1} and \ref{n4kp1}. Let $M$ be the tensorial product of all the matrices $H$ and $M_{p_j}$ (where $M_{p_j}$ is taken with the same multiplicity as $p_j$ in $n$). As $r\le m$ for any $n\in \NN_m$, the errors do not accumulate, and we still have $|[M]_{i,j}|=(1+o_m(1))\frac{1}{\sqrt{n}}$.

\medskip

(ii) The matrices $M_{p_j}$ are circulant by the constructions of Propositions \ref{had4km1} and \ref{n4kp1}. Let $\vx_j$ denote the first row of $M_{p_j}$.  If $s=0$, $n=p_1 \dots p_r$ with all $p_j$ distinct, then $\ZZ_n\equiv \ZZ_{p_1}\dots \ZZ_{p_r}$. Let $\vk=(k_1, \dots, k_r)\in \ZZ_{p_1}\dots \ZZ_{p_r}$ and let $\vy(\vk)=\prod_{j=1}^r \vx_j(k_j)$. Then $\vy$ generates a circulant matrix which is orthogonal (because all $\hat \vx_j$ are unimodular, and hence so is $\hat \vy$), and the entries of $\vy$ are of absolute value $(1+o_m(1))\frac{1}{\sqrt{n}}$. When $s=2$ we can incorporate the $4\times 4$ real circulant Hadamard matrix in the same manner.
\end{proof}

\begin{remark}\rm
Problem \ref{probcircfn} concerning the circulant case has an interesting connection to {\it ultraflat polynomials}. It is well-known that $\vx=(x_1, \dots, x_n)$ generates a circulant orthogonal matrix if and only if the Fourier transform $\hat \vx=(w_1, \dots, w_n)$ is unimodular on $\hat \ZZ_n$, i.e. $|w_j|=|w_k|$ for all $j, k$. If $\vx$ is real, then $\vw$ is conjugate symmetric, i.e. $w_j=\overline w_{n-j}$. We also want that all $|x_j|\approx \frac{1}{\sqrt{n}}$. Considering $w_1, \dots, w_n$ as variables we are led to the problem of constructing a polynomial $P(z)=\sum_{j=1}^n w_jz^j$ where $w_j=\overline w_{n-j}$ and $|w_j|=1$ (after re-normalization), such that $P(z)$ is ``flat'' at the $n$th roots of unity, i.e. $|P(\omega^j)|=(1+o(1))\sqrt{n}$, where $\omega=e^{2i\pi/n}$ and $j=0, \dots n-1$. Dropping the restriction $w_j=\overline w_{n-j}$ one can require $P(z)$ to be flat all over the unit circle (i.e. $|P(z)|=(1+o(1))\sqrt{n}$ for all $|z|=1$), and such polynomials are called {\it ultraflat}. The existence of ultraflat polynomials was proven by Kahane \cite{kahane}. However, the extra condition $w_j=\overline w_{n-j}$ prevents $P(z)$ from being ultraflat as shown by Remark 5.1 in \cite{erdelyi}: the restriction $w_j=\overline w_{n-j}$ implies $\max_{|z|=1} |P(z)|\ge (1+\e)\sqrt{n}$ with $\e=\sqrt{4/3}-1$. However, this does not answer Problem \ref{probcircfn} because we require $P(z)$ to be flat only at the $n$th roots of unity. Problem \ref{probcircfn} can therefore be regarded as a discretized version (at the $n$th roots of unity) of the question of existence of ultraflat polynomials with the self-conjugacy restriction $w_j=\overline w_{n-j}$.

Notice also that Corollary \ref{cor1} gives an affirmative answer to Problem \ref{probcircfn} for some dimensions $n$ (e.g. when $n=4p$, $p$ being a prime). Therefore, trying to prove Ryser's conjecture on the non-existence of circulant Hadamard matrices by giving a negative answer to Problem \ref{probcircfn} cannot possibly work. In the other direction, hoping to construct a circulant Hadamard matrix by first constructing a circulant flat orthogonal matrix and then modifying its entries is also rather naive and hopeless in our opinion. \hfill $\square$
\end{remark}

We saw in Proposition \ref{had4km1} that given any Hadamard matrix $H$ one can reduce the dimension by 1, and construct a flat orthogonal matrix. It is natural to try also to increase the dimension by 1.
The general construction given in Proposition \ref{smallerror} allows to do so, but will introduce an entry of size $1/2$, thus destroying
the flatness. Our last specific construction concerns the increase of dimension by $1$ without destroying flatness, but it only works under some restrition on the dimension. We recall that a Hadamard matrix is called {\it regular} if the row sums and column sums of $H$ are all equal.

\begin{proposition}\label{regular}
Assume the dimension $n$ is such that a regular Hadamard matrix $H$ of order $n$ exists (this implies $n=4k^2$ for some $k$). Then there exists an orthogonal matrix $M$ of size $(n+1)\times (n+1)$ with all entries having modulus $(1+o(1))\frac{1}{\sqrt{n+1}}$.
\end{proposition}

\begin{proof}
Let $b=\frac{1-(n+1)^{-1/2}}{n}$ and replace the positive entries of $H$ by $\frac{1}{\sqrt{n}}-b$, while the negative entries with $-\frac{1}{\sqrt{n}}-b$ (the sign of $b$ is purposefully negative in both cases). Next, extend this matrix by a new row and column filled with entries $a=\frac{-1}{\sqrt{n+1}}$, and let $M$ be the arising matrix. The entries of $M$ are of modulus $(1+o(1))\frac{1}{\sqrt{n+1}}$ (in fact, the error term is of the order $\frac{1}{\sqrt{n}}$), and an easy calculation shows that $M$ is orthogonal.
\end{proof}

\medskip

We end this note by emphasizing that Problems \ref{probln}, \ref{probfn}, \ref{probcircfn} all remain open for dimensions $n$ with small prime factors, i.e. for dimensions $n$ which are not covered by Corollary \ref{cor1}.