The lattice Landau gauge photon propagator for 4D compact QED

In this work we report on the Landau gauge photon propagator computed for pure gauge 4D compact QED in the confined and deconfined phases and for large lattices volumes: $32^4$, $48^4$ and $96^4$. In the confined phase, compact QED develops mass scales that render the propagator finite at all momentum scales and no volume dependence is observed for the simulations performed. Furthermore, for the confined phase the propagator is compatible with a Yukawa massive type functional form. For the deconfined phase the photon propagator seems to approach a free field propagator as the lattice volume is increased. In both cases, we also investigate the static potential and the average value of the number of Dirac strings in the gauge configurations $m$. In the confined phase the mass gap translates into a linearly growing static potential, while in the deconfined phase the static potential approaches a constant at large separations. Results shows that $m$ is, at least, one order of magnitude larger in the confined phase and confirm that the appearance of a confined phase is connected with the topology of the gauge group.


I. INTRODUCTION AND MOTIVATION
The quest to understand quark confinement has, long ago, lead to the formulation of QED on a hypercubic lattice [1] by Wilson. He was able to show that the compact formulation of QED has two phases. In the strong coupling limit, i.e. at low β = 1/e 2 where e is the bare coupling constant, the theory is confining in the sense that the static potential between fermions grows linearly with the distance. In the weak coupling regime, i.e. at large β values, it was argued that the results of perturbation theory for the continuum formulation should be recovered. These predictions for the phase diagram of QED have been confirmed with numerical simulations and by various theoretical analyses of the theory [2][3][4][5][6][7][8][9][10][11][12][13]. In numerical calculations, they typically use a mapping into its dual formulation associated with the Z n symmetry group [14,15] instead of performing a sampling of the compact formulation of QED directly.
For compact QED the transition between the confining and the non-confining phase occurs at β = β c ≈ 1, see [7,16] and references therein. For β < β c the static potential V (R) is compatible with a Cornell type of potential with the dimensionless string tension σ being a decreasing function of β as its critical value is approached from below [12]. Numerical simulations for β > β c suggest that the theory contains a massless photon and reproduces the behaviour associated with a free field theory. * loverilc@piercecollege.edu † orlando@uc.pt ‡ psilva@uc.pt However, in the numerical simulations the reproduction of a free field theory for β > β c is never perfect [9,17,18].
Compact QED, having a confined and a deconfined phase, can also be viewed as a laboratory to try to understand the differences between the two phases. In 3D, the confining mechanism is associated with the presence of monopole configurations [19,20] at low β values; see also [21,22] for numerical simulations. However, in 4D these classical configurations or their equivalent are not able to generate a linearly rising potential at large separations and, therefore, cannot be at the origin of the confining mechanism. It has been argued that a Cornell potential type for the static potential is related to the topology of the gauge group, see e.g. [23] and references therein.
The interest in compact QED is not limited to its phase diagram. For example, it is not yet clear if the continuum limit of compact QED is a sensible theory. If the fermion sector of the theory is ignored, its continuum theory is a free field theory that one expects to recover for sufficiently high β values. The simulations performed so far suggest that this is the case, but full agreement with the results of a free field theory has not been achieved. The transition of the confined phase to the deconfined phase seems to be first order which, once more, raises the question of how to take the continuum limit for compact QED. Besides the question of the continuum limit, U (1) gauge theories are relevant to understand the Standard Model, the comprehension of the Higgs sector calls for simulations of Abelian Higgs models, and to describe many properties in Condensed Matter Physics. Furthermore, recently the lattice QCD community has started to compute QED corrections to QCD and, certainly, a good understanding of lattice Abelian models is of paramount importance. U (1) gauge theories are also a laboratory towards building simulations of gauge theories on quantum computers.
In the current work we are mainly interested in identifying the differences between the confined and deconfined phase by looking at the Landau gauge photon propagator in momentum space for the two phases of compact QED in the quenched approximation, i.e. for the compact U(1) theory. The simulations described in this work do not take into account the contribution of the fermionic degrees of freedom to the dynamics. For the confined phase the theory develops a mass scale and the photon propagator is finite and, at least for the simulations discussed, the presence of the mass scale prevents or reduces the finite volume effects. On the other hand, in the deconfined phase the photon propagator is compatible with a divergent 1/k 2 or a higher power of k behaviour at low momenta k. Finite volume effects are sizeable and the propagator approaches that of a free field theory as the lattice volume is increased. We hope that by studying the photon propagator one can arrive at a better understanding of the confinement mechanism for QCD.
The Landau gauge photon propagator for the confined phase qualitatively follows the same type of behaviour as the Landau gauge gluon propagator in QCD for pure Yang-Mills theories; see, for example, [24][25][26][27][28][29][30] for lattice simulations of the Landau gauge gluon propagator and [31][32][33][34][35][36][37][38] for continuum estimations of the same correlation function (see also the references therein). This suggests that the confinement mechanisms for the gluon and the photon have some similarities and, in particular, that the confining theory is associated with dynamically generated mass scales that make the propagator finite in the full momentum range for compact QED and for QCD. For QED a photon mass has been related to a meaningful physical value of the expectation value of the vector potential squared, that is connected with the existence of topological structures for the theory [13].
Our work also aims to fill the gap in the literature [18,39] and provide a large volume 4D lattice simulation of the Landau gauge photon propagator in momentum space for β below and above β c . The two phases are distinguished not only by the qualitative behaviour of the propagators but also by their topological structures and the static potential as computed from the Wilson loops. Our findings for the topological structures and static potential confirm previous results that can be found in the literature.
This work is organised as follows. In Sec. II we provide the definitions used in our calculation ranging from the Wilson action, the gauge fixing procedure, the definition of the electromagnetic potential, the propagator and number of Dirac strings crossing each plaquette. In Sec. III the propagator and the static potential are discussed for the confined phase, while in Sec. IV we discuss these properties for the deconfined phase. Finally, in Sec. V we summarise our results and discuss the differences between the two phases.

