An efficient algorithm for numerical computations of continuous densities of states

In Wang–Landau type algorithms, Monte-Carlo updates are performed with respect to the density of states, which is iteratively refined during simulations. The partition function and thermodynamic observables are then obtained by standard integration. In this work, our recently introduced method in this class (the LLR approach) is analysed and further developed. Our approach is a histogram free method particularly suited for systems with continuous degrees of freedom giving rise to a continuum density of states, as it is commonly found in lattice gauge theories and in some statistical mechanics systems. We show that the method possesses an exponential error suppression that allows us to estimate the density of states over several orders of magnitude with nearly constant relative precision. We explain how ergodicity issues can be avoided and how expectation values of arbitrary observables can be obtained within this framework. We then demonstrate the method using compact U(1) lattice gauge theory as a show case. A thorough study of the algorithm parameter dependence of the results is performed and compared with the analytically expected behaviour. We obtain high precision values for the critical coupling for the phase transition and for the peak value of the specific heat for lattice sizes ranging from 84\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8^4$$\end{document} to 204\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20^4$$\end{document}. Our results perfectly agree with the reference values reported in the literature, which covers lattice sizes up to 184\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18^4$$\end{document}. Robust results for the 204\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20^4$$\end{document} volume are obtained for the first time. This latter investigation, which, due to strong metastabilities developed at the pseudo-critical coupling of the system, so far has been out of reach even on supercomputers with importance sampling approaches, has been performed to high accuracy with modest computational resources. This shows the potential of the method for studies of first order phase transitions. Other situations where the method is expected to be superior to importance sampling techniques are pointed out.

4 Discussion, conclusions and future plans 24 A Reference scale and volume scaling 27 1

Introduction and motivations
Monte-Carlo methods are widely used in Theoretical Physics, Statistical Mechanics and Condensed Matter (for an overview, see e.g.[1]).Since the inception of the field [2], most of the applications have relied on importance sampling, which allows us to evaluate stochastically with a controllable error multi-dimensional integrals of localised functions.These methods have immediate applications when one needs to compute thermodynamic properties, since statistical averages of (most) observables can be computed efficiently with importance sampling techniques.Similarly, in Lattice Gauge Theories, most quantities of interest can be expressed in the path integral formalism as ensemble averages over a positive-definite (and sharply peaked) measure, which, once again, provide an ideal scenario for applying importance sampling methods.However, there are noticeable cases in which Monte-Carlo importance sampling methods are either very inefficient or produce inherently wrong results for well understood reasons.Among those cases, some of the most relevant situations include systems with a sign problem (see [3] for a recent review), direct computations of free energies (comprising the study of properties of interfaces), systems with strong metastabilities (for instance, a system with a first order phase transition in the region in which the phases coexist) and systems with a rough free energy landscape.Alternatives to importance sampling techniques do exist, but generally they are less efficient in standard cases and hence their use is limited to ad-hoc situations in which more standard methods are inapplicable.Noticeable exceptions are micro-canonical methods, which have experienced a surge in interest in the past fifteen years.Most of the growing popularity of those methods is due to the work of Wang and Landau [4], which provided an efficient algorithm to access the density of states in a statistical system with a discrete spectrum.Once the density of states is known, the partition function (and from it all thermodynamic properties of the system) can be reconstructed by performing one-dimensional numerical integrals.The generalisation of the Wang-Landau algorithm to systems with a continuum spectrum is far from straightforward [5,6].To overcome this limitation, a very promising method, here referred to as the Logarithmic Linear Relaxation (LLR) algorithm, was introduced in [7].The potentialities of the method were demonstrated in subsequent studies of systems afflicted by a sign problem [8,9], in the computation of the Polyakov loop probability distribution function in two-colour QCD with heavy quarks at finite density [10] and -rather unexpectedly -even in the determination of thermodynamic properties of systems with a discrete energy spectrum [11].
The main purpose of this work is to discuss in detail some improvements of the original LLR algorithm and to formally prove that expectation values of observables computed with this method converge to the correct result, which fills a gap in the current literature.In addition, we apply the algorithm to the study of Compact U(1) Lattice Gauge Theory, a system with severe metastabilities at its first order phase transition point that make the determination of observables near the transition very difficult from a numerical point of view.We find that in the LLR approach correlation times near criticality grow at most quadratically with the volume, as opposed to the exponential growth that one expects with importance sampling methods.This investigation shows the efficiency of the LLR method when dealing with systems having a first order phase transition.These results suggest that the LLR method can be efficient at overcoming numerical metastabilities in other classes of systems with a multi-peaked probability distribution, such as those with rough free energy landscapes (as commonly found, for instance, in models of protein folding or spin glasses).The rest of the paper is organised as follows.In Sect. 2 we cover the formal general aspects of the algorithm.The investigation of Compact U(1) Lattice Gauge Theory is reported in Sect.3. A critical analysis of our findings, our conclusions and our future plans are presented in Sect. 4. Finally, some technical material is discussed in the appendix.Some preliminary results of this study have already been presented in [12].
2 Numerical determination of the density of states

