# HLRZ-1998-56 WUB-98/32 Multicanonical Hybrid Monte Carlo: Boosting Simulations of Compact QED

## 1 Introduction

It appears exceedingly intriguing to define variants of QED by studying vacuum states other than the usual perturbative vacuum. Lattice techniques have the potential power to deal with this situation, whenever they provide us with phase transition points of second and higher orders.

It is embarrassing that lattice simulations of compact QED still have not succeeded to clarify the order of the phase transition near $`\beta =1`$, the existence of which was established in the classical paper of Guth . This is mainly due to the failure of standard updating algorithms, like metropolis, heatbath or metropolis with reflections , to move the system at sufficient rate between the observed metastable states near its phase transition. The tunneling rates decrease exponentially in $`L^3`$ and exclude the use of lattices large enough to make contact with the thermodynamic limit by finite size scaling techniques (FSS) .

In this paper, we propose to make use of the multicanonical (MUCA) algorithm  within the hybrid Monte Carlo (HMC) updating scheme  in order to boost the tunneling rates. Since both algorithms are inherently of global nature, their combination will facilitate the parallelization of MUCA which could not be achieved otherwise.

In the early days of simulations on the hyper-torus, the U(1) phase transition was claimed to be second-order , however, with increasing lattice sizes, metastabilities and double peak action distributions became manifest, strongly hinting at its first-order character . This picture is in accord with various renormalization group studies .

However, the latent heat is found to decrease with the lattice size and the critical exponent $`\nu `$ is neither 0.25 (first order) nor $`0.5`$ (trivially second-order) . These facts allow for two possible propositions:

1. the double peak structure is a finite size effect and might vanish in the thermodynamic limit, leading to the signature of a second-order phase transition;
2. the phase transition is weakly first-order, i. e., the correlation length $`\xi `$ is finite, but large in terms of the available lattice extent $`L`$; this would fake, on small lattices, the signature of a second-order transition, and a stabilized value of the latent heat would only become visible in the thermodynamic limit $`L>\xi `$ .

In search for a lattice formulation of QED with a second-order transition point, the action was generalized to include a piece in the adjoint representation with coupling $`\gamma `$ , in the expectation that the phase transition would be driven towards second-order, at sufficiently small negative values of $`\gamma `$. However, simulations on the hyper-torus, up to $`\gamma =0.4`$, revealed the reappearance of a double peak on large enough lattices . Thus, the hypothesis of a second-order phase transition at some finite negative value of $`\gamma `$  is again doubted; furthermore, renormalization group investigations indicated that the second-order phase transition is located at $`\gamma \mathrm{}`$ .

Ref.  has speculated that the mechanism behind the lattice heuristics of metastabilities is driven by monopole loops that wrap around the hyper-torus. According to this scenario, the inefficiency of local updating algorithms to create and annihilate such monopole constellations causes their slowing down, in agreement with the earlier proposition (1). Results were presented in support of this view by switching to spherical lattices with trivial homotopy group where such wrapping loops are no more topologically stabilized . But on spherical lattices equivalent to $`L=26`$ at $`\gamma =0.2`$, double peak structures have recently been reported to reappear , corroborating earlier observations with periodic boundary conditions at $`\gamma =0`$: the suppression of monopole loop penetration through the lattice surface turned out to be incapable of preventing the incriminated double peak signal to show up on large lattices, to say $`L=32`$ .

It appears that a clarification of the situation on the Wilson line is mandatory for further progress in the understanding of compact lattice QED! This challenge requires the design of more powerful updating algorithms. A promising method is based on simulated tempering , enlarging the Lagrangian by a monopole term whose coupling is treated as an additional dynamical variable. Multi-scale update schemes in principle can alleviate critical slowing down (CSD) which is associated with the increase of the correlation length $`\xi `$ (as measured in a non-mixed phase) near the critical coupling, $`\beta _c`$ . However, the exponential supercritical slowing down (SCSD) which is a consequence of the surface tension at first-order phase transitions cannot be overcome by such type of scale-adapted methods. In that instance, one expects the autocorrelation times to grow exponentially with the system size due to the occurrence of two 3-dimensional interfaces, leading to