II. COMPACT QED: DEFINITIONS, LATTICE SETUP AND DETAILS OF THE SIMULATION
In the current work we will consider the compact version of QED defined on an hypercubic lattice that is described by the Wilson action which, in Euclidean space, is given by where β = 1/e 2 , with e being the bare coupling constant, and the plaquette operator is written in terms of the link fields with a being the lattice spacing,ê µ the unit vector along direction µ and A µ is the bare photon field. In the continuum limit the plaquette operator (3) can also be written where C is any closed curve around point x and infinitesimally close to it. On an hypercubic lattice, the term that appears in the exponential is the change of the photon field around a plaquette centered at x + a (ê µ +ê ν )/2 and we write It follows from the definitions that everywhere in the lattice −π e a A µ π and −π e a ∆A µν (x) π, i.e. the quantities e a A µ and e a ∆A µν take values on compact spaces. Note that, in general, e a ∆A µν is not given by the sum of e a A µ over each of the links that define the plaquette but, instead, where m µν (x) is an integer number that measures the number of Dirac strings that cross the plaquette 1 . The integer field m µν (x) can be measured by combining information on the links and plaquettes [40]. We will not study m µν (x) in detail but will report on its mean value over the lattice 2 , i.e.
where V is the number of lattice points. Indeed, as discussed below m can be used to distinguish between the confined and deconfined phases, with the configurations in the confined phase having a much larger mean number of Dirac strings crossing each plaquette. The generating functional of the compact QED Green's functions is given by where S W is defined in Eq. (2). For the importance sampling we rely on an implementation of the Hybrid Monte Carlo method [41] based on QDP++ and Chroma libraries [42]. The rotation of the links obtained by importance sampling to the Landau gauge is formulated as an optimization problem, over the gauge orbits. Setting the optimization problem depends on the definition of the photon field.
The gauge field can be computed with a linear definition that, for the U (1) theory, reads If one relies on this definition the gauge fixing is performed by maximizing the functional where V is the total number of lattice points and D the Euclidean spacetime dimension. It can be shown, see e.g. [43], that in this way the continuum Landau gauge condition is reproduced up to corrections O(a 2 ). For the U (1) gauge theory there is no clear way to set the lattice spacing and, therefore, this procedure can introduce significant deviations of the continuum Landau gauge when applied to the confined (β 1) or to the deconfined (β 1) phase. The setup just described, that uses a linear definition of the gauge field, is used for non-Abelian gauge theories defined on a lattice from simulations that are close to continuum physics. In a first step to define the Landau gauge for the photon field, we used the above procedure relying on a steepest descent algorithm with Fourier acceleration, see [43] for details and references, and controlling the approach to the Landau gauge with the quantity a lattice version of − ∂ · A(x). The maximisation was stopped when On the other hand, the Euclidean photon field can be computed using a logarithmic definition This is an exact definition, up to machine precision, that does not call for the use of a small lattice spacing. Then, following [44] adapted to the Abelian theory, the Landau gauge condition is achieved by maximizing the functional over the gauge orbits. In the Eq. (15) the field e a A (g) is the photon field given by Eq. (14) after the links U µ (x) have been gauge transformed by g(x). The approach towards the Landau gauge can be monitored using once more a lattice version of − ∂ · A(x). In our computation of the Landau gauge propagator, after the maximization problem associated with the functional given in Eq. (11), a maximization of the functional (15) is also performed. In this way one aims to reduce possible deviations of the continuum Landau gauge for both phases of the theory. In this second maximization we use again a steepest descent algorithm with Fourier acceleration and the gauge fixing was stopped when The definition of gauge field potential from the link variables, its implications on the choice of the optimization functional to define the Landau gauge on the lattice, how the different definitions impact the gluon correlation functions as been the subject of several studies, see e.g. [54][55][56][57][58][59] and references therein. These studies address the problem for non-Abelian gauge theories and explore essentially the linear definition given in Eq. (10) and its variants. For the non-Abelian gauge theories the different choices for the definition of the gauge potential lead to the same propagators. However, as already stated, the linear definition is possible due to asymptotic freedom of SU (N ) gauge theories. For a U (1) gauge theory on the lattice there is not a clear way of setting the physical scale and, therefore, to define the lattice spacing that would allow (or not) to the use for the linear definition of the photon field.
In both stages the maximisation of the gauge fixing functional is done with a Fourier accelerated steepest descent method that calls for the PFFT library [45] to do the required fast Fourier transformations. The complete numerical simulation, i.e. the importance sampling, the gauge fixing and the computation of all quantities, were performed in the Navigator cluster [46] of the University of Coimbra.
The use of the Hybrid Monte Carlo to perform the importance sampling can be questionable as, long ago [49,50], very long-living metastable states were observed in the Monte Carlo history of the system. In our study, the Monte Carlo time evolution of the mean value of the plaquette does not reveal metastable states. Indeed, at least for the two values of β considered, that are far from the transition between the two known phases, the Monte Carlo history for the mean value of the plaquette show a initial thermalization phase followed by a fluctuation pattern. No sign of metastable states are found and we have considered large Monte Carlo time histories. The time histories are larger that those mentioned in [49] but our lattice sizes are much larger than those considered there. However, the good agreement between all the simulations for each of the β values and the non-observation of metastable states gives confidence on the results of the simulation.
Another issue of concern connected with the computation of the photon propagator is the presence of Gribov copies. In the lattice formulation of the compact U (1) gauge theory, the Gribov copies are associated with different extrema of the maximizating functional that defines the lattice Landau gauge, see Eqs (11) and (15). The presence of the Gribov in the Landau gauge was also studied long ago, see e.g. [51][52][53], with observed measurable effects in the photon propagator. There are two important differences between our implementation of the Landau gauge and that used in the those earlier studies. The first one is the precision on the definition of the lattice version of ∂A = 0 that differs by several orders of magnitude. It well known that from the analysis of the non-Abelian theories that a small value of θ of 10 −12 or smaller should be used to have a good implementation of the Landau gauge on the lattice. The relaxation of θ towards higher values affects mainly the low momentum region propagator 3 . The second difference is related to the definition of the gauge field. The previous studies use the linear definition for the gauge field (10), while we use the logarithmic definition (14). This does not eliminate the problem of the Gribov copies but it is not obvious that the observations reported in [51][52][53] apply in our formulation. The studies of pure gauge non-Abelian gauge theories, see e.g. [43], show that when θ ∼ 10 −12 or smaller and for statistics with similar number of gauge configurations, typically, the effects of Gribov copies are diluted within the statistical error. We hope the same applies to compact U (1) theory and, in the following, we will assume that this is the case.
From the definition (14) for the Euclidean spacetime photon field, the momentum space photon field is given by and the Landau gauge propagator reads where · · · stands for the vacuum expectation value. In a lattice simulation, the vacuum expectation values are accessed via the generation of a set of configurations sampled accordingly with the probability distribution (9) and taking averages of the products of gauge fields, such as those in Eq. (19), over the full set of gauge configurations. For the analysis of the propagator, it will be assumed that the propagator has the same tensor structure as the continuum theory, i.e.
where the function D(p), named propagator below, is a function of the tree level improved momentâ p = 2 a sin π L n µ , where L is the number of lattice points in each side of the hypercubic lattice. The rationale to usep instead of the naive lattice momenta comes from lattice perturbation theory that requiresp instead of p. In the lattice evaluation of the gluon propagator the improved momentum also helps to suppress finite spacing effects in the propagator [47]. In order to further suppress the effects due to the use of finite lattice spacing we perform the conical and cylindrical cuts introduced in [47] for momenta ap > Λ IR and, following the procedure devised in [24], below this threshold we consider all the momenta available to get information on the infrared region. The choice of the infrared threshold Λ IR is a compromise between taking into account extra data, allowing larger fluctuations, and resulting in a smooth curve for D(p 2 ). The choice of this threshold does not change the overall behaviour of the lattice data and Λ IR will be chosen differently for each simulation. In the following we use Λ IR = 0.2 for the largest lattice volume and Λ IR = 0.4 for the two smallest lattices. The description of the lattice propagator with the continuum tensor structure as given by Eq. (20) is questionable, especially concerning the confined phase. Similar studies for the gluon propagator show that the lattice propagator has other tensor structures not considered in Eq. (20) and, in principle, they should also be considered here. However, given that the definition of the lattice Landau gauge returns a transverse gauge field, one expects a gauge propagator that should also be transverse. Furthermore, the studies performed for the gluon propagator suggest that the introduction of momentum cuts selects the set of momenta where the finite lattice effects are minimised. This gives us confidence that the same should apply to the photon propagator.
If one assumes a tensor structure as given by Eq. (20), then the type of momentum considered in the projector is irrelevant as long as one measures the propagator form factor using We remind the reader that for zero momentum the propagator is given by δ µν D(0) and, therefore, the computation of D(0) requires a different normalisation factor. In the current work, we aim to see how the photon propagator behaves in the confined and deconfined phases. To achieve such a goal we perform Monte Carlo simulations of the theory at β = 0.8 (confined phase) and at β = 1.2 (deconfined phase). In order to check for finite volume effects in both cases we perform simulations on 32 4 , 48 4 and 96 4 hypercubic lattices. For each β value and lattice volume, the propagators were computed using the last (in the Markov chain) 200 gauge configurations. The configurations used in the calculation of propagator have a separation of 10 trajectories for the smaller lattice volume, for both β values considered herein, and also for the 48 4 simulation in the confined phase (β = 0.8). In the remaining simulations we used a separation of 100 trajectories in the corresponding Markov chain.
In Fig. 1 the mean values of the plaquette over the lattice are shown for each of the Markov chains. In the simulations the value of the plaquette at the end of each trajectory was not always kept and, therefore, in the reconstruction of the plaquette history we lost some of the data points. As the Fig. shows the mean value of the plaquette seems to be independent of the lattice volume for each β and in the deconfined region, i.e. for the simulation with β = 1.2, the plaquette is significantly larger. This result suggests that the U (1) links approach unity as β is increased. For the computation of statistical errors for all the quantities reported here, i.e. propagators, Wilson loops and monopole densities, we rely on the bootstrap method with a 67.5% confidence level. The quoted errors associated with the fits assume Gaussian error propagation.

