Gauge invariant determination of charged hadron masses

In this paper we show, for the first time, that charged-hadron masses can be calculated on the lattice without relying on gauge fixing at any stage of the calculations. In our simulations we follow a recent proposal and formulate full QCD+QED on a finite volume, without spoiling locality, by imposing C-periodic boundary conditions in the spatial directions. Electrically charged states are interpolated with a class of operators, originally suggested by Dirac and built as functionals of the photon field, that are invariant under local gauge transformations. We show that the quality of the numerical signal of charged-hadron masses is the same as in the neutral sector and that charged-neutral mass splittings can be calculated with satisfactory accuracy in this setup. We also discuss how to describe states of charged hadrons with real photons in a fully gauge-invariant way by providing a first evidence that the proposed strategy can be numerically viable.


Introduction
QED radiative corrections to hadronic observables are generally rather small but they become phenomenologically relevant when the target precision is at the percent level. For example, hadron masses and leptonic decay rates of light pseudoscalar mesons are among the best measured hadronic observables and they have to be calculated at the same level of precision. Presently, these quantities can be calculated with percent accuracy by performing lattice simulations of QCD+QED, see e.g. refs. [1][2][3][4][5][6] for a selection of recent papers on the subject and refs. [7][8][9] for recent reviews. All these calculations have been performed by using non-compact gauge-fixed lattice formulations of QED in a finite box, see ref. [9].
In ref. [10] it has been argued that charged-hadron masses can be calculated on the lattice from first principles, in a completely gauge-invariant setup, without spoiling basic QFT principles in finite volume, in particular locality. This result is far from obvious. The construction is possible thanks to two crucial ingredients: a slightly unconventional compact formulation of lattice QED, and properly chosen boundary conditions in the spatial directions.
In a gauge theory physical states are invariant under local gauge transformations. Therefore, in order to avoid gauge fixing, physical states have to be probed by using interpolating operators that are invariant under local gauge transformations. Building these operators is trivial in the neutral sector of the theory. For example, in order to compute the mass of a neutral kaon one can usesγ 5 d as the interpolating operator. Since the down and strange quarks have the same electric charge the operator is electrically neutral and invariant under both local and global U(1) gauge transformations.
Remarkably, in infinite volume, one can build interpolating operators that are invariant under local gauge transformations also in the charged sector of the theory. The existence of these operators was first pointed out by Dirac in an illuminating classic paper [11] (see section 2). In principle, Dirac's interpolating operators can be used to calculate observables associated with charged particles, e.g. the mass of the electron or of a charged kaon, in a fully gauge-invariant way. In practice, in order to obtain a fully gauge-invariant formulation of QCD+QED one has to provide a regularisation of the theory where Dirac's construction can be implemented without any theoretical ambiguity.
Dirac's construction cannot be implemented on the periodic torus. In operatorial formalism, the generator of local gauge transformations is ∂ k E k − j 0 , where E k is the electric field and j 0 is the charge density, such that Q = d 3 x j 0 . Identifying physical states, |Ψ , with gauge-invariant states is equivalent to requiring that physical states must satisfy the Gauss law. In particular this implies that Q|Ψ = d 3 x ∂ k E k |Ψ . Therefore, with periodic boundary conditions in space the global constraint imposed by the Gauss law forbids states with non-zero charge. Equivalently, no interpolating operator exists on the periodic torus which is electrically charged and invariant under local gauge transformations.
In ref. [10] it has been proposed to discretise QCD+QED on a finite lattice by using the compact formulation and, as first suggested in refs. [12][13][14][15], with C-periodic (or C ) boundary conditions in space. A detailed theoretical analysis of the theory, called QCD+QED C , has shown that the Gauss law implies a less restrictive global constraint in this case. Some electrically charged states can be probed by implementing Dirac's original construction in a fully consistent theoretical setup (see section 3), i.e. by using charged interpolating operators which are invariant under local gauge transformations.
While the theoretical analysis of ref. [10] opens the attractive possibility to perform firstprinciples non-perturbative lattice simulations of QCD+QED in a fully gauge-invariant setup, no evidence was provided concerning the numerical viability of the proposal 1 . In this paper we make a first step in the direction of filling this gap. We provide clear numerical evidence that charged-hadron masses can be effectively calculated in QCD+QED C from the gauge-invariant interpolating operators with the same signal-to-noise ratio as their neutral almost-degenerate counterparts. We also discuss how to describe states of charged hadrons with real photons in a fully gauge-invariant way. On the other hand, the cost of the generation of configurations will be analysed in future work.
The paper is organised as follows. In section 2 we review Dirac's original construction of gauge-invariant interpolating operators for charged states. In section 3 we recall the finite-volume formulation of QCD+QED with C boundary conditions of ref. [10] and the lattice construction of gauge-invariant electrically-charged operators. In section 4 we present our numerical results for charged and neutral meson masses both in the vector and pseudoscalar channel. In particular, in subsection 4.1 we discuss the implementation of a strategy to probe charged-hadron states with real photons. We draw our conclusions in section 5. Finally, in appendix A we discuss some of the subtleties arising in the charged sector when the U(1) gauge is fixed, and in appendix B we provide some technical details concerning the numerical evaluation of the correlators used in this study. 2 Dirac's interpolating operator Dirac [11] has shown that charged states in infinite-volume QED can be described in a fully gaugeinvariant setup in terms of physical degrees of freedom. In Dirac's original construction the state of an electron can be interpolated by means of the operator where k is a spatial index, A µ (x) and ψ e (x) are the photon and electron fields while Φ(x) is the electrostatic potential satisfying Under a gauge transformation λ(x) the fundamental fields transform as If λ(x) has compact support, the integral appearing in the definition of Ψ c e (x) transforms as The operator Ψ c e (x) is invariant under local gauge transformations, but transforms non-trivially under global gauge transformations. When acting on the vacuum, Ψ c e (x) generates a physical state (i.e. invariant under local gauge transformations) with total charge different from zero.
An important observation concerning this construction is that in Coulomb gauge ∂ k A k (x) = 0 the interpolating operator is identically equal to ψ e (x). On the one hand, this means that Dirac's construction can be circumvented and that the mass of the electron can be calculated in Coulomb gauge by using ψ e (x) as interpolating operator. This is presumably the reason why Dirac's paper went almost forgotten. On the other hand, Dirac's construction explains why gauge-invariant physical quantities can be conveniently extracted by working at fixed gauge.
The gauge-invariant language is very useful in order to identify and clarify some of the subtleties arising with commonly used gauge-fixing conditions. For instance the Landau-gauge elementary field ψ e (x) is identical to the following generalisation of Dirac's original operator, This implies for the two-point function Since Ψ e (x) is non-local in time, a standard interpretation as an interpolating operator is not possible. The phase in eq. (2.5) should rather be viewed as a term in the action. Since the term is linear in the electromagnetic field, this is in fact the coupling to a non-real external electromagnetic current. This mechanism is quite general. As discussed in appendix A, gauge fixing introduces (except special cases, of which Coulomb gauge is the most notable one) a violation of the Gauss law in the sector of non-zero charge, which can be interpreted as the effect of coupling the physical system to an external electromagnetic four-current. This current, and consequently the Hamiltonian, is generally time dependent. In Euclidean spacetime, as an effect of the Wick rotation, the external charge density is real while the external current density is imaginary and the Hamiltonian turns out to be non-hermitean. This implies that a spectral decomposition of two-point functions as a sum of exponentials of the form n a n exp(−tE n ) is simply incorrect. For reasonable enough gauges (e.g. covariant gauges) the external four-current vanishes asymptotically far away from the interpolating fields in the two-point function, and the long-distance behaviour of the two-point function is dictated by the ground state of the physical Hamiltonian, i.e. in absence of the external four-current. However, in a setup in which observables are not expanded in powers of α em , it is not obvious at all how to extract excited physical states, such as the finite-volume counterparts of states of charged hadrons with real photons. A gauge-invariant construction of n-point functions becomes of utmost relevance precisely when excited states are of interest. Because the gauge-invariant Hamiltonian is hermitean and time independent, standard spectral theory applies, and gauge invariance ensures that only physical states (i.e. states that satisfy the Gauss law) propagate at any intermediate time.
In ref.
[10] Dirac's construction has been used to provide a theoretically consistent definition of electrically charged states in a finite volume within the framework of local field theory, as we will review in the next section.

