%-----------------------------------------------------------------------
% Method
%-----------------------------------------------------------------------

\section{Simulation method}
\label{sec:method}

Our study employs a path-integral Monte Carlo (PIMC) method for simulating
bose systems.  Details on the application 
of this method to bulk helium can be found
in the review by Ceperley,\cite{ceprev} and our modifications for the
simulation of adsorbed helium films may be found in a previous 
publication.\cite{pierce99a} We briefly review the procedure now
in order to explain the calculations presented in this paper.

PIMC evaluates properties of an $N$-body quantum system at 
the inverse temperature $\tau$ by sampling the partition function 
$Z$.  $Z$ is expanded as a path-integral 
by inserting $M$ intermediate configurations:
%%%%%%%%%%%%%%%%
\begin{eqnarray}
Z \ = & & \frac{1}{N!} 
	\sum_P \int...\int d^3 R_1 ... d^3 R_M d^3 R \nonumber \\
  & & \times \rho ({\bf R}_1,{\bf R}_2; \tau) \rho({\bf R}_2, {\bf R}_3; \tau) 
  	\ldots \rho ({\bf R}_M,P {\bf R}_1; \tau),
\label{eq:pathint}
\end{eqnarray}
%%%%%%%%%%%%%%%%
where $\rho$ is the density matrix, ${\bf R}_i$ is a configuration of 
$N$ particles, and $\tau=\beta/M$.  By taking $M$ large enough, the
density matrices at $\tau$ can be accurately approximated.  The particular
PIMC method that we employ ergodically samples both particle positions
and particle permutations.

The key element needed in the first layer simulation is an accurate
approximation for the high temperature density matrices.  For bulk
helium and helium on a smooth substrate, the starting approximation
for the density matrix can be taken with a temperature as low as
40 K.  In this paper, results for the first layer on a smooth
substrate were obtained using the same starting approximation
for the density matrix that was used in previous calculations
for the second layer.\cite{pierce99a}

However, as discussed in our previous publication,\cite{pierce99a}
a helium-graphite interaction that includes substrate corrugations 
makes the starting approximation used for smooth substrates impractical,
and so we must use a simpler form.
For these calculations,
we instead use the high temperature
expansion of the density matrix.  At sufficiently small values of the
inverse temperature, the density matrix is given by
the Trotter approximation:
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
\label{eq:trotter}
\exp(-\tau \hat{H}) \approx \exp(-\tau \hat{T}) \exp(-\tau \hat{V}),
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
where $\hat{H}$, $\hat{T}$ and $\hat{V}$ are the Hamiltonian and the 
kinetic and potential energy operators, respectively, and
$\tau$ is the inverse temperature.  This allows us
to approximate the density matrices at $\tau$ as
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
\label{eq:hitemp}
\rho({\bf R},{\bf R}';\tau) \propto & &\exp[-\pi ({\bf R}_i
	-{\bf R}_{i+1})/\lambda_T^2 \nonumber \\
        & & -\tau (V({\bf R}_i)+V({\bf R}_{i+1}))/2],
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
where $R$ and $R'$ are the position vectors for $N$ particle configurations, 
and $\lambda_T$ is the thermal wavelength.  This is sometimes referred
to as the semiclassical approximation, and is accurate for sufficiently
high starting temperatures.
Averaging over the potential energy terms is referred to as the
endpoint approximation.  The potential energy term $V({\bf R}_i)$ is
the sum of all helium-helium and helium-graphite interactions.
The exponent of Eq.\ \ref{eq:hitemp} is the first term in an expansion
in powers of $\tau$.\cite{imada84}

In the semiclassical calculations, we have used 200 K as the starting 
temperature, meaning that 200 inverse-temperature slices are 
required to reach 1 K.  We have verified that this temperature
is sufficiently high for the approximation to be accurate by comparing
the energy calculated at 4 K using 200 K and 320 K as the starting 
temperatures.  The energy values obtained agreed within error bars.  We have
also adopted a three-level bisection, $l=3$, for sampling the positions.
We verified that this level gave a lower energy at 4 K than
calculations using $l=2$ and $l=4$.  The acceptance rate using
$l=3$ is approximately 40\%.  We have further verified that the 
density matrix is well approximated at 200 K by Eq.\ \ref{eq:hitemp} by
including the next term in the series expansion in powers of $\tau$.  
Energy values calculated with this starting approximation agreed within error
bars with those that used Eq.\ \ref{eq:hitemp}.  

We did find that higher order terms in $\tau$ must be included in the
energy estimator.  The values we report in Sec.\ \ref{sec:sqr3} are
obtained using 
%%%%%%%%%%%%%%%%%
\begin{eqnarray}
\label{engest}
<E> = \frac{3 N}{2 \tau} + <T> + <V> + \tau^2 \hbar^2/(8 m) <(\nabla V)^2>
\end{eqnarray}
%%%%%%%%%%%%%%%%%
for the energy expectation value,
where $N$ is the number of particles, $<T>$ and $<V>$ are the kinetic and
potential energy expectation values, and $m$ is the mass of the helium
atom.  Note that in PIMC the effective action, defined as
the natural log of the density matrix, not the energy estimator, is
used to choose between configurations in the Monte Carlo procedure.

Once the machinery for the Monte Carlo method is in place, helium-helium
and helium-substrate potentials are required for input.  We use
the Aziz\cite{aziz92} potential for the helium-helium interaction.
For the helium-graphite interaction, we have adopted the 
anisotropic Lennard-Jones potential proposed by
Carlos and Cole.\cite{cole80}  This
potential can be expressed as a Fourier series in the reciprocal lattice
vectors, {\bf G}, of the graphite substrate.  In cylindrical coordinates
($\rho$,z) the expansion is
%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
\label{eq:vgr}
V({\bf r}) = V_0(z) + \sum_{{\bf G}} V_{{\bf G}}(z) \exp(i{\bf G 
		\cdot \rho}), 
\end{eqnarray}
%%%%%%%%%%%%%%%%%%%%
where $V_0(z)$ is the laterally averaged potential, and $V_{{\bf G}}(z)$
gives the corrugation strength.  The mathematical forms for these
terms are given elsewhere, \cite{cole80} and the series
converges rapidly.  For the smooth substrate, we
use only $V_0(z)$.  For calculations that included corrugations, we
kept the six lowest, equivalent values of {\bf G} in the expansion. 

The limitation of our method is thus related to the accuracy of the
potentials available to us.  It is possible, for instance, that
the substrate may substantially mediate the helium-helium 
interaction.\cite{bruch82}
This is the so-called McLachlan interaction \cite{mclachlan64}.  We have
performed calculations both with and without this term.  Another possible
concern is the helium-graphite potential, Eq.\ \ref{eq:vgr}, which
may overestimate the corrugation strength.  Lowering the corrugation
height will effect the properties of the first layer, with the favored
phase becoming liquid instead of solid at sufficiently small corrugations.
We have examined this effect by repeating some of the calculations with
the corrugation strength set to 50 \% and to 75 \% of the value obtained from
the Carlos-Cole model.  Details of these tests of the limits of our
interaction model are given below.

Particle permutations were also included in the calculations and were
observed in the film on the smooth substrate at low densities.  
However, we have found that permutations do not play a role in the 
first layer on the corrugated substrate.  We have allowed for 
permutations at intermediate densities for temperatures as low as
0.571 K with the level of the path bisectioning taken as large as
5 (32 beads updated in one move) but never observed any particle exchanges.  
This has been checked in calculations both with and without the McLachlan
interaction.

