Free energy of the self-interacting relativistic lattice Bose gas at finite density

The density of state approach has recently been proposed as a potential route to circumvent the sign problem in systems at finite density. In this study, using the Linear Logarithmic Relaxation (LLR) algorithm, we extract the generalised density of states, which is defined in terms of the imaginary part of the action, for the self-interacting relativistic lattice Bose gas at finite density. After discussing the implementation and testing the reliability of our approach, we focus on the determination of the free energy difference between the full system and its phase-quenched counterpart. Using a set of lattices ranging from $4^4$ to $16^4$ , we show that in the low density phase, this overlap free energy can be reliably extrapolated to the thermodynamic limit. The numerical precision we obtain with the LLR method allows us to determine with sufficient accuracy the expectation value of the phase factor, which is used in the calculation of the overlap free energy, down to values of ${\cal O}(10^{-480})$. When phase factor measurements are extended to the dense phase, a change of behaviour of the overlap free energy is clearly visible as the chemical potential crosses a critical value. Using fits inspired by the approximate validity of mean-field theory, which is confirmed by our simulations, we extract the critical chemical potential as the non-analyticity point in the overlap free energy, obtaining a value that is in agreement with other determinations. Implications of our findings and potential improvements of our methodology are also discussed.


I. INTRODUCTION
Monte Carlo simulations of the system regularised on an Euclidean spacetime lattice provide the most efficient method for extracting quantitative information from non-supersymmetric non-Abelian gauge theories at zero density. The associated general methodology consists in generating configurations according to the Boltzmann weight W (S) = e −S , with S the Euclidean action of the system, and then computing averages of observables over the generated sample. In order for the method to work, the action S needs to be real. This guarantees that the Euclidean path integral be positive definite, which allows us to interpret it as the partition function of an equivalent statistical system. However, there are physically relevant cases in which the action is complex and the Euclidean path integral becomes non-positive-definite. This gives rise to large numerical cancellations that generate noise overcoming by several orders of magnitude the typical signal one would like to observe, resulting in the inapplicability of importance sampling Monte Carlo for extracting physical observables. This cancellation, known in the literature as the sign problem (see e.g. [1] for a recent review) characterises, among others, finite density systems in Quantum Field Theory and strongly correlated electron systems in Condensed Matter Physics.
Resolving the difficulties caused by the sign problem would enable us to make substantial progress for systems such as finite density QCD, which at the moment can not be reliably studied either numerically or analytically. While a general algorithm that solves the sign problem for any system has been shown not to exist [2], it is possible to devise numerical approaches that make the problem tractable in specific cases. 1 Recent examples include the Complex Langevin approach [3], thimble regularisation [4] and the density of states route, which, following the original proposal of [5] and further refinements [6,7], has been recently revisited in [8,9]. Key to the latter two studies is the introduction of a restricted sampling [10] in terms of the independent variable used to to define the density of states. This allows us to determine the logarithm of the density of states with exponential error reduction, hence enabling us to perform extremely accurate measurements. If the precision of the determination of the density of states is high enough, one may eventually overcome cancellations that arise when computing numerical integrals. Examples of successful applications of the density of state method along these lines to systems affected by the sign problem have been provided in [8,9].
Among theories used to test techniques to tame the sign problem, the self-interacting Bose gas at non-zero chemical potential is amongst the most widely studied. Here, we shall investigate this system at finite density as a function of the chemical potential using the density of state approach. Recently, this system has been studies with Complex Langevin [3] and dualization approaches [11,12], which, together with the good agreement with mean-field theory [13], can be used to validate our results and hence to assess the viability of our proposal for the self-interacting Bose gas. These studies provide a scan in the parameters' space of the system and extract the associated phase structure. Using those results, we will fix the other action parameters to a set of values such that system is known to undergo a phase transition from zero particle net content to a dense phase for a critical value of the chemical potential. The density of states will be determined using the Linear Logarithmic Relaxation (LLR) algorithm [10,14], which has been shown to provide exponential error reduction in a range of Lattice Gauge Theories applications involving e.g., tunnelling suppression at a first order phase transition [14], the determination of the free energy for a system with shifted boundary conditions [15] and the de-correlation of the topological charge near the continuum limit [16]. Earlier studies of complex action systems with the LLR and the closely related FFA method can be found respectively in [8,[17][18][19][20][21] and [9,[22][23][24] (see also [25]). In this work, we will focus on the full -phase quenched overlap free energy. This quantity is defined as the free energy difference between the original system and the system obtained by setting to zero the imaginary part of the action, the latter being referred to as phase quenched system. The reason for choosing this observable is twofold: (1) as it will be shown below, the overlap free energy controls the severity of the sign problem; (2) in the formulation of the density of states used in this work, this free energy is a central quantity to determine the density of particles (see e.g. [8]). Given its characteristic exponential error reduction, the LLR algorithm provides an efficient method to compute this free energy, or equivalently the ratio of the two partition functions that defines this observable. The purpose of this work is to understand whether the method we are introducing can be used for characterising the phases of the system and if the numerical measurements are precise enough for studying the phase transition that occurs at a critical value of the chemical potential.
The rest of the paper is organised as follows. In Sect. II we review the density of state method and we discuss the main observables targeted in our study. Sect. III describes the self-interacting Bose gas at finite chemical potential on a spacetime lattice. A description of the numerical methodology used in our work is reported in Sects. IV to VIII. Numerical results are presented and discussed in Sect. IX. Finally, Sect. X contains a critical discussion of our findings and their implications. Some earlier results related to the current study have been reported in [26] and, more recently, in [27].