The density of states
Owing to formal similarities between the two fields, the approach we are proposing can be applied to both Statistical Mechanics and Lattice Field Theory systems.In order to keep the discussion as general as possible, we shall introduce notations and conventions that can describe simultaneously both cases.We shall consider a system described by the set of dynamical variables φ, which could represent a set of spin or field variables and are assumed to be continuous.The action (in the field theory case) or the Hamiltonian (for the statistical system) is indicated by S and the coupling (or inverse temperature) by β.
Since the product βS is dimensionless, without loss of generality we will take both S and β dimensionless.We consider a system with a finite volume V , which will be sent to infinity in the final step of our calculations.The finiteness of V in the intermediate steps allows us to define naturally a measure over the variables φ, which we shall call Dφ.Properties of the system can be derived from the function which defines the canonical partition function for the statistical system or the path integral in the field theory case.The density of state (which is a function of the value of S[φ] = E) is formally defined by the integral In terms of ρ(E), Z takes the form The vacuum expectation value (or ensemble average) of an observable O which is function of E can be written as Hence, a numerical determination of ρ(E) would enable us to express Z and O as numerical integrals of known functions in the single variable E. This approach is inherently different from conventional Monte-Carlo calculations, which relie on the concept of importance sampling, i.e. the configurations contributing to the integral are generated with probability P β (E) = ρ(E) e βE /Z(β) .
Owing to this conceptual difference, the method we are proposing can overcome notorious drawbacks of importance sampling techniques.

The LLR method
We will now detail our approach to the evaluation of the density of states by means of a lattice simulations.Our initial assumption is that the density of states is a regular function of the energy that can be always approximated in a finite interval by a suitable functional expansion.If we consider the energy interval [E k , E k + δ E ], under the physically motivated assumption that the density of states is a smooth function in this interval, the logarithm of the latter quantity can be written, using Taylor's theorem, as Thereby, for a given action E, the integer k is chosen such that Our goal will be to devise a numerical method to calculate the Taylor coefficients and to reconstruct from these an approximation for the density of states ρ(E).By introducing the intrinsic thermodynamic quantities, T k (temperature) and c k (specific heat) by we expose the important feature that the target coefficients a k are independent of the volume while the correction R k (E) is of order δ 2 E /V .In all practical applications, R k will be numerically much smaller than a k δ E .For a certain parameter range (i.e., for the correlation length smaller than the lattice size), we can analytically derive this particular volume dependence of the density derivatives.Details are left to the appendix.Using the trapezium rule for integration, we find in particular Using this equation recursively, we find Note that N δ E = O(1).Exponentiating (2.3) and using (2.7), we obtain where we have defined an overall multiplicative constant by We are now in the position to introduce the piecewise-linear and continuous approximation of the density of states by i.e., N is chosen in such a way that E N ≤ E < E N + δ E for a given E. With this definition, we obtain the remarkable identity which we will extensively use below.We will observe that ρ(E) spans many orders of magnitude.The key observation is that our approximation implements exponential error suppression, meaning that ρ(E) can be approximated with nearly-constant relative error despite it may reach over thousands of orders of magnitude: (2.12) We will now present our method to calculate the coefficients a k .To this aim, we introduce the action restricted and re-weighted expectation values [7] with a being an external variable: where we have used (2.1) to express N k as an ordinary integral.We also introduced the modified Heaviside function If the observable only depends on the action, i.e., W [φ] = O(S[φ]), (2.13) simplifies to Let us now consider the specific action observable and the solution a of the non-linear equation ∆E k (a) = 0 . (2.17) Inserting ρ(E) from (2.8) into (2.15) and defining ∆a = a k − a, we obtain: Let us consider for the moment the function .
It is easy to check that F is monotonic and vanishing for ∆a = 0: We therefore conclude for any δ E that if (2.18) does have a solution, this solution is unique.
For sufficiently small δ E there is a solution, and, hence, the only solution is given by: The later equation is at the heart of the LLR algorithm: it details how we can obtain the log-rho derivative by calculating the Monte-Carlo average ∆E k (a) (using (2.13)) and solving a non-linear equation, i.e., (2.17).
In the following, we will discuss the practical implementation by addressing two questions: (i) How do we solve the non-linear equation?(ii) How do we deal with the statistical uncertainty since the Monte-Carlo method only provides stochastic estimates for the expectation value ∆E k (a)?
Let us start with the standard Newton-Raphson method to answer question (i).Starting from an initial guess a (0) for the solution, this method produces a sequence a (0) → a (1) → a (2) → . . .→ a (n) → a (n+1) . . ., which converges to the true solution a k .Starting from a (n) for the solution, we would like to derive an equation that generates a value a (n+1) that is even closer to the true solution: Using the definition of ∆E k a (n+1) in (2.18) with reference to (2.16) and (2.15), we find: We thus find for the improved solution: We can convert the Newton-Raphson recursion into a simpler fixed point iteration if we assume that the choice a (n) is sufficiently close to the true value a k such that Without affecting the precision with which the solution a of (2.18) can be obtained, we replace Hence, the Newton-Raphson iteration is given by We point out that one fixed point of the above iteration, i.e., a which, indeed, is the correct solution.We have already shown that the above equation has only one solution.Hence, if the iteration converges at all, it necessarily converges to the true solution.Note that convergence can always be achieved by suitable choice of under-relaxation.We here point out that the solution to question (ii) above will involve a particular type of under-relaxation.
Let us address the question (ii) now.We have already pointed out that we have only a stochastic estimate for the expectation value ∆E k (a) and the convergence of the Newton-Raphson method is necessarily hampered by the inevitable statistical error of the estimator.This problem, however, has been already solved by Robbins and Monroe [13].
For completeness, we shall now give a brief presentation of the algorithm.The starting point is the function M (x), and a constant α, such that the equation M (x) = α has a unique root at x = θ.M (x) is only available by stochastic estimation using the random variable N (x): with E[N (x)] being the ensemble average of N (x).The iterative root finding problem is of the type where c n is a sequence of positive numbers sizes satisfying the requirements It is possible to prove that under certain assumptions [13] on the function M (x) the lim n→∞ x n converges in L 2 and hence in probability to the true value θ.A major advance in understanding the asymptotic properties of this algorithm was the main result of [13].If we restrict ourselves to the case where σ 2 ξ is the variance of the noise.Hence, the optimal value of the constant c, which minimises the variance is given by Adapting the Robbins-Monro approach to our root finding iteration in (2.24), we finally obtain an under-relaxed Newton-Raphson iteration which is optimal with respect to the statistical noise during iteration.