III. PHOTON IN THE CONFINED PHASE
The Landau gauge photon propagator for compact QED in the confined phase with β = 0.8 and for the various lattice volumes can be seen in Fig. 2. The data does not follow the behaviour of a free particle propagator and deviations from a 1/p 2 functional form are clearly seen. Indeed, the various data sets seem to be closer to the qualitatively behaviour of the QCD gluon propagator [24,25,27]. Moreover, the propagator being finite over the full range of momenta suggests that 4D com- pact QED generates a mass gap dynamically, as is also observed in 3D simulations [21,22]. The data for the various volumes is compatible within one standard deviation and, therefore, shows no volume dependence.
It seems that the presence of the mass gap is sufficient to reduce the volume dependence of D(p 2 ). This contrasts with what is observed for the propagator in the deconfined phase; more on this topic later.
A possible way to identify the mass gap is by fitting the lattice data to a given functional form. We found that, for all volumes, the lattice data is well described by a Yukawa type propagator where z 0 , p 2 and m 2 are dimensionless quantities.  Fig. 2 the solid black line represents the functional form given in Eq. (24) with z 0 and m 2 given by the estimation of the fit to the lattice data from the largest volume. Similar curves using the other two sets of parameters could be drawn but the curves are indistinguishable to the naked eye from the curve shown. The low β phase of compact QED was investigated by Wilson in [1], where he computed the static potential from Wilson loops. Indeed, it was shown that, at low β, the Wilson loop follows an area law and, therefore, the as-  sociated static potential grows linearly with the distance between sources. It is in this sense that compact QED at low β values is a confining theory. This observation motivated us to compute Wilson loops and we took only those loops whose spatial part is along one of the lattice axis to measure the static potential V (R). The Wilson loop can be seen in Fig. 3. Note that we use no trick to improve the signal to noise ratio, the noise for W (R, T ) is large for some cases and increases with R. This is an indication that the static potential grows with R given that Further, it is clear from Fig. 3 that exponential behaviour sets in for quite small T . Then, from the data for and in Fig. 4 we show V (R) computed from taking T = 2 and T = 3. The data in Fig. 4 should be regarded as an upper bound on V (R). The results summarised in Figs. 3 and 4 confirms that V (R) grows with R and suggest that the data is compatible with linear behaviour at large R.
In this sense the simulation confirms that compact QED is a confining theory at low β values. The static potential for 4D compact QED was computed in [12,60,61] from Polyakov loops, exploring duality transformations, and it was found that in the confined phase V (R) grows linearly with the distance for sufficiently large R as also found in our simulations.