II. GENERALISED DENSITY OF STATES
The generalised density of states method provides a straightforward approach to the problem of simulating systems with a complex action where we have explicitly separated the real and imaginary part. For convenience, we have assumed a linear dependence of the latter on the chemical potential, and suppressed any parameter dependence of S R and S I , which, in particular, may also depend on µ. Without loss of generality, any partition function can be written as a functional integral over the fields, Introducing a generalised density of states (DoS) function, the partition function can be obtained from a 1-dimensional integration This reformulation suggests to split up the problem of evaluating the partition function of systems with complex action in two separate steps: first, to evaluate ρ(s) numerically to a high level of precision, and then tackle the influence of the imaginary part of the action separately by performing the remaining one dimensional integral. Although we still expect a sign problem manifesting from the need of cancellations over multiple orders of magnitude in the oscillatory integral, we have transformed a multidimensional oscillatory integration to a softer variant where the resulting one dimensional Fourier transform is separated from the Monte Carlo integration. The theory obtained by neglecting the imaginary part of the action is usually referred to as "phase quenched" and the associated partition function is given by It is worth noting that the phase-quenched system can be studied via standard importance sampling techniques, as it is a sign problem-free system. However, the physics described by it is not a good representation of the physics of the full system. To evaluate the hardness of the sign problem it is possible to evaluate the overlap factor between the full and phase quenched theory defined as the ratio of the two partition functions We can interpret this quantity as the expectation value of the phase in the phase quenched theory, from here on defined as e iϕ pq . Thanks to the symmetry ρ(s) = ρ(−s) of the DoS, the phase factor is obtained from the real part of the Fourier transform of the DoS Physically, e iϕ pq is related to the free energy difference between the full and phase quenched systems. Specifically, writing with F (F pq ) being the free energy per unit of volume V of the full (phase quenched) system, we have that Since e iϕ pq ≤ 1, the phase quenched model provides a lower bound for the free energy. To provide the expected finite ∆F in the thermodynamic limit, | log e iϕ pq | ∝ V , hence e iϕ pq has to be exponentially small in V , implying that the oscillatory integral (7) that defines it must provide cancellations over many orders of magnitude.
In our density of states approach for complex action systems, we will show that high precision data for the discretised DoS can be obtained by specialised Monte Carlo methods as provided by the linear logarithmic relaxation (LLR) algorithm. However, in order to obtain a precise and numerically stable evaluation of the Fourier integral, the full continuous DoS must be reconstructed. The relativistic Bose gas studied in this work provides a concrete frame of application for our methodology and at the same time a probing benchmark to assess how effective it can be.

