Simulations of QCD and QED with C* boundary conditions

We present exploratory results from dynamical simulations of QCD in isolation, as well as QCD coupled to QED, with C* boundary conditions. In finite volume, the use of C* boundary conditions allows for a gauge invariant and local formulation of QED without zero modes. In particular we show that the simulations reproduce known results and that masses of charged mesons can be extracted in a completely gauge invariant way. For the simulations we use a modified version of the HiRep code. The primary features of the simulation code are presented and we discuss some details regarding the implementation of C* boundary conditions and the simulated lattice action.


Introduction
Calculating electromagnetic corrections to hadronic observables via lattice simulations requires a consistent formulation of QCD+QED in finite volume. Such a formulation, based on C boundary conditions, was proposed and thoroughly discussed in [1]. Here we present the first exploratory simulations showing that the proposed setup works in practice. For the simulations we use a modified version of the HiRep code, as discussed in section 2, and in section 3 we present our exploratory results. In particular, we show that masses of charged hadrons can be extracted in a gauge invariant way, with a signal-to-noise ratio equivalent to that of neutral hadrons.

Simulation code
In this section we present a modified version of the HiRep code [2] used for simulating QCD and QED with C boundary conditions. In particular we outline the primary features of the code and discuss the approach used to implement the boundary conditions.

Overview
The HiRep code was developed as a flexible code for studying BSM models and for this reason it has native support for different gauge groups and higher dimensional fermion representations. Especially the latter feature was one of the main reasons for choosing to modify this code, as discussed later on.

Oct 2017
The modified code presented here can simulate QCD, QED and QCD+QED with either periodic or C boundary conditions in space. The main features of the code include: • Wilson fermions with O(a) improvement • Lüscher-Weisz gauge action for SU(3) field • Plaquette gauge action for U(1) field • Compact QED with optional Fourier acceleration • Support for rational approximations • Hierarchical OMF integrators • Several inverters, such as MINRES, BiCGstab, and CG with multishift support • Gradient flow observables for both SU(3) and U(1) field • Measurements of charged and neutral meson correlators The approach used for implementing the C boundary conditions has been described in the appendix of [1]. For the discussion at hand we recall that C boundary conditions correspond to performing a charge conjugation when wrapping around the torus. Denoting the matter fields by ψ f and the U(1) and SU(3) gauge potentials by A µ and B µ , respectively, the boundary conditions read: Here C is the charge conjugation matrix. These boundary conditions pose no additional complications for the gauge fields. In terms of the link variables, for both gauge fields, the boundary conditions correspond to complex conjugating the links when wrapping around the torus.
In practice this constraint is implemented by having an extended lattice with a static border around the actual lattice. The relevant field variables are then copied into the static border and the boundary conditions are applied. This operation must be repeated every time the field changes.