Charged states in finite volume
The formulation of QCD+QED C has been discussed in ref. [10] together with a detailed analysis of its symmetries and an analytical calculation of the leading finite volume effects on the masses of charged hadrons. Here, in order to make the paper self-contained, we briefly discuss the compact lattice formulation of the theory.
Gauge degrees of freedom are encoded in the link variables U µ (x) ∈ U(1) and V µ (x) ∈ SU(3). All the fields obey C boundary conditions along the spatial directions, namely where ψ f are the quark fields, f is the flavour index and C is the charge-conjugation matrix 2 . We have simulated the theory by imposing periodic boundary conditions in time.
The spatial boundary conditions for the gauge fields are imposed in a completely straightforward way. However, since C boundary conditions mix ψ andψ, the Dirac operator D f cannot be defined as an operator acting on the space of the fields ψ only, but it has to be thought as an operator acting on the quark-antiquark doublet which satisfies the following boundary condition where the Pauli matrix σ 1 acts on the quark-antiquark components. An explicit expression for the Dirac operator D f will be given at the end of this section. Once the fermions are integrated out, the lattice-discretised path-integral measure turns out to be where S g and S γ are the SU(3) and U(1) gauge actions respectively, Pf denotes the Pfaffian, and we choose to have three dynamical quarks for definiteness. The Pfaffian is proven to be real at finite lattice spacing, and positive in the continuum limit (see appendix D in [10]). The probability to find a negative value is expected to be negligible in our simulations with fairly heavy quarks. Therefore, we have simulated the absolute value of the Pfaffian, and monitored that the lowest eigenvalue stays significantly away from zero. For the SU(3) gauge action S g we use the Lüscher-Weisz discretisation [17], while the U(1) gauge action is defined as where e 0 is the bare electric charge of the positron and U µν (x) is the U(1) gauge plaquette, i.e.
The point to be noticed in previous formulae is the unconventional normalisation of the U(1) gauge action, namely the factor 18/e 2 0 instead of 1/2e 2 0 . The canonically-normalised continuum action is obtained by setting To be consistent with this normalisation, the covariant derivatives acting on the quark fields are defined with the 6q f -power of the U(1) gauge links, where q f is the charge of ψ f in units of e 0 . For example the forward covariant derivative acting on the flavour f is given by The peculiar normalisation of S γ is due to the fact that quarks have fractional electric charges, q u,c,t = 2/3 and q d,s,b = −1/3, and to the fact that with this choice Dirac's interpolating operators can be discretised using analytical functions of the link variables. In the lattice formulation one can choose The k-th term in the sum above is the unique U(1) gauge-invariant extension of the quark field in axial gauge 3 U k (x) = 1. The corresponding expression in the finite-volume continuum theory is A discretization of Dirac's original interpolating operator, i.e. the one corresponding to Coulomb gauge, can be obtained by considering where ∇ k and ∇ * k are the free forward and backward lattice derivatives,∇ k = (∇ k + ∇ * k )/2, ∆ = ∇ k ∇ * k , and F µν is a discretisation of the U(1) field tensor. In this work we have used the standard clover discretisation for the field tensor. Notice that A c µ is a gauge-invariant discretisation of the photon field in Coulomb gauge,∇ k A c k = 0. In the formal continuum limit where Φ(x) is the unique electrostatic potential on the finite volume with antiperiodic boundary conditions. Therefore, is a consistent discretisation of Dirac's interpolating operator. In our numerical calculations, we used both the string operator Ψ s f and the Coulomb operator Ψ c f . Fully gauge-invariant interpolating operators for charged hadrons can be obtained by starting from the usual expressions, e.g.sγ 5 u, and by replacing the quark fields with the chosen Dirac's interpolating operator, e.g.S c γ 5 U c .
Before closing this section we give the explicit expression of the O(a)-improved Wilson-Dirac operator used in our simulations 14) The forward derivative acts on the quark-antiquark doublet η f as 15) and is defined at the boundary by means of the relation (3.3). The backward derivative ∇ f * µ is defined analogously. G µν and F µν are the clover discretisations of the SU(3) and U(1) field tensors respectively, and σ µν = i[γ µ , γ ν ]/2. The field tensors are normalised in such a way that tree-level improvement is achieved by choosing c QCD sw,f = c QED sw,f = 1.