III. MODEL: RELATIVISTIC BOSE GAS
In this paper we will concentrate on the relativistic Bose gas. This model has been extensively studied in the context of the sign problem of lattice field theories via Complex Langevin dynamics and mean-field approximation [13] as well as complete dualization [12]. Therefore, it is a good candidate for a benchmark study to test whether or not the DoS approach can be used to mitigate the sign problem in lattice field theories.
The lattice discretised action can be expressed as: Splitting the field into its real and imaginary part, φ x = φ 1,x +iφ 2,x , we can separate real and imaginary part of the action, It is worth noting that the phase quenched system differs from the system at zero chemical potential due to the presence of the cosh(µ) term in the real part of the action. Throughout all our studies we have set the action parameters λ = m = 1.0. This choice is motivated by the existence of extensive literature results for this set of parameters. Moreover, it has been shown that, in this setup, numerical investigations results are well described by a mean-field evaluation.

IV. LLR ALGORITHM
The LLR (Linear Logarithmic Relaxation) algorithm provides a way to estimate the slope of the density of states by solving a non-linear stochastic equation. In the following we briefly review the method and our implementation.
We define a restricted and reweighed expectation value of a general operator O as following for a given reweighing parameter a, and with ρ(s) defined as in (3). The normalisation factor N is defined as The heart of the LLR algorithm is the dynamical tuning of a, such that the reweighing factor e −a(s−S I k ) counterbalances the intrinsic density of states distribution of the system, resulting in an uniform sampling in a interval around S I k . To achieve such a result we consider the specific observable O(s) = s − S I k . As it has been shown in [14], in the limit of vanishing ∆ the expectation value of this observable has a monotonous behaviour in a, and, more importantly, it vanishes when the corresponding value of a coincides with the derivative of the DoS logarithm Corrections of order ∆ are not present due to the symmetry of the integrand function.
To solve this implicit equation for a, we use two different techniques. Initially we use a Newton-Raphson method, generating a chain of reweighing factor a (n) k according to Approximating the variance of the distribution by our actual update step can be written as As shown in Fig. 1, Newton-Raphson manages to approach the root extremely rapidly. However, due to the stochastic nature of Eq. (14), the statistical uncertainty intrinsic to the evaluation of ∆S I k (a k ) eventually prevents the Newton-Raphson method to converge to high level of precision. To overcome this issue, we employ the Robbins-Monro procedure [28], applied to the determination of the a k once the once the Newton-Raphson method starts to oscillate around the solution. The Robbins-Monro method is based on iterative procedure where the conditions on the parameters c n ensure convergence to the correct root in the limit of N RM → ∞, where by N RM we indicate the total number of Robbins-Monro steps, even in presence of non white noise in the iteration estimator. To maximise the speed of convergence, we choose c n = 1/(n+1) to maximises the damping while respecting the bounds of the Robbins-Monro procedure leading to the sequence Such a procedure converges in L 2 norm and hence in probability to the exact value, meaning that for each interval S I k − ∆/2, S I k + ∆/2 the estimator a (n) k is normally distributed around a k with variance scaling asymptotically as 1/N RM . We report in Figs. 1 and 2 a detailed study of the convergence properties of the mean and standard deviation of the Robbins-Monro and Newton-Raphson algorithms, performed over various independent replicas reported as coloured lines in the plots.
By applying this combined root finding procedure to all intervals, we can determine d log ρ/ds for an uniformly distributed set of S I k values in the imaginary phase domain. In Sect. VI we will discuss how to extend these results of the LLR procedure to the full domain of the DoS to fully reconstruct ρ(s).

V. LLR INTRINSIC BIAS
As it has been discussed in the previous section the LLR algorithm is exact for ∆ → 0. However, this regime is unfavourable in numerical simulations as the Robbins-Monro step size scales as ∆ −2 leading to huge jumps in the root finding procedure and a consequent really long convergence time. For this reason we are interested in studying the behaviour of the LLR algorithm when ∆ is small, but not so much that the higher order correction to the DoS are negligible compared to the linear relaxation in the interval S I k − ∆/2, S I k + ∆/2 . To do so we consider ∆S I k (a) with a = a k = d log ρ/ds | s=S I k , writing for ease of notation ρ(s) = exp(f (s)) and including also higher order corrections In the above the first order (O ∆ 2 ) term vanishes for a = a k = f (S I k ), the second order (O ∆ 3 ), linked to f (S I k ), vanishes for symmetry of the integral, making the term of third order in the derivative the leading term. This term has the important characteristic of not depending on a, the reweighing parameter, meaning that it will introduce the same systematics in the Robbins-Monro procedure regardless of the distance from the root. In particular we can treat this term as an additive shift to Eq. (14), We can now evaluate what is the impact of this additive term on the reweighing parameter a by solving the previous equation with the lhs set to zero. Obtaining, Therefore the bias will depend on two parameters: f (3) (S I k ), the third derivative of the logarithm of the DoS that is system specific thus impossible to control a priori, and, as expected, ∆, the interval width. As shown in Fig. 3  to define a region for which the bias influence is negligible compared to the statistical uncertainty of independent simulations.