IV. PHOTON IN THE DECONFINED PHASE
The nature of the photon propagator at large β is expected to be rather different than that observed in Fig. 2. Indeed, as can be seen in Fig. 5, for the deconfined phase with β = 1.2 the photon propagator seems to diverge at zero momentum. Furthermore, if at low β the propagator is blind to the finite volume effects, the data for the various volumes in Fig. 5 are not compatible with each other within one standard deviation. The propagator data for the smallest volume 32 4 is above the other two sets of propagator data in the mid range momenta, while the data associated with the 48 4 lattice is between the data computed with the smallest and the largest lattice volumes. However, at zero momenta the largest a 2 D(0) is associated with the largest volume, followed by the 48 4 data and by the 32 4 data in decreasing order of values. For momenta such that a p 2 all the data sets seems to be compatible within one standard deviation, see the inner plot in Fig. 5.
The data in Fig. 5 suggest that the photon propagator diverges as momentum approaches zero. If the data is to be associated with a free field theory, it should reproduce the behaviour of a free field propagator. However, in a finite volume Monte Carlo simulation deviations from the continuum free field theory are expected as the simulation is performed on a finite lattice. The approach to the continuum behaviour can be tested by fitting the lattice data to the functional form If the theory reproduces a free field theory a Z 1 = 0 is a manifestation of finite volume effects and one expects Z 1 to become smaller as the lattice volume is increased. The direct fit using the full range of momenta and taking into account the statistical errors of the lattice data returns values of the χ 2 /d.o.f. 12. For the smallest volume, for ap 0.6 the fit has a χ 2 /d.o.f. = 1.74 with Z 0 = 2.425 (14) and Z 1 = 0.986(46). For the 48 4 data and for ap 0.4 it follows that χ 2 /d.o.f. = 1.01 with Z 0 = 2.4057(76) and Z 1 = 0.355 (13). On the other hand for the largest lattice volume, due to large fluctuations that are observed at larger momenta, one can never achieve a reasonable χ 2 /d.o.f. However, by doubling the statistical errors on the definition of the χ 2 , the data becomes compatible with (27) for ap 0.15. In this case the fit has χ 2 /d.o.f. = 1.93 with Z 0 = 2.374(13) and Z 1 = 0.2958(68). Note that in all cases one has Z 0 ≈ 2.4, while Z 1 decreases with the lattice volume. 4 .
The lattice data together with the fits can be seen in Fig. 6. In general and for the corresponding fitting ranges, the curves overlap the Monte Carlo data. Further, the coefficient Z 0 ≈ 2.374 (13), we are quoting the value of the fit to the largest lattice volume, seems to be nearly independent of the volume. The data in Fig. 6 suggest that Z 0 is independent of L, while Z 1 is sensitive to L. Indeed this coefficient goes from Z 1 = 0.986(46) for the smallest volume to Z 1 = 0.2958(68) for the largest volume, which is about 1/3 of the value for the smallest volume; note that the inverse of the ratio of the lattice sizes is 1/3. This result for Z 1 suggests that the data for the propagator seems to converge to the propagator of a free field theory in the infinite volume limit. This statement has to be read with care due to the use of 2 σ in the definition of the minimising χ 2 for the largest volume. That the fitting range does not start at the smallest non vanishing momentum for each volume is not surprising, as finite volume effects, that should appear at the smallest momenta, are to be expected. We have also tried fitting the data with an almost free field propagator, i.e. assuming D(p 2 ) = Z 0 /(p 2 ) α and leaving Z 0 and α as free parameters. The fits to this last functional form have the same problems as those mentioned before but it turns out that α ≈ 1, i.e. the lattice data for the propagator follows closely the behaviour of a free field theory 5 .
The above analysis suggests that the Monte Carlo propagator data almost reproduce a free field theory propagator. Let us check the results for the static potential, computed from Wilson loops as was done for the confined phase. The Wilson loop for various values of R 4 The data for Z 0 and Z 1 in Fig. 6 suggest that Z 0 approaches its infinite volume limit linearly. A linear extrapolation towardsL → ∞ results in Z 0 = 2.352 (10). On the other hand, Z 1 clearly does not follow a linear behavior towards the infinite volume limit. Given that our simulations only considered three different volumes it is not possible to estimate Z 1 (L → ∞). 5 In order to quantify the typical values of α let us report on its values given by fitting to the propagator data replacing σ by 2σ in the definition of the minimising χ 2 . Demanding that the χ 2 /d.o.f. 2, it follows that for the smallest lattice volume the fitting range starts at ap = 0.6 and has α = 1.095/10), the fitting range for the 48 4 data starts at ap = 0.4 and has α = 1.0529 (60), while for the largest volume the fitting range starts at ap = 1 and has α = 0.9977(74).  is given in Fig. 7 for the largest lattice volume and it looks rather different from the Wilson loop for the confined phase reported in Fig. 3. If for the confined phase the slope increases with R, for the deconfined phase the slope of the log W (R, T ) seems to be the same for all R. Indeed, measuring V (R) from the effective mass and taking its value for T = 9, one gets the bottom plot of Fig. 7. The large distance behaviour of V (R) is sensitive to finite volume effects and the slope of V (R) for large R becomes smaller as L is increased. In Fig. 8 we show the string tension measured by fitting V (R) in the range R = 9−12 for the various volumes. The corresponding χ 2 /d.o.f. for the various fits are always below 0.5. The dashed blue line connects the origin with the result for the largest volume, while the shaded region takes into account the one standard deviation on σa 2 for L = 96. Our results seems to be compatible with a vanishing string tension in the infinite volume limit. The short distance behaviour of V (R) is difficult to understand from the computed Wilson loop directly. The Monte Carlo data for the photon propagator is compatible with free field propagator behaviour at high p and approaching 1/p 2 as the volume is increased and, therefore, one expects to have, in the infinite volume limit, V (R) ∝ 1/R at short distances. We have tried to disentangle the short distance behaviour from the V (R) Monte Carlo data but the results were inconclusive.
Our simulations for the deconfined phase of compact QED suggest that for the β considered here, the finite volume effects are still not negligible even for a lattice volume as large as 96 4 .

V. SUMMARY AND CONCLUSION
In the current work the Landau gauge photon propagator is investigated for compact QED in the strong coupling (confining) and weak coupling (free field theory) regimes and for various lattice volumes. By computing the static potential, our simulation confirms that at low β the theory is confining and the behaviour of the photon propagator in momentum space follows closely a Yukawa type of propagator, i.e. 4D compact QED has a mass gap. Moreover, in the confining phase the theory develops a mass scale that makes the photon propagator finite in the full momentum range.
For the deconfined phase, the photon propagator seems to approach a free field type of propagator as the infinite volume is approached. We have observed that the matching with a free field theory is not perfect with both the photon propagator and the static potential showing some deviations from the expected behaviour, that we interpreted as being due to finite volume effects. Indeed, the deviations from a free field theory are reduced, in all the computed quantities, as the lattice volume is increased. Given that at low momenta the propagator of a free field theory diverges and the lattice regularizes both the UV and IR divergences, in a sense the deviations from the free field theory results are not unexpected.
Comparing the confined and deconfined phase results, it seems that it is the generation of a mass gap that occurs for the confined phase that turns the propagator essentially independent of the lattice volume. This is a situation that is also seen in the simulations for QCD.
If the phase diagram for compact QED as a function of β has two different phases, one can ask how can one distinguish them. According to [19] the appearance of the confining phase for the 3D theory is due to the presence of monopole configurations. The monopoles are connected with the topology of the gauge group and they becomes irrelevant for the dynamics at large β values. The 4D equivalent to the monopole configurations are Dirac strings that should be seen on a finer analysis of the gauge configurations. In Fig. 9 we report on the average number of Dirac strings over the lattice, computed with the definitions (7) and (8), for all the lattices. The plots show that m is independent of L for the confined phase (β = 0.8), with m having larger fluctuations in the smaller lattices. For the deconfined phase (β = 1.2) m decreases with the lattice volume and the results suggest, once more, that in the simulations performed in the deconfined phase the infinite volume limit is not yet realised, despite the large volume considered. Furthermore, m is about a factor of fifty larger in the confined phase when compared to its value in the deconfined phase. This result suggests that, indeed, the Dirac strings are responsible for the confined phase in 4D compact QED in agreement with the suggestion of [19]. However, further studies are required to draw firm conclusions.

ACKNOWLEDGMENTS
This work was partly supported by the FCT -Fundação para a Ciência e a Tecnologia, I.P., under Projects Nos. UIDB/04564/2020 and UIDP/04564/2020. P. J. S. acknowledges financial support from FCT (Portugal) under Contract No. CEECIND/00488/2017. The authors acknowledge the Laboratory for Advanced Computing at the University of Coimbra (http://www.uc.pt/lca) for providing access to the HPC resource Navigator.