\section{Introduction}

Observations of type Ia supernovae, cosmic microwave background (CMB)
radiation anisotropy, and large-scale structure of galaxies have revealed that
our present universe is dominated by two dark components \cite{Recent-Obs}.
Based on the Friedmann world model, the energy content
of the present universe is occupied by 73\% of dark energy and 23\% of
dark matter, whose nature is the fundamental mystery of the present day
theoretical physics and cosmology
(see Ref.\ \cite{Komatsu-etal-2011} for recent cosmological interpretation
of astronomical observations).

The dark matter is often modeled as pressureless, non-relativistic particles,
which is known as cold dark matter (CDM).
Although the CDM model has been greatly successful in explaining many
cosmological observations, there remain certain conflicts at the galactic
scales. First, the CDM-based simulations of galactic halos show density
profiles with a central cusp at kpc scale, while observations indicate
constant density cores \cite{vanEymeren-etal-2009}.
Secondly, the CDM model predicts overpopulation of low mass halos within
a galaxy or a galaxy group, which is an order of magnitude larger than
the present observations \cite{Klypin-etal-1999}.

These CDM problems at galactic scales may be overcome by introducing
a suppression of small-scale power in the density perturbation.
There have been efforts to explain the CDM conflicts using the
dark matter models with the Bose-Einstein condensate (BEC) and
the ultra-light scalar field.

In the BEC dark matter model the bosonic particles occupy the lowest
quantum state of the external potential and can be described as a
coherent matter wave of macroscopic size due to Bose enhancement
\cite{Sin-1994,Lee-1996-2010,BEC-refs,Rodriguez-Montoya-2010,
Kain-2012,Chavanis-2012}.
Along this line there is fuzzy cold dark matter model proposed
by Hu et al.\ \cite{Hu-2000} which is a case of free particles
ignoring the self-interaction terms in the BEC based on the generalized
dark matter approach \cite{Hu-1998}.
In this model the coherent wave of ultra-light dark matter with
$m \sim 10^{-22}~\textrm{eV}$ can suppress the kpc-scale cusp in the dark matter
halos and reduce the abundance of low mass halos.
Recent investigations on the growth of perturbation in the universe
with BEC dark matter can be found in Refs.\ \cite{Kain-2012,Chavanis-2012}.
The BEC dark matter has a close connection to the axion dark matter model.
Recently, Sikivie and Yang \cite{Sikivie-2009} showed that
the cold dark matter axions form a BEC through their self-interactions
and gravitational interactions (see also Ref.\ \cite{Erken-2012}).

Another possibility is the ultra-light scalar field dark matter model
\cite{Peebles-1999,Matos-1999,Sahni-2000,Matos-2000-2001,Arbey-2001,
Matos-2009,Marsh-2010,Marsh-2011,Suarez-2011,Magana-2012}
(see Ref.\ \cite{Magana-2012-review} for a review).
Based on a specific scalar field potential with appropriate initial conditions
the evolution of scalar field mimics the property of the CDM fluid.
For example, authors of Refs.\ \cite{Sahni-2000,Matos-2000-2001}
proposed $\cosh$-like potential $V=V_0 (\cosh \lambda\phi-1)$
as a good candidate of dark matter (see also Ref.\ \cite{Matos-2009}).
Authors of Refs.\ \cite{Suarez-2011,Magana-2012} investigated the structure
formation in the scalar-field dark matter model with
$V=\frac{1}{2}m^2\phi^2+\frac{1}{4}\lambda\phi^4$ by using the fluid and
the field approaches. Authors of Refs.\ \cite{Marsh-2010,Marsh-2011}
studied the effect of ultra-light scalar field with $V=\frac{1}{2}m^2\phi^2$
on the growth of structure of the universe,
in a situation where only a substantial fraction of ultra-light scalar field
dark matter constitutes the total dark matter including CDM.

It is well known that the axion as a coherently oscillating scalar field
behaves as CDM in cosmologically relevant scales \cite{Axion-refs}.
Khlopov et al.\ \cite{Khlopov-etal-1985} investigated the growth of
perturbations due to gravitational instability in a universe dominated
by a scalar axion field.
The quantum description of cosmological axion perturbations was given
by Nambu and Sasaki \cite{Nambu-1990}, while the full relativistic description
based on the linear perturbation theory was presented in Refs.\
\cite{Ratra-1991,Hwang-1997,Hwang-2009}. For a typical mass with
$10^{-5}~\textrm{eV}$ the behavior of the axion is equivalent to that of
CDM on the cosmological scales including the superhorizon scale \cite{Hwang-2009}.

In this work, we consider an axion dark matter model
with extremely low mass ($m \le 10^{-22}~\textrm{eV}$) and investigate
the effect of ultra-light axion on the baryon matter density and the CMB
anisotropy power spectra and the perturbation growth.
Our calculations are based on the full relativistic linear cosmological
perturbation theory. Note that some of previous studies on BEC dark matter
models relies on Newtonian cosmology and the relativistic description
was given only in a limited fashion.
Due to the common property of coherent oscillations,
our description for the low-mass axion is closely related with that from
the BEC or ultra-light scalar field dark matter models.

The structure of this paper is as follows.
Section \ref{sec:basic_equations} reviews the basic equations
used to describe background and perturbation evolutions based on the
relativistic linear perturbation theory, dealing with the multi-component
fluids together with the minimally coupled scalar field.
In Sec. \ref{sec:axion} the basic properties of the axion dark matter are
briefly described. A complete set of perturbation equations for the axion,
baryon, and radiation fluids is presented, and initial conditions for
the perturbation variables in the early radiation era are derived for
two asymptotic limits. Our primary results are summarized in Sec.\
\ref{sec:results}, including baryon matter density and CMB power spectra
and the growth of baryon density perturbation for various values of
extremely light axion mass. We also compare the cosmological effects
of light axion with the CDM model mixed with the massive neutrinos
as a hot dark matter component. We discuss our results in Sec.\
\ref{sec:discussion}. Throughout this paper, we set $c\equiv 1 \equiv \hbar$.

\section{Basic equations}
\label{sec:basic_equations}

We consider scalar-type perturbations with the metric \cite{Bardeen-1988}
\bea
   & & ds^2 = -(1+2\alpha)dt^2-2a\beta_{,\alpha}dtdx^\alpha
   \nonumber \\
   & & \qquad
         + a^2 \left[(1+2\varphi)g_{\alpha\beta}^{(3)}
         + 2 \gamma_{,\alpha|\beta} \right] dx^\alpha dx^\beta,
\eea
and introduce
\bea
   \chi\equiv a(\beta+a\dot\gamma), \quad
       \kappa \equiv 3H\alpha -3\dot\varphi-{\Delta \over a^2} \chi.
\eea
Here $a(t)$ is the cosmic scale factor, $\alpha$, $\beta$,
$\varphi$, $\gamma$ are metric perturbations, a dot and a comma indicate
the time and spatial derivatives, respectively, and a vertical bar
represents the covariant derivative based on the comoving three-space metric
$g_{\alpha\beta}^{(3)}$.
The energy-momentum tensor is
\bea
   & & T^0_0 = - \mu - \delta \mu, \quad
       T^0_\alpha = - {1 \over k} \left( \mu + p \right) v_{,\alpha},
   \nonumber \\
   & & T^\alpha_\beta = \left( p + \delta p \right) \delta^\alpha_\beta,
\eea where we ignore the anisotropic stress;
$k$ is the wave number with $\Delta \equiv - k^2$ in Fourier space.

The background evolution is described by
\begin{eqnarray}
   && H^2
       = \frac{8\pi G}{3} \mu -\frac{K}{a^2} + \frac{\Lambda}{3}, \nonumber \\
   && \dot\mu + 3H(\mu+p)=0,
\end{eqnarray}
where $H \equiv \dot a/a$ is the Hubble parameter, $K$ is the sign of spatial curvature, and
$\Lambda$ is the cosmological constant.
In the multiple component case we have
\begin{eqnarray}
  && \mu=\sum_j{\mu_j}, \quad p=\sum_j{p_j},
\end{eqnarray}
and
\bea
  && \dot\mu_i + 3H(\mu_i+p_i)=0,
\eea
where we ignore direct interactions among fluids.
For a minimally coupled scalar field we have the equation of motion
\bea
  & & \ddot\phi+3H\dot\phi+V_{,\phi}=0,
\eea
and the fluid quantities
\bea
  & & \mu_\phi=\frac{1}{2}\dot\phi^2+V,\quad
       p_\phi  =\frac{1}{2}\dot\phi^2-V.
\eea
We need the following equations for the perturbations \cite{Bardeen-1988,Hwang-1991}
\bea
   & & \dot\kappa+2H\kappa-4\pi G(\delta\mu+3\delta p)
         +\left(3\dot{H}-\frac{k^2}{a^2}\right)\alpha = 0,
   \label{eq:kappa_dot} \\
   & & \delta\dot\mu_i+3H(\delta\mu_i+\delta p_i)
       -(\mu_i+p_i)\left(\kappa-3H\alpha-\frac{k}{a}v_i\right)
   \nonumber \\
   & & \qquad
         =0,
   \label{eq:mu_dot} \\
   & & \frac{1}{a^4} \left[a^4 (\mu_i+p_i)v_i \right]^{\displaystyle\cdot}
        = \frac{k}{a} \left[ \left( \mu_i + p_i \right) \alpha
        + \delta p_i \right],
   \label{eq:v_dot}
\eea
where
\bea
   & & \delta\mu=\sum_j \delta\mu_j, \quad \delta p=\sum_j \delta p_j,
   \nonumber \\
   & & (\mu+p)v=\sum_j (\mu_j+p_j) v_j.
\eea
In the presence of a scalar field we have the equation of motion \cite{Hwang-1991}
\bea
   & & \delta\ddot\phi+3H\delta\dot\phi +\frac{k^2}{a^2}\delta\phi
         +V_{,\phi\phi}\delta\phi
   \nonumber \\
   & & \qquad
         =2\ddot\phi\alpha
         +\dot\phi \left(\dot\alpha+6H\alpha+\frac{k^2}{a^2}\chi
         -3\dot\varphi\right),
    \label{eq:field_eom}
\eea
and the fluid quantities
\bea
   & &  (\mu_\phi+p_\phi)v_\phi=\frac{k}{a}\dot\phi\delta\phi, \quad
          \delta\mu_\phi=\dot\phi\delta\dot\phi-\dot\phi^2\alpha
                       +V_{,\phi}\delta\phi,
   \nonumber \\
   & & \delta p_\phi=\dot\phi\delta\dot\phi-\dot\phi^2\alpha-V_{,\phi}\delta\phi,
   \label{eq:field_fluid}
\eea
with vanishing anisotropic stress.

\section{Axion}
\label{sec:axion}

We consider the axion in the presence of baryonic dust and radiation (photons and neutrinos).
The axion is considered as a minimally coupled scalar field with a potential
$V=\frac{1}{2}m^2\phi^2$, where $m$ is the axion mass.
We strictly ignore $H/m$ higher order term in our analysis.
At the current epoch it is
\begin{equation}
   \frac{H_0}{m} = 2.133\times 10^{-28} h
               \left(\frac{m}{10^{-5}~\textrm{eV}}\right)^{-1},