VI. DOS RECONSTRUCTION TECHNIQUES
We aim at a faithful reconstruction of the DoS of the system over the whole domain. Thanks to the knowledge of d log ρ(s)/ds many different reconstruction strategies can be formulated. In the following we will present two different choices of reconstruction and we will highlight the biases associated with each choice.
Assuming the logarithmic derivative to be constant in each interval leads to the piecewise definition ρ pw (s) = kρ k (s) witĥ andρ k (s) = 0 for s outside the interval. The parameters C k are chosen to ensure continuity The LLR method achieves exponential error suppression, meaning that the relative error of ρ(s) stays constant throughout the entire range of the imaginary action. However, the piecewise approximation introduces a finite number of second order discontinuity in the imaginary action domain where the neighbouring exponentials are linked at the edge of the intervals. Such discontinuities will lead to precision issues in the evaluation of the oscillatory integral, Eq. (7).
To overcome this limitation, we introduce our second reconstruction technique, the polynomial fitting [8] to substantially improve on the piecewise approximation.
In the polynomial fit approach the LLR results are fitted to a polynomial p l (s) = l i=0 c i s i . An analytic integration of Eq. (14) allows to directly evaluate the a k . Due to the symmetry properties of the DoS ρ(s) = ρ(−s), only odd powers of s enter into the polynomial, p l (s) = l i=1 c (2i−1) s 2i−1 and c i are determined by fitting our LLR results for a k . The density of state resulting from the polynomial fitting can be expressed as where we are normalising the DoS to have ρ fit(l) (0) = 1 as p l (0) = 0. As displayed in Fig. 4, this approximation provides much smoother behaviour than the piecewise one, from which it shows bounded relative deviations. The jagged appearance of the plotted quantity results from the artefacts of the linear approximation involved in the reconstruction of the piecewise DoS. To compare the values of the phase factor obtained with both approximations, we define the quantity This function can be evaluated with both the piecewise or fitted definition of the DoS. In addition, I(S I max ) is related to the expectation value of the phase factor by where M is the size of an ensemble of gaussianly distributed realisations of the a k . In Fig. 5 we compare the absolute values of the partially integrated phase factor as a function of S I max for both approximations for two different volumes V = 6 4 , 10 4 at µ = 0.8 and a particular realisation of the a k . Our results show that using the piecewise approximation generates much larger fluctuations than the polynomial interpolation: for the 10 4 the relative amplitude of the fluctuations is of around 45 orders of magnitude. While, when averaged over multiple realizations the a k coefficients, both definitions give compatible results, only the polynomial interpolation provides a value that is accurately different from zero within the precision of the calculation and hence allows us to detect the severe cancellations generated by the sign problem. Indeed the piecewise approximation fails to achieve sufficient precision as the integration generates an intrinsic error of O ∆ 2 for each interval (due to the correction to the linear approximation neglected in this procedure). When the sign problem gets exponentially hard, an exponentially large number of small intervals should be taken into account to achieve the required precision in order to suppress the intrinsic error. On the other hand, the polynomial approximated DoS seems to show no difficulty at obtaining an accurate result broadly compatible with the mean-field calculation also for the harder V = 10 4 case. The different precision obtained by the two methods can be analysed further by computing the phase average. In order to check the convergence of the result, we study the latter quantity for different coarse-graining of the a k , obtained by taking subsets with different spacing between consecutive values (subsampling). In Fig. 6, we contrast the level of precision on e iϕ pq obtained with the two methods as the spacing between two central values of the S I used to calculate the a k varies. For our finest determinations (i.e., for N/N tot approaching the value of one, which means that all the values we have determined are used in the reconstruction), the data converge to an asymptotic value, from which they deviate for coarser spacing (N/N tot 1). However, while the polynomial fit provides a reliably accurate determination of the average sign, the use of the piecewise interpolation generates a statistical error that makes the result compatible with zero. Fig. 7 shows the quality of the determination of the free energy difference ∆F corresponding to the polynomial fit.
In addition to the polynomial fit, we have performed other interpolations to understand the regime of validity of the results. In particular, we have interpolated the a k using • Expansions in an L 2 basis (in particular, using Hermite functions); • Continuous local fits (loess/lowess), whereby a local low-order polynomial fit is convoluted with localised weight functions; • Gaussian processes, which use a multi-variate Gaussian a priori ansatz for modelling the distribution of correlations among n-tuples of observations.
Somehow surprisingly, all these methods produced results that were less accurate than the simple (and a priori simplistic) polynomial interpolator, basically failing at disentangling a non-zero average phase from the noise when a hard sign problem is present. The different reasons for the observed failures are instructive: • The L 2 expansion converges slowly, hence requiring a high number of terms or equivalently a high number of fitted parameters, which results in detectable overfitting (we will discuss overfitting in the context of the polynomial interpolation below); • Continuous local fits are too sensitive to the locally projecting functions, generating noise at a frequency that is roughly the inverse of the amplitude of the window on which one performs the projection; • Gaussian processes presented localised high-frequency oscillations that also resulted in noisy measurements for the Fourier transform.
Amongst those methods, perhaps the most surprising failure is associated to Padé approximants, which in general are expected to converge faster than polynomial interpolations. The better outcomes obtained with the latter may indicate that the variation of the a k with S I is indeed described by a (near-)polynomial function. We believe that this information is physically relevant for understanding the system. A possible explanation of the success of the polynomial interpolation may be inferred from the behaviour of the a k as a function of the relevant observable (e.g., the energy) for systems at zero density, whereby higher power contributions are suppressed by powers of the volume [14]. It is possible that this local property holds on a wider scale.
Having shown that, unlike other choices, the polynomial fitting approach of the DoS allows us to achieves a high level of precision for the determination of the average sign, we shall now address the numerical stability of the approach with respect to the order of the chosen polynomial and estimate possible systematics related to the determination of the maximum power of s appearing in the polynomial.