Fermion representation
While the boundary conditions for the gauge fields are trivially implemented, special care is needed for the matter fields because of the mixing between fermions and antifermions at the boundary. Due to this mixing it is advantageous to write the fermion action in terms of a new doublet χ containing both fermion fields.
In this notation the C boundary conditions simply swap the two components of the doublet.
The boundary conditions can be diagonalized via a unitary transformation in which case the basis is the two eigenstates ψ ± of the charge conjugation operator.
In this basis the boundary conditions read and the link variables appearing the Dirac operator are 6 × 6 real matrices constructed via the map The relation D[U] T = CD[U * ]C −1 can now be used to rewrite the action.
Because the boundary conditions have been diagonalized we can evaluate the associated path integral. This results in the Pfaffian of the Dirac operator instead of the usual determinant. [ Using the fact that CD[U η ] is an antisymmetric matrix, the absolute value of the Pfaffian can be related to the determinant via the identification Due to this relation, on the lattice we can use the pseudofermion method to represent the Pfaffian, but it requires the use of rational approximations in the action.
Because the HiRep code already supported higher dimensional representations, defining the new representation U η used in the Dirac operator, amounted to writing a new function that implements the map defined in Eq. (10).

Lattice action
Because the quarks are fractionally charged, when using the compact formulation of QED it is necessary to rescale the elementary charge to ensure gauge covariance. We define the U(1) action as where the bare coupling β em now depends on the electromagnetic coupling constant α and the elementary charge q el . The latter is a tunable parameter chosen to be q el = 1/6 in our simulations.
The Dirac operator used for QCD+QED simulations is written in terms of the U(3) links defined by where V µ (x) is the U(1) link and U µ (x) is the SU (3) link. The integer exponentq = q/q el is the fermion charge in units of the elementary charge. With these definitions, the O(a) improved Dirac operator can be written as where we use the clover definition of the field tensors. For the simulations presented in the next section, we use the tree-level coefficient for the U(1) clover term given by c QED sw =q.

Interpolating operators
The construction of gauge-invariant interpolating operators for charged states was discussed in [1]. In the next section we use the so-called "Coulomb operator" for measurements of charged meson states.
The electrostatic potential Φ(x) must be anti-periodic and satisfy ∂ k ∂ k Φ(x) = δ 3 (x). This operator is invariant under local gauge transformations, but transforms under global gauge transformations, and as such it defines a charged state. In the Coulomb gauge Ψ(x) = ψ(x), and the operator Ψ(x) is therefore the unique gauge-invariant extension of the quark field defined in the Coulomb gauge. We refer again to [1] for a definition of the discretized operator implemented in the code.

Exploratory studies
In this section we discuss some exploratory simulations of QCD with C boundary conditions as well as dynamical QCD+QED simulations. For the QCD simulations we use parameters from the CLS collaboration to allow for a comparison with their results. For the QCD+QED simulations we start from the parameters of one of the QCD simulations, while adding the QED interactions in the Dirac operator. Introducing the QED interactions will, most notably, shift the value of the critical masses, and hence affect the mass spectrum. We study how the mass spectrum changes when adding the QED interactions, and in particular, we show that the masses of charged hadrons can be extracted in a completely gauge invariant way.

QCD
We have performed several QCD simulations to show that our setup with C boundary conditions works in practice. For the simulations we have chosen N f = 3 dynamical fermions at the isospin symmetric point with the bare parameters taken from the CLS collaboration [3]. The simulated ensembles are listed in Table 1 together with the bare parameters and the number of thermalized configurations used in the analysis. Except from ensemble A1 we use C boundary conditions in the spatial directions, and periodic boundary conditions in time. Ensemble A1 was simulated with periodic boundary conditions in all directions to allow for a direct comparison between the two cases. For all ensembles we have measured the PCAC mass, the pion/kaon mass and decay constant and the gradient flow observable t 0 . The results are listed in Table 2, where the PCAC mass is unrenormalized, but O(a) improved using c A from [4], and the decay constant has been renormalized using Z A from [5].
The results for the A (B) ensembles should be compared with CLS ensemble H101 (H200) in [6]. Because our volumes are smaller than the corresponding CLS ensembles we have significant finite volume effects, especially on the B1 ensemble. On the two largest volumes A3 and B2 our results agree with the CLS values within 5% in all cases. This seems reasonable given the smaller volume and smaller amount of statistics.

QCD+QED
For our first test simulations of dynamical QCD+QED we took the parameters for the B1 ensemble and added the QED interactions. While, in these simulations, the bare masses are still degenerate for all three quarks, in the Dirac operator we use the physical charges i.e. we have one up-type quark with charge q = 2/3 and two down-type quarks with q = −1/3. Because the charges are different, the quark masses are also renormalized differently, and for this reason it makes little sense to keep the bare masses degenerate. However, for these exploratory simulations we just wanted a simple setup for studying the two-point functions.
In Table 3 we list the two simulated QCD+QED ensembles, with the only difference being the value of electromagnetic coupling constant. In the first ensemble we use an unphysically large value of the coupling constant α ≈ 7α phys and in the second ensemble we use the physical value.
As previously stated, the primary goal of these exploratory QCD+QED simulations is to prove that we can extract the masses of charged hadrons in a gauge invariant way using the operator defined in section 2.4. In particular we have studied the mass of the charged and neutral kaon. In Fig. 1 Table 3. Ensembles and parameters for the QCD+QED simulations.
the Q2 ensemble. In the first row we show the correlators for the two states, in the second row the corresponding effective masses, and in the last row we show the effective mass splitting. This effective splitting, denoted ∆ K in the plots, is defined as the mass difference divided by the average mass, i.e.
For the Q1 ensemble we observe a large mass splitting with ∆ K = 0.2378(46) due to the unphysically large value of the electromagnetic coupling. As expected, we also observe that the masses increase compared to the QCD simulation (ensemble B1), due to the QED interactions shifting the critical bare masses. On the Q1 ensemble the mass splitting is naturally smaller with ∆ K = 0.065(15). While these simulations are affected by finite volume effects, we are clearly able to extract a statistically significant signal even for the physical value of the electromagnetic coupling constant.

Optimizations
As argued in section 2.2, the RHMC algorithm is necessary for simulations with C boundary conditions. This will naturally make the simulations significantly more expensive, but there are several ways to reduce the cost, such as pole-splitting for the rational approximation or the introduction of multiple pseudofermions [7]. We explored the possibility of using multiple pseudofermions for our QCD simulations and we were eventually able to reduce the cost by more than a factor of two 1 . In our setup we can simulate N f = 3 degenerate fermions with a single pseudofermion when using the fraction 3/4 in the rational approximation. This was, however, quite expensive because of the relatively ill-conditioned Dirac operator. For this reason we changed the action to include two pseudofermions, each with a fraction of 3/8, which resulted in a much better conditioned setup. In practice it allowed us to reduce the number of integration steps by roughly a factor of five without loosing acceptance.
In Fig. 2 we show how the mean value and standard deviation of the numerical MD forces decrease when using two pseudofermions for the A3 ensemble, which explains why the number of integration steps can be decreased. On the plot |F| 2 is the squared norm of the force vector, and the average is over all positions and directions. Naturally, we also tried using three pseudofermions, but no futher gain was observed in this case. For our QCD+QED simulations we naturally need two pseudofermions because of the different charge assignments for the quarks, which automatically leads to a reasonably well-conditioned setup. What makes a difference for these simulations is the use of Fourier acceleration for the U(1) field. This was originally implemented to decrease the autocorrelation, but unexpectedly, the number of integration steps could also be decreased by a factor of two, without loosing acceptance, compared to the same simulation without Fourier acceleration. This means that, while the use of Fourier acceleration is slightly more expensive, it is well compensated by the possible decrease in the number of integration steps. 1 We use a simple two-level integration scheme with fermions on the outer level and gauge on the inner level.

Conclusion
We presented a modified version of the HiRep code suitable for simulations of QCD and QCD+QED with C boundary conditions. The primary features of the code have been discussed together with some implementation specific details. The code has been used to perform the first exploratory simulations and we have shown that we are able to reproduce known results and that masses of charges hadrons can be extracted in a gauge invariant way, with a signal-to-noise ratio equivalent to that of neutral hadrons. Moreover, we have discussed a few possible ways of optimizing the cost of the simulations.