\end{equation}
where $H_0 \equiv 100 h~\textrm{km}\textrm{s}^{-1}\textrm{Mpc}^{-1}$ is
the Hubble expansion rate at present.

In the axion model the equation of motion leads to the background solution
as \cite{Ratra-1991}
\begin{equation}
   \phi(t)=a^{-3/2} [\phi_{+0}\sin(mt)+\phi_{-0}\cos(mt)],
\end{equation}
where $\phi_{+0}$ and $\phi_{-0}$ are constant coefficients.
We use the fluid formulation of axion.
We have
\begin{equation}
   \mu_a=\frac{m^2}{2a^3} (\phi_{+0}^2+\phi_{-0}^2), \quad p_a=0,
\end{equation}
thus, the axion behaves as a zero-pressure fluid; for the fluid quantities
the time-averaging has been applied for highly oscillating scalar field.
For baryon, radiation (as a fluid), and axion, we have
\begin{equation}
   \mu=\mu_b + \mu_r + \mu_a, \quad p=\frac{1}{3}\mu_r.
\end{equation}

To perturbed order, with an ansatz \cite{Ratra-1991}
\begin{equation}
   \delta\phi(\bik,t)=\delta\phi_{+}(\bik,t)\sin(mt)
                     +\delta\phi_{-}(\bik,t)\cos(mt),
\end{equation}
the fluid quantities in Eq.\ (\ref{eq:field_fluid}) give \cite{Hwang-2009}
\bea
  & & \mskip-36mu \delta\mu_a =
        a^{-3/2} m [ m(\phi_{+0}\delta\phi_{+}+\phi_{-0}\delta\phi_{-}),
  \nonumber \\
  & & \mskip-36mu \qquad
        +\frac{1}{2}(\phi_{+0}\delta\dot\phi_{-}
        -\phi_{-0}\delta\dot\phi_+)]-\mu_a\alpha,
  \nonumber \\
  & & \mskip-36mu \delta p_a = \frac{1}{2} a^{-3/2} m
        (\phi_{+0}\delta\dot\phi_{-}-\phi_{-0}\delta\dot\phi_{+})
        -\mu_a \alpha,
  \nonumber \\
  & & \mskip-36mu {a \over k} (\mu_a+p_a)v_a=\frac{1}{2}a^{-3/2}m
        (\phi_{+0}\delta\phi_{-}-\phi_{-0}\delta\phi_{+}),
\label{eq:field_fluid_avg}
\eea
where for the fluid quantities we have taken the time average.
Up to this point the perturbed quantities are spatially gauge invariant and
the temporal gauge condition has not been taken yet \cite{Bardeen-1988}.

\subsection{Axion-comoving gauge}

As the temporal gauge (hypersurface or slicing) condition we
take the axion-comoving gauge \cite{Hwang-2009}
\begin{equation}
   \left( \mu_a + p_a \right) v_a \equiv 0.
\end{equation}
Equation (\ref{eq:field_fluid_avg}) gives
\begin{equation}
   \frac{\delta\phi_+}{\phi_{+0}}=\frac{\delta\phi_{-}}{\phi_{-0}}
\end{equation}
and
\begin{equation}
   \frac{\delta\mu_a}{\mu_a}
     =2 a^{3/2} \frac{\delta\phi_+}{\phi_{+0}}-\alpha,
   \quad
   \frac{\delta p_a}{\mu_a}=-\alpha.
\label{eq:field_fluid_v}
\end{equation}
From Eqs.\ (\ref{eq:kappa_dot})
and (\ref{eq:mu_dot}) for the axion, we have
\bea
    & & \dot\kappa+2H\kappa
          =4\pi G(\mu_b\delta_b+2\mu_r\delta_r+\mu_a\delta_a)
    \nonumber \\
    & & \qquad
          + \left[12\pi G \left(\mu_b+\frac{4}{3}\mu_r \right)
          +\frac{k^2-3K}{a^2} \right]\alpha,
    \label{eq:kappa_dot2} \\
    & & \dot\delta_a=\kappa,
\label{eq:delta_dot2}
\eea
where $\delta_i\equiv \delta\mu_i/\mu_i$.
For the baryon and radiation, under the fluid formulation,
Eqs.\ (\ref{eq:mu_dot}) and (\ref{eq:v_dot}) give
\begin{eqnarray}
   && \dot\delta_b = \kappa-3H\alpha-\frac{k}{a}v_b,
   \\
   && \dot{v}_b+Hv_b=\frac{k}{a}\alpha,
   \\
   && \dot\delta_r=\frac{4}{3}\left(\kappa-3H\alpha-\frac{k}{a}v_r\right),
   \\
   && \dot{v}_r = \frac{k}{a}\left(\alpha+\frac{1}{4}\delta_r\right).
   \label{eq:v_dot2}
\end{eqnarray}
For a zero-pressure fluid (like a CDM) we have $\alpha=-\delta p/\mu=0$.
Now, in the case of axion, $\alpha$ can be determined in terms of $\delta_a$
using Eq.\ (\ref{eq:field_eom}). From Eq.\ (\ref{eq:field_fluid}) we have
\begin{equation}
   \delta_a=2 a^{3/2}\frac{\delta\phi_+}{\phi_{+0}}-\alpha, \quad
   \alpha=-\frac{a^{3/2}}{2m^2}\frac{k^2}{a^2}\frac{\delta\phi_+}{\phi_{+0}}.
\end{equation}
From these two relations, we have \cite{Hwang-2009}
\begin{equation}
   \alpha=-\frac{k^2}{4m^2a^2}
           \frac{1}{1+\frac{k^2}{4m^2a^2}}\delta_a.
\label{eq:alpha}
\end{equation}
Equations (\ref{eq:kappa_dot2})--(\ref{eq:v_dot2}) and
(\ref{eq:alpha}) constitutes the complete set for axion, baryon, and radiation.
In a realistic situation we have to use the Boltzmann equations
for radiation (photons and neutrinos) or the tight-coupling approximation
for photon-baryon fluid \cite{Ma-Bertschinger-1995,Hwang-2001}. We note that the above equations are valid
in the presence of $K$ and $\Lambda$ in the background.

From Eqs.\ (\ref{eq:kappa_dot2}), (\ref{eq:delta_dot2})
and (\ref{eq:alpha}) we have
\bea
   & & \ddot \delta_a + 2 H \dot\delta_a
         - 4 \pi G \left( \mu_a \delta_a
         + \mu_b \delta_b + 2 \mu_r \delta_r \right)
   \nonumber \\
   & & \qquad
         + \left[ 12 \pi G \left( \mu_b + {4 \over 3} \mu_r \right)
         + \frac{k^2-3K}{a^2} \right]
   \nonumber \\
   & & \qquad
         \times
         \frac{k^2}{4m^2a^2}
         \frac{1}{1+\frac{k^2}{4m^2a^2}}\delta_a
         = 0.
\eea
The equation is valid at {\it all} scales.
Ignoring $K$, and for axion only, we obtain the axion Jeans scale,
\begin{equation}
   \lambda_J \equiv \frac{2\pi a}{k_J}
      = \frac{2\pi}{\sqrt{2\pi G\mu_a
                         +\sqrt{4\pi^2 G^2 \mu_a^2+16\pi G\mu_a m^2}}}.
\label{eq:Jeans_scale_general}
\end{equation}
In the limit of $k^2 / (4m^2a^2) \ll 1$, with present
$\Omega_{a} \equiv 8 \pi G \mu_a / (3H^2) \simeq 0.27$,
the current Jeans scale becomes
\cite{Khlopov-etal-1985,Nambu-1990,Sikivie-2009,Hwang-2009}
\begin{equation}
\begin{split}
   \lambda_J
       =\left(\frac{\pi^3}{G\mu_{a0} m^2}\right)^{1/4}
       = 50 h^{-1/2} \left(\frac{m}{10^{-5}~\textrm{eV}}\right)^{-1/2}
            ~\textrm{AU},
\end{split}
\label{eq:Jeans_scale}
\end{equation}
which is about the solar system size for $m\sim 10^{-5}~\textrm{eV}$.
For extremely low mass axion, we have
\begin{equation}
   \lambda_J = 2.4 h^{-1/2} \left(\frac{m}{10^{-25}~\textrm{eV}}\right)^{-1/2}
            ~\textrm{Mpc},
\label{eq:Jeans_scale2}
\end{equation}
which is a cosmologically significant scale.
On scales larger than the axion Jeans scale the axion behaves as the CDM.
In the following we are interested in the observable consequences of such
a small mass axion as a dark matter candidate.

\subsection{Initial conditions}

In order to have initial conditions for perturbation variables in the early
radiation dominated era,
we consider the baryon, radiation, and axion in a flat background.
In the radiation dominated era ($a\propto t^{1/2}$),
Eqs.\ (\ref{eq:kappa_dot2})-(\ref{eq:v_dot2}) give
\begin{eqnarray}
 && \ddot\delta_a+\frac{1}{t}\dot\delta_a
     =\frac{3}{4t^2}\delta_r+\left(\frac{3}{2t^2}+\frac{k^2}{a^2}\right)\alpha,
     \nonumber \\
 && \ddot\delta_r+\frac{1}{2t}\dot\delta_r
     +\left(-\frac{1}{t^2}+\frac{k^2}{3a^2}\right)\delta_r
     =-\frac{2}{3t}\dot\delta_a-\frac{2}{t}\dot\alpha+\frac{3}{t^2}\alpha,
     \nonumber \\
 && \ddot\delta_b+\frac{1}{t}\dot\delta_b
     =\frac{3}{4t^2}\delta_r-\frac{3}{2t}\dot\alpha+\frac{3}{2t^2}\alpha.
\end{eqnarray}
We consider asymptotic solutions in the large-scale limit, thus
$k/(aH) \ll 1$.

(I) For $k^2/(m^2 a^2)
\ll 1$, we have $\alpha=0$; this is
the same as CDM. For $\delta_a \propto t^n$ we have
$(n-1)(2n-1)n(n+1)=0$, thus $n=1$, $\frac{1}{2}$, $0$, $-1$, and
the growing solution is
\begin{equation}
   \delta_a = \delta_b = \frac{3}{4} \delta_r \propto t.
\label{eq:ic1}
\end{equation}
This is the same as the adiabatic solution in the case of CDM.

(II) For $k^2/(m^2 a^2)
\gg 1$, we have $\alpha=-\delta$.
For $\delta_a \propto t^n$, we have $(n-1)(n-\frac{1}{2})(n^2+n+\frac{3}{2})=0$,
and the growing mode is
\begin{equation}
   \delta_a = \frac{2}{5}\delta_b=\frac{3}{10}\delta_r \propto t.