VII. FIT VALIDATION
Concerning the stability of the polynomial fit, the choice of the polynomial order l is of course of crucial importance. In this section we will illustrate how to ensure that the functional form choice avoids the two most likely source of systematics, under -fitting and over -fitting.

Under-fitting
Under-fitting happens when the proposed polynomial is too simple (the order is too low) to represent all the features of the data. In this case a χ 2 analysis of the fit residuals is able to pin down the minimum number of polynomial coefficients needed to describe the LLR results. As shown in Fig. 8, the χ 2 value decreases while increasing the order of the fitted polynomial. For the data considered at l = 7 we start to see a plateau forming. After the onset of the plateau, due to the statistical uncertainty of our data, higher order polynomials will start to pick up the statistical noise rather than improve the approximation. For this reason one could be tempted to choose simply the smallest order in accordance with the Occam's razor principle. Instead, we evaluate (7) for a various choices of the polynomial fit order. If we see a plateau also in the expectation values of the phase factor we will accept the results, otherwise we will reject the result and proceed to increase the precision of the simulation.

Over-fitting
The other way in which the fitting process could introduce a systematic error is over-fitting, when the polynomial order is so high that the fitted function will start to introduce noise not related to the statistical uncertainty of the fitted data. To control the over-fitting we are going to study the expectation values of the second order derivative of log ρ with the intent of comparing the numerical values with the derivative of the fitted polynomial.
The second derivative of log ρ can be determined numerically by evaluating restricted and reweighed expectation values (12). We start by determining (∆S I ) 2 k with a = a k , where we are writing ρ(s) = exp{f (s)}, from which it is possible to obtain the value of the second derivative as This quantity is measurable to an acceptable level of statistical relevance in our simulations, as shown in Fig. 9, despite coming from an evaluation of the second moment of the distribution. Rather than  Figure 9. Second logarithmic derivative f (S I k ) obtained from the evaluation of (∆S I ) 2 from a simulation at lattice volume of 10 4 at µ = 0.8. using the second derivative directly in the fitting procedure we look at how well the polynomial fit of the a k describes this quantity (i.e. we perform an a posterior validation of the functional form of fit). We compare f (s) with the derivative of the polynomial fit p l by defining a χ 2 -like function To illustrate this principle we consider a sub-sampling of our data, so that the effects of the over-fitting are evident even at lower order of the fitting polynomial. The results of this analysis are shown in Fig. 10, where three different behaviours are visible: for small orders of the polynomial fit(l = 1 ∼ 5) the discrepancies are big as the polynomial is not a good approximation of the original data; in the middle section (l = 7 ∼ 11) we see a plateau indicating that in this region the fit not only approximates well f (s), but also its derivative; for high orders (l 13) we see a clear indication of over-fitting, as, while the χ 2 of the fit of the a k would still be on the plateau, χ 2 f shows that the fitted function does not approximate well the numerical derivative. This gives us a quantitative indication of whether the chosen functional form is over-fitting the data. When using the full set of data no sign of over-fitting has been found in any of our analysis.