Observables and convergence with δ E
We have already pointed out that expectation values of observables depending on the action only can be obtained by a simple integral over the density of states (see (2.2)).
Here we develop a prescription for determining the values of expectations of more general observables by folding with the numerical density of states and analyse the dependence of the estimate on δ E .
Let us denote a generic observable by B(φ).Its expectation value is defined by In order to relate to the LLR approach, we break up the latter integration into energy intervals: (2.32) Note that B[φ] does not depend on δ E .We can express B[φ] in terms of a sum over double-bracket expectation values by choosing ).Without any approximation, we find: where 14).The above result can be further simplified by using (2.11):

.35)
We now define the approximation to B[φ] by Since the double-bracket expectation values do not produce a singularity if δ E → 0, i.e., lim The latter formula together with (2.36) provides access to all types of observables using the LLR method with little more computational resources: Once the Robbins-Monro iteration (2.30) has settled for an estimate of the coefficient a k , the Monte-Carlo simulation simply continues to derive estimators for the double-bracket expectation values in (2.36) and (2.37).
With the further assumption that the double-bracket expectation values are (semi-)positive, an even better error estimate is produced by our approach: This implies that the observable B[φ] can be calculated with an relative error of order δ 2 E .Indeed, we find from (2.33,2.34,2.35)that Thereby, we have used The assumption of (semi-)positive double expectation values is true for many action observables, and possibly also for Wilson loops, whose re-weighted and action restricted double expectation values might turn out to be positive (as it is the case for their standard expectation values).In this case, our method would provide an efficient determination of those quantities.This is important in particular for large Wilson loop expectation values, since they are notoriously difficult to measure with importance sampling methods (see e.g.[14]).We also note that, in order to have an accurate determination of a generic observable, any Monte-Carlo estimate of the double expectation values must be obtained to good precision dictated by the size of δ E .A detailed numerical investigation of these and related issues is left to future work.
For the specific case that the observable B[φ] only depends on the action S[φ], we circumvent this problem and evaluate the double-expectation values exactly.To this aim, we introduce for the general case W [φ] k the generalised density w k (E) by (2.41) We then point out that if W [φ] is depending on the action only, i.e., W [φ] = f (S[φ]), we obtain: With the definition of the double expectation value (2.13), we find: Rather than calculating W [φ] k by Monte-Carlo methods, we can analytically evaluate this quantity (up to order O(δ 2 E ) ).Using the observation that for any smooth (C 2 ) function g and using this equation for both, numerator and denominator of (2.42), we conclude that Let us now specialise to the case that is relevant for (2.39) with B depending on the action only: This leaves us with (2.45) Inserting (2.43) together with (2.44) into (2.36),we find: (2.47) Below, we will numerically test the quality of expectation values obtained by the LLR approach using action observables only, i.e., B[φ] = O(S[φ]).We will find that we indeed achieve the predicted precision in δ 2 E for this type of observables (see below Fig. 6).