\label{eq:ic2}
\end{equation}
The condition $k/(ma) \gg 1$ together with $k/(aH) \ll 1$ demands
$H/m \gg 1$, thus {\it violating} our basic assumption of the axion as
a coherently oscillating scalar field. Although the
numerical integration of the set of equations is possible, the equations
are not valid in this case.
In the small axion mass we are interested in, such a violation is inevitable
at some epoch in the early universe (Sec.\ \ref{sec:violation}).
To avoid this difficulty, as a simple solution, though {\it ad hoc},
we can introduce an evolving axion mass model
(Sec.\ \ref{sec:evolving}).

\section{Observational effects}
\label{sec:results}

In this section, we present our numerical calculations of background
and perturbation evolutions in the axion dark matter model with
extremely low mass.
As a fiducial model, we use the flat $\Lambda\textrm{CDM}$ model consistent
with the Wilkinson Microwave Anisotropy Probe
7-year observation
($\Omega_{b0}h^2=0.02260$, $\Omega_{c0}h^2=0.1123$, $\Omega_\Lambda=0.728$,
$n_s=0.963$, $\tau=0.087$) including massless neutrinos
with $N_\nu=3.04$ (see Table 14 of Ref.\ \cite{Komatsu-etal-2011}).
As the dark matter component we consider the axion
with $\Omega_{a0}=\Omega_{c0}$.
In order to calculate the matter and CMB power spectra, we solve
a system composed of radiation, baryon, and axion in the axion-comoving gauge.
The radiation components are handled by the Boltzmann equations
while the photon-baryon fluid by the tight coupling approximation.
The detailed numerical methods are presented in Ref.\ \cite{Hwang-2001}.

\subsection{Violation of axion-fluid nature in the early era}
\label{sec:violation}