VIII. BIAS OPTIMISED SIMULATIONS
The following scheme ensures a bias free and performance optimised simulation: • Run a low precision simulation (fewer Monte Carlo samplings (N M C ) as well as Robbins-Monro steps (N RM )) with a small and constant ∆ for each interval, extract the values of the a k , and use those to estimate the bias over the complex action range taken into consideration.
• Scale the simulation parameters (N M C , N RM and ∆) so that bias σ a k . We use the known scalings bias ∝ ∆ 2 and σ a k ∝ (∆ · √ N M C · N RM ) −1 , and the fact that the simulation runtime is proportional to N M C · N RM .
• With the scaled parameters run a high precision simulation, the results of which will be used to rebuild the DoS.
• Finally, using the high precision results double check that the bias is negligible in comparison to the statistical noise of the results.

IX. RESULTS
Following the scheme described in the previous sections we have been able to obtain the a k estimates for a wide range of values in the chemical potential, ranging from µ = 0 to 2.0, and volumes ranging from 4 4 to 16 4 . The typical values of the simulation parameters are reported in Tab I.
A representative set of results of this evaluation has been reported in Fig. 11. A general feature of the a k as a function of S I is the appearance of a sharp change of behaviour for large S I when µ is close to a critical value, the location of the inflection point decreasing for larger volumes. If the inflection point is present in the interval of imaginary action relevant for our numerical integration the fitting procedure defined in the previous section fails to converge. As a consequence, we can estimate the free energy only if the above change of behaviour does not occur, hence for large volumes we can do the integral only outside a region around the critical µ. Phase structure away from criticality In Fig. 12 we show the results of our simulations up to V = 10 4 (reported in Tab. II) and for chemical potential values ranging from zero to µ = 2.0 and λ = m = 1.0. We stress that our procedure fails to converge for large volumes in a window close to the critical µ, while no issues are found for values of µ sufficiently far from the critical µ.
In the region µ c 1.15 (as predicted by mean-field analysis) we expect, and observe, a phase transition. A clear difference in the behaviour of the free energy is visible in the two phases, distinctly for 4 4 and 6 4 , and reasonably clearly also for 8 4 and 10 4 . In the region µ < µ c the free energy difference has been fitted to the functional form ∆F (µ) = aµ 2 + bµ 4 + cµ 6 , while in the µ > µ c a linear fit is enough to describe the behaviour of the data. By intersecting the fits in the two regions we have been able to give an estimate for the critical value of the chemical potential as well as its error via the confidence intervals of the fits. Our data are reported in Tab. III. As shown in Fig. 13, for volumes 4 4 and 6 4 the extrapolation obtained has a good level of precision (respectively .6% and .4% relative error), while the results for the larger volumes suffer from the lack of points close to the phase transition, resulting in relative errors of 1% for 8 4 and 4% for 10 4 .
Since our ability to study values near µ c decreases as the volume increases, we have not performed an extrapolation to the thermodynamic limit. However, our calculation shows that our results are compatible with the mean-field calculations [13] (µ c 1.15) as well as with the value obtained in [12] (µ c = 1.146 ± 0.001) with a dual formulation of the same theory in a work more focused to the study of the phase transition than the present one. The statistical uncertainty in our result is of the order of a few percent. A careful determination of the systematic error would require an improvement of our method in order for us to be able to simulate closer to µ c on larger lattices.  [12] (the error band is also indicated, but it is barely visible on the scale of the plot).