The numerical algorithm
So far, we have shown that a piecewise continuous approximation of the density of states that is linear in intervals of sufficiently small amplitude δ E allows us to obtain a controlled estimate of averages of observables and that the angular coefficients a i of the linear approximations can be computed in each interval i using the Robbins-Monro recursion (2.30).
Imposing the continuity of log ρ(E), one can then determine the latter quantity up to an additive constant, which does not play any role in cases in which observables are standard ensemble averages.The Robbins-Monro recursion can be easily implemented in a numerical algorithm.Ideally, the recurrence would be stopped when a tolerance for a i is reached, i.e. when a with (for instance) set to the precision of the computation.When this condition is fulfilled, we can set a i = a (n+1) i . However, one has to keep into account the fact that the computation of ∆E i requires an averaging over Monte-Carlo configurations.This brings into play considerations about thermalisation (which has to be taken into account each time we send a ), the number of measurements used for determining ∆E i at fixed a (n) i and -last but not least -fluctuations of the a (n) i themselves.Following those considerations, an algorithm based on the Robbins-Monro recursion relation should depend on the following input (tunable) parameters: • N TH , the number of Monte-Carlo updates in the restricted energy interval before starting to measure expectation values; • N SW , the number of iterations used for computing expectation values; • N RM , the number of Robbins-Monro iterations for determining a i ; • N B , number of final values from the Robbins-Monro iteration subjected to a subsequent bootstrap analysis.
The version of the LLR method proposed and implemented in this paper is reported in an algorithmic fashion in the box Algorithm 1.This implementation differs from that provided in [7,8] by the replacement of the originally proposed root-finding procedure based on a deterministic Newton-Raphson like recursion with the Robbins-Monro recursion, which is better suited to the problem of finding zeros of stochastic equations.Since the a i are Algorithm 1: The LLR method as implemented in this work.
Evolve the whole system with an importance sampling algorithm for one sweep according to the probability distribution accepting only configuration such that Compute E (j) , the value of the energy in the current configuration j; Repeat N B times to produce N B candidates a i for a subsequent bootstrap analysis determined stochastically, a different reiteration of the algorithm with different starting conditions and different random seeds would produce a different value for the same a i .The stochastic nature of the process implies that the distribution of the a i found in different runs is Gaussian.The generated ensemble of the a i can then be used to determine the error of the estimate of observables using analysis techniques such as jackknife and bootstrap.
The parameters E min and E max depend on the system and on the phenomenon under investigation.In particular, standard thermodynamic considerations on the infinite volume limit imply that if one is interested in a specific range of temperatures and the studied observables can be written as statistical averages with Gaussian fluctuations, it is possible Figure 1: Left: For contiguous energy intervals if a transition between configurations with energy in the same interval requires going through configurations with energy that are outside that interval, the simulation might get trapped in one of the allowed regions (in green).Right: For overlapping energy intervals with replica exchange, the simulation can travel from one allowed region to the other through excursions to the upper interval.
to restrict the range of energies between the energy that is typical of the smallest considered temperature and the energy that is typical of the highest considered temperature.Determining a reasonable value for the amplitude of the energy interval δ E and the other tunable parameters N SW , N TH , N RM and N A requires a modest amount of experimenting with trial values.In our applications we found that the results were very stable for wide ranges of values of those parameters.Likewise, āi , the initial value for the Robbins-Monro recursion in interval i, does not play a crucial role; when required and possible, an initial value close to the expected result can be inferred inverting E(β) , which can be obtained with a quick study using conventional techniques.The average . . .imposes an update that restricts configurations to those with energies in a specific range.In most of our studies, we have imposed the constraint analytically at the level of the generation of the newly proposed variables, which results in a performance that is comparable with that of the unconstrained system.Using a simple-minded more direct approach, in which one imposes the constraint after the generation of the proposed new variable, we found that in most cases the efficiency of Monte-Carlo algorithms did not drop drastically as a consequence of the restriction, and even for systems like SU(3) (see Ref. [7]) we were able to keep an efficiency of at least 30% and in most cases no less than 50% with respect to the unconstrained system.