\begin{figure
\mbox{\epsfig{file=jfig1.ps,width=85mm,clip=}}
\caption{Top panels: evolution of $H/m$ (left) and  $k/(ma)$
         (right panel) as a function of scale factor $a(t)$
         for different choices of axion mass $m$ and comoving wave number $k$.
         The left panel shows that for $m = 10^{-24}~\textrm{eV}$,
         our basic assumption of $H/m \ll 1$ is {\it violated} at
         $a \lesssim 10^{-5.5}$. On the right panel showing the
         evolution of $k/(ma)$ for axion mass $m=10^{-24}~\textrm{eV}$,
         four different comoving scales from $10^{-4}$ to
         $0.1~h\textrm{Mpc}^{-1}$ have been chosen.
         Bottom panels: typical evolution of density perturbations $\delta_i$
         ($i=a$, $b$, $c$, $\gamma$) in the axion-comoving gauge.
         We considered the $\Lambda\textrm{CDM}$ model (solid curves)
         and the $\Lambda$-axion-DM model with different choices of axion mass
         (dashed curves).
         Here we consider two cases of axion mass, $m=10^{-24}~\textrm{eV}$
         (left) and $10^{-22}~\textrm{eV}$ (right panel).
         In both cases, compared with the cold dark matter in the
         $\Lambda\textrm{CDM}$ model ($\delta_c$; blue solid curve),
         the axion density perturbation ($\delta_a$; blue dashed curve)
         starts to evolve with lower amplitude; see Eq.\ (\ref{eq:ic2}).
         The behavior of $\delta_a$ slightly changes at a certain transition
         epoch
         which corresponds to when $k/(ma)$ becomes unity.
         }
\label{fig:Hm_kma_and_delta_imassv0}
\end{figure}

Due to the time dependence of $H$, our basic assumption of $H/m \ll 1$
for axion as a coherently oscillating scalar field is {\it violated}
in the early era especially for a low-mass axion.
Top panels of Fig.\ \ref{fig:Hm_kma_and_delta_imassv0}
shows evolutions of $H/m$ and $k/(ma)$ for different values of axion mass
and comoving scales.
Despite such a violation of the basic axion-fluid assumption,
here we present the baryon matter density and CMB anisotropy
power spectra for axion with a low mass.  In section \ref{sec:evolving}
we will introduce models which avoid such a violation,
and will show that the results are similar.

Examples of evolution of individual density perturbation variables are
shown in the bottom panels of Fig.\ \ref{fig:Hm_kma_and_delta_imassv0}
(see the caption for detailed descriptions), where the initial conditions for
$k/(ma) \gg 1$ limit have been used for the axion mass and the comoving scales
considered here.
The violation of $H/m \ll 1$ together with the large-scale assumption leads to
the transition between $k/(ma) \ll 1$ and $\gg 1$ at the initial epoch
(see Fig.\ \ref{fig:Hm_kma_and_delta_imassv0}, top-right panel).
In our calculation we have imposed initial conditions for $k/(ma) \gg 1$
presented in Eq.\ (\ref{eq:ic2}) because in most
cases of comoving scale and axion mass considered the values of $k/(ma)$
are larger than unity.
However, we have checked that the differences in power spectra
obtained with different sets of initial conditions in Eqs.\ (\ref{eq:ic1})
and (\ref{eq:ic2}) are negligibly small.
We also note that if we apply the initial conditions depending on the value of
$k/(ma)$ it induces a significant discontinuity in the final power spectra.

\begin{figure
\mbox{\epsfig{file=jfig2.ps,width=85mm,clip=}}
\caption{Baryonic matter (top-left) and CMB temperature, $E$-mode polarization,
         $TE$-cross power spectra (other panels)
         for axion mass $3\times 10^{-26}~\textrm{eV} \le m
         \le 10^{-24}~\textrm{eV}$ (colored curves).
         On the top-panels, the ratio of the power spectrum amplitude
         relative to the fiducial $\Lambda\textrm{CDM}$ prediction
         (black curves) are shown to indicate the differences.
         As the mass decreases below $10^{-24}~\textrm{eV}$
         we have strong small-scale suppression of the matter power spectrum
         due to the increase of axion Jeans scale.
         The CMB temperature power spectrum show small changes in higher
         multipoles corresponding to the scales below the axion Jeans scale.
         Compared with the baryon power spectrum, the change in the CMB power
         spectrum is minor.
         No noticeable changes appear in the polarization power spectra.
         }
\label{fig:matter_CMB_PS}
\end{figure}

The matter and CMB power spectra expected in the low-mass axion model
with a mass range $3\times 10^{-26}~\textrm{eV} \le m \le 10^{-24}~\textrm{eV}$
are shown in Fig.\ \ref{fig:matter_CMB_PS}.
As expected the small axion mass induces a significant damping
in the baryon matter density power spectrum.
As the axion mass decreases, the axion Jeans scale increases.
From Eq.\ (\ref{eq:Jeans_scale2}) we can see that for
$m \sim 10^{-25}~\textrm{eV}$ we have the current Jeans scale
$\lambda_J \sim 2.4~\textrm{Mpc}$.
However, the damping of the matter power spectrum appears
at the comoving scale larger than the current Jeans scale.
Note that the critical scale where the damping occurs
is determined by the Jeans scale at the radiation-matter equality.
Since the comoving Jeans wave number scales with the scale factor as
$k_{J}\propto a^{1/4}$ during the matter era, we expect
$k_{J\textrm{eq}} \approx 0.15~h\textrm{Mpc}^{-1}$
($\lambda_{J\textrm{eq}} \approx 40~h^{-1}\textrm{Mpc}$) for
$\Omega_{a0}=0.27$ and $a_{\textrm{eq}} \simeq 3\times 10^{-4}$
(see Fig.\ \ref{fig:matter_CMB_PS}), see also \cite{Hu-2000}.
Similarly, the CMB power spectra shows the deviation at
scales smaller than the axion Jeans scale, thus at very high multipoles;
for $m=3\times 10^{-26}~\textrm{eV}$, the damping of matter
power spectrum occurs around at $k_d \gtrsim 0.1~h\textrm{Mpc}^{-1}$,
which translates into the angular scales
$\ell=r_{\textrm{dec}} k_d \gtrsim 1000$
(with $r_{\textrm{dec}} \approx 10^{4}~h^{-1}\textrm{Mpc}$ the distance
to the decoupling surface); the deviation increases up to
15\% at $\ell \sim 3000$ compared to the $\Lambda\textrm{CDM}$ model
(see Refs.\ \cite{Rodriguez-Montoya-2010,Marsh-2011}
for predictions of the CMB anisotropy power spectra based on the BEC and
the ultra-light scalar field dark matter models).

\subsection{Evolving axion mass model}
\label{sec:evolving}

\begin{figure
\mbox{\epsfig{file=jfig3.ps,width=85mm,clip=}}
\caption{
         Evolution of $H/m$ (left) and $k/(ma)$ (right panels)
         as a function of scale factor $a(t)$ for evolving axion mass model
         [Eq.\ (\ref{eq:evolving_mass})].
         Here $a_\textrm{tr}=7\times 10^{-5}$ has been adopted and
         the axion mass ranges over $3\times 10^{-26}~\textrm{eV} \le m_0 \le
         10^{-5}~\textrm{eV}$.
         On the right panel the variation of
         $k/(ma)$ is shown for a fixed comoving scale
         $k=0.1~h\textrm{Mpc}^{-1}$
         and for axion mass $3\times 10^{-26} \le m_0 \le 10^{-20}~\textrm{eV}$.
         For comparison, the case of constant axion mass with
         $m=3\times 10^{-26}~\textrm{eV}$ is indicated by red dashed curves.
         }
\label{fig:Hm_kma_imassv3}
\end{figure}

We have observed that the condition for axion as a rapidly oscillating scalar
field is violated in the early universe.
As a simple solution to avoid this problem, we can introduce the evolving
axion mass in the early universe.
In this way we can investigate the role of constant-low-mass axion in the
later evolution of cosmic structures within the axion-fluid formalism.
One simple mass-evolution model is to assume that the axion mass varies
in proportional to the Hubble
parameter before a transition epoch $a_\textrm{tr}$ and becomes
constant after the transition.
In the early epoch of radiation dominated era, the Hubble parameter
varies with $a^{-2}$. Thus, assuming the transition occurs in the radiation era,
the evolving axion mass model can be written as
\begin{equation}
m(a)= \left\{
   \begin{array}{cl}
         m_0 (a_\textrm{tr}/a)^2 & a<a_\textrm{tr},  \\
         m_0                     & a \ge a_\textrm{tr},
   \end{array} \right.
\label{eq:evolving_mass}
\end{equation}
where $m_0$ is the constant axion mass after the transition $a_\textrm{tr}$.
Alternatively, we may introduce a evolving axion mass model with smooth
variation with e.g., $m(a)=m_0 [1+(a_\textrm{tr}/a)^2]$.
Although not shown here, the results are quite similar to those from the simple
evolving mass model adopted here.
For the results shown in this section, we adopt
$a_\textrm{tr}=7\times 10^{-5}$, which has been set from the condition
that $H/m \approx 0.1$ at $a< a_\textrm{tr}$ for axion mass
$m_0=3\times 10^{-26}~\textrm{eV}$, the lowest axion mass considered
in this work.
Figure \ref{fig:Hm_kma_imassv3} shows evolution of $H/m$ and $k/(ma)$
for different axion mass that varies with this functional form.

\subsection{Effects of small-scale damping on power spectra and growth rate}

Here we present the baryonic matter and the CMB anisotropy
power spectra, evolution of perturbation variables, and the growth factor
based on the simple evolving axion mass model in Eq.\ (\ref{eq:evolving_mass}).
The results for the power
spectra are shown in Fig.\ \ref{fig:ps_mvar_imassv3_v2} (solid curves),
where the power spectra for the constant axion mass (dotted curves;
the same as in Fig.\ \ref{fig:matter_CMB_PS}) are shown for comparison.

\begin{figure*
\mbox{\epsfig{file=jfig4.ps,width=150mm,clip=}}
\caption{Baryonic matter (left) and CMB temperature anisotropy (right) power
         spectra for evolving axion mass model
         with $a_\textrm{tr}=7\times 10^{-5}$ for different axion masses
         (solid curves). Dotted curves represent the similar power spectra
         obtained when the constant axion mass is assumed which are
         the same as those in top panels of Fig.\ \ref{fig:matter_CMB_PS}.
         The bottom panels show power ratios relative to
         the fiducial $\Lambda\textrm{CDM}$ model prediction.
         In the CMB power spectrum the deviations at high
         multipoles ($\ell \gtrsim 1000$) show
         opposite trends depending on the constant mass and the evolving mass models.
         }
\label{fig:ps_mvar_imassv3_v2}
\end{figure*}

Compared to the cases of constant axion mass, the power spectra for
varying axion mass shows a small difference: the damping of the power
spectrum is slightly less significant as expected.
The damping in the matter power spectrum depends on the transition
epoch $a_\textrm{tr}$; more recent transition gives smaller damping.
Compared with the $\Lambda\textrm{CDM}$ fiducial model, the deviation
in the CMB anisotropy power spectrum at high multipole $\ell$ is now
to the negative direction. Note that in the case of constant axion mass
the deviation appears in the positive direction.
In both cases, the magnitude of the maximum deviations
are very similar to each other.

In order to determine the low-mass limit of axion we present the baryon
density power spectrum in Fig.\ \ref{fig:ps_highk}, which shows the substantial
small-scale cut-off in the power spectrum for $m < 10^{-22}~\textrm{eV}$.
For $m \ge 10^{-22}~\textrm{eV}$ the matter power spectrum is almost identical
to the CDM case on the scales considered
($\lambda \lesssim 3~h^{-1}\textrm{Mpc}$).

\begin{figure
\mbox{\epsfig{file=jfig5.ps,width=75mm,clip=}}
\caption{Baryonic matter density power spectra expected in the evolving
         axion mass model with
         $10^{-25}~\textrm{eV} \le m_0 \le 10^{-22}~\textrm{eV}$
         (solid curves).
         Now the range of comoving wave number has been extended to the
         higher values (smaller scales) around $k=2~h\textrm{Mpc}^{-1}$.
         Similarly to Fig.\ \ref{fig:ps_mvar_imassv3_v2}, the power spectra
         for the constant axion mass are shown as dotted curves.
         The bottom panel shows the ratio of power relative to the
         $\Lambda\textrm{CDM}$ prediction, indicating the power spectrum
         damping on small scales.
         }
\label{fig:ps_highk}
\end{figure}

\begin{figure
\mbox{\epsfig{file=jfig6.ps,width=85mm,clip=}}
\caption{Evolution of density perturbation variables for photon, baryon,
         dark matter for evolving axion mass model, Eq.\
         (\ref{eq:evolving_mass}) with $a_\textrm{tr}=7\times 10^{-5}$.
         Shown are typical evolutions of density perturbations $\delta_i$
         ($i=a$, $b$, $c$, $\gamma$) in the axion-comoving gauge
         for $m_0=2\times 10^{-25}$, $3\times 10^{-26}~\textrm{eV}$
         and $k=0.2$, $0.5~h\textrm{Mpc}^{-1}$ (dashed curves).
         We also considered the $\Lambda\textrm{CDM}$ model whose density
         perturbations are shown as solid curves in each panel.
         Compare the results with those in Fig.\
         \ref{fig:Hm_kma_and_delta_imassv0} (bottom panels)
         }
\label{fig:delta_i_imassv3_v2}
\end{figure}
Figure \ref{fig:delta_i_imassv3_v2} shows the evolutions of density
perturbation variables for photon, baryon, and dark matter in the evolving
axion mass model at small comoving scales ($k=0.2$ and $0.5~h\textrm{Mpc}^{-1}$)
and for two extremely small axion mass ($m_0=3\times 10^{-26}$ and
$2\times 10^{-25}~\textrm{eV}$).
In all cases, the initial conditions for $k/(ma) \ll 1$ limit have been used.

\begin{figure
\mbox{\epsfig{file=jfig7.ps,width=85mm,clip=}}
\caption{Evolution of the normalized perturbation growth factor
         $g\equiv (\delta_b /a)$ at comoving scales of $k=0.1$,
         $0.5~h\textrm{Mpc}^{-1}$ for different axion masses
         (colored curves). Here we use an evolving-axion-mass model
         with $a_\textrm{tr}=7\times 10^{-5}$.
         Black curves represent the growth factor of the fiducial
         $\Lambda\textrm{CDM}$ model.
         }
\label{fig:ggg_imassv3}
\end{figure}

In Fig.\ \ref{fig:ggg_imassv3} we show the growth factor
$g \equiv \delta_b /a$ (normalized to unity at present)
at comoving scales $k=0.1$, $0.5~h\textrm{Mpc}^{-1}$ and
for various axion mass values. At $k=0.5~h\textrm{Mpc}^{-1}$,
noticeable positive deviations from the $\Lambda\textrm{CDM}$ prediction
are seen.
Despite the significant suppression of overall perturbation growth,
the recent perturbation growth factor is larger in the low-mass axion
model, especially at the small comoving scales.
This can be interpreted as the more growth in the axion mass model
(see Refs.\
\cite{Suarez-2011,Kain-2012,Chavanis-2012,Magana-2012-review}).

\subsection{Mixed cold and hot dark matter model}
\label{sec:HDM}

\begin{figure
\mbox{\epsfig{file=jfig8.ps,width=85mm,clip=}}
\caption{Matter and CMB temperature anisotropy power spectra
         expected in a model containing small fraction of massive neutrino
         component, where we have used the CAMB software \cite{CAMB}.
         Here we consider massive neutrinos ($N_\nu=3.04$)
         whose mass ranges from
         $m_\nu=0.154~\textrm{eV}$ ($\Omega_{\nu 0}=0.01$; red)
         to $0.769~\textrm{eV}$ ($\Omega_{\nu 0}=0.05$; violet curves),
         with a relation $(\Omega_{\nu 0}+\Omega_{c0})h^2=0.1123$.
         Black curves represent the power spectrum of the fiducial
         $\Lambda\textrm{CDM}$ model with massless neutrinos.
        The curves in the bottom panels indicate the ratios of powers
         relative to the $\Lambda\textrm{CDM}$ model prediction.
         }
\label{fig:pow_massive_neutrino}
\end{figure}

The warm dark matter model or the mixture of the cold and hot dark matter model
are the alternative candidates to the
CDM model with small-scale suppression in the matter power spectra.
It would be interesting to compare the small-scale suppression of our
low-mass axion model with those models.
In Fig.\ \ref{fig:pow_massive_neutrino} we present the case of massive
neutrinos contributing as the hot dark matter added to the dominant CDM model.
The results are similar to those found in Ref. \cite{Lesgourgues-2006}.

Although the low-mass axion dark matter model shows a sharp damping
in the baryonic matter power spectrum, the CDM model with a small fraction
of massive neutrinos gives a damping of power with almost constant factor at
$k > 0.1~h\textrm{Mpc}^{-1}$. Besides, the CMB anisotropy power spectrum
is very sensitive to the massive neutrino contribution; the massive neutrinos
affect all the acoustic oscillatory features at intermediate and
high angular scales ($\ell \gtrsim 100$). On the other hand, the behavior of
the axion dark matter model is quite different, and is almost insensitive
to the axion mass except at the higher multipoles.
In the low-mass axion case the changes from the CDM in both power spectra
occur only at scales smaller than the axion Jeans scale.

\section{Discussion}
\label{sec:discussion}

In this work, we have studied effects of extremely low-mass
($m \le 10^{-24}~\textrm{eV}$) axion on the baryon matter density
and the CMB anisotropy power spectra, and the perturbation
growth based on the full relativistic linear perturbation analysis.
With a low mass, however, the basic assumption about the coherently
oscillating scalar field is inevitably broken in the early universe
($H/m \gg 1$).
We have introduced the simple evolving-axion-mass-in-the-early-universe model
to avoid this problem.

We showed that axion mass smaller than $10^{-24}~\textrm{eV}$ induces the
characteristic significant damping in the baryon density power spectrum
on scales smaller than the axion Jeans scale, and changes in the higher
multipole in the CMB anisotropy power spectra.
Except for small changes in the higher multipoles ($\ell \gtrsim 1000$)
corresponding to the scales smaller than the axion Jeans scale, the CMB
power spectrum remains the same as the CDM case.
The CDM nature is also preserved in the baryon matter power spectrum above
the axion Jeans scale. We showed that the small-scale damping nature of
our low-mass axion model differs from the one expected in the CDM model
mixed with the massive neutrinos as a hot dark matter component.

Whether the small-scale damping of the baryon matter density
power spectrum can help alleviating the excess small-scale clustering problem
of the CDM model requires further studies in the nonlinear clustering
properties of the light mass axion.

\noindent{\bf Acknowledgments:}
J.H.\ was supported by KRF Grant funded by the Korean Government
(KRF-2008-341-C00022). H.N.\ was supported by grant No.\
2010-0000302 from KOSEF funded by the Korean Government (MEST).

\def\and{{and }}


\section{Introduction}

\subsection{High-temperature expansions}

Calculating exact finite-temperature properties of quantum lattice models can
be very challenging. One general approach to achieve that goal is to devise
series expansions for the lattice model in question in the
thermodynamic limit~\cite{HTE,linked}. A common example of these series expansions
are high-temperature expansions (HTEs), in which the partition function
$\mathcal{Z}$ and other extensive properties of the system are expanded
in powers of the inverse temperature $\beta=(k_BT)^{-1}$ (in what follows we set the
Boltzman constant $k_B$ to unity). For example, consider the
Ising Hamiltonian with nearest-neighbor interactions:
\begin{equation}
\label{ising}
\hat{H}=-J\sum_{\left<i,j\right>} \sigma_i \sigma_j,
\end{equation}
where $J$ is the strength of the interaction, $\left<..\right>$ denotes
nearest neighbors, and $\sigma_i~(=\pm1)$ is the Ising
spin on site $i$. The partition function can be written as
\begin{equation}
 \mathcal{Z}=\sum_{\{\sigma\}} e^{-\beta
\hat{H}}=\sum_{\{\sigma\}} e^{\beta J \sum_{\left<i,j\right>}
\sigma_i \sigma_j},
\label{eq:z}
\end{equation}
where the sum runs over all possible configurations of the spins.
We define $K=\beta J$, which serves as a small parameter in the expansion:
\begin{equation}
\label{z2}
 \mathcal{Z}=\sum_{\{\sigma\}} \prod_{\left<i,j\right>} e^{K
\sigma_i \sigma_j}=\sum_{\{\sigma\}} \prod_{\left<i,j\right>} \sum_{l=0}^{\infty}
\frac{K^l}{l !} (\sigma_i \sigma_j)^l
\end{equation}
If we associate $(\sigma_i \sigma_j)^l$ to an $l$-fold bond between sites
$i$ and $j$, a typical term in the above expansion can be represented
graphically as the lattice with various number of lines (including no
line) connecting every two nearest-neighbor sites. Therefore,
one can write the partition function in terms of contributions from all
possible graphs, up to a certain size, that can be embedded on the
lattice. (In the Ising model, the fact that $\sigma=\pm1$ makes
the calculations much easier than in quantum models such as the Heisenberg
model.) The order in $K$ that each graph contributes to is equal to
the number of bonds it has (see Ref.~\cite{linked}
for more details on this derivation). Equation~\eqref{z2} implies that the
series converges only at high temperatures, above or of the order of $J$.
In what follows, we will see how this type of expansion is
related to linked-cluster expansions.

\subsection{Low-temperature expansions}

Similar to HTEs, low-temperature expansions (LTEs) can be devised to describe
properties of a system with discrete excited states close to its ground
state, i.e., for $\beta\rightarrow\infty$. In this approach, the partition function is written as
\begin{equation}
\mathcal{Z}=e^{-\beta E_{0}} \left[ 1+\sum_{n \neq 0}
e^{-\beta(E_{n}-E_{0})}\right],
\end{equation}
where $E_0$ ($E_n$) is the ground state ($n^{\text{th}}$ excited state).
If there is an energy increment, $\epsilon$, satisfying $E_{n}-E_{0}=m_n\epsilon$
with $m_n$ being integer for all $n$, any thermodynamic property
can be expressed as an expansion in powers of $e^{-\beta \epsilon}$. Then,
cluster expansions similar to the ones discussed above for HTEs follow~\cite{linked}.

\subsection{Linked-cluster expansions}
\label{sec:lce}

The idea behind linked-cluster expansions (LCEs)~\cite{linked,sykes} is that for any extensive property $P$
of a lattice model (such as the logarithm of the partition function or the internal energy)
one can compute its value per lattice site $P(\mathcal{L})/N$ in the thermodynamic limit
in terms of a sum over contributions from all clusters $c$ that can be embedded
on the lattice:
\begin{equation}
\label{NLCE1}
P(\mathcal{L})/N=\sum_{c}L(c)\times W_{P}(c),
\end{equation}
where $L(c)$ is the multiplicity of $c$, namely, the number of ways per site in which
cluster $c$ can be embedded on the lattice, and $W_{P}(c)$ is the
weight of that cluster for the property $P$. $W_{P}(c)$ is defined according to the
inclusion-exclusion principle:
\begin{equation}
\label{NLCE2}
 W_{P}(c)=P(c)-\sum_{s \subset c} W_{P}(s),
\end{equation}
where
\begin{equation}
\label{eq:p}
P(c)=\frac{\textrm{Tr} \left[\hat{P}(c)e^{-\beta \hat{H_c}}\right]} {\textrm{Tr} \left[e^{-\beta \hat{H_c}}\right]}
\end{equation}
is the property calculated for the finite cluster $c$ and
the sum on $s$ runs over all sub-clusters of $c$. In Eq.~\eqref{eq:p}, $\hat{H_c}$
is the Hamiltonian of cluster $c$. One can check that the
weight of a disconnected cluster vanishes because $P(c)$ can be written as
the sum of its parts (see Sec.~\ref{subsec:connect}), hence,
the name linked-cluster expansions.

The convergence of the series in Eq.~\eqref{NLCE1}, when all the terms are
considered, is guaranteed by the definition of weights in Eq.~\eqref{NLCE2}.
In fact, by swapping the place of $P(c)$ and $W_{P}(c)$ in Eq.~\eqref{NLCE2},
one can write the property of a {\em finite} cluster $c$ as the sum of
its weight and the weights of its sub-clusters. Taking the thermodynamic
limit $c \to \mathcal{L}$ brings one back to Eq.~\eqref{NLCE1}. However, in that
limit, only a finite number of terms can be accounted for, and the series
has to be truncated.

Because of the inclusion-exclusion principle [Eq.~\eqref{NLCE2}],
the weight of every cluster contains only the contribution to the property
that results from correlations that involve all the sites in the cluster,
and in a unique fashion in accord with its specific geometry. At low
temperature, when correlations grow beyond the size of the largest
clusters considered in the series, the results show a divergent behavior
(due to the missing contributions of clusters in higher orders of the expansion).
In most of the two-dimensional quantum models of interest, this occurs
near or at zero temperature, e.g., for the nearest-neighbor
antiferromagnetic (AF) Heisenberg model on a bipartite lattice.

The clusters in the sum \eqref{NLCE1} are usually grouped together based
on common characteristics to represent different
orders~\cite{M_rigol_06,m_rigol_07a}. For instance, in the {\em site expansion},
where sites are used as building blocks to generate the clusters,
the order of the expansion is determined by the number of sites of the
largest clusters. All the clusters with $n$ sites are said to
belong to the $n^{\text{th}}$ order. In LCEs, one has the freedom to
devise an expansion (with a certain building block for generating the
clusters in different orders) that best suites the particular model of
interest. Some of these include the site, bond, or square expansions on
the square lattice, and site, bond, or triangle expansions on the triangular and Kagome
lattices, and so on~\cite{m_rigol_07a}.

Despite the simple form of the LCE equations above, its computational
implementation can be challenging, as one has to (i) generate all the
linked clusters that can be embedded on the lattice, (ii) identify their
symmetries and topologies to compute multiplicities [this step dramatically
reduces the computational effort as, for any given model, many
different clusters have identical values of $P(c)$], (iii) identify the
sub-clusters of each cluster to calculate weights, and (iv)
calculate the property of each cluster and perform the sums. LCEs
are computationally very demanding as the number of embedded clusters,
and their sub-clusters, grow exponentially with increasing the order
of the expansion. Below, we explain all those steps in the context
of an example (the site expansion on a finite square lattice).
We should stress that the HTE explained above for the Ising model can be
seen as a LCE in which the property for each cluster is calculated using a
perturbative expansion of Eq.~\eqref{eq:p} in terms of $K$.

\subsection{Numerical linked-cluster expansions}
\label{sec:nlce}

In this work, we present a pedagogical overview of the {\em numerical}
linked-cluster expansions (NLCEs) introduced in Ref.~\cite{M_rigol_06}. NLCEs
use the basis of LCEs, but employ exact diagonalization (ED), instead of
perturbation theory as done in HTEs, to calculate the properties of finite
clusters in the series. This means that the properties of each cluster are calculated
to all orders in $\beta$. The main advantage of NLCEs over HTEs is that, for models
with short-range correlations, it is possible to access arbitrarily low temperatures
because of the lack of a small parameter in the series. Furthermore,
for models in which correlations grow slowly as the temperature is lowered,
NLCEs can converge well below the temperature at which HTEs diverge.
In NLCEs, as opposed to HTEs, the convergence temperature is controlled by the
correlations in the model, and by the highest order
in the series that can be calculated.

The basics of NLCEs, and results for various spin and itinerant models
in the square, triangular and kagome lattices, were presented in
Refs.~\cite{M_rigol_06,m_rigol_07a,m_rigol_07b}. Recent applications
of this method exploring properties of frustrated magnetic systems
can be found in Refs.~\cite{singh_oitmaa_1,singh_oitmaa_2,applegate}. NLCE studies of the
Hubbard model in the square and honeycomb lattices were reported in
Refs.~\cite{E_khatami_11b,khatami_rigol_12,tang_paiva_12}. How to deal with
Hamiltonians and observables that break some of the symmetries
of the underlying lattice was discussed in
Refs.~\cite{E_khatami_11,khatami_helton_12_69}. Finally, how to generalize
NLCEs to calculate ground-state as well as low-temperature properties
of lattice Hamiltonians with dimerized ground states was the subject
of Ref.~\cite{khatami_singh_11_66}, while direct ground state calculations
in the thermodynamic limit were done in Ref.~\cite{irving} and dynamics
were explored in Ref.~\cite{dynamics}.

Here, we introduce NLCEs for finite clusters, in order to show how they converge
to the exact result with increasing the order in the expansion, and, later,
discuss NLCEs in the thermodynamic limit. The exposition is organized as follows.
In Sec.~\ref{sec:NLCE}, we present the algorithmic details and the
numerical implementation of NLCE in two dimensions for a finite $4\times 4$
lattice. In Sec.~\ref{sec:res}, we report results
of the expansion for the AF Heisenberg model on this cluster and compare
them, in each order, to exact results that can be obtained by means of full
exact diagonalization of the $4\times 4$ lattice. In Sec.~\ref{sec:thermo}, we
discuss how to extend NLCEs to the thermodynamic limit, and compare NLCE
results for the Heisenberg model to those from large-scale quantum Monte
Carlo (QMC) simulations. In Sec.~\ref{sec:resum}, we describe two resummation
techniques that have been found to be widely applicable to accelerate
the convergence of the NLCEs for thermodynamic properties
of lattice models of interest.

\section{Implementation of NLCEs for Finite Systems}
\label{sec:NLCE}

In this section, we sketch the NLCE algorithm for an arbitrary Hamiltonian that
only connects nearest-neighbor sites on a $N=16$ site square lattice
with open boundary conditions, which is shown in Fig.~\ref{fig:lattice16}. We
have chosen this small lattice because it allows us to carry out the NLCE
to all orders in the site expansion. It also allows us to compare the properties
in each order of the expansion to exact results calculated directly by ED of
the entire 16-site system.

\begin{figure}[t]
\centerline {\includegraphics*[width=1.2in]{16-site-cluster1.eps}}
\caption{The $4\times 4$ finite lattice with open boundary conditions used as an
example in the derivation of the NLCE.}
\label{fig:lattice16}
\end{figure}

\subsection{Embedded clusters}
\label{subsec:gener}

\begin{table}[!b]
\caption{Total number of embedded clusters (second column), number of linked-clusters (third column),
number of linked-clusters that are not related by point-group symmetries (fourth column), and
number of linked-clusters that are topologically distinct (fifth column), in each order of the site
expansion for a 16-site lattice with nearest-neighbor interactions (Fig.~\ref{fig:lattice16}).}

\vspace{0.05in}

\begin{tabular}{|c|c|c|c|c|}
\hline
\color{blue} $n$ & \color{blue} $\frac{N!}{n!(N-n)!}$ & \color{blue}Connected &
\color{blue}Sym. distinct & \color{blue} Topo. \\
\hline
\color{blue}1 & 16 & 16 & 1 & 1 \\
\color{blue}2 & 120 & 24 & 1 & 1  \\
\color{blue}3 & 560 & 52 & 2 & 1 \\
\color{blue}4 & 1820 & 113 & 5 & 3 \\
\color{blue}5 & 4368 & 244 & 14 & 4 \\
\color{blue}6 & 8008 & 496 & 43 & 10 \\
\color{blue}7 & 11440 & 912 & 94 & 19 \\
\color{blue}8 & 12870 & 1474 & 197 & 49 \\
\color{blue}9 & 11440 & 2032 & 296 & 92  \\
\color{blue}10 & 8008 & 2286 & 327 & 167 \\
\color{blue}11 & 4368 & 2000 & 265 & 190 \\
\color{blue}12 & 1820 & 1236 & 169 & 152 \\
\color{blue}13 & 560 & 488 & 66 & 65 \\
\color{blue}14 & 120 & 116 & 20 & 20  \\
\color{blue}15 & 16  & 16 & 3 & 3 \\
\color{blue}16 & 1   & 1 & 1 & 1 \\
\hline
\end{tabular}
\label{tab:generation}
\end{table}

There exists $\binom{N}{n}=\frac{N!}{n!(N-n)!}$ clusters in the \nth order
on a lattice with $N$ sites. To identify them in a computer, one can
number the lattice sites in an arbitrary fashion. Then, construct an array
of size $N$ whose $i^{\textrm{th}}$ element is a binary number representing
site number $i$ on the lattice. A 0 as an element of this array indicates that
the corresponding site is not part of the generated cluster whereas a 1
indicates that it is part of the cluster. Therefore, all clusters in
the \nth order can be generated by exploring all possible configurations of
1 and 0 as elements of the above array, while keeping the total number of
nonzero elements fixed at $n$.

In the second column of Table~\ref{tab:generation},
we list $\binom{N}{n}$ for our example of the 16-site lattice in
Fig.~\ref{fig:lattice16}. Once we have all the site clusters, we need to
connect the sites present on them with bonds. This is done following the
Hamiltonian of interest. For models with nearest-neighbor interactions, the
case of interest here, bonds are added to every pair of nearest-neighbor sites
\cite{M_rigol_06}. Hence, each site in our clusters is connected
by bonds with up to 4 other sites. Examples of clusters within the site expansion
are given in Fig.~\ref{fig:generation}.

\subsection{Connected clusters}
\label{subsec:connect}

Many of the site clusters generated in the step above contain disconnected parts.
Those clusters should be discarded because their weight $W_{P}(c)$, as given by
Eq.~\eqref{NLCE2}, is zero. This can be easily seen if one assumes that cluster
$c$ consists of two disconnected sub-clusters $c_{1}$ and $c_{2}$. Then, since
$P(c)$ is extensive, it can be written as the sum of the properties of the two
sub-clusters:
\begin{equation}
\label{start1}
 P(c)=P(c_{1})+P(c_{2}).
\end{equation}
But, in addition to all of their subclusters, $c_1$ and $c_2$ themselves
are among the subclusters of $c$. Therefore,
\begin{eqnarray}
W_{P}(c)&=& P(c)-\sum_{s \subset c} W_{P}(s)\nonumber \\
&=&P(c)- \left[W_{P}(c_{1})+\sum_{s \subset c_{1}}
W_{P}(s)\right]\nonumber \\
&\ &\ \ \ \ \ \ \ \ -\left[W_{P}(c_{2})+\sum_{s \subset c_{2}} W_{P}(s)
\right]\nonumber\\ &=&P(c)-P(c_{1})-P(c_{2})=0.
\end{eqnarray}

\begin{figure}[!b]
\centerline {\includegraphics*[width=3.in]{sample.eps}}
\caption{A sample of linked clusters that can be embedded on our finite 16-site
lattice. Clusters 3 and 5 are related by point-group symmetry and
are topologically the same as cluster 4.}
\label{fig:generation}
\end{figure}

One then needs to find a way to distinguish between connected and disconnected
clusters. A cluster is connected if starting from any site
one can reach any other site moving along the bonds (only nearest-neighbor bonds
in our case). We have implemented this idea by reconstructing the
cluster from scratch. In our codes, starting from one of the sites in
the cluster, we add bonds between that site and all of its neighboring
sites, provided that they are part of the cluster.
The same process is repeated for each of the new sites and this continues
until there are no more options for adding bonds. The resulting cluster
is then compared to the original one. If they are not identical,
the cluster must have disconnected parts and is discarded.

In the third column in Table~\ref{tab:generation}, we show the number of
clusters in each order of the site expansion after the disconnected ones
have been removed. Figure~\ref{fig:generation}, depicts some of the linked
clusters generated for our 16-site lattice.

\subsection{Point-group symmetries}
\label{subsec:symm}

If the Hamiltonian preserves the symmetries of the underlying lattice,
clusters related by point-group symmetries have the same weight
for all observables. Therefore, they can be grouped together in order
to avoid repeating the calculation of the weights for all of them.
Depending on the lattice, one may have different number
of point-group symmetries. For the square lattice, which is the case in
our example, one has the following eight point group symmetries:
\\

\begin{tabular}{lll}
$\bullet$ Identity & \ \ & $\bullet$ Reflection about $x=0$\\
$\bullet$ Rotation by 90$^\circ$ & \ & $\bullet$ Reflection about $y=0$\\
$\bullet$ Rotation by 180$^\circ$ &\  & $\bullet$ Reflection about $x=y$ \\
$\bullet$ Rotation by 270$^\circ$ & \ & $\bullet$ Reflection about $x=-y$
\end{tabular}
\vspace{0.2in}

Hence, at this point in the implementation of the NLCEs, the goal is to
identify and keep only those clusters that are not related by point-group
symmetries. We achieve this through the following steps, which need to be
followed independently for each order of the expansion:
\begin{itemize}

\item[(i)] Create a list of symmetrically distinct clusters and their
multiplicity. This list begins with the first linked cluster in the
order under consideration with multiplicity one.

\item[(ii)] Take the next connected cluster and generate all clusters that
are symmetrically related to it by applying the symmetry operations
mentioned above.

\item[(iii)] For each cluster generated in (ii), check whether
it is present in the list created in (i). This is achieved by finding out
whether there is a translation vector that converts one cluster into the other.
If yes, increase the multiplicity of the stored cluster that is a match
by one, and go to (ii).

\item[(iv)] If none of the clusters generated in (ii) is in the list created
in (i), store the connected cluster in the list, set its multiplicity to
one, and go to (ii).

\end{itemize}

This step significantly reduces the number of clusters that one has to consider.
Among the clusters shown in Fig.~\ref{fig:generation}, clusters 3 and 5 are
related by a 90$^\circ$ rotation and only one of them needs to be stored.
The fourth column in Table \ref{tab:generation} shows the number of clusters that
are not related by point-group symmetries in each order of the site expansion.

\subsection{Topological clusters}
\label{subsec:topol}

Furthermore, even if some clusters are not related by point-group symmetries,
their Hamiltonian may be identical. For example, in models with only
nearest-neighbor terms, clusters 3 and 4 in Fig.~\ref{fig:generation}
have the same Hamiltonian. We then say that these two clusters are topologically
equivalent, and only one of them needs to be diagonalized. It is important to note
that topologically equivalent clusters share the same thermodynamic properties,
but they may lead to different results for correlation functions. For example,
the distance between the two extreme sites in clusters 3 and 4 in
Fig.~\ref{fig:generation} is $2a$ and $\sqrt{2}a$ ($a$ is the lattice spacing),
respectively, and that difference needs to be taken into account when calculating
correlation functions. Still, since the full exact diagonalization of each cluster
is the most time consuming part of the NLCE calculations, the fact that one only
needs to diagonalize topologically different clusters (or simply,
{\em topological clusters}) reduces the computation time dramatically.

At this point, we need to generate a list of all
topological clusters. It can be easily verified that they share the same
topologically distinct sub-clusters, i.e., only diagonalizing topological
clusters allows one to carry out the entire NLCE calculations.

The identification of the cluster topologies can be done through adjacency
matrices ($M$). An adjacency matrix contains all the information about
the spatial relations between sites in the cluster. In the simplest
form, an adjacency matrix of an $n$-site cluster for a model with only
nearest-neighbors interaction can be a $n\times n$ matrix whose
rows/columns represent different sites and the matrix elements,
$M_{ij}$, are either 1, if sites $i$ and $j$ are nearest neighbors
(are connected by a bond), or 0 otherwise. Such matrix will be symmetric
with zeros as diagonal elements. A generalized version of such matrices
can be used for models with bond anisotropy or correlations beyond
nearest neighbors.

If two clusters of the same size are topologically equivalent, then there
exists a permutation of sites that transforms the adjacency matrix of one
onto the other. Therefore, similar to the approach taken in the previous step,
we make a list of topologically distinct clusters for each order as follows:

\begin{itemize}

\item[(i)] Create a list of topological clusters and their
multiplicity. This list begins with the first symmetrically
distinct cluster in the order under consideration with its
multiplicity.

\item[(ii)] Take the next symmetrically distinct cluster and
construct its $n!$ possible adjacency matrices (given the
$n!$ site permutations).

\item[(iii)] For each adjacency matrix generated in (ii), check whether
it coincides with the adjacency matrix of any of the clusters in the
list created in (i). If yes, increase the multiplicity of the stored
cluster that is a match by the multiplicity of the new cluster,
and go to (ii).

\item[(iv)] If none of the adjacency matrices generated in (ii)
matches the one of a cluster in the list created in (i), store
the newly found topological cluster in the list, set its multiplicity
to the one it had as a symmetrically distinct cluster, and go to (ii).

\end{itemize}

In practice, this part of the algorithm can be very time consuming not
only because the number of permutations can be very large, but also
because of the matrix comparisons. So, alternatively, one can generate a
{\em key} (an individualized number) for each adjacency matrix and use
that key for comparison purposes. To do that, one has to pick only one,
out of the $n!$ adjacency matrices to represent a cluster. This can be
accomplished by labeling sites according to a unique
criterion for the adjacency matrix. For example, if we think of each
adjacency matrix as a large binary number by attaching its columns/rows
one after another, one can pick the permutation of site numbers
that yields the largest (or smallest) binary number. In Fig.~\ref{fig:key},
we show one of the adjacency matrices of a four-site cluster and the
same matrix after exchanging site labels 1 and 2. The latter, which
yields the smallest binary number that consists of columns of the matrix,
is the one that will represent this cluster and that will be stored.

\begin{figure*}[ht]
\parbox[l]{2in}{
\centerline {\includegraphics*[width=0.6in]{AdjMat.eps}}}
\parbox[r]{3.5in}{
\begin{equation}
M=\left( \begin{array}{cccc}
       0 & 1 & 0 & 0  \\
       1 & 0 & 1 & 1  \\
       0 & 1 & 0 & 0  \\
       0 & 1 & 0 & 0  \\
       \end{array}
\right) \stackrel{\stackrel{\rm SORT}{(1\leftrightarrow 2)}}{\vector(1,0){30}}
\left( \begin{array}{cccc}
       0 & 1 & 1 & 1  \\
       1 & 0 & 0 & 0  \\
       1 & 0 & 0 & 0  \\
       1 & 0 & 0 & 0  \\
       \end{array}
\right) \vector(1,0){10} \ \ {\rm key}= 1113 \nonumber
\end{equation}}
\caption{An example of the adjacency matrix that represents a 4-site
cluster (left). Rows and columns are representatives of the sites
on the cluster. A 1 (0) as the $(i,j)$ element of the matrix indicates
that sites $i$ and $j$ are connected (not connected) by a nearest-neighbor
bond. The matrix on the right is sorted by exchanging sites 1 and 2.
A key is associated to the sorted matrix according to an arbitrary
formula [here, we use Eq.~\eqref{eq:key}].}
\label{fig:key}
\end{figure*}

The recipe for generating the keys is arbitrary and, for large cluster
sizes, one may not be able to find a recipe that guarantees a unique key
for each cluster, i.e., multiple adjacency matrices may end up sharing
the same key. The following is an example of one such recipe for creating
the keys:
\begin{equation}
\label{eq:key}
\textrm{key}=\sum_{i,j} M_{ij} (i\times j)^{i+j}.
\end{equation}
If the key is not unique, then when the key of a new cluster matches the
one of a stored topological cluster, one has to directly compare
the adjacency matrices in order to determine whether the two clusters are
topologically the same. We emphasize that the use of keys accelerates
the process of comparing adjacency matrices tremendously and has been
implemented in all our NLCEs.

Up to this point, we have identified and stored the topological clusters
of each size, their multiplicities, and their adjacency matrices with their
respective keys. The number of topological clusters in each order of our
example is shown in the last column in Table.~\ref{tab:generation}. This table
now makes apparent the significant reduction of the number of clusters that
need to be fully diagonalized from the initial $\binom{N}{n}$. In the following
step, we explain how to enumerate the sub-clusters of each topological cluster
to be able to ultimately calculate its weight.

We should stress that the approach described here to identify topologically
distinct clusters works so long as the cluster sizes are small enough that
all permutations of the introduced adjacency matrices can be explored. For
larger cluster sizes, one needs to generate more sophisticated adjacency
matrices whose description is beyond the scope of this
introduction to NLCEs.

\subsection{Sub-clusters}
\label{subsec:subcluster}

Given the description so far for finding topological clusters,
identifying all subclusters of a topological cluster is a relatively
straightforward task. Note that finding the topological clusters of our
finite 16-site lattice can, in itself, be interpreted as finding all
subclusters of a $4\times 4$ cluster in the list of topological
clusters of a larger lattice. Therefore, following the same procedure
as in Sec.~\ref{subsec:gener} and Sec.~\ref{subsec:connect}, and replacing
the $4\times 4$ system with any of the topological clusters, we first
generate all possible connected subclusters. This time, there is no need
for checking the point-group symmetries of the subclusters
(Sec.~\ref{subsec:symm}). This is because all clusters that are related
by symmetry share the same topology and have the same adjacency matrix. So,
after generating a new connected subcluster, it is sufficient to construct
its adjacency matrix and compare it to those of clusters with the same size
that are stored in the list of topological clusters. Since
there will always be a match, one only needs to keep track of how many
subclusters of a certain topology a topological cluster has. The latter can be achieved
by considering a multiplicity matrix, $Y$, whose $(i,j)$ element gives
the number of times the $j$th topological cluster can be embedded in the $i$th
topological cluster.

The steps described from Sec.~\ref{subsec:gener} to Sec.~\ref{subsec:subcluster}
need to be carried out just once for all Hamiltonians involving only
nearest-neighbor terms in the square lattice. The table with all topological
clusters and subclusters can be stored and used for different nearest
neighbor Hamiltonians and, for any given Hamiltonian, for different
microscopic parameters.

\subsection{Weights}
\label{subsec:weights}

Now, we have all the necessary tools to account for the subcluster
subtractions in order to calculate the weights in Eq.~\eqref{NLCE2}.
The steps that follow need to be carried out every time that a new
Hamiltonian or Hamiltonian parameter is explored. One starts with
the smallest cluster in the first order. That cluster has no subclusters
and, therefore, the weight of the property $P$ in the cluster is simply
equal to the property $P(1)$. Then, to compute the weight of
the next topological cluster (a single bond in the second order of our
site expansion, see Fig.~\ref{fig:generation}), we have to subtract
the weight of the cluster in the first order from its property,
according to the multiplicities given by the matrix $Y$:
\begin{eqnarray}
W_P(1) &=& P(1) \nonumber \\
W_P(2) &=& P(2) - Y_{21}W_P(1) = P(2) - Y_{21}P(1)\nonumber\\
\vdots
\end{eqnarray}
In the above equations, indices 1 and 2 refer to the cluster number, $c$,
which is not necessarily the order number. As mentioned in Sec.~\ref{sec:lce},
in NLCEs, the property $P(c)$ is calculated using ED.
For instance, the internal energy is:
\begin{equation}
\label{eq:E}
E(c)=\langle \hat{H_c} \rangle=\frac{\sum_{n=1}^{M_c}\varepsilon_{n}(c)
\exp[-\beta \varepsilon_{n}(c)]} {\sum_{n=1}^{M_c} \exp[-\beta
\varepsilon_{n}(c)]},
\end{equation}
where $\hat{H_c}$ is the Hamiltonian for the finite cluster $c$ with
eigenvalues $\varepsilon_n(c)$, and $M_c$ is the size of the Hilbert
space on that cluster.

We define partial sums, $S$, to group together contributions from topological
clusters with a common characteristic in one order of the NLCE.
For example, in the site expansion, the most natural characteristic is the
number of sites. Therefore, $S_n$, which is the sum of the contributions
to property $P$ from all $n$-site topological clusters in the series
($c_n$) reads
\begin{equation}
S_n = \sum_{c_n} L(c_n)W_P(c_n).
\end{equation}
The property in the $m$th order of NLCE is then
\begin{equation}
\label{eq:order}
P_m(\mathcal{L})/N = \sum_{n=1}^{m} S_n.
\end{equation}

\subsection{Results for the AF Heisenberg model}
\label{sec:res}

Now that we have sketched the site expansion for a 16-site lattice and for
generic Hamiltonians that connect only nearest-neighbor sites, let us
consider a specific example of one such Hamiltonians, namely, the
AF Heisenberg model
\begin{equation}
\label{eq:ham}
\hat{H}=J\sum_{\left<i,j\right>} \hat{\bf S}_{i} \cdot \hat{\bf S}_{j},
\end{equation}
where $\hat{\bf S}_i$ is the spin operator at site $i$. For simplicity, we set
$J=1$ so that, in what follows, all energies are given in units of $J$.

In Fig.~\ref{fig:finite}, we show results for the energy per site (a), the entropy per site (b),
\begin{equation}
S=\frac{1}{N}\left( \ln Z +
\frac{ \left\langle{\hat{H}}\right\rangle}{T}
\right),
\end{equation}
and the specific heat per site (c),
\begin{equation}
C_v= \frac{1}{N}\frac{\langle {\hat{H}}^2 \rangle-\langle {\hat{H}} \rangle^2}{T^2},
\end{equation}
on the 16-site lattice. The exact results (labeled ``Exact'' in Fig.~\ref{fig:finite}) were
obtained by means of full diagonalization of the Hamiltonian, and are compared
to those obtained in different orders of the site expansion [Eq.~\eqref{eq:order}].
Two things are apparent in those plots: (i) with increasing order, the NLCE results converge
to the exact ones at lower temperatures. (ii) At any given order, quantities that can be
represented by lower order derivatives of the partition function converge to lower temperature.
Note the contrast
between the results for the energy and $C_v$. Interestingly, even in the 15th order,
when only the contribution from the 16-site cluster is missing, the energy has a
visible discrepancy with the exact result at temperatures as high as $T=0.3$.

\begin{figure}[!t]
\centerline {\includegraphics[width=0.4\textwidth]{NLC_16_E_S_CV_Final.eps}}
\caption{(a) Energy, (b) entropy, and (c) specific heat per site vs
temperature for the AF Heisenberg model on the 16-site cluster in
Fig.~\ref{fig:lattice16}. Exact results are obtained by full exact
diagonalization of the 16-site cluster with open boundary conditions and are
compared to those in different orders of the site expansion explained
throughout this work.}
\label{fig:finite}
\end{figure}

\section{Thermodynamic Limit}
\label{sec:thermo}

The generalization of the NLCE presented above to the
thermodynamic limit is straightforward. Only certain steps in the algorithm
change when dealing with the infinite-size system. In this section, we
discuss those changes and show NLCE results for the Heisenberg model
[Eq.~\eqref{eq:ham}] in the thermodynamic limit.

First, for the generation of the clusters, we start with the smallest one
(a site in the site expansion) and systematically add more building blocks
(sites in the site expansion). Hence, to generate clusters with $n+1$ sites,
we consider all symmetrically distinct clusters that have $n$ sites and add
a nearest-neighbor site to every site at the edge of the cluster. Note that
this way of building site clusters works because of  translational
invariance in the infinite lattice, and it is not appropriate for the finite
system in Fig.~\ref{fig:lattice16} where not all sites are equivalent.
This approach not only guarantees the generation of all possible clusters
that can be embedded on the infinite lattice, but also produces only connected
clusters, i.e., the step presented in Sec.~\ref{subsec:connect} (examining
the connectivity of the clusters) can be skipped. The second column in Table
\ref{tab:generationTL} show the number of connected clusters per site,
with up to 17 sites,
in the site expansion of the square lattice in the thermodynamic limit.

For finite systems without translational symmetry, all the different
embeddings of a particular cluster are automatically
generated in the approach discussed in Sec.~\ref{subsec:gener}.
Thus, the multiplicity of a symmetrically
distinct cluster is calculated by counting the times it appears
in the process of generating all clusters. For the infinite system,
because of translational symmetry, the multiplicity is simply equal
to the number of point-group symmetries that transform a cluster to
another cluster that is not related to the first one by a translation.
Hence, a process similar to the one carried out for identifying
the symmetries of the clusters on a finite lattice needs to be implemented
for the clusters of the infinite system. However, in the latter,
the multiplicity of symmetrically distinct clusters are determined immediately
after building them by just examining their point-group symmetries.
Moreover, if by adding sites to a cluster in the previous order one finds
a new cluster that is related by a symmetry operation to one of the
clusters stored in the new list, we simply dismiss it. The third column in
Table \ref{tab:generationTL} show the number of symmetrically distinct clusters per site,
with up to 17 sites, in the site expansion of the square lattice in the
thermodynamic limit.

\begin{figure}[!t]
\centerline {\includegraphics[width=0.4\textwidth]{NLC_thermo_E_S_CV_Final.eps}}
\caption{(a) Energy, (b) entropy, and (c) specific heat per site of the AF
Heisenberg model on the square lattice in the thermodynamic limit as a
function of temperature. NLCE bare sums up to 15th order are compared with QMC
results for a 256 $\times$ 256 lattice.}
\label{fig:E_Thermo_Bare}
\end{figure}

\begin{table}[!b]
\caption{Total number of linked-clusters (second column),
number of linked-clusters that are not related by point-group symmetries (third column), and
number of linked-clusters that are topologically distinct (fourth column) per site, up to the 17th order
of the site expansion of the square lattice (in the thermodynamic limit) with nearest-neighbor
interactions.}

\vspace{0.05in}

\begin{tabular}{|c|c|c|c|c|}
\hline
\color{blue} $n$ & \color{blue}Connected &
\color{blue}Sym. distinct & \color{blue} Topo. \\
\hline
\color{blue}1 & 1 & 1 & 1 \\
\color{blue}2 & 2 & 1 &  1  \\
\color{blue}3 & 6 & 2 &  1 \\
\color{blue}4 & 19 & 5 & 3 \\
\color{blue}5 & 63 & 12 & 4 \\
\color{blue}6 & 216 & 35 & 10 \\
\color{blue}7 & 760 & 108 & 19 \\
\color{blue}8 & 2725 & 369 & 51 \\
\color{blue}9 & 9910  & 1285 & 112 \\
\color{blue}10 & 36446  & 4655 & 300 \\
\color{blue}11 & 135268 & 17073 & 746 \\
\color{blue}12 & 505861 & 63600 & 2042 \\
\color{blue}13 & 1903890 & 238591 & 5450 \\
\color{blue}14 & 7204874 & 901971 & 15197  \\
\color{blue}15 & 27394666 & 3426576 & 42192 \\
\color{blue}16 & 104592937 & 13079255 & 119561 \\
\color{blue}17 & 400795844 & 50107909 & 339594 \\
\hline
\end{tabular}
\label{tab:generationTL}
\end{table}

After finding all symmetrically distinct clusters, we follow exactly the
same procedure described in Secs.~\ref{subsec:topol} and \ref{subsec:subcluster} to determine all
topological clusters and their subclusters. For the site expansion in the square lattice, the last
column in Table \ref{tab:generationTL} show the number of topological clusters per site
in the thermodynamic limit.
We should stress that all the steps described above need to be carried out
just once for all lattice Hamiltonians that only connect nearest-neighbor
sites in the square lattice.

Once the topological clusters and subclusters are determined, one can
study any nearest-neighbor model of interest. In Fig.~\ref{fig:E_Thermo_Bare},
we show results for the energy (a), the entropy (b), and the specific heat (c)
per site for the AF Heisenberg model on the square lattice, and up to the 15th
order of the site expansion. NLCE results are compared to those obtained
by means of QMC simulations, previously reported in Ref.~\cite{E_khatami_11},
for a very large periodic system with $256\times256$ sites. Those plots
exhibit qualitatively similar features to the ones in Fig.~\ref{fig:finite}
for a finite cluster. Namely, as the order increases, NLCE results
converge to lower temperature. In addition, the energy
and the entropy can be seen to converge to lower temperature than
the specific heat. We note that, only in the 15th order, there are 42,192
topological clusters that need to be diagonalized (see Table
\ref{tab:generationTL}). Because of the very large number of clusters,
one can use an embarrassingly parallel code that distributes groups of
clusters to different processors so that one diagonalizes
many of them at once in every order of the NLCEs.

\section{Resummation Algorithms}
\label{sec:resum}

As seen from the results in Fig.~\ref{fig:E_Thermo_Bare}, the bare sums of the NLCE
convergence down to a temperature that depends on the order up to which the expansion
can be carried out and, on a more fundamental level, it depends on the build up of
correlations in the system as the temperature is lowered. The longer the correlations
the larger the cluster sizes that need to be included in the series to
achieve convergence. Fortunately, even if one cannot calculate more orders
of the NLCE, because of the exponential increase of the number of clusters and
of the Hilbert space of each cluster, several numerical resummation algorithms
have been developed that accelerate the convergence of NLCEs \cite{m_rigol_07a}.

The goal of those algorithms is to estimate
$P(\mathcal{L})/N=\text{lim}_{m\rightarrow\infty}P_m$ in Eq.~\eqref{eq:order} from
a finite set $\{P_m\}$. Resummation techniques then provide results for the observables of
interest in regions where the bare sums [Eq.~\eqref{eq:order}] do not converge.
Here we will focus our discussion on two such methods:
Wynn's algorithm~\cite{Wynn} and Euler's transformation~\cite{Euler}.
They have been shown to be extremely useful in accelerating the convergence
of NLCEs for thermodynamic properties in several models of interest
\cite{M_rigol_06,m_rigol_07a,m_rigol_07b,E_khatami_11b,khatami_rigol_12,tang_paiva_12,
E_khatami_11,khatami_helton_12_69,khatami_singh_11_66}.

In Wynn's algorithm, one defines
\begin{eqnarray}
\label{eq:wynn}
 \varepsilon_{n}^{(-1)}=0,~~ \varepsilon_{n}^{(0)}=P_n, \nonumber \\
 \varepsilon_{n}^{(k)}=\varepsilon_{n+1}^{(k-2)}+\frac{1}{\Delta\varepsilon_{n}^{(k-1)}},
\end{eqnarray}
where the differentiating operator $\Delta$ is applied to subscripts and
is expressed as
\begin{equation}
  \Delta\varepsilon_{n}^{(k-1)}=\varepsilon_{n+1}^{(k-1)}-\varepsilon_{n}^{(k-1)}.
\end{equation}

\begin{figure}[!t]
\centerline {\includegraphics[width=0.4\textwidth]{Thermo_resum_Final.eps}}
\caption{(a) Energy, (b) entropy, and (c) specific heat per site of the AF
Heisenberg model on the square lattice as a function of temperature. Here we show results of
last two orders (thick lines are used for the last order and thin lines for the one to last order)
for the NLCE bare sums (up to the 15th order), Wynn's algorithm with 6 cycles of improvement,
and Euler's transformation (for the last 11 orders). Results from a large-scale QMC and exact
diagonalization of a 16-site cluster with periodic boundary conditions are also shown.}
\label{fig:Thermo_limit}
\end{figure}

It is expected that the even entries $\varepsilon_{n}^{(2l)}$ converge to
$P(\mathcal{L})/N$, while the odd ones $\varepsilon_{n}^{(2l+1)}$ usually diverge,
where $l$ is defined as the number of cycles of improvement. As an example,
assume that the highest order in the NLCE that can be considered is $m$, i.e., the
bare sum yields $P_m$ for the property in the last order. Now, assume that
we use the Wynn algorithm with $l=1$ (one cycle of improvement) instead.
In that case,
\begin{equation}
P_m\rightarrow \varepsilon_{m-2}^{(2)}=P_{m}-\frac{S_{m}}{1-\frac{S_{m-1}}{S_{m}}},
\end{equation}
where the second term in the right hand side is the first order correction to
the result of the bare sum. Yet, it can often significantly improve the convergence.
Note that, for every cycle of improvement, there will be two terms less in
the new series determined by means of Eq.~\eqref{eq:wynn} ($P_m\rightarrow \varepsilon_{m-2l}^{(2l)}$).

We have found that this algorithm is one of the best for accelerating the
convergence of NLCEs. In Fig.~\ref{fig:Thermo_limit}(a),
we show the same energy as in Fig.~\ref{fig:E_Thermo_Bare}(a) but include the results
obtained after six cycles of improvement using Wynn's algorithm. One can see that the
last two orders of the Wynn sum (red solid lines) agree with each other down
to much lower temperatures ($T\sim 0.2$) in comparison to the bare sums. In addition,
the last order (thick solid line) is in perfect agreement with the QMC results in the entire
temperature region shown. Very similar results can be seen for the entropy in
Fig.~\ref{fig:Thermo_limit}(b). Furthermore, as depicted in Fig.~\ref{fig:Thermo_limit}(c),
and unlike the bare sums, the Wynn algorithm can precisely capture the position
and the height of the peak in the specific heat of this model, and even converges
to temperatures lower than those accessible to the QMC calculations depicted
in the figure.

Note that the maximum cluster sizes considered in NLCE shown there
(15 sites) are far smaller than the size of the lattice in the QMC simulations
($256\times 265$). Hence, we see that a thorough exploration of all topologies in
NLCEs, up to those cluster sizes, completely eliminates finite-size effects. In addition,
in combination with resummation algorithms, one can compute observables up to quite
low temperatures. NLCE results can also be contrasted against those obtained from ED of a
16-site cluster with periodic boundary conditions. As seen in Fig.~\ref{fig:Thermo_limit},
the finite cluster results depart from QMC at temperatures higher than the
convergence temperature of NLCE even with the bare sum. Using slightly larger
clusters, which can still be solved by means of full exact diagonalization,
does not improve this trend either~\cite{E_khatami_11}.

We should point out that the AF Heisenberg model on the square lattice is one of the
most challenging models that one can attempt to solve with NLCEs. This is because
the antiferromagnetic correlations grow exponentially with decreasing temperature,
and quickly exceed the linear size of the largest clusters that can be treated
within ED. For this reason, frustrated magnetic systems, for which QMC techniques
run into the infamous sign problem, and where correlation lengths are small
and grow slowly with decreasing temperature, are ideal to be explored by means of
NLCEs.

Another very useful resummation technique for alternating series, i.e.,
series in which $S_n$ alternates in sign, is the Euler transformation. Within
this approach, $S_n$ is replaced by $u_n=(-1)^n S_n$ and $P(\mathcal{L})/N$ is
approximated by the following sum:
\begin{equation} \label{Euler}
u_0-u_1+u_2+\dots-u_{n-1}+\sum_{l=0}^{m-n}\frac{(-1)^l}{2^{l+1}}\Delta^{l}u_{n},
\end{equation}
where $\Delta$ is defined as forward differencing operator
\begin{eqnarray}
\Delta^0 u_{n} &=& u_{n}, \nonumber \\
\Delta^1 u_{n} &=& u_{n+1} - u_{n}, \nonumber \\
\Delta^2 u_{n} &=& u_{n+2} - 2u_{n+1} + u_{n}, \nonumber \\
\Delta^3 u_{n} &=& u_{n+3} - 3u_{n+2} + 3u_{n+1} - u_{n}, \\
\vdots \nonumber
\end{eqnarray}
and $n-1$ is the number of terms for which a bare sum is performed
before the Euler transformation sets in for higher order terms. The first
few terms of the approximation can be written as
\begin{equation}
P_m\rightarrow P_{n-1}+(-1)^n\left[\frac{1}{2}S_{n}+\frac{1}{4}(S_{n}+S_{n+1})+\cdots\right].
\end{equation}

It is evident from the results in Fig.~\ref{fig:Thermo_limit} after implementing
the Euler transformation for the last 11 terms ($n=5$), that the latter
algorithm also dramatically improves the convergence from that of the bare sum.

\section*{Acknowledgments}

This research was supported by the National Science
Foundation (NSF) under Grant No. OCI-0904597 and enabled
by an allocation of advanced computing resources, supported
by the NSF. The computations in Sec.~\ref{sec:thermo}
were performed on Kraken at the National Institute for
Computational Science under Account No. TG-DMR100026.