$$\tau _{\text{SCSD}}\mathrm{exp}\left(2\sigma L^3\right).$$
(1)   
Torrie and Valleau  have shown how to generate arbitrary “non-physical” sampling distributions. Their method, termed “umbrella sampling”, has been introduced to span large regions of phase diagrams. The method is capable to improve the efficiency of stochastic sampling for situations when dynamically nearly disconnected parts of phase space occur by biassing the system to frequent the dynamically depleted, connecting regions of configuration space. They interpreted their method as sampling “a whole range of temperatures” .

In recent years the idea of “umbrella sampling” has been popularized and extensively applied under the name “multi-canonical algorithm” (MUCA) by Berg and Neuhaus  to the simulation of a variety of systems exhibiting first-order phase transitions . In this procedure, the biassing weight $`w(S)`$ of a configuration with action $`S`$ is dynamically adjusted (bootstrapped) such as to achieve a near-constant overall frequency distribution over a wide range of $`S`$ within a single simulation.

Obviously, MUCA in principle offers a powerful handle to deal with SCSD. It remains then a practical question whether one can indeed proceed to large lattices by boosting tunneling rates from the SCSD behaviour (eq. (1)) to the peak efficiency of local Monte Carlo methods (characterized by $`O((L^4)^2)`$ complexity). This leads us to the key point of this paper: it is a severe shortcoming of the multi-canonical algorithm that its implementation is seemingly restricted to sequential computers, as it requires knowledge of the global action, even during local updating. We will show that HMC is from the very outset able to implement MUCA in a parallel manner.

In sections 2.1 and 2.2, we will give a short review of the MUCA and the HMC algorithms. In section 2.3, we merge MUCA with HMC. From our ongoing simulation project of U(1) theory on the Wilson line , we determine the tunneling efficiency compared to the standard metropolis algorithm which in our case is complemented by three reflection steps. In section 3, we shall present our results for lattice sizes up to $`16^4`$ and predict the tunnelling rates for lattice sizes up to $`24^4`$, as would be required for a proper FSS.

## 2 Multicanonical Hybrid Monte Carlo

The hybrid Monte Carlo (HMC) algorithm  produces a global trial configuration by carrying out a molecular dynamics evolution of the field configuration very close to the surface of constant action. Subsequently, a Monte Carlo decision is imposed which is based on the global action difference $`\mathrm{\Delta }S`$, being small enough to be frequently accepted. Within HMC, all degrees-of-freedom can be changed simultaneously and hence in parallel. This then provides the straightforward path to implement MUCA as part of HMC on parallel machines<sup>1</sup><sup>1</sup>1A first attempt in this direction has been made in Ref.  in the framework of the Higgs-Yukawa model.: one just uses the values of the global action, as provided by HMC, to compute the bias function $`w(S)`$ for the MUCA algorithm.

### 2.1 Multicanonical Algorithm

“Canonical” Monte Carlo generates a sample of field configurations, $`\left\{\varphi \right\}`$, within a Markov process, according to the Boltzmann weight,

$$P_\beta (S)=\frac{1}{Z_{\text{can}}}e^{\beta S(\varphi )},$$
(2)   
which follows from maximizing the entropy with respect to all possible probability distributions $`P[\varphi ]`$. The partition function $`Z`$ normalizes the total probability to 1,

$$Z_\beta =[d\varphi ]e^{\beta S[\varphi ]}.$$
(3)   
$`S`$ is the action (the energy in the case of statistical mechanics) and $`\beta `$ the coupling (or inverse temperature $`\frac{1}{kT}`$).

The canonical action density which in general exhibits a double peak structure at a first-order phase transition, can be rewritten as

$$N_{\text{can}}(S,\beta )=\rho (S)e^{\beta S},$$
(4)   
with the spectral density $`\rho (S)`$ being independent of $`\beta `$. Usually, $`N_{\text{can}}(S,\beta )𝑑S`$ is set to 1.

The multicanonical approach aims at generating a flat action density

$$N_{\text{MUCA}}(S,\beta )=\text{const.,       for }S_{\text{min}}SS_{\text{max}},$$
(5)   
in a range of $`S`$ that covers the double peaks at the first-order phase transition.

Therefore, instead of sampling canonically according to $`e^{\beta S}`$, one modifies the sampling by a weight factor $`W_{\text{MUCA}}`$:

$$W_{\text{MUCA}}(S,\beta )=\{\begin{array}{ccc}\hfill \frac{1}{\rho (S_{\text{min}})}& e^{\beta S_{\text{min}}}\hfill & \text{for}S<S_{\text{min}}\hfill \\ \hfill \frac{1}{\rho (S)}& e^{\beta S}\hfill & \text{for}S_{\text{min}}SS_{\text{max}}\hfill \\ \hfill \frac{1}{\rho (S_{\text{max}})}& e^{\beta S_{\text{max}}}\hfill & \text{for}S>S_{\text{max}},\hfill \end{array}$$
(6)   
which is constant outside the relevant action range. Such ‘corrigez la fortune’ is equivalent to a net sampling according to the yet unknown probability distribution

$$w(S)=\frac{1}{\rho (S)},\text{for}S_{\text{min}}SS_{\text{max}}.$$
(7)   
Since $`W(S,\beta )`$ is unknown at the begin of the simulation, it is instrumental for MUCA to follow Münchhausen and bootstrap from good guesstimates . We shall do so by starting from an observed histogram of the canonical action density, $`\widehat{N}_{\text{can}}(S,\widehat{\beta }_c)`$, see eq. (4), at the supposed location of the phase transition<sup>2</sup><sup>2</sup>2Hatted quantities refer to stochastic estimates., $`\widehat{\beta }_c`$. From the action density, we compute $`\widehat{W}_{\text{MUCA}}(S,\widehat{\beta }_c)`$ according to eq. (6). The sampling then proceeds with the full MUCA weight,

$`\widehat{P}_{\text{MUCA}}(S)`$ $``$ $`e^{\widehat{\beta }_cS}\widehat{W}_{\text{MUCA}}(S,\widehat{\beta }_c)`$ (8)   
$``$ $`e^{\left(\widehat{\beta }_c+\widehat{\beta }(S)\right)S\widehat{\alpha }(S)}.`$   
This latter formulation can be interpreted as a simulation proceeding at varying couplings (temperatures), hence the naming ‘multicanonical’.

In order to compute expectation values of observables $`𝒪`$, one has to reweight the resulting action density at the end of the day by the factor $`\widehat{W}_{\text{MUCA}}(S,\widehat{\beta }_c)`$, which reconstitutes the proper canonical density:

$$\widehat{𝒪}_{\widehat{\beta }_c}=\frac{_i𝒪_{\widehat{\beta }_c}^i\frac{1}{\widehat{W}_{\text{MUCA}}(S_i,\widehat{\beta }_c)}}{_i\frac{1}{\widehat{W}_{\text{MUCA}}(S_i,\widehat{\beta }_c)}}.$$
(9)   
Additionally, $`\widehat{𝒪}_{\widehat{\beta }_c}`$ simulated at $`\widehat{\beta }_c`$ can be reweighted to any desired $`\beta `$ (following ), given that the corresponding region of phase space has been covered by the MUCA simulation. We emphasize that eq. (9) is only useful complemented by a proper error analysis. The canonical error computed from the multicanonical ensemble has been elaborated in Ref. .

Note that there are many possible choices for the form of the multicanonical weight. Just for technical reasons we require it to be continuous in $`S`$. One can either guess an analytic function, or choose a polygonal approximation such as given in eq. (8). In this case, the multicanonical weight is expressed in terms of the functions $`\widehat{\beta }(S)`$ and $`\widehat{\alpha }(S)`$ which are actually characteristic functions of the bins. $`\widehat{\beta }(S)`$ can be considered as an effective temperature .

The computation of the weights requires the knowledge of the global and not just the local change in action for each single update step. For this reason, even for a local action, one cannot perform local updating moves in parallel, such as the well-known checkerboard pattern. As a consequence, MUCA is not parallelizable for local update algorithms. For remedy, we propose here to go global and utilize the HMC updating procedure.

### 2.2 Hybrid Monte Carlo

The HMC consists of two parts: the hybrid molecular dynamics algorithm (HMD) evolves the degrees of freedom by means of molecular dynamics (MD) which is followed by a global Metropolis decision to render the algorithm exact.

In addition to the gauge fields $`\varphi _\mu (x)`$ one introduces a set of statistically independent canonical momenta $`\pi _\mu (x)`$, chosen at random according to a Gaussian distribution $`\mathrm{exp}(\frac{\pi ^2}{2})`$. The action $`S[\varphi ]`$ is extended to a guidance Hamiltonian

$$[\varphi ,\pi ]=\frac{1}{2}\underset{\mu ,x}{}\pi _\mu ^2(x)+\beta S[\varphi ].$$
(10)   
Starting with a configuration $`(\varphi ,\pi )`$ at MD time $`t=0`$, the system moves through phase space according to the equations of motion

$`\dot{\varphi }_\mu `$ $`=`$ $`{\displaystyle \frac{}{\pi _\mu }}=\pi _\mu ,`$   
$`\dot{\pi }_\mu `$ $`=`$ $`{\displaystyle \frac{}{\varphi _\mu }}={\displaystyle \frac{}{\varphi _\mu }}[\beta S],`$ (11)   
leading to a proposal configuration $`(\varphi ^{},\pi ^{})`$ at time $`t=\tau `$. Finally this proposal is accepted in a global Metropolis step with probability

$$P_{acc}=\mathrm{min}(1,e^\mathrm{\Delta }),\text{with}\mathrm{\Delta }=[\varphi ^{},\pi ^{}][\varphi ,\pi ].$$
(12)   
The equations of motion are integrated numerically with finite step size $`\mathrm{\Delta }t`$ along the trajectory from $`t=0`$ up to $`t=N_{\text{md}}`$. Using the leap-frog scheme as symplectic integrator the discretized version of eq. (2.2) reads:

$`\varphi ^{n+1}`$ $`=`$ $`\varphi ^n+\mathrm{\Delta }t\pi ^n{\displaystyle \frac{\mathrm{\Delta }t^2}{2}}{\displaystyle \frac{}{\varphi }}\left(\beta S[\varphi ^n]\right)`$   
$`\pi ^{n+1}`$ $`=`$ $`\pi ^n{\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \frac{}{\varphi }}\left(\beta S[\varphi ^n]\right){\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \frac{}{\varphi }}\left(\beta S[\varphi ^{n+1}]\right).`$ (13)   
Here we have presented the scheme with both the momenta and the gauge fields defined at full time steps<sup>3</sup><sup>3</sup>3Note that the actual implementation computes the momenta at half time steps according to the sequence    $`\pi (t+\mathrm{\Delta }t/2)`$ $`=`$ $`\pi (t\mathrm{\Delta }t/2)\mathrm{\Delta }t{\displaystyle \frac{}{\varphi }}\left(\beta S[\varphi (t)]\right),`$    $`\varphi (t+\mathrm{\Delta }t)`$ $`=`$ $`\varphi (t)+\mathrm{\Delta }t\pi (t+\mathrm{\Delta }t/2),`$  initialized and finished by a half-step in $`\pi `$ . Each sequence approximates the correct $``$ with an error of $`O(\mathrm{\Delta }t^3)`$., $`t=n\mathrm{\Delta }t`$.

To ensure that the Markov chain of gauge field configurations reaches a unique fixed point distribution $`\mathrm{exp}(S[\varphi ])`$ one must require the updating procedure to fulfil detailed balance, which is guaranteed by the iterative map of eq. (2.2) $`f:(\varphi ,\pi )(\varphi ^{},\pi ^{})`$ being

* time reversible: $`f(\varphi ^{},\pi ^{})=(\varphi ,\pi )`$
* and measure preserving: $`[d\varphi ^{}][d\pi ^{}]=[d\varphi ][d\pi ]`$.

It is easy to prove that these two conditions also hold for the multicanonical action. Note that the guidance Hamiltonian, eq. (10), defining the MD may differ from the acceptance Hamiltonian in eq. (12), which produces the equilibrium distribution proper. In the following, we shall exploit this freedom to develop two variant mergers of MUCA and HMC.

### 2.3 Merging MUCA and HMC for Compact QED (MHMC)

We consider a multicanonical HMC for pure 4-dimensional $`U(1)`$ gauge theory with standard Wilson action defined as

$$S[\varphi ]=\underset{x,\nu >\mu }{}\left[1\mathrm{cos}\left(\theta _{\mu \nu }(x)\right)\right],$$
(14)   
where

$$\theta _{\mu \nu }(x)=\varphi _\mu (x)+\varphi _\nu (x+\widehat{\mu })\varphi _\mu (x+\widehat{\nu })\varphi _\nu (x)$$

is the sum of link angles that contribute to one of six plaquettes interacting with the link angle $`\varphi _\mu (x)`$.

Eq. 8 suggests to consider an effective action $`\widehat{S}`$ including the “multicanonical potential” $`V_{\text{MUCA}}`$

$`\widehat{S}`$ $`=`$ $`\widehat{\beta }_cS+\widehat{V}_{\text{MUCA}}(S,\widehat{\beta }_c),`$ (15)   
with

$`\widehat{V}_{\text{MUCA}}(S,\widehat{\beta }_c)`$ $`=`$ $`\mathrm{log}\left({\displaystyle \frac{1}{\widehat{W}_{\text{MUCA}}(S,\widehat{\beta }_c)}}\right)`$ (16)   
$`=`$ $`\widehat{\beta }(S)S+\widehat{\alpha }(S).`$   
There are two natural options to proceed from here:

performs molecular dynamics using the canonical guidance Hamiltonian

$$_{md}=\frac{1}{2}\pi ^2+\widehat{\beta }_cS,$$

with standard action $`S`$. The resulting gluonic force is given by

$`F(x)`$ $`=`$ $`\dot{\pi }_\mu (x)`$ (17)   
$`=`$ $`\widehat{\beta }_c{\displaystyle \underset{\nu \mu }{}}\left[\mathrm{sin}\theta _{\mu \nu }(x\widehat{\nu })\mathrm{sin}\theta _{\mu \nu }(x)\right].`$   
makes use of the multicanonical potential as a driving term within the Hamiltonian,

$$_{md}=\frac{1}{2}\pi ^2+\widehat{S},$$

inducing an additional drift term

$`\dot{\pi }_\mu (x)`$ $`=`$ $`F(x){\displaystyle \frac{}{\varphi _\mu (x)}}\widehat{V}_{\text{MUCA}}(S\left[\varphi (x)\right],\widehat{\beta }_c)`$ (18)   
$`=`$ $`\left(\widehat{\beta }_c+\widehat{\beta }(S)\right){\displaystyle \underset{\nu \mu }{}}\left[\mathrm{sin}\theta _{\mu \nu }(x\widehat{\nu })\mathrm{sin}\theta _{\mu \nu }(x)\right],`$   
with the effective $`\widehat{\beta }`$ as defined in eq. (8).

For both options, the Hamiltonian governing the accept/reject decision, eq. (12), reads equally:

$$_{\text{acc}}=\frac{1}{2}\pi ^2+\widehat{S}.$$
(19)   
The latter method is governed by the dynamics underlying the very two peak structure: as one can see from Fig. 5, $`V_{\text{MUCA}}`$ is repelling the system out of the hot (cold) phase towards the cold (hot) phase, thus increasing its mobility and enhancing flip-flop activity.

We comment that the implementation of method 2 requires the computation of the global action (to adjust the correct multicanonical weight, eq. (8)) at each integration step along the trajectory of molecular dynamics to guarantee reversibility. In the polygon approximation, this amounts to a determination of the effective coupling $`\widehat{\beta }`$ at each time step in MD. Note that $`\widehat{\alpha }`$ does not influence MD but enters into the global Metropolis decision, eq. (12). Method 1 is much simpler, running at fixed trial coupling $`\widehat{\beta }_c`$ and avoiding the effort of computing $`N_{\text{md}}2`$ global sums while travelling along the trajectories. It turned out that of both versions of MHMC, method 2 performs better than method 1, when autocorrelation and difference in computational effort are taken into account. Thus, we continue our analysis by investigation of method 2.

## 3 Results

In order to evaluate the efficiency of the MHMC and to set the stage for a proper extrapolation, we generated time series of the U(1) action on lattices of size $`6^4`$ up to $`16^4`$. The runs are summarized in Tab. 1. We applied method 2 as being the more promising one on real large lattices.

To arrive at an action density $`N_{\text{MUCA}}(S,\beta )`$ which is approximately flat in the desired region between the two peaks it is crucial to find a good estimate $`\widehat{N}_{\text{can}}(S,\widehat{\beta }_c)`$ of the canonical action density. However, it becomes more and more delicate for large volumes to find the proper multicanonical weight, $`\widehat{P}_{\text{MUCA}}^L(S)`$ (eq. (8)). Fig. 1 shows the evolution of the action density as the lattice size increases.

For lattices $`<16^4`$ it was sufficient to perform a short canonical run to generate an action density $`\widehat{N}_{\text{MUCA}}`$ suitable to compute a proper weight factor $`\widehat{W}_{\text{MUCA}}`$. On the $`16^4`$ system, however, the resulting multicanonical distribution becomes quite sensitive to the choice of the weight. Therefore, in the case of large lattices ($`L16`$) we cannot rely on canonical simulations to start with. Even if we perform $`O(10^6)`$ sweeps using standard Metropolis update with 3 reflection steps<sup>4</sup><sup>4</sup>4The metropolis algorithm with reflection steps (MRS) is considered as a very effective local update algorithm for U(1) ., SCSD prevents a sufficiently accurate determination of the phase weight. Therefore, we install a recursive procedure: from a previous guess $`\widehat{N}_\text{i}(S,\widehat{\beta }_c^i)`$ we go through MHMC and arrive at $`\widehat{N}_{\text{i+1}}(S,\widehat{\beta }_c^{i+1})`$. This computational scheme is initialized by a standard canonical “short run”. We found that one such learning cycle is sufficient. Fig. 2 illustrates the evolution of the multicanonical action density on the $`16^4`$ lattice.

On larger volumes the determination of a good guess can be considerably boosted by a crank-up extrapolation that starts from smaller systems . In Fig. 3, we display the quality of “flatness” of $`\widehat{N}_{\text{MUCA}}(S,\beta )`$ achieved in our investigations for the various lattice sizes.

### 3.1 Tunneling Behaviour

With our estimate for $`\widehat{N}_{\text{MUCA}}^{16}(S)`$ at $`\beta =1.010753`$, depicted in Fig. 2, we have generated the time history of the action per site, $`s=S/6V`$, as shown in the upper part of Fig. 4.

For reference, we have included the time history from the MRS algorithm on the same lattice. The figure demonstrates the success of MHMC: the method provides us with a gain factor in tunneling rate of about one order of magnitude on the $`16^4`$ lattice.

In order to quantify this achievement, we introduce the average flip time, $`\tau _{\text{flip}}`$, a quantity that is readily measurable. $`\tau _{\text{flip}}`$ is defined as follows: we histogram the time series of $`s`$ using $`N`$ bins, as illustrated in Fig. 5.

A suitable number is $`N=8`$. The total binning range is adjusted such that $`99.9`$ % of the events are covered symmetrically by the 8 bins. A flip (flop) is given when the system travels from $`I_8`$ to $`I_1`$ (and vice versa). $`\tau _{\text{flip}}`$ is defined as the inverse number of the sum of flips and flops multiplied by the total number of trajectories. In Tab. 2, $`\tau _{\text{flip}}`$ is given for the various lattices. The error in $`\tau _{\text{flip}}`$ has been computed by a jackknife error analysis.

### 3.2 Scaling Behaviour

With the results for $`\tau _{\text{flip}}`$ on lattices up to $`16^4`$ we are in the position to estimate the scaling behaviour of MHMC in comparison to standard MRS updates. Fig. 6 shows $`\tau _{\text{flip}}`$ both for MHMC and MRS as a function of the lattice size $`L`$ at their respective $`\beta `$-values, as listed in Tab. 2.

According to the expected exponential behaviour of $`\tau _{\text{flip}}^{\text{MRS}}`$ which, in the asymptotic regime $`L\mathrm{}`$, should be given by eq. (1), we perform a $`\chi ^2`$-fit with the ansatz:

$$\tau _{\text{flip}}^{\text{MRS}}=aL^be^{cL^3}.$$
(20)   
It yields the following parameter values:

$`a`$ $`=`$ $`11.9(3.7)`$   
$`b`$ $`=`$ $`2.01(18)`$   
$`c`$ $`=`$ $`6.7(8)10^4,`$ (21)   
with $`\chi _{\text{per d.o.f.}}^2=0.897`$. As a result, we find a clear exponential SCSD behaviour for the MRS algorihm<sup>5</sup><sup>5</sup>5One is tempted to extract the interfacial surface tension $`\sigma `$ from the fit to the MRS data. We find    
$$\sigma =3.35(39)10^4.$$
.

On the other hand, for the tunneling times of the MHMC, we expect a monomial dependence in $`L`$:

$$\tau _{\text{flip}}^{\text{MHMC}}=pL^q,$$
(22)   
We obtain for the fit parameters:

$`p`$ $`=`$ $`11.6(1.6)`$ (23)   
$`q`$ $`=`$ $`2.238(68).`$ (24)   
The power law ansatz is well confirmed by the fit quality with $`\chi _{\text{per d.o.f.}}^2=0.795`$.

We also took the pessimistic ansatz and tried to detect a potentially exponential increase of $`\tau _{\text{flip}}^{\text{MHMC}}`$. The exponential fit gives $`\chi _{\text{per d.o.f.}}^2=0.975`$. As can be seen in Fig. 6, the exponential contribution remains suppressed in the extrapolation. A potentially dominating exponential behaviour for MHMC can only be detected in future MHMC simulations on larger lattices. In other words, parallel MHMC is capable to overcome SCSD in compact QED in practical simulations, at least up to lattices sizes $`24^4`$.

## 4 Cost Estimates for a FSS Study.

Finally, we try to assess the compute effort required to perform a FSS study on a series of lattices ranging up to $`24^4`$.

$`\tau _{\text{flip}}`$ being readily accessible, we relate this quantity to the effective integrated autocorrelation time, $`\tau _{\text{eff}}`$, defined in Ref.  by

$$\sigma _{\text{MUCA}}^2=\sigma _{\text{can}}^2\frac{2\tau _{\text{eff}}}{N_{\text{ts}}}.$$
(25)   
$`\sigma _{\text{MUCA}}^2`$ is the squared error of the observable $`𝒪`$ computed from the multicanonical ensemble, see eq. (9). $`\sigma _{\text{can}}^2`$ denotes the canonical variance of $`𝒪`$ (computed from the reweighted canonical ensemble) and $`N_{\text{ts}}`$ is the length of the multicanonical time series. We can determine $`\sigma _{\text{MUCA}}^2`$ in a numerically quite stable way from jackknife blocking. From the time series of $`S`$ on the $`14^4`$ and a $`16^4`$ MUCA system, with about 550000 and 300000 entries, respectively, we have determined $`\tau _{\text{eff}}`$ to be

$`\tau _{\text{eff}}^{14}=35(10)10^3\times \tau _{\text{flip}}^{14},`$   
$`\tau _{\text{eff}}^{16}=41(12)10^3\times \tau _{\text{flip}}^{16}.`$ (26)   
As a result, we found $`\tau _{\text{flip}}0.038(8)\times \tau _{\text{eff}}`$<sup>6</sup><sup>6</sup>6Note that $`\tau _{\text{flip}}`$ strongly depends on the difference between $`I_8`$ and $`I_1`$ in Fig. 5. It remains to be confirmed that the ratio between $`\tau _{\text{eff}}`$ and $`\tau _{\text{flip}}`$ does not vary too much going to larger lattices.. From hereon we can estimate the number of de-correlated subsamples (independent measurements) out of a time series of length $`N_{\text{ts}}`$ to be roughly

$$N_{\text{indep}}=\frac{N_{\text{ts}}}{2\tau _{\text{eff}}}=\frac{N_{\text{ts}}}{79(16)10^3\tau _{\text{flip}}}.$$
(27)   
Assuming the inverse square-root of $`N_{\text{indep}}`$ to be an upper bound to the relative error $`r`$ of an observable $`𝒪`$, we arrive at

$$K=\frac{79(16)10^3}{r^2},$$
(28)   
with $`K`$ being the required number of flips to achieve a relative error $`<r`$.

We thus conclude that $`O(100)`$ flip-flops might allow to determine quantities like the specific heat and the Binder-Landau cumulant with a relative error of 3 %.

Obviously, the costs of MHMC and MRS simulation increase with the volume of the lattice, $`V=L^4`$. Additionally, for MHMC, we want to keep the average acceptance probability of the leapfrog scheme constant. To this end, we have to lower the step size according to $`\mathrm{\Delta }tV^{1/4}`$. In a detailed tuning investigation we have confirmed that the scaling rule of constant acceptance probability  leads to optimal performance. From a spectral analysis of the molecular dynamics we can find an optimized trajectory length, $`N_{md}`$, (in the sense that the average acceptance probability is maximized) fulfilling $`\mathrm{\Delta }tN_{md}=\text{const}(\beta )`$, with only a slight dependence on $`\beta `$ near the phase transition. We choose step-size $`\mathrm{\Delta }t`$ according to

$$P_{\text{acc}}=\text{erfc}\left(c(\beta )V^{\frac{1}{2}}\mathrm{\Delta }t^2\right)=\text{const.},$$
(29)   
adjusted such that the product $`\tau _{\text{int}}N_{md}`$ is minimized finally. In our case, the optimal acceptance probability is 65 % .

Fig. 7 confirms the scaling of MHMC: the ratio of measured CPU-times for a sweep, $`\frac{t_{\text{MUCA}}}{t_{\text{MRS}}}`$ increases quite linear with $`L=V^{1/4}`$. Minor deviations from this behaviour stem from the use of suboptimal run parameters $`\mathrm{\Delta }t`$ and $`N_{\text{md}}`$.

The product $`\tau _{\text{flip}}\times t_{\text{MUCA}}`$ reflects the efficiency of MHMC. Tab. 3 lists the effective gain factors achieved taking MHMC instead of MRS. The two columns refer to the power law and exponential extrapolations for MHMC.

Let us finally translate these factors into real costs: in Fig. 8, we extrapolate the sustained CPU time in Tflop-hours required to generate one flip.

We conclude, that the integrated CPU time to generate the required 100 flips with MHMC on a $`24^4`$ lattice amounts to about 3 Tflop-hours sustained CPU time.

## 5 Summary and Outlook

We have demonstrated that the fully parallel MHMC algorithm is a very effective tool which is able to overcome SCSD as present in the pronounced metastabilities of 4-dimensional $`U(1)`$ gauge theory. A FSS study up to a lattice size of $`L=24`$ with about 100 flip events for each lattice is feasible within half a year runtime, given a sustained performance of about 2 Gflops, due to the improvements achieved by MHMC. These performance figures should be obtainable on a 32 node partition of a Cray T3E-600.

Less well known is the influence of the delicate part of MHMC, i. e. the determination of a suitable estimate for $`V_{\text{MUCA}}`$, which is carried out in an iterative manner. So far, we have encouraging experiences on the $`16^4`$ lattice. The success of the crank-up procedures described in Ref.  gives us hope that the $`P_{\text{MUCA}}`$ determination will carry through with only marginal deterioration of the improvement factors estimated here.

The investigations presented form part of an ongoing study that aims at a conclusive FSS analysis of compact QED on the Wilson line .

## Acknowledgements

The present investigation is based on the diploma thesis of G. A. carried out on the Connection machine CM5 at Wuppertal university. We thank P. Ueberholz and P. Fiebach for their friendly support. We are indebted to Thomas Riechmann and Dr. Claus-Uwe Linster of the Institut für Mathematische Maschinen at the computer center of Erlangen university, Germany, for providing us a substantial amount of computer time on their 32 node connection machine CM5. Without this support, the present investigation would not have been possible.