Ergodicity
Our implementation of the energy restricted average • • • assumes that the update algorithm is able to generate all configurations with energy in the relevant interval starting from configurations that have energy in the same interval.This assumption might be too strong when the update is local 2 in the energy (i.e. each elementary update step changes the energy by a quantity of order one for a system with total energy of order V ) and there are topological excitations that can create regions with the same energy that are separated by high energy barriers.In these cases, which are rather common in gauge theories and statistical mechanics3 , generally in order to go from one acceptable region to the other one has to travel through a region of energies that is forbidden by an energy-restricted update method such as the LLR.Hence, by construction, in such a scenario our algorithm will get trapped in one of the allowed regions.Therefore, the update will not be ergodic.In order to solve this problem, one can use an adaptation of the replica exchange method [15], as first proposed in [16].The idea is that instead of dividing the whole energy interval in contiguous sub-intervals overlapping only in one point (in the following simply referred to as contiguous intervals), one can divide it in sub-intervals overlapping in a finite energy region (this case will be referred to as overlapping intervals).With the latter prescription, after a fixed number of iterations of the Robbins-Monro procedure, we can check whether in any pairs of overlapping intervals (I 1 , I 2 ) the energy of both the corresponding configurations is in the common region.For pairs fulfilling this condition, we can propose an exchange of the configurations with a Metropolis probability I 1 and a I 2 are the values of the parameter a at the current n-th iterations of the Robbins-Monro procedure respectively in intervals I 1 and I 2 and E C 1 (E C 2 ) is the value of the energy of the current configuration C 1 (C 2 ) of the replica in the interval I 1 (I 2 ).If the proposed exchange is accepted, C 1 → C 2 and C 2 → C 1 .With repeated exchanges of configurations from neighbour intervals, the system can now travel through all configuration space.A schematic illustration of how this mechanism works is provided in Fig. 1.As already noticed in [16], the replica exchange step is amenable to parallelisation and hence can be conveniently deployed in calculations on massively parallel computers.Note that the replica exchange step adds another tunable parameter to the algorithm, which is the number N SWAP of configurations swaps during the Monte-Carlo simulation at a given Monte-Carlo step.A modification of the LLR algorithm that incorporates this step can be easily implemented.

Reweighting with the numerical density of states
In order to screen our approach outlined in subsections 2.2 and 2.3 for ergodicity violations and to propose an efficient procedure to calculate any observable once an estimate for the density of states has been obtained, as an alternative to the replica exchange method discussed in the previous section, we here introduce an importance sampling algorithm with reweighting with respect to the estimate ρ.This algorithm features short correlation times even near critical points.Consider for instance a system described by the canonical ensemble.We define a modified Boltzmann weight W B (E) as follows: (2.50)Here E min and E max are two values of the energy that are far from the typical energy of interest E: E min E E max . (2.51) If conventional Monte-Carlo simulations can be used for numerical studies of the given system, we can chose β 1 and β 2 from the conditions (2.52) If importance sampling methods are inefficient or unreliable, β 1 and β 2 can be chosen to be the micro-canonical β µ corresponding respectively to the density of states centred in E min and E max .These β µ are outputs of our numerical determination ρ(E).The two constants c 1 and c 2 are determined by requiring continuity of W B (E) at E min and at E max : lim Let ρ(E) be the correct density of state of the system.If ρ(E) = ρ(E), then for and a Monte-Carlo update with weights W B (E) drives the system in configuration space following a random walk in the energy.In practice, since ρ(E) is determined numerically, upon normalisation and the random walk is only approximate.However, if ρ(E) is a good approximation of ρ(E), possible free energy barriers and metastabilities of the canonical system can be successfully overcome with the weights (2.50).Values of observables for the canonical ensemble at temperature T = 1/β can be obtained using reweighting: These criteria can be used to confirm that the numerically determined ρ(E) is a good approximation of ρ(E).The flatness of the histogram is not influenced by the β of interest in the original multi-canonical simulation.This is particularly important for first order phase transitions, where traditional Monte-Carlo algorithms have a tunnelling time that is exponentially suppressed with the volume of the system.Since the modified ensemble relies on a random walk in energy, the tunnelling time between two fixed energy densities is expected to grow only as the square root of the volume.This procedure of using a modified ensemble followed by reweighting is inspired by the multi-canonical method [18], the only substantial difference being the recursion relation for determining the weights.Indeed for U(1) lattice gauge theory a multi-canonical update for which the weights are determined starting from a Wang-Landau recursion is discussed in [19].We also note that the procedure used here to restrict ergodically the energy interval between E min and E max can be easily implemented also in the replica exchange method analysed in the previous subsection.
3 Application to Compact U(1) Lattice Gauge Theory

The model
Compact U(1) Lattice Gauge Theory is the simplest gauge theory based on a Lie group.Its action is given by where β = 1/g 2 , with g 2 the gauge coupling, x is a point of a d-dimensional lattice of size L d and µ and ν indicate two lattice directions, indicised from 1 to d (for simplicity, in this work we shall consider only the case d = 4), θ µν plays the role of the electromagnetic field tensor: if we associate the compact angular variable θ µ (x) ∈ [−π; π[ with the link stemming from i in direction μ, The path integral of the theory is given by the latter identity defining the Haar measure of the U(1) group.
The connection with the framework of SU(N) lattice gauge theories is better elucidated if we introduce the link variable U µ (x) = e iθµ(x) . (3.4) With this definition, S can be rewritten as with the plaquette variable and U µ (x) is the complex conjugate of U µ (x).Working with the variables U µ (x) allows us to show immediately that S is invariant under U(1) gauge transformations, which act as with λ(x) ∈ [−π; π[ a function defined on lattice points.The connection with U(1) gauge theory in the continuum can be shown by introducing the lattice spacing a and the non-compact gauge field a A µ (x) = θ µ (x)/g, so that Taking a small and expanding the cosine leads us to x,µ,ν with ∆ µ the forward difference operator.In the limit a → 0, we finally find with F µν being the usual field strength tensor.This shows that in the classical a → 0 limit S becomes the Euclidean action of a free gas of photons, with interactions being related to the neglected lattice corrections.It is worth to remark that this classical continuum limit is not the continuum limit of the full theory.In fact, this classical continuum limit is spoiled by quantum fluctuations.These prevent the system from developing a second order transition point in the a → 0 limit, which is a necessary condition to be able to remove the ultraviolet cutoff introduced with the lattice discretisation.The lack of a continuum limit is related to the fact that the theory is strongly coupled in the ultraviolet.Despite the non-existence of a continuum limit for Compact U(1) Lattice Gauge Theory, this lattice model is still interesting, since it provides a simple realisation of a weakly first order phase transition.This bulk phase transition separates a confining phase at low β (whose existence was pointed out by Wilson [20] in his seminal work on Lattice Gauge Theory) from a deconfined phase at high β, with the transition itself occurring at a critical value of the coupling β c 1. Rather unexpectedly at first side, importance-sampling Monte-Carlo studies of this phase transitions turned out to be demanding and not immediate to interpret, with the order of the transition having been debated for a long time (see e.g.[21][22][23][24][25][26][27][28][29][30]).The issue was cleared only relatively recently, with investigations that made a crucial use of supercomputers [31,32].What makes the transition difficult to observe numerically is the role played in the deconfinement phase transition by magnetic monopoles [33], which condense in the confined phase [33,34].The existence of topological sectors and the presence of a transition with exponentially suppressed tunnelling times can provide robust tests for the efficiency and the ergodicity of our algorithm.This motivates our choice of Compact U(1) for the numerical investigation presented in this paper.

Simulation details
The study of the critical properties of U(1) lattice gauge theory is presented in this section.In order to test our algorithm, we investigated the behaviour of specific heath as function of the volume.This quantity has been carefully investigated in previous studies, and as such provides a stringent test of our procedure.In order to compare data across different sizes, our results will be often provided normalised to the number of plaquette 6L 4 = 6V .We studied lattices sizes ranging from 8 4 to 20 4 and for each lattice size we computed the density of states ρ(E) in the entire interval E min ≤ E ≤ E max (see Tab. 1).The rational behind the choice of the energy region is that it must be centred around the critical energy and it has to be large enough to study all the critical properties of the theory, i.e. every observable evaluated has to have support in this region and have virtually no correction coming from the choice of the energy boundaries.We divided the energy interval in steps of δ E and for each of the sub-interval we have repeated the entire generation of the log-linear density of states function and evaluation of the observables N B = 20 times to create the bootstrap samples for the estimate of the errors.The values of the other tunable parameters of the algorithm used in our study are reported in Tab. 1.An example determination of one of the a i is reported in Fig. 3.The plot shows the rapid convergence to the asymptotic value and the negligible amplitude of residual fluctuations.Concerning the cost of the simulations, we found that accurate determinations of observables can be obtained with modest computational resources compared to those needed in investigations of the system with importance sampling methods.For instance, the most costly simulation presented here, the investigation of the 20 4 lattice, was performed on 512 cores of Intel Westmere processors in about five days.This needs to be contrasted with the fact that in the early 2000's only lattices up to 18 4 could be reliably investigated with importance sampling methods, with the largest sizes requiring supercomputers [31,32].One of our first analyses was a screening for potential ergodicity violations with the LLR approach.As detailed in subsection 2.5, these can emerge for LLR simulations using contiguous intervals as it is the case for the U(1) study reported in this paper.To this aim, we calculated the action expectation value E for a 12 4 lattice for several values using the LLR method and using the re-weighting with respect to the estimate ρ.Since the latter approach is conceptually free of ergodicity issues, any violations by the LLR method would be flagged by discrepancy.Our findings are summarised in Fig. 2 and the corresponding table.We find good agreement for the results from both methods.This suggests that topological objects do not generate energy barriers that trap our algorithm in a restricted section of configuration space.Said in other words, for this system the LLR method using contiguous interval seems to be ergodic.

Volume dependence of log ρ and computational cost of the algorithm
As first investigation we have performed a study of the scaling properties of the a i as function of the volume.In Fig. 4 we show the behaviour of the a i with the lattice volume.
The estimates are done for a fixed δ E /V , where the chosen value for the ratio fulfils the request that within the errors all our observables are not varying for δ E → 0 (we report on the study of δ E → 0 in section 3.5).As it is clearly visible from the plot the data are scaling toward a infinite volume estimate of the a i for fixed energy density.As mentioned before, the issue facing importance sampling studies at first order phase transitions are connected with tunnelling times that grow exponentially with the volume.With the LLR method, the algorithmic cost is expected to grow with the size of the system as V 2 , where one factor of V comes from the increase of the size and the other factor of V comes from the fact that one needs to keep the energy interval per unit of volume δ E /V fixed, as in the large-volume limit only intensive quantities are expected to determine the physics.One might wonder whether this apparently simplistic argument fails at the first order phase transition point.This might happen if the dynamics is such that a slowing down takes place at criticality.In the case of Compact U(1), for the range of lattice sizes studied here, we have found that the computational cost of the algorithm is compatible with a quadratic increase with the volume.

Numerical investigation of the phase transition
Using the density of states it is straightforward to evaluate, by direct integration (see subsection 2.3), the expectation values of any power of the energy and evaluate thermodynamical quantities like the specific heat As usual we define the pseudo-critical coupling β c (L) such as the coupling at which the peak of the specific heat occurs for a fixed volume.The peak of the specific heat has been located using our numerical procedure and the error bars are computed using the bootstrap method.Our results are summarised in Tab. 2 with a comparison with the values in [31].
Once again, the agreement testify the good ergodic properties of the algorithm.Table 2: β c (L) evaluated with the LLR algorithm and reference data from [31].
data it is possible to make a precise estimate of the infinite volume critical beta by means of a finite size scaling analysis.The finite size scaling of the pseudo-critical coupling is given by where β c is the critical coupling.We fit our data with the function in Eq. (3.11), the results are reported in Tab.
Table 4: C V (β c (L)) evaluated with the LLR algorithm and reference data from [31].
Results for a 20 4 lattice have never been reported before in the literature.
Another quantity easily accessible is the latent heat, this quantity can be related to the height of the peak of the specific heat at the critical temperature through: where G is the latent heat.The results for this observable are reported in Tab. 4. We fit the result with Eq. (3.12), see Tab. 5.
The latent heat can be obtained also from the knowledge of the location of the peaks of the probability density at β c (of infinite volume), indeed in this case the latent heat is equal to energy gap between the peaks.This direct measure can be used as crosscheck of the previous analysis.In the language of the density of states the probability density is simply given by

.13)
We have performed the study of the location in energy of the two peaks of P βc (E) and we have reported them in Tab. 6.Also in this case we have performed a finite size scaling analysis to extract the infinite volume behaviour:  A fit of the values in Tab.6 yields χ 2 red,1 = 0.67, 1 = 0.6279(9) and χ 2 red,2 = 0.2, 2 = 0.65485(4).The latent heat can be evaluated as G = 2 − 1 = 0.0270 (9) which is in perfect agreement with the estimates obtained by studying the scaling of the specific heath.partition function to performing an ordinary integral.Wang-Landau type algorithms perform Markov chain Monte-Carlo updates with respect to ρ while improving the estimate for ρ during simulations.The LLR approach, firstly introduced in [7], uses a non-linear stochastic equation (see (2.17)) for this task and is particularly suited for systems with continuous degrees of freedom.To date, the LLR method has been applied to gauge theories in several publications, e.g.[8][9][10]12], and has turned out in practice to be a reliable and robust method.In the present paper, we have thoroughly investigated the foundations of the method and have presented high-precision results for the U(1) gauge theory to illustrate the excellent performance of the approach.
Two key features of the LLR approach are: (i) It solves an overlap problem in the sense that the method can specifically target the action range that is of particular importance for an observable.This range might easily be outside the regime for which standard MC methods would not be able to produce statistics.
(ii) It features exponential error suppression: although the density of states ρ spans many orders of magnitude, its linear approximation ρ has a nearly-constant relative error (see subsection 2.2) and the numerical determination of ρ preserves this level of accuracy.
We point out that feature (i) is not exclusive of the LLR method, but is quite generic for multi-canonical techniques [18], Wang-Landau type updates [4] or hybrids thereof [19].
Key ingredient for the LLR approach is the double-bracket expectation value [7] (see (2.13)).
It appears as a standard Monte-Carlo expectation value over a finite action interval of size δ E and with the density of states as a re-weighting factor.The derivative of the density of states a(E) emerges from an iteration involving these Monte-Carlo expectation values.This implies that their statistical error interferes with the convergence of the iteration.This might introduce a bias preventing the iteration to converge to the true derivative a(E).We resolved this issue by using the Robbins-Monro formalism [13]: we showed that a particular type of under-relaxation produces a normal distribution of potential values a(E) with the mean of this distribution coinciding with the correct answer (see subsection 2.2).
In this paper, we also addressed two concerns, which were raised in the wake of the publication [7]: (1) The LLR simulations restrict the Monte-Carlo updates to a finite action interval and might therefore be prone to ergodicity violations.
(2) The LLR approach seems to be limited to the calculation of action dependent observables only.
To address the first issue, we have proposed in subsections 2.5 and 2.6 two procedures that are conceptually free of ergodicity violations.The first method is based upon the replica exchange method [15,16]: using overlapping action ranges during the calculation of the double-bracket expectation values offers the possibility to exchange the configurations of neighbouring action intervals with appropriate probability (see subsection 2.5 for details).
The second method is a standard Monte-Carlo simulation but with the inverse of the estimated density of states, i.e., ρ−1 (E), as re-weighting factor.The latter approach falls into the class of ergodic Monte-Carlo update techniques and is not limited by a potential overlap problem: if the estimate ρ is close to the true density ρ, the Monte-Carlo simulation is essentially a random walk in configuration space sweeping the action range of interest.
To address issue (2), we firstly point out that the latter re-weighting approach produces a sequence of gauge field configurations that can be used to calculate any observable by averaging with the correct weight.Secondly, we have developed in subsection 2.2 the formalism to calculate any observable by a suitable sum over a combination of the density of states and double-bracket expectation values involving the observable of interest.We were able to show that the order of convergence (with the size δ E of the action interval) for these observables is the same as for ρ itself (i.e., O(δ 2 E )).
In view of the features of the density of states approach, our future plans naturally involve investigations that either are enhanced by the direct access to the partition function (such as the calculation of thermodynamical quantities) or that are otherwise hampered by the overlap problem.These, most notably, include complex action systems such as cold and dense quantum matter.The LLR method is very well equipped for this task since it is based upon Monte-Carlo updates with respect to the positive (and real) estimate of the density of states and features exponential error suppression which might beat the resulting overlap problem.Indeed, a strong sign problem was solved by LLR techniques using the original degrees of freedom of the Z 3 spin model [8,9].We are currently extending these investigations to other finite density gauge theories.QCD at finite densities for heavy quarks (HDQCD) is work in progress.We have plans to extend the studies to finite density QCD with moderate quark masses.

Figure 3 :
Figure 3: Estimated a i as a function of the Robbins-Monro iteration, on a 20 4 lattice and for action E/(6V ) = 0.59009548 at the centre of the interval with δ E /V = 1.91 × 10 −4 .

Figure 4 :
Figure 4: Estimate of a i as function of the energy density for various volume, the right panel is a zoom of the interesting region.

Figure 5 :
Figure 5: Probability density for L = 20 at β c .The probability is plotted at β c of infinite volume hence the peaks are not of equal height.

Figure 6 :
Figure 6: The peak of the C V (β C (L)) as function δ E .
[17] where denotes average over the canonical ensemble and W average over the modified ensemble defined in (2.50).The weights W B (E) guarantees ergodic sampling with small auto-correlation time for the configurations with energies E such that E min ≤ E ≤ E max , .56)doesnotpresent any overlap problem.The role of E min and E max is to restrict the approximate random walk only to energies that are physically interesting, in order to save computer time.Hence, the choice of E min , E max and of the corresponding β 1 , β 2 do not need to be fine-tuned, the only requirement being that Eqs.(2.57) hold.These conditions can be verified a posteriori.Obviously, choosing the smallest interval E max − E min where the conditions (2.57) hold optimises the computational time required by the algorithm.The weights (2.56) can be easily imposed using a metropolis or a biased metropolis[17].Again, due to the absence of free energy barriers, no ergodicity problems are expected to arise.This can be checked by verifying that in the simulation there are various tunnellings (i.e.round trips) between E min and E max and that the frequency histogram of the energy is approximately flat between E min and E max .Reasonable requirements are to have O(100 − 1000) tunnellings and an histogram that is flat within 15-20%.
Comparison between the plaquette computed with the LLR algorithm (see subsection 2.2) and via re-weighting with respect to the estimate ρ (see subsection 2.6) for a L = 12 lattice.

Table 1 :
Values of the tunable parameters of the LLR algorithm used in our numerical investigation.

Table 3 :
3. Estimates of β c for various choices of the fit parameters.In bold the best fits.

Table 5 :
Estimates of G for various choices of the fit parameters.In bold the best fits.

Table 8 :
The coefficient b dis for different lattice sizes.