Low density region
Far from the phase transition the integration procedure poses no threat. Hence, we could study more precisely the low chemical potential region (µ = 0 ∼ 0.9), extending the results to higher volumes where the sign problem get exponentially harder. The minimum polynomial order required to describe the a k data in this region ranges from 5 to 9, and for all the volumes and values of the chemical potential at least three subsequent polynomial orders (i.e. {7, 9, 11}) managed to integrate to statistically comparable values.
We report in Tab. IV our results for ∆F as a function of µ for volumes up to V = 16 4 . In the same table, we show also the thermodynamic extrapolation of ∆F , obtained with the ansatz which is a good description of our data (see Fig. 14 for some representative examples showing the fit quality and Fig. 15, top, for a zoomed out picture of the extrapolation in the whole range of µ). The behaviour of ∆F as a function of µ for µ < µ c is displayed in Fig. 15, bottom.

X. DISCUSSION AND CONCLUSIONS
In this work, we have further refined the LLR method for complex action systems, studying the main sources of systematic errors in an application to the Bose gas at finite density. Using the expected scaling with the size of the imaginary action intervals for restricted sampling, we were able to eliminate for all practical purposes the error related to this discretisation. In addition, we have further investigated the necessity of interpolating the a k in order to obtain a robust result for the oscillating integral. We expect that these lessons are generalisable to other studies of complex action systems with the LLR method. Concerning the a k interpolations, we studied several possibilities. Somewhat unexpectedly, the data support the necessity of a polynomial interpolation, which has been shown to be the only one in the set of those we analysed that is able to produce a controlled result for the highly oscillating integral. In particular, we have shown that a polynomial fitting approach  produces numerically stable and reliable results for phase factors down to O 10 −480 occurring in our model for the scenarios with the hardest sign problem we have explored. Reasons for the failure of the other studied interpolating methods have been analysed. However, at the moment it is unclear whether the polynomial interpolation would have the same degree of success on other systems. Other studies in the literature (e.g. [8,20]) also find the polynomial interpolator sufficiently accurate, although no alternative methods have been considered in these works. In our study, we also established a criterium that allows us to assess the requested order of the fitting polynomial.
Armed with this machinery, we have then performed a numerical investigation of the free energy difference ∆F between the full and phase-quenched system. In the low density phase, we have been able to determine this observable and to extrapolate it to the thermodynamic limit for up to volume V = 16 4 , for chemical potential values for which the sign problem is indeed hard (as already mentioned, we have successfully resolved and compared with other methods phase factors of O 10 −480 ). Our results are compatible with those obtained with Complex Langevin, mean-field calculations and dual methods. Our method also allows us to determine ∆F for µ values that put the system in the dense phase, although the method fails in a region around µ c whose upper bound seems to increase with the volume. For the maximum volume we have simulated in the dense phase, V = 10 4 , we have been able to extract ∆F only for µ ≥ 1.6 (while the critical value is µ c 1.15). The failure of the approach in the proximity of µ c is explained by the observation that, despite the a k turn out to be very well determined, the polynomial interpolation is not able to account for a sudden change of behaviour of these coefficients as S I increases in the region that gives non-negligible contributions to the integral. We leave to future studies to understand whether a more suitable ansatz can enable us to make progress in the currently inaccessible region. Similarly, we defer to further investigations the question of whether the upper bound of the currently inaccessible region keeps increasing with µ or stabilises.
Despite those difficulties, we have shown that the intersection of two simple interpolation ansatz for ∆F defines a critical value of µ that is accurate at the order of the percent over the range of the volumes we have simulated and compares well with the current literature. Our analysis for the determination of µ c implicitly assumes the validity of mean-field (in fact, since we can not determine ∆F sufficiently close to µ c , we are insensitive to any potential critical behaviour beyond mean-field). In four dimensions, one expects that the critical region grows logarithmically with the volume. Hence, the critical domain possibly is very small and hidden in the region we can not access with our numerical simulations. However, it is an interesting question whether much larger lattices than those used currently in the literature (including also studies based on dual and Langevin methods) would be necessary in order to pin down the correct critical behaviour of the system.