Numerical explorations
In this section we discuss some exploratory simulations of QCD+QED with C boundary conditions. The main goal of this study is to show that the masses of charged mesons can be extracted in a completely gauge invariant way, with the same quality of the numerical signal as for neutral mesons. A preliminary calculation of excited states that would correspond to states of charged mesons with one real photon at α em = 0 is also shortly presented. The simulations have been performed by using a modified version of the HiRep code [18] (see ref. [19] for more details concerning the implementation) and we have checked our results by performing dedicated runs with the publicly-available openQ*D code [20] developed independently within the RC collaboration (see ref. [21]). While the HiRep code has been preferred in this exploratory work because of its simplicity, the optimized openQ*D code is currently used by the RC collaboration to perform realistic QCD+QED simulations.
We have generated two SU(3)×U(1) ensembles which differ only for the electromagnetic coupling, one with α em = 1/137 and one with α em = 0.05 = 6.85/137. The lattice is 48 × 24 3 with periodic boundary conditions in time and C boundary conditions in all spatial directions. The Lüscher-Weisz action and the action in eq. (3.5) have been used for the SU(3) and U(1) gauge fields respectively. Three dynamical Wilson fermions with Dirac operator given in eq. (3.14) have been simulated, one up-type quark with charge q = 2/3 and two down-type quarks with q = −1/3. The QCD bare parameters have been taken from one of the N f = 2 + 1 CLS ensembles at the symmetric point, i.e. the H200 ensemble in ref. [22] with β = 3.55, κ = 0.137, c QCD sw, = 1.824865, and complemented with the tree-level value c QED sw, = 1. The values in physical units of the lattice spacing and of the pseudoscalar meson masses are given in table 1.
In order to obtain a similar physics in the QCD and QCD+QED ensembles, the bare parameters would need to be retuned. In particular the bare masses of the up and down quarks should be retuned separately. However for sake of simplicity, in these exploratory simulations we chose to keep  [22]. The three ensembles share the same value of β = 3.55 and κu = κ d = κs = 0.137. For our simulations at αem = 0 we use qu = 2/3 and q d = qs = −1/3. The lattice spacing has been estimated by rescaling the CLS value with our measured a/ √ t0, and the error is estimated to be of order 10 −3 fm. The error on our pseudoscalar masses is estimated to be of order 15 MeV.
the bare parameters fixed and to measure the QED effects on the physical quantities. In particular we observe that QED corrections on the lattice spacing are fairly small even at the larger value of α em . The effect on the critical bare mass is general larger, as expected since this is an ultraviolet divergent quantity. Nevertheless we observe that in our ensemble with α em = 1/137 the pseudoscalar mesons have reasonable masses, of the order of the physical kaon mass.
Our simulations use a volume that is smaller than the original CLS ensemble. This is potentially an issue since masses in QCD+QED have finite volume corrections that decay as inverse powers of L rather than exponentially. An estimate of the finite-volume effects can be obtained by calculating the universal 1/L and 1/L 2 corrections (see sec. 5 in [10]), which turn out to be well below 1% for both values of α em .

Charged and neutral mesons
With C boundary conditions the eigenstates of the momentum are also eigenstates of charge conjugation. In particular zero-momentum states are also even under charge conjugation. The boundary conditions break the U(1) global gauge symmetry down to its Z 2 subgroup. As a consequence, if Q is the electric charge operator, then Q is not conserved, but (−1) Q is. When we talk about neutral states we really talk about states with (−1) Q = +1, and when we talk about charged states we really talk about states with (−1) Q = −1.
We consider the following C-even, zero-momentum, neutral interpolating operators and the following C-even, zero-momentum, charged interpolating operators where the non-local operatorsS I and U I are constructed as in eqs.   The effective masses are shown in fig. 1 for the P states and in fig. 2 for the V states for both values of α em . For comparison, in fig. 1 we also report the effective mass calculated on a QCD-only ensemble, generated with the CLS H200 bare parameters on a 48 × 24 3 lattice with C boundary conditions. In all cases we have used 500 configurations and 8 stochastic sources per configuration. As expected, we observe that the pseudoscalar masses are larger with respect to the ones quoted in ref. [22] for the H200 ensemble, because of the mass shift due to the electromagnetic interactions.
The most important result of this paper is the fact that effective masses can be extracted with similar errors in the neutral and charged channels. In fact, the introduction of the non-local gauge-invariant operators for charged states does not affect much the quality of the signal in correlators and effective masses. We also observe that in these channels, the string and Coulomb operators behave very similarly. Moreover, we observe that the statistical errors in the QCD+QED pseudoscalar effective mass are very similar to their QCD-only counterparts.
While these simulations are performed at unphysical values of the quark masses, the chargedneutral mass splittings can clearly be extracted with a statistically significant accuracy for both the pseudoscalar and vector states at α em = 0.05. Remarkably, the mass splitting in the pseudoscalar channel is statistically significant even at α em = 1/137.

Charged mesons with real photons
The goal of this subsection is to sketch a strategy to extract states of charged mesons with real photons. Let us focus on the charged vector channel. In finite volume, the spectral decomposition can   In these formulae we assume that the T → ∞ limit has been taken already. Since the full QCD+QED Hamiltonian does not conserve the photon number, the states |n, r are not eigenstates of the photon number operator. However at the leading order in α 1/2 em , the state |n, r is nothing but the tensor product of a QCD state with r free real photons, and its energy is given by the energy of the QCD state plus the energy of the free photons. Therefore it makes sense to refer to |n, r as a state with r real photons, as long as α em is small enough. Notice that these states are gauge invariant by construction, therefore only physical polarizations of the photon contribute.
If the volume is large enough, the ground state of the C I V (t) correlator is a state with one real photon. At the leading order in α 1/2 em this state contains a charged P particle and a real photon in a kinematic configuration with zero momentum and zero angular momentum. Some tedious but standard group theory reveals that, in order to be able to construct a state in the vector (T − 1 ) representation of the cubic group O h , the minimum-norm momentum allowed for the photon is up to isometries of the cube 4 . This state has energy equal to 9) and is created at the leading order in α 1/2 em by the following interpolating operator whereP I andÃ c are defined as and A c k is the gauge-invariant representation of the Coulomb-gauge photon field defined in eq. (3.11). Notice that the operatorP I (t, p) is C-odd, contrarily to the analogous operator defined in the previous subsection. This is due to the fact that states with momentump are antiperiodic, i.e. they are odd under a translation by a distance L in any of the spatial directions, and therefore odd under charge conjugation.
If the volume is large enough and α em is small enough, then the inequality E 0,1 < E 0,0 comes from the observation that M P is always smaller than M V , and in particular this is true at α em = 0. However as the volume goes to zero, the relative momentum of the two particles in the |0, 1 state diverge and so does E 0,1 . Therefore, if the volume is small enough, then E 0,1 > E 0,0 . It will turn out that this is the kinematic region of our simulations.
One can set up a generalised eigenvalue problem with two operators: V I k and W I k . If α em is small enough, V I k has maximal overlap with the state |0, 0 and W I k has maximal overlap with the state |0, 1 . At moderate value of α em , or in the regime in which the P+γ state is almost degenerate with a P+P state (which is in fact the case in our simulations), a larger operator basis may be necessary. In this exploratory calculation we will ignore these subtleties and proceed with the simple two-operator setup. If C I (t) is the 2 × 2 matrix of correlators constructed with the operators V I k and W I k , we solve the generalised eigenvalue problem given by We have extracted the ground state, λ I 0 (t, t 0 ), and the excited state, λ I 1 (t, t 0 ), eigenvalues by using both the string and Coulomb interpolating operators by obtaining statistically consistent results with essentially the same quality of the signal-to-noise ratio. In fig. 3 we plot the effective masses extracted from λ n (t, t 0 ) = λ s n (t, t 0 ) + λ c n (t, t 0 ) 2 (4.14) for n = 0, 1, corresponding to α em = 1/137 and α em = 0.05 respectively. The presented results are obtained with t 0 = 8, but we have checked the stability of our results in the range t 0 ∈ [4, 10].
On the one hand, from a quantitative analysis of the excited-state energy it turns out that (as anticipated in the discussion above) we cannot discriminate between a P+γ and a P+P state within the present statistical uncertainties. Since this may be due to the unphysical values of the bare parameters used in this study, we postpone a more detailed numerical analysis to future work on this subject. This will certainly require more statistics and, possibly, an extended basis of interpolating operators.
On the other hand, some qualitative information can be drawn from the plots in fig. 3. In our opinion, the quality of the numerical signals makes us pretty confident of the possibility to probe charged states containing real photons by using a fully non-perturbative gauge-invariant strategy along the lines of the one sketched in this section.

Conclusions
We have performed numerical lattice simulations of the compact formulation of QCD+QED with C-periodic boundary conditions in the spatial directions. In this setup, following ref. [10], chargedhadron masses can be calculated from first principles without relying on gauge fixing at any stage of the calculation.
Our simulations are performed at unphysical values of the bare parameters, with pseudoscalar meson masses of the order of the physical kaon at α em = 1/137. For this reason our results do not have phenomenological relevance but do have, in our opinion, deep theoretical implications. We provide a clear evidence that the strategy of ref. [10] is numerically viable and that charged states can be efficiently probed in a gauge-invariant way.
In particular, we show in section 4 that the masses of charged hadrons can be extracted with the same numerical accuracy as their almost-degenerate neutral counterparts. This is true both in the pseudoscalar and in the vector meson channels. At the values of the bare parameters used in our study, the pseudoscalar-meson charged-neutral mass splitting can be extracted with statistical significance even in the simulation performed at α em = 1/137.
We have also sketched a strategy to probe states of charged mesons with real photons. The proposal consists of using gauge-invariant interpolating operators that, at leading order in α em , have maximal overlap with states having a fixed number of real photons. Although much more work is certainly needed in this direction, the results of subsection 4.2 represent a promising indication on the numerical validity of this approach. CERN, managed by the HPC team in the IT Department, and on the Marconi system at CINECA under the initiative INFN-LQCD123.

A Gauge-fixed two-point functions
The goal of this appendix is to illustrate some of the subtleties that arise in the charged sector, when the U(1) gauge is fixed. For definiteness we work here with the familiar case of covariant gauge, in continuum notation. In order to avoid potential issues with IR divergences, we consider QCD+QED in a spatial box with size L 3 and C boundary conditions for all fields. For simplicity we consider an infinite time extent. In Euclidean spacetime, the action in covariant gauge is where S 0 is the gauge-invariant part of the action, A µ and B µ are the photon and gluon fields, while ψ andψ are the quark fields, and the scalar product is defined as Let h(x) be some local operator which interpolates a hadron with electric charge q h , and leth(x) the interpolating operator with the corresponding antiparticle. We are interested in the two-point function The integrands do not depend on λ, therefore the auxiliary integral over λ gives an infinite constant which simplifies in the ratio. We change variables in the two integrals to the gauge-transformed fields The interpolating operator and action transform as After this change of variables, the integral over λ is Gaussian and can be calculated analytically, yielding the following gauge-invariant representation where the current J µ (z) is defined by the equation Because of C boundary conditions, the Laplacian = ∂ µ ∂ µ is defined with antiperiodic boundary conditions in space and is therefore invertible. The expectation value in eq. The integration by part (J µ , ∂ µ λ) = −(∂ µ J µ , λ) does not generate boundary terms since the product J µ λ satisfies periodic boundary conditions. The factor e iq h [λ(y)−λ(z)] in the above equation cancels the phase generated by the gauge transformation of h(z)h(y). As a consequence, the observable in eq. (A.6) is invariant under local gauge transformations. It is tempting to interpret the gaugeinvariant observable as a possible interpolating operator for the charged hadron h. In fact this operator is formally very similar to Dirac's interpolating operator. However H(x) is non-local in time and a standard interpretation as an interpolating operator is not possible. The Hamiltonian representation of the expectation value in the r.h.s. of eq. (A.6) is obtained by interpreting the phase as a term of the action. As in the case of J = 0, the action S 0 + i(J µ , A µ ) defines a constrained Hamiltonian system. States propagating in the gauge-invariant two-point function satisfy the Gauss law in presence of the charge density j 0 (x) = f q f ψ † f ψ f (x) of the dynamical degrees of freedom, and the external time-dependent charge density J 0 (t, x), i.e.
The evolution of states is governed by a time-dependent non-hermitean Hamiltonian where H 0 is the standard gauge-invariant Hamiltonian without external current. Notice that for z 0 t y 0 , the four-current vanishes exponentially, i.e.
On the one hand, this is a way to see that the leading exponential behaviour of the two-point function is determined by the ground state in the charged sector of the gauge-invariant Hamiltonian H 0 . Therefore the mass defined by means of the two-point function in covariant gauge is the correct one. On the other hand, the unphysical exponentials in the external current mimic the contribution of excited states in the long-distance behaviour of the two-point function. For this reason the covariant gauge is not a suitable choice for the extraction of excited states from two-point functions.

B Explicit expressions for two-point functions
In this appendix we provide explicit expressions for the two-point functions used in this work, in which fermions have been integrated out. Because of C boundary conditions, the fermion Wick contractions are not the usual ones in terms of the original fields ψ f andψ f . For instance, the ψψ Wick contraction does not vanish. For this reason, we find more convenient to work with the quark-antiquark doublet η f defined in eq. (3.2). The neutral meson operators considered in this work can be easily written in terms of the η f field,s We rewrite the interpolating operator for a P+γ state in the V channel as 14) The two new correlators used for the generalised-eigenvalue problem in section 4.