Lattice simulations with G-parity Boundary Conditions

We discuss G-parity lattice boundary conditions as a means to impose momentum on the pion ground state without breaking isospin symmetry. This technique is expected to be critical for the precision measurement of $K\rightarrow(\pi\pi)_{I=0}$ matrix elements where physical kinematics demands moving pions in the final state and the statistical noise caused by disconnected contributions will make it difficult to use multi-exponential fits to isolate this as an excited state. We present a formalism for computing hadronic Green's functions with G-parity boundary conditions, derive the discretized action and its symmetries, discuss how the strange quark can be introduced and detail techniques for the numerical implementation of these boundary conditions. We demonstrate and test these methods using several $16^3\times 32$ dynamical domain wall ensembles with a $420$ MeV pion mass and G-parity boundary conditions in one and two spatial directions.


I. INTRODUCTION
Important theoretical advances [1,2] have opened a path for the determination of physical decay amplitudes comprising multi-particle final states using lattice QCD. In particular, calculations of the magnitude of direct CP-violation in the decays of kaons into two pions are now possible [3][4][5], which offers a novel and exciting test of the Standard Model. The study of such decays with physical kinematics invariably requires moving final-state particles, for which it may be difficult to isolate the signal over the typically much larger contribution involving stationary particles. This can be avoided by modifying the lattice boundary conditions to remove the state with stationary pions from the particle spectrum. For charged pions this can be achieved by imposing antiperiodic boundary conditions on just the down or up quark propagators, which results in pion momenta that are odd-integer multiples of π/L. This technique has been used successfully in the calculation of the ∆I = 3/2 K → ππ amplitude [3,4], but is limited by the fact that it applies only to charged pions and also explicitly breaks the isospin symmetry. These issues were avoided in the aforementioned case by utilizing the Wigner-Eckart theorem to relate the desired K + → π + π 0 decay to the unphysical decay K + → π + π + , for which the final state contains only charged pions and is protected from mixing with other isospin states by virtue of being the only charge-2 state with those quantum numbers. This is not the case for the calculation of the ∆I = 1/2 K → ππ amplitude, where the final state necessarily contains neutral pions and the isospin breaking cannot be circumvented. G-parity boundary conditions (GPBC) [6][7][8] provide an alternative approach that results in antiperiodic charged and neutral pions, and also preserves isospin symmetry.
The RBC & UKQCD collaborations have successfully employed GPBC to compute the ∆I = 1/2 K → ππ amplitude [5]. This document is intended as a partner to that work in which we lay the theoretical groundwork for employing G-parity boundary conditions in zero temperature lattice simulations and perform thorough numerical investigations of a number of its key features and difficulties using three 16 3 ×32×16, 2+1 flavor dynamical domain wall ensembles with GPBC in 0, 1 and 2 directions, respectively.
We begin in Section II by outlining the transformation of the pion and quark fields under Gparity, and introduce notation for later use. In Section III we derive the appropriate discretized action in our notation that makes explicit the boundary field mixing. We also discuss the numerical implementation of GPBC. The symmetries of the lattice action are investigated in Section IV, and we detail the appropriate operators for light and heavy meson states in Sections V and VI. In the latter we also discuss how the strange quark can be introduced into this framework. In Sections VII and VIII we present numerical demonstrations of the G-parity technique on the aforementioned ensembles. Our conclusions are presented in Section IX.
A related development of charge conjugation boundary conditions, introduced for a different purpose, has been carried out by Lucini, Patella, Ramos and Tantalo [9].

II. G-PARITY BOUNDARY CONDITIONS
G-parity is conventionally defined as a product of charge conjugation and an isospin rotation by π radians around the y-axis:Ĝ where the hat-symbol is used to denote operators. Here the choice of the y-axis is dictated by the convention that the representation matrices for I x and I z are real while those for I y are imaginary and impliesĜ commutes with isospin rotations. The charged and neutral pions are G-parity odd, therefore if imposed as a spatial boundary condition the pions become antiperiodic and have discretized momenta that are odd-integer multiples of π/L, where L is the lattice spatial extent.
For QCD we are concerned with the action of G-parity upon the quarks. In this section we derive the appropriate transformation and introduce a convenient notation for studying its implications that we will use throughout the remainder of this document.

A. Action upon quark fields
The isospin rotation around the y-axis transforms the flavor-doublet of light-quark annihilation operators as follows: To describe charge conjugation we recall the action ofĈ on the fermion fields q andq [10]: where C is the 4 × 4 charge conjugation matrix (not to be confused with the operatorĈ) which obeys (4) C −1 γ µ C = −γ T µ .
The usual requirement thatĈ 2 = 1 implies that C T = −C as well as the important property that when acting on states with isospin I:Ĝ 2 = (−1) I .
While we will not adopt specific conventions for the Euclidean γ matrices used in this paper, we find it convenient to choose a chiral basis where γ 5 is real and diagonal which, together with Eq. (3), implies that [C, γ 5 ] = 1. In addition we will assume that the matrix C is unitary and real, conditions consistent with the conventions, for example, of Ref. [10]. We summarize the properties of the matrix C: Combining the two operations we find that G-parity has the following action on the quark flavor Notice that this involves an explicit mixing of the flavors and also of the spinors and conjugate spinors. We will see that this leads to a number of complications.

B. Notation
For this document we adopt a convenient notation in which we place the field operators u, d We will refer to the indices of these vectors as 'flavor' indices. Note that ψ and ψ transform in the same way as the quark fields under charge-conjugation (Eq. (3)). The benefit of this notation is that the action of the G-parity operator upon ψ takes on a simple form: where σ 2 is the second Pauli matrix.
Using Eq. (7) we can write q = ū ,d = ψ T CF 21 + ψF 12 , and the inverse transformations, ψ =qF 21 + q T CF 12 , where These relations allow us to transform back and forth between the two notations. For later use we also define C

. Translations of fermion fields
In this document we define our coordinates on a torus (i.e. modulo the lattice size), such that x µ = L µ is equivalent to x µ = 0; x µ = −1 is equivalent to x µ = L µ − 1; and so on. The corresponding field variables are treated in the same way: ψ(x µ = L µ ) ≡ ψ(x µ = 0); ψ(x µ = −1) ≡ ψ(x µ = L µ − 1); etc. We will assume that the spatial box dimensions are all equal to L for convenience. The action of the boundary condition upon the field is treated explicitly in the context of a spatial translation via the introduction of a coordinate-dependent coefficient given below. In the following section we will rewrite the usual covariant derivative of the fermion action in terms of the translation operator such that the symmetry of the action under these modified translations can be imposed.
The boundary condition implies that as we move from site to site along a G-parity direction from the origin we encounter the fermion fields in the following order: where || indicates the location of the lattice boundary, and we have suppressed the coordinates along the other three directions for clarity. This implies the action of the translation operatorT µ on the fermion field in a G-parity direction µ is as follows: where we have again suppressed coordinates along the other three directions. These equations can be neatly summarized by introducing unitary flavor matrices for which G µ is unity in directions with GPBC and zero otherwise. These objects are related as follows: The translations then become: Note the four operatorsT µ obey the commutation relations required for elements of the group of translations, [T µ ,T ν ] = 0.

III. DISCRETIZED ACTION
In this section we derive the appropriate discretized action for two light quark flavors with G-parity boundary conditions.
A. Lattice QCD action with periodic BCs in our notation

Fermion action
We begin with the usual four-dimensional Euclidean lattice fermion action for two-flavor QCD with periodic boundary conditions in the spatial directions, where the spin matrices Γ ± µ = 1 2 (1 ∓ γ µ ) for the Wilson/domain wall actions or Γ ± µ = ± 1 2 γ µ for the naïve action. Here and below the sum over each of the four components of the coordinate x runs from x µ = 0 to x µ = L µ − 1 where x µ and L µ are expressed in lattice units. As mentioned previously we will assume the spatial dimensions are all equal in extent: L x = L y = L z = L for convenience.
Note that in these expressions we interpret the quantities q andq as Grassmann variables that would appear in a path integral. We will use these operator and Grassmann interpretations interchangeably and specify a particular choice only when it is necessary.
The above can be rewritten in terms of the quark field vectors ψ and ψ defined in Eq. (7) using Eqs. (9) and (10): (22) This expression can be simplified by introducing the matrix which has both flavor and color indices, and with which the action becomes For use below it is convenient to rewrite the covariant derivative in terms of the translation operators:

Gauge action
As the fermion action is expressed in terms of the flavor-matrix gauge links, it is convenient to do the same for the gauge action. We assume the Wilson action, although it it straightforward to generalize this result. The action is (26) is the usual plaquette. We can construct a similar "plaquette" from the flavored gauge linksŨ µ : for which (29) where it is understood that the trace on the right-hand side includes the flavor indices. The Wilson action then becomes As with the covariant derivative, it is useful to express the plaquette in terms of translation operators: B. Lattice QCD action with G-parity BCs

Fermion action
The covariant derivative for G-parity BCs can be obtained by inserting Eqs. (20a) and (20b) into Eq. (25): The translational properties of the flavor-matrix gauge links can be found by imposing gauge invariance upon the derivative term of the action, Under a gauge transformation V the field ψ and its conjugate transform as The forwards component of S ∇ then transforms as, whereŨ ′ µ is the gauge transformation ofŨ µ . Invariance of this term under the gauge transformation then implies (36) For the backwards component of S ∇ the termT −1 µŨ † µ (x)T µ enters. If we assume that with α and β flavor matrices, the backwards component of S ∇ transforms under a gauge transformation as follows: We can now write down the covariant derivative with G-parity BCs: the corresponding (complete) fermion action and the Dirac matrix, Note that throughout this document we consistently ignore the fact that domain wall fermions have an additional discrete index associated with their coordinate in the fifth dimension, and instead treat them identically to Wilson fermions. In Appendix A we demonstrate that in our (suitably extended) ψ-field notation the G-parity boundary condition does not affect the fifth dimensional coordinate despite the reflection in this dimension induced by charge conjugation, further evidencing the power of this notation and justifying us dropping this index.
It is important to recognize that the introduction of G-parity boundary conditions does not alter the usual γ 5 -hermiticity of the Dirac propagator, This can be seen from the expression given in Eq. (42) for the Dirac operator, which has the structure of the usual Euclidean lattice Dirac operator with the exception of the appearance of the 2 × 2 matrices B ± µ . We can then deduce γ 5 -hermiticity following the usual steps using Eqs. (19a) and (19b) and [γ 5 , B ± µ ] = 0. The quark propagators obtained by inverting this Dirac matrix are 2 × 2 matrices in flavor space as well as being matrices in spin and color space. In practise this leads to additional diagrams that must be evaluated when computing hadronic observables. In addition, we must typically compute the inverse with separate sources for each flavor, doubling the number of matrix inversions required to compute the full propagator. However in many cases this requirement can be circumvented by taking advantage of the isospin symmetry, as we demonstrate in Section. IV A.

Gauge action
Expressing the Wilson gauge action Eq. (30) in terms of translation operators acting on the gauge links via Eq. (31), we have The action of the translation operator upon the links has thus far been derived only for links in the same direction as the translation (Eq. (39)). The corresponding relation for links that are orthogonal to the translation can be derived by imposing gauge invariance on S W . Assuminĝ where α − δ are flavor matrices, and then applying Eq. (36) it is straightforward to show that the following action is gauge invariant: for which α = β † = B + µ (x µ ) and γ = δ † = B + ν (x ν ). Combining the above results for α −δ with Eq. (39) we obtain the general action of translations on the flavored gauge links: Note that usingT µŨν (x)Ũ † ν (x)T −1 µ = 1 it is easy to see that the translation operators act onŨ † ν in the same way as they do onŨ ν .
Eq. (46a) implies that as we again move from site to site along a G-parity direction µ from the origin, the gauge fields are encountered in the following order: This implies the links transform asŨ µ → (iσ 2 )Ũ µ (−iσ 2 ) =Ũ * µ under the action of the boundary, i.e. the links must obey complex conjugate (or equivalently, charge conjugation) boundary conditions. This of course requires new gauge configurations to be generated for GPBC.

C. One-flavor equivalence
Consider the upper component of the flavor doublet field ψ for G-parity in one direction. As we move from site to site in the G-parity direction we encounter the fields in the following order: where || again indicates the position of the lattice boundary and we suppress coordinates other than in the G-parity direction. This infinite series is antiperiodic in 2L, and the subset between 0 and 2L contains all of the fermionic degrees of freedom in the G-parity setup. We can therefore define a field Ψ on a lattice of size 2L (denoting the corresponding coordinates with capital letters) with the following mapping: , which contains all the fermionic degrees of freedom and obeys antiperiodic boundary conditions in 2L. Similarly, the gauge links U ν defined on this doubled lattice are mapped as follows: Consider the forwards component of the action, Eq. (41), in the G-parity direction µ, and suppress the coordinates in the other directions: where B + µ (X µ ) = exp(iπδ X µ ,2L−1 ) imposes the appropriate sign for the term crossing the boundary in 2L. For the backwards component we likewise obtain where B − µ (X µ ) = exp(iπδ X µ ,0 ). The two terms together comprise the action for a single flavor quark field residing on a lattice of size 2L with antiperiodic boundary conditions in the G-parity direction µ, with the only additional condition being that the gauge links on the second half of the lattice are the complex conjugates of those on the first. This establishes a direct equivalence between the two-flavor theory and a one-flavor theory on a doubled lattice that proves very useful when it comes to implementing these boundary conditions on a computer.
Let us consider extending this formalism to GPBC in a second direction, y (with the original doubling in the x-direction). This can be achieved by doubling the lattice again in the y-direction and imposing APBC in this second direction. However in this approach some care is required to recognize that the fermion fields on the four quadrants of the resulting lattice are not independent but are related according to Figure 1. This has implications for example in the construction of propagator sources: a down-quark source on timeslice τ with spatial smearing function Θ centered at position x 0 in the two flavor setup corresponds to the following in the twice-doubled single-flavor setup: where D LL ( X) and D UR ( X) are functions that are unity on the lower-left and upper-right quadrants, respectively, and zero elsewhere (referring to Figure 1). Here the minus sign between the terms arises because the fermion field in the upper-right quadrant has the opposite sign to that of the lower-left quadrant. Note also that this second doubling of the lattice doubles the cost of performing the Dirac matrix inversion relative to the two-flavor approach, and therefore the twice-doubled one-flavor approach is of limited practical use for GPBCs in multiple directions.
An alternative treatment for GPBC in two directions, easily generalized to three, is to modify the boundary conditions in the y-direction such that passing through the boundary is accompanied by a translation by L in the x-direction; this is illustrated in Figure 2. This approach avoids further doubling for GPBC in more than one direction and therefore has the same computational cost as the two-flavor approach, but is somewhat complicated to implement numerically, and the non-nearest neighbor communications pattern in the y-direction is suboptimal for most parallel machines.
We therefore conclude that this one-flavor mapping, while a useful cross-check, is not practical for high precision lattice calculations with GPBC in more than one direction. The alternative, which we employ in practise, is to implement the full two-flavor theory with the appropriate mixing of the flavors at the boundary directly.

IV. LATTICE SYMMETRIES
The use of G-parity boundary conditions has a number of symmetry implications that we detail in this section. Given that the action away from the boundary is simply a rewrite of the usual lattice action, we need only consider the effects of global symmetry operations upon the boundary terms. To do so it is convenient to rewrite the matrices B ± defined in Eqs. (18a) and (18b) as such that the contribution to action from the G-parity boundary conditions is contained in the following expression: (56) The total action then becomes (57) where S std is the 'standard' action (58) for a periodic field. Here we see that S GPBC acts to subtract the boundary term of the standard action and add the correct G-parity flavor 'twist'.
We can simplify Eq. (56) somewhat by noting that for coordinates defined on a torus (in one We can then write (60) where for the second term we have transposed the Grassman variables and used the following identity: In terms of the usual flavor doublets we can rewrite Eq. (60) as which proves convenient in some cases.

A. Isospin
G-parity commutes with isospin rotations around the y-axis by construction. This leads to a very useful identity for the quark propagators: The action of such a rotation by π radians on the quark fields is given in Eq. (2). In our G-parity notation the fields ψ and ψ therefore transform as (63a) e −Î y π ψeÎ y π = (iσ 2 )Cψ T (63b) e −Î y π ψeÎ y π = ψ T C(−iσ 2 ) .
We write the path integral in which we have integrated over only the fermion fields on a fixed gauge background as with which the propagator -the inverse of the Dirac matrix on a given configuration -can be written as where α,β are combined spin/color/flavor indices.
As the action is invariant under isospin rotation we can write where we have absorbed the sign from commuting the fields using C −1 = −C. γ 5 -hermiticity then implies Examining the flavor structure explicitly we find   G 00 (y, z) G 01 (y, z) Notice that this implies the second column of the flavor-matrix propagator (source flavor index 1) can be obtained entirely from the first column (source flavor index 0), hence we can calculate the full propagator using only the matrix computed from sites of a single flavor. We will exploit this property in Section VIII.

B. Baryon number
Since our G-parity boundary conditions change quarks to anti-quarks, baryon number symmetry is violated. This is dramatically illustrated by the transformation of a proton (uud), with baryon number B = 1, which becomes an anti-neutron (ddu) with baryon number -1 at the boundary.
The baryon number violation has an additional manifestation at the quark level: The mixing of quark flavor at the G-parity boundary allows for the Wick contraction of up and down field operators: As a result there are typically additional diagrams involving propagators that cross the boundary.
In the first of the above contractions, quark flavor flows towards the boundary on both sides.
Likewise, quark flavor flows away from the boundary in the second contraction. We may interpret this as the boundary respectively destroying and creating flavor.
Other than the mixing of the quark flavors at the boundary which we handle explicitly in our two-flavor formulation, the breaking of the baryon number symmetry is not important for calculations involving only mesonic states.

C. Flavor non-singlet axial vector transformations
The action is invariant under the flavor non-singlet vector (isospin) transformations by con-

struction. However for an axial transformation
For Wilson/domain wall fermions, the Wilson term in Eq. (62) explicitly breaks the axial symmetry, therefore to separate the effects of G-parity we consider naïve fermions, for which Γ ± µ = ±γ µ . For the second and third terms of Eq. (62), we have Aγ µ = γ µ A † , such that the terms are invariant.
For the other two terms we note (71a) Therefore if we write S GPBC → S GPBC + ∆S GPBC (where ∆S GPBC = 0 would imply invariance of the action), we find The action is therefore not invariant under the flavor non-singlet axial transformations.
This explicit breaking of chiral symmetry gives rise to a non-zero chiral condensate even for zero quark mass. This may provide a useful tool when studying finite temperature QCD and was one of the original motivations for the study of G-parity boundary conditions by Wiese [6].
For our purposes we are interested only in the low-temperature behavior of the massive theory, where the chiral symmetry is also broken spontaneously by the dynamics and explicitly by the lattice fermion formulation. In this regime the effects of the G-parity boundary conditions enter in two ways: The first are due simply to the different sets of allowed momenta for G-parity even (integer multiples of 2π/L) and odd (odd-integer multiples of π/L) states. Such effects are straightforward to take into account in our measurements. There are also more subtle effects which, in the context of the low-energy theory of interactive massive pions, enter as a change in the momentum discretization of pion loop diagrams. These effects are exponentially suppressed in m π L according to the Poisson summation formula and are therefore comparable in size to other, more conventional finite volume effects. As a result we need not be concerned that this boundaryinduced chiral symmetry breaking will have a significant effect upon our measurements.

D. Parity
It is convenient to define the parity transformation P acting on the points in our finite lattice thus: x i → x P i = L − 1 − x i for i in spatial directions, such that the lattice coordinates are inverted: 0, 1, 2 . . .(L − 1) → (L − 1) . . . 2, 1, 0. This is equivalent to reflecting about the midpoint (L − 1)/2 of each spatial direction.
Under parity, the gauge links on a standard periodic lattice transform as where P(i) = −i for spatial directions i = 1, 2, 3 and P(4) = 4 for the time direction, and The flavor-matrix gauge linksŨ µ = diag(U µ ,U * µ ) therefore transform as where the analog to the right-hand side of Eq. (74) should take into account the non-trivial boundary condition on the links. To determine the appropriate form it is convenient to rewrite this equation in terms of the translation operators, then use Eq. (46b): We will also require the corresponding action of parity onŨ i (x −î) for a spatial direction i, which can be found usingPŨ µ (x)Ũ † µ (x)P −1 = 1 and Eq. (75), resulting in and thus The action of parity on the down quark field d = ψ 0 is: For the second component, ψ 1 = Cū T , we have The action of parity upon the full G-parity fermion doublet is therefore For determining the transformation properties of the action under parity it is more convenient to return to the full action, Eq. (41). Applying the parity transformation the first term of the action in the spatial direction i: Similarly, the second term transforms as These match the second and first terms of Eq. (41) respectively, written in terms of the transformed coordinates. The temporal components are trivially invariant because Γ ± 4 commutes with γ 4 and B ± 4 ≡ 1. We therefore see that the action is invariant under the parity transformation.
In the above we have used (84) 4 are unit matrices hence this relation is trivially applicable.) We also recognize that for any flavor matrix F, which follows from the fact that B ± µ are unitary and are either 1 or ±iσ 2 depending on the coordinate, and that (±iσ 2 ) 2 = −1.

E. Translational symmetry
Because the gauge and fermion actions are constructed as a sum of local terms containing differences between neighboring sites which are obtained using the translation operatorsT µ , we can write the complete action as the sum: where s(0) and s W (0) are the terms in the fermion and gauge actions corresponding to the point x = 0. Since the operatorT µ L is a symmetry of both the fermion and gauge actions, the summand in Eq. (87) depends on the integer summation variables n µ only through (n µ mod L). This implies that if we conjugate S + S W with a translation operatorT κ we will increase the summation variable n κ by one which simply permutes the terms in the sum and leaves the complete action unchanged.

F. Translational covariance of field operators
In Eqs. (20a) and (20b) we observe that the naïve translational covariance of the quark field is broken at the boundary, where it picks up an additional matrix structure ±iσ 2 . This implies that the quantity (88) is not an not an eigenstate of translation, i.e. (89) for a G-parity direction µ. This is problematic as we ultimately wish to construct states of definite momentum.
We can easily form eigenstates of iσ 2 that simply pick up coefficients at the boundaries: The factor of 1 2 in Eq. (90) is arbitrary; we choose it such that ψ can be interpreted as the sum of the two eigenvectors. The fields ψ ± have the following boundary conditions: where µ is a G-parity direction.

G. Rotational symmetry
Examining Eq. (95) more closely, we observe that a non-zero field operator with a definite momentum can only be created if the flavor projection operators for each momentum direction have the same sign, i.e. e in p i π is the same for all G-parity directions. With GPBC in two directions for example, this then implies where m is an integer. In other words, the momentum components are constrained to differ only by integer multiples of 2π/L. We can therefore simplify the expression for the translationally covariant field in Eq. (95) to The implications of this observation can be seen by considering the two momenta 2π L ( 1 4 , 1 4 , 0) and 2π L ( 1 4 , − 1 4 , 0), which are related by a cubic rotation. In the former the momentum components are identical, hence this is an allowed momentum. However, in the latter the two momentum components in the G-parity directions differ only by π/L, and therefore this is not an allowed quark momentum. This implies that the rotational symmetry has been broken at the quark level by GPBC in multiple directions. For example if we imposed GPBC in all three spatial directions in a cubic box, the naïve cubic symmetry of this choice is broken by this fixed relation between the three components of the allowed quark momenta.
The breaking of the rotational symmetry can be shown in a different manner by considering the structure of the Brillouin lattice: In Section III C we demonstrated that the two-flavor theory with GPBC in one direction of size L is equivalent to a single-flavor theory on a lattice of size 2L with the quarks obeying antiperiodic boundary conditions. In Figure 3(a) we plot the allowed momenta for this doubled-lattice setup. In Figure 3(b) we plot those obtained if we double the lattice again and impose antiperiodic BCs in the second direction. In each case the number of points corresponds to the number of fermionic degrees of freedom, which is double in the latter compared to the former and is therefore not equivalent to GPBC in two-directions, for which the number of degrees of freedom is independent of the number of G-parity directions. (In Section III C this binding of field variables between the quadrants required us to carefully construct propagator sources in order to correctly reproduce the two-flavor theory.) We must therefore eliminate half of these points to reproduce the two-flavor theory.
Given that parity symmetry demands the Brillouin lattice be symmetric about the diagonals, and that the sites are equally spaced along both axes, there are only two possible configurations for the remaining Brillouin lattice sites: The first has all of the momentum points distributed along the positive diagonal as we plot in Figure 3(c), and matches the allowed momenta we identified above.
The second is that in which the momentum sites reside in the center of the currently-unoccupied grid squares of this figure such that they lie instead along the negative diagonal. Note that the diagram is not invariant under 90 • rotations, thus showing that the cubic symmetry is broken.
It is interesting to consider why the allowed momenta lie along the positive and not the negative diagonal. The reason is due to our conventions: As mentioned in Section II A, the direction around the y-axis that we perform the isospin rotation is arbitrary. We choose to perform the rotation in the anticlockwise direction at every upper boundary of the lattice, which results in the favored direction being along the positive diagonal; we could just as easily have chosen to perform the rotation in the opposite direction at one or more boundaries, resulting in a change in the favored direction.
It is important to recognize that the choice of convention does not affect the boundary conditions on the pion wavefunction, which obeys antiperiodic boundary conditions in all cases. As a result, the pion energy remains invariant under cubic rotations, as we will demonstrate numerically in Section VIII C. We do however find that the rotational symmetry breaking has a measurable im- pact on the amplitudes is a bilinear operator. Such an operator can be constructed as where p 1 + p 2 = p and the ellipses denote spin and flavor structure that will be elaborated in the coming section. Note that both the ψ − ψ + and ψ + ψ − forms can be used providing p 1 and p 2 are chosen from the appropriate set of allowed momenta. In the aforementioned section we determine that the discrepancy in amplitudes can be substantially reduced by averaging the O + π ( p) and O − π ( p) forms. This observation is vital to constructing s-wave ππ states in our K → ππ calculation.

A. Local light-quark bilinear operators
As a result of the flavor mixing, many typical hadronic states (for example the proton) are no longer eigenstates of the Hamiltonian. For this work and in the measurement of the K → (ππ) I=0 amplitude we are only concerned with meson states. On the lattice we form bilinear operators in the quark fields that, when applied to the vaccuum state, create a linear combination of all mesonic states with the quantum numbers specified by the operator. In this section we consider bilinear operators in which both quark fields act at the same space-time point, and in Section V B we consider operators involving quark fields operating at different positions, for which the forms are further restricted by the requirement that the operators project onto eigenstates of the translation operator (and hence have definite momentum) in the context in which either of the quarks can independently cross the boundary, and hence change flavor.
In anticipation of including the strange quark, we henceforth refer to the fermion doublet comprising the light fields with subscript l.
There are three generic forms for a point bilinear operator involving only the light quark fields: where Γ is a generic spin-color matrix, a-c are c-number coefficients indexed with i ∈ {0..3}, and we have exploited the fact that any 2 × 2 complex matrix can be written as a linear combination of the Pauli matrices and the unit matrix (written here as σ 0 ).
From the color structure of the operators in Eq. (102) one can recognize that local, gaugeinvariant operators of the O 1 ll form can be obtained only if the flavor matrices σ 0 or σ 3 appear, while for O 2 ll and O 3 ll the matrices σ 1 or σ 2 must be used. Using Eq. (8) we can easily see that the G-parity odd operators are those for which the flavor structure anticommutes with σ 2 , i.e. σ 1 and σ 3 . Similarly, G-parity even operators are those that contain σ 2 or the unit matrix. In Table I we have compiled a list of local bilinear operators that are invariant under the above transformations and project onto states of definite G-parity eigenvalue.
From the table we can easily read off the operators that create states with the quantum numbers of the pion: It is illustrative to consider here the Wick contractions of the π 0 two-point correlation function π 0 (x)π 0 (y) . In addition to the usual connected diagram, this will also contain a diagram comprising two disconnected quark loops each of the form written is the two-flavor formalism as well as the standard form. Here Γ is a generic spin-color matrix..
In the limit of isospin symmetry such diagrams must vanish, which becomes clear if we return to standard notation but is not so obvious in this form. We can see that it vanishes in the notation of Eq. (104) as follows: where on the second line we have used the γ 5 Hermiticity, Eq. (43), of the Euclidean propagator and on the third line we have used Eq. (67).

B. Non-local bilinear operators
In this section we restrict our attention to operators of the form O 1 ll , although the concepts can be easily generalized to the other forms.
In practice we typically achieve better overlap with a chosen state of a particular momentum p + q using spatially smeared operators of the form where φ is some smearing function and Γ and Σ are arbitrary spin and flavor matrices respectively.
Here and below it is assumed that the spatial links at the time t are gauge-fixed such that this operator is gauge invariant.
To be useful, an operator O( p + q,t) specified by Eq. (106) must add the momentum p + q to the state to which it is applied. That is, it must be an eigenstate of the translation operatorT j in the spatial j th direction defined in Eq. (20): We will discuss how this can be done in the next two subsections.

Point operator
Let us first consider the special case of a local operator: φ (| x − y|) = δ x, y . Under a spatial translation this becomes (108) where x ′ = x +ĵ modulo the lattice size. We see that for Σ = σ 0 , B + j commutes with Σ and the above is translationally covariant providing p j + q j = n 2π L for integer n; this is just the operator we identified as being a G-parity even eigenstate in the previous section. Similarly for Σ = σ 3 , B + µ anticommutes and we require p j + q j = (n + 1 2 ) 2π L ; this corresponds to the G-parity odd operator. As a result, providing the correct momenta are chosen, the local operator is an eigenstate of the translation operator.

Non-local operator
For a non-local operator, there are additional terms where only one quark crosses the boundary that render the operator non-translationally covariant: the boundary term enters for one quark but not the other. In Section IV F we identified quark field operators,ψ( p,t), that are explicitly translationally covariant by virtue of their flavor-projector structure. Inserting these into Eq. (106), we obtain an operator that has a well defined momentum: where the integers n p and n q are associated with the momenta p and q respectively, via Eq. (100).
As B ± j commute with σ 2 , the restrictions on the allowed momenta for different choices of Σ are the same as for the local operator. The 1 2 (1 ± σ 2 ) projection operators enforce these same restrictions: For Σ = σ 3 , on commuting the first projector through Σ we find the projectors will cancel unless e in p π = e in q π and thus n q = n p + 2m for integer m. This implies p i = q i + 2mπ/L in G-parity directions, and thus i.e. this just restricts the total momentum to odd integer multiples of π/L, the allowed momenta of G-parity odd states. The corresponding condition for Σ = σ 0 requires −e in p π = e in q π and so n q = n p + (2m + 1), which in turn implies p i + q i = 2π L (n p + m + 1), the allowed momenta of G-parity even states.

Non-local neutral pion operator
For a neutral pion the local pseudoscalar operator is ψγ 5 σ 3 ψ and the corresponding non-local operator has the form: where translational covariance and the projection operators restrict the total momentum in each G-parity direction j ∈ G to p j + q j = (2m + 1) π L for integer m as desired. Let us examine the quark content of this operator more carefully: Expanding the parentheses in Eq. (111) we find two independent flavor structures, The second has unphysical flavor structure, ψγ 5 σ 1 ψ =dγ 5 Cū T + u T Cγ 5 d, and arises in cases where one of the quark fields crosses the boundary when the other does not. Naïvely this looks incorrect, but the form of the second term is, in a sense, artificial; given the translational invariance of the action, we are free to shift the boundary such that, for example, in one dimension, where stands for the ensemble average. Thus the apparently unphysical form of the second term is merely an artifact of our arbitrary choice as to where to place the boundary.

VI. THE STRANGE QUARK
We wish to simulate with a single strange quark whose discretized action is consistent with the charge conjugation boundary conditions on the gauge fields. The most obvious choice is to also impose charge conjugation boundary conditions on the strange quark, i.e. on crossing the boundary (115) s → Cs T ands → s T C .
Unfortunately, with this choice of strange quark boundary condition it is impossible to form a pseudoscalar operator that projects onto the K 0 state and that is invariant under translations such that it will create a K meson at rest. To see this, consider the two operators,sγ 5 d andūγ 5 s, that transform into each other under the boundary conditions: A linear combination of these operators cannot be chosen that is invariant under translation as, under G-parity, (117) which would simultaneously require α = β and α = −β . As a result we cannot create a kaon-like state of zero momentum. This is ultimately related to the fact that the strange quark with charge conjugation BCs is periodic in 2L, whereas the up/down quarks with GPBCs are anti-periodic in 2L and periodic in 4L.
A solution is to introduce a fictional degenerate partner to the strange quark, which we refer to as the s ′ quark, and to impose GPBC within that pair: and we can form eigenstates, with eigenvalue ±1, i.e. that obey either periodic or antiperiodic BCs. The former can be used to produce a stationary state whose physical component projects onto the neutral kaon and is therefore suitable for a K → ππ calculation.

A. Local heavy-light bilinear operators
In this section we consider local bilinear operators containing the strange quark and its fictional partner, s ′ . As with the light quarks, we write which is distinguished from the light-quark flavor doublet by the subscript h.
Below we use the following conventions for operators that create states with the quantum num-bers of the kaon: These transform under G-parity as follows: where we have denoted the fictional G-parity partners to the physical kaons with a prime ( ′ ) superscript.
For operators comprising only the heavy quark fields, the results obtained in Section V A also apply. For heavy-light bilinears we have four operator forms: As before, gauge invariance restricts O 1 hl and O 4 hl to the choices σ 0 and σ 3 , and likewise O 2 hl and O 3 hl are restricted to σ 1 and σ 2 . However, here the lack of symmetry under interchange of the fields implies that there is no restriction on the spin-color matrices Γ. We have again compiled the operators along with their standard forms and G-parity eigenvalues in Table II.
Consider the first line of Table II with Γ = γ 5 . We have This operator creates a stationary state whose physical component corresponds to theK 0 . Similarly the equivalent operator for the K + can be obtained from the sixth line of the table with Γ = Cγ 5 : .
From the table we can therefore read off the operators that project onto stationary states whose physical components correspond to the full set of charged and neutral kaons. We denote operators of definite G-parity quantum number with a tilde (∼), and display the quantum number in the subscript:K We also have moving states with the opposite G-parity quantum number: B. Operators acting on the physical kaon When measuring the K → ππ amplitudes or B K in the G-parity framework we are concerned with operators that act only on the physical kaon state and not the fictional partner. However, in order to determine matrix elements between physical states with known momenta we must work with the G-parity eigenstates that contain the fictional partner. As this transforms into the physical state when it crosses the boundary, we might expect it to make a non-trivial unphysical contribution to the measured amplitude. In order to analyze the size of this effect, consider the infinite-volume where O phys is chosen to be an operator that acts on the physical kaon, does not involve the unphysical s ′ quark operator and induces a mixing/decay to some final state |φ . Now we introduce the fictional partner to the strange quark, s ′ , which introduces a new state, |K ′ 0 , that is degenerate with |K 0 but has different flavor quantum numbers. (Here the infinite-volume states |K 0 and |K ′ 0 are QCD energy and momentum eigenstates.) This infinite-volume setup permits no mixing between these two states, hence We can therefore combine these equations in infinite volume, giving the result If we restrict ourselves to momenta whose components are integer multiples of 2π/L, the linear combination that we label |K 0 + is a valid eigenstate of the system in infinite volume and, if viewed as a two-component wave function, (K 0 , K ′ 0 ), will also obey G-parity boundary conditions for a finite volume of size V = L 3 . The infinite-volume composite particles K 0 and K ′ 0 receive only exponentially suppressed corrections when confined to move in a finite volume [11]. Thus, up to terms O(e −m K L ), we can obtain the infinite volume matrix element φ |O phys |K 0 that appears on the right-hand side of Eq. (131) from the matrix element on the left-hand side evaluated in finite volume. Thus, where the matrix element on the left-hand side is evaluated in infinite volume while that on the right in the finite volume V .
This is an instructive, highly simplified example of the more familiar relation between infinitevolume two-particle scattering states and finite-volume energies and matrix elements discovered by Lüscher [12] and by Lellouch and Lüscher [2]. If we consider the case where only the swave two-particle scattering phase shift is non-zero, the analysis in Ref. [12] can be understood as adjusting the relative momentum of the two-pions until the s-wave two-particle scattering state combines with the non-s-wave component of Lüscher's Helmholtz function to produce a function obeying the finite-volume boundary conditions. In our simpler case the momenta of the |K 0 and |K ′ 0 states are adjusted so that the combined state (|K 0 + |K ′ 0 )/ √ 2 obeys the required G-parity boundary conditions. The Lellouch-Luscher factor ∝ ∂ (δ 0 + φ )/∂ k of Ref. [2] corrects for the higher partial waves that are present in the finite volume ππ state and affect its normalization but do not contribute directly to the s-wave decay matrix element. The factor √ 2 in Eq. 132 plays a similar role, correcting for the shift in the normalization of the state caused by the otherwise irrelevant |K ′ 0 component of the finite-volume eigenstate, |K 0 + . In the discussion above, we have focused on the treatment of the kaon state that appears on the right-hand side of the matrix element in Eq. (131) when G-parity boundary conditions are present. Of course, the left-hand state φ | in that matrix element is also affected by these boundary conditions. In the case of the mixing parameter B K the K 0 | state appears on the left-hand side of that matrix element and can be treated in the same way as proposed for the |K 0 state on the righthand side. In the case that the left-hand state φ | is an interacting multi-particle final state there are additional contributions arising from the scattering, the form of which depends on the particles in question. In practice we are primarily interested in K → ππ decays, where the pions interact in the finite volume and individually obey antiperiodic boundary conditions. Here the contribution is completely described by the Lellouch-Lüscher formula [2,13] generalized to the antiperiodic case [14].

C. Strange quark determinant
With the introduction of the fictional s ′ quark, the theory now contains four flavors: two degenerate pairs, one light and one heavy. In addition to the effect of the fictional s ′ as a valence quark that must appear in states which are eigenstates of spatial momenta, that was discussed in the previous section, the s ′ quark will also appear as a sea quark through the fermion determinant of the heavy-quark Dirac operator which by definition acts on the s/s ′ quark doublet. However, in practice we wish to simulate a 2+1 flavor theory with a single strange quark species so the contribution of the s ′ quark must be removed.
We can represent the Dirac matrix in this heavy flavor space schematically as The determinant of this block matrix is In infinite volume the flavor mixing components vanish and the determinant reduces to the product of the determinants of the elements connecting individual quark flavors: where in the last equality we have used the degeneracy of the s and s ′ quarks. In this limit the contribution of a single flavor can be obtained simply by taking the square root of the two-flavor determinant.
At finite volume the Dirac matrix cannot be factored and the rooting prescription does not result in the determinant of a local operator. The non-locality of the resulting effective action then leaves no guarantee that the continuum limit of the rooted determinant defines a theory that lies in the correct universality class. This issue has received much attention in the context of staggered fermions, where each quark flavor corresponds to four 'tastes' which are coupled at high momenta and where a rooting prescription is used to reduce the number of quark tastes from four to one.
Fortunately our situation is easier to analyze.
We will consider two quite similar determinants. The first is the case of direct interest: an s/s ′ doublet with the s and s ′ species connected at a boundary by G-parity boundary conditions.
In this case the determinant does not factorize and taking the square root of that determinant produces an effective action which does not precisely correspond to that of a local field theory of fermions. The second case also involves a degenerate q 1 /q 2 doublet but each obeys independent charge-conjugation boundary conditions and the resulting determinant is a simple square whose square root corresponds to the Pfaffian of a local field theory of a single fermion species obeying charge-conjugation boundary conditions [8].
We will demonstrate that these two determinants differ by terms which fall exponentially as the system size grows so that the square root of the G-parity determinant differs from that of For simplicity, we have assumed that G-parity boundary conditions are imposed only in the x direction.
For the case of the doublet of fermions obeying charge-conjugation boundary conditions, the corresponding boundary terms are similar: We can now compare these two theories by expanding their respective determinants in powers of these boundary terms. These expansions will take the form of a sum of products of closed loops of lines will have the gamma and link matrix structure Γ + x CU x (L − 1) while those at which two fermion lines are absorbed will contain CΓ + x U x (L − 1) * . As suggested in Fig. 4 the leading term when L is large will involve fermion propagators which connect operators in two boundary terms at the same side of the volume. That is either x = 0 is connected with x = 0 or x = L − 1 is connected with x = L − 1. A propagator joining operators located at x = 0 and x = L − 1 will be exponentially suppressed. The connections that are shown in Fig. 4 correspond to a term which is not exponentially suppressed.
Inspecting Eqs. (136) and (137), we can recognize that the graphs describing terms of the form shown in Fig. 4, which are not exponentially suppressed, are identical for these two theories.
For the theory with charge-conjugation boundary conditions a graph such as that in Fig. 4 will correspond to a specific Feynman amplitude with a prefactor of 1/3 because of the symmetry of the graph under cyclic shifts of the vertices by two positions along the ring and an additional factor of 2 to account for the two independent flavors.
Exactly the same Feynman amplitude will appear in the G-parity theory with the same factor of 1/3 arising from the cyclic symmetry. The factor of 2 also appears because for the G-parity theory the propagators must alternate between that of an s and that of an s ′ quark as one moves around the ring. The exchange of s and s ′ everywhere yields the factor of two. The two minus signs which appear in Eq. (136) do not appear in Eq. (137) . However, as can be seen from the pattern of contractions which appears in Fig. 4 these minus signs always appear in pairs for the G-parity case and hence have no effect.
We conclude that if we neglect terms which are exponentially suppressed, our rooted s/s ′ determinant equals the Pfaffian corresponding to a local theory of fermions obeying charge conjugation boundary conditions. We therefore expect no subtle difficulties to arise from our use of the square root of this determinant.

VII. 2 + 1F DOMAIN WALL FERMION ENSEMBLES
In this and the following section we present results for full QCD simulations on a 16 3 × 32 lattice volume with G-parity boundary conditions in zero, one and two directions. We begin by discussing the generation of the G-parity gauge configurations and then present and interpret the results of a variety of measurements made on these ensembles.

A. Simulation Parameters and Generation
In Section III C we discussed the 'one-flavor' implementation of G-parity boundary conditions in a single direction whereby one simulates with a lattice of doubled extent and antiperiodic quark boundary conditions. This technique is very easy to implement in any lattice code library, but suffers from an additional factor of 2 n−1 increase in computational cost when applied in n > 1 G-parity directions. A much cleaner approach is simply to simulate with two separate quark fields that mix at the boundary; we refer to this as the 'two-flavor' approach. In practise this requires significant code modifications to handle the increased memory size of the two-flavor field and its unusual boundary condition, as well as careful implementation of the complex conjugate boundary conditions on the gauge fields. Nevertheless, using the one-flavor method as a cross-check, we implemented the two-flavor technique in the CPS++ library using the Bagel/BFM library [15] for optimized fermion inversions on IBM Blue Gene/Q machines [16]. An independent implementation in the Grid framework [17] has since been developed and cross-verified.
For our first fully dynamical simulations, we generated two 2+1 flavor domain wall ensembles with the Iwasaki gauge action at β = 2.13 (a −1 = 1.73(3) GeV [18]) and lattice size 16 3 × 32 × 16 with G-parity boundary conditions in one and two directions. We henceforth refer to these as the GP1 and GP2 ensembles respectively. For both ensembles we used an input light quark mass of m u/d = 0.01 and a strange quark mass of m s = 0.032. These parameters were chosen to match those of a previously generated ensemble with periodic boundary conditions, details of which were originally published in Ref. [19], although there a heavier strange mass of 0.04 was used; the ensembles were later extended with the value of m s = 0.032 used for this study, which is closer to the physical strange mass. For comparison with the G-parity ensembles we make use of the latter configurations and refer to this ensemble by the label GP0. With G-parity boundary conditions the Dirac matrix is intrinsically two-flavor and therefore the square of the determinant describes four flavors. As a result we must also use RHMC in the light-quark sector in order to obtain the two-flavor determinant; as this is an exact square there are no systematics introduced by this procedure. For the strange quark we must use the fourth root to obtain the contribution of a single flavor, and the implications of this rooting were discussed in Section VI C. In practise we found the use of RHMC in the light sector required more precise rational approximations than is is typically required for the strange quark in order to obtain good Metropolis acceptance. This was not a severe hindrance for these cheap ensembles, but for our physical-point ensembles that we generated to measure the K → (ππ) I=0 amplitudes [5] we found the associated linear algebra overheads were significant and limited the value of introducing multiple Hasenbusch masses to reduce the inversion cost as well as making it difficult to take advantage of mixed-precision techniques.
Note that, although it was not used for this analysis, the RHMC algorithm for the light quarks can be replaced by the "exact one-flavor action" developed by TWQCD [20,21], whereby a Hermitian and positive-definite operator is constructed describing a single flavor (or two flavors with GPBC) of Wilson/domain wall fermions. In Ref. [22] we describe how to implement this technique in a highly efficient manner and demonstrate a factor of 4.2× improvement in the speed of generating the G-parity ensembles for our K → ππ analysis.
The ensembles used for the analysis in this document were generated using three layers of

C. Ensemble properties
For each of the three ensembles we measured the plaquette, chiral and pseudoscalar condensates after every trajectory. We also measured the topological charge after every fifth trajectory starting from 500 MD time units of evolution. We plot the Monte Carlo time histories of these quantities in Figure 5. From the plaquette, chiral and pseudoscalar condensates we measure the accumulated integrated autocorrelation times as a function of the MD time separation, the results of which are shown in Figure 6. Here the errors were obtained by binning the correlations C(t,t + ∆) with a fixed separation ∆ over neighboring configurations, using the first technique described in Section II.D of Ref. [23]. We also measured the autocorrelation using the topological charge measurements, but found that these measurements, performed every 5 MD time units, were sufficiently separated that no observable correlations were found within the statistical errors; we therefore excluded these data from the figures. The figures suggest integrated autocorrelation times of ∼15 MD time units for the G-parity ensembles, and ∼10 for the GP0 ensemble, such that there are 30 and 20 MD time units between independent measurements, respectively.
The expectation values of the plaquette, chiral condensate and pseudoscalar density are given in Table III. Here the errors are determined after binning uniformly over 40 MD time units. We observe excellent consistency between the plaquettes, and the pseudoscalar density is consistent with zero indicating good topological sampling. The chiral condensate agrees well between the GP0 and GP1 ensembles but differs between the GP0 and GP2 ensembles by 1.9(7)% level. While this may be simply statistics, as supported by the agreement between the GP0 and GP1 ensembles, this difference may also be attributed to the explicit breaking of the flavor non-singlet axial symmetry by the boundary conditions that we identified in Section IV C. We would expect that any change in the chiral condensate would also manifest directly in the pion energy, which is related to the chiral condensate at leading order in chiral perturbation theory (χPT). In Section VIII A below, we find that the pion energy on this ensemble is indeed different than expected at the 1.5% level, but the value obtained is lower than expected, in the opposite direction to that suggested by χPT for an upwards shift in the chiral condensate. We therefore conclude that the observed effect is likely statistical in nature or that a more elaborate analysis using finite-volume χPT is needed.

VIII. RESULTS
In this section we present measurements of the residual mass, the axial-current renormalization factor (via the Ward-Takahashi identity), the pion and kaon energies and decay constants, and the neutral kaon mixing parameter, B K . These quantities are computed on each of the three ensembles discussed above, and compared to examine the effects of the boundary. semble we use momenta p = ± π 2L (1, 0, 0) and on the GP2 ensemble we use p = ± π 2L (1, 1, 0). The sources are placed on time-slice 0 on all configurations, and we generate propagators with both periodic (p) and antiperiodic (a) temporal boundary conditions. From these we take linear combinations F = p + a and B = p − a that project out the forwards and backwards moving components of the propagator respectively, effectively doubling the temporal periodicity and significantly reducing round-the-world finite-time-extent effects.
For the G-parity ensembles we must construct the 2 × 2 flavor matrix propagators, with elements (138) where f -i are flavor indices,Ṽ is the gauge-fixing matrix, and η is a generic, flavor-matrix-valued source. For a wall source, η is just the unit matrix for all sites. In order to obtain Θ one can simply invert propagators for both source flavor indices separately, i.e. construct Θ f 1 and Θ f 2 (for all f ) with two inversions of the two-flavor Dirac matrix. With this approach, one must perform four inversions in order to both compute Θ( p) and its parity partner Θ(− p).
In most cases we can use the exact relation between the propagator and its complex conjugate, Eq. (67), to obtain both Θ( p) and Θ(− p) with only two inversions of the Dirac matrix as follows: Taking the complex conjugate of the solution vector, For any source that has the property Cγ 5 η * ( y) = η( y)Cγ 5 (which trivially applies to any source proportional to the unit matrix, such as the wall source in question), the above reduces to Multiplying from both sides by σ 2 and examining the equation in component form, we find Notice that we can obtain the first column of the left-hand side, i.e. those elements with source index 2, from the first column of the right-hand side, which contains only elements with source index 1. Using this relation we can obtain both Θ( p) and Θ(− p) by computing only Θ f 1 ( p) and Θ f 1 (− p), from which we obtain Θ f 2 (− p) and Θ f 2 ( p), respectively.

A. Pion two-point function
Using the momentum-eigenstate field operators obtained in Section IV we derived the appropriate operator for a neutral pion state in Eq. (111). Taking p i = π/2L for each G-parity direction (i.e. n p = 0 in Eq. (99)), and utilizing a wall source and point sink, we obtain the following form for the pion two-point function: where Θ l are the light-quark solution vectors defined above. The LW superscript indicates that the correlation function has a local (point) sink and wall source. At the sink location we have not included the (1 − σ 2 ) projector as it is not necessary for a local operator. We remind the reader that Wick contractions resulting in disconnected loops at the sink location vanish because of the γ 5 -hermiticity of the propagator and its relation to its complex conjugate (cf. Section V A).
It is straightforward to show that the last line of Eq. (142) can also be obtained by considering the charged pion operators given in Eq. (103), which is necessary due to the exact isospin symmetry of the formalism. We therefore generically assign the source/sink operator with the label P for pseudoscalar.
The corresponding pion two-point function for the GP0 ensemble is, For all ensembles we compute the correlation functions with both the forwards (F) and backwards (B) moving propagators and average the results (after the appropriate temporal reflection) to improve statistics. As the amount of data is large, we were unable to accurately resolve a correlation matrix and therefore performed uncorrelated fits (i.e. with a diagonal correlation matrix) here and for the other results presented in this document.
For use below we also fit the correlation functions with axial-vector sink operators appropriate for neutral pions, is real and positive. With this convention the correlators with spatial axial-vector sinks are pureimaginary, which we will show is necessary to fulfill the axial Ward-Takahashi identity.
We also compute the wall-source, wall-sink PP correlation function, which takes the following form on the G-parity ensembles: where again we have suppressed the gauge fixing matrices at source and sink.
On each ensemble we simultaneously fit the various correlation functions with a common exponent. We also include the point-sink two-point function with the A 4 operator at the sink and the Hermitian conjugate of this operator (ensuring a positive sign) at the source, in order to improve the signal for the exponent. The data were fit to the following forms, where O 1 and O 2 are the sink and source operators respectively, T = 32 is the lattice temporal length, and N are the amplitudes. We use the notation P for the pseudoscalar operator and A µ for the axial operators. The sign of the backwards propagating component (the second term in the brackets) depends on the how the operators transform under time reflection: The PP and A i P, for i in a spatial direction, require the plus (cosh-like) form and the A 4 P the minus (sinh-like) form. Note the temporal length is doubled due to our use of the F and B combinations of temporal boundary condition.
The fitted masses and amplitudes that we obtain are given in Table IV, alongside the fit ranges chosen by eye based on effective mass plots (used uniformly for all correlators on a given ensemble) and the uncorrelated χ 2 /dof. A far more precise measurement of m π on the GP0 ensemble was performed in Ref. [24], which we also include in this table. In Figure 7 we show effective mass plots for the PP point-sink channel.
In Table IV we also list the predicted values of the energy that one obtains by applying the continuum dispersion relation, where n is the number of G-parity directions, L = 16 is the spatial box size, and √ nπ/L are the magnitudes of the expected pion momenta. Here m π is the fitted mass obtained from the GP0 data in Ref. [24]. From the table we observe good agreement between the predicted and measured energies on the GP1 ensemble. The value on the GP2 ensemble is 1.8(1.0)% smaller than the predicted value, although if we replace the continuum dispersion relation with the following lattice expression, E π = m 2 π + n sin 2 (π/L), the discrepancy drops to 1.4(1.0)% which is statistically compatible with zero. On the other hand it is possible that this effect is related to the larger value for the chiral condensate observed on this ensemble that we discussed in Section VII C.

B. Signal-to-noise ratio of pion two-point functions
An interesting observation from Figure 7 is that the signal-to-noise ratio for the G-parity ensembles appears to reduce with time, and degenerates more quickly as we increase the number of G-parity directions. This will unfortunately hamper our ability to measure the ππ energy and K → ππ amplitudes with G-parity boundary conditions. The origin of the falling signal-to-noise ratio can be understood via Lepage's argument [25]   and wall-source-wall-sink correlators respectively. For the A 1 P and A 2 P correlators, the amplitude N is of the imaginary part of the correlator, and for the remainder it is of the real part. The value of E π on the fifth line is that measured in Ref. [24]. E pred π is the predicted pion energy obtained by applying the continuum dispersion relation to this, more precise value.
on a gauge configuration i in an ensemble of size N, the standard error on the mean, σ , of the can be obtained as follows: hence the error must grow as the square-root of the expectation value of the square of the Green's function. The signal to noise ratio S/N, where S = | O | and N = σ is Here E 2 and A 2 are the energy and amplitude of the lightest state in the squared correlation function OO † , and E 1 and A 1 are those of the state that we are measuring (in this case the pion). We can see that the time dependence of the signal-to-noise ratio drops out when E 2 = 2E 1 , and that the signal-to-noise ratio will decrease with time when E 2 < 2E 1 . For the two-point functions in the previous section, the Green's function OO † describes the propagation of two quarks and two antiquarks. On the GP0 ensemble, the ground state of this system is just two pions at rest, hence the S/N is constant. On the G-parity ensembles this state does not exist: instead the lightest two-pion state contains only moving pions, and we again have We might therefore expect that S/N would also remain constant for these ensembles.
However, the four-quark Green's function also contains a contribution from two pseudoscalar SU(2) flavor singlet particles with quark contentūu +dd. As only the connected components are present in the noise, this state has the same energy as the stationary pion on the GP0 ensemble, and the signal-to-noise will remain constant. However, on the G-parity ensembles this state is Gparity even and will therefore have the energy of the stationary pion rather than the moving pion, such that the states entering the noise are lighter than the moving pions in the signal. To confirm this we measured the energy of the flavor-singlet, whose operator isψγ 5 ψ in our formalism (cf.  Note that if we replace one of the light quark propagators with that of the strange/strangeprime quark pair (which differ only in their input mass), the contractions for the flavor singlet are identical to those of the G-parity kaonic state,K 0 + (see below). As a result the flavor singlet might also be considered a "light kaon" in this sense.
To demonstrate that it is indeed this flavor singlet state that dominates the noise we performed single-exponential fits to the signal-to-noise ratio of the PP wall-source-point-sink correlation function on all three ensembles; the results are given in Table V and the data in Figure 8. We observe excellent consistency between the predicted and measured exponential falloff of the signalto-noise ratio on both the GP1 and GP2 ensembles. On the GP0 ensemble however the LePage argument predicts a constant behavior but instead we observe a significant exponential falloff in S/N suggesting a state of energy 0.370(32) appears in the noise. This discrepancy has not been understood but may result from the presence of multiple heavier states contributing with different signs or the absence of disconnected diagrams.

C. Effects of quark-level rotational symmetry breaking
In Section IV G we observed that on a lattice with GPBC, the quark's momentum components in orthogonal directions must differ only by integer multiples of 2π/L. As a result the momentum (1, 1, 0) π 2L for GPBC in two directions is allowed but (1, −1, 0) π 2L is not. Given that these two momenta are related by a cubic rotation, the above restriction constitutes a breaking of the cubic rotational symmetry at the quark level.
In this subsection we will describe a numerical investigation of these rotational symmetry and the Iwasaki+DSDR gauge action with β = 1.75 that was generated as part of the preparation for our K → ππ project [5]. Aside from the lattice size, Möbius parameter and a heavier pion mass of 440 MeV, this ensemble is otherwise identical to the 32 3 × 64 × 12 ensemble used for the K → ππ calculation, and has G-parity BCs in three directions. Instead of the Coulomb gauge fixed wall sources used elsewhere in this document, this investigation was performed using hydrogen wavefunction smeared sources and the correlation functions were computed using the all-to-all propagator technique [26] used for the K → ππ calculation. In this technique an approximation to the all-to-all propagator is generated by combining exact low eigenmodes (we used 100 eigen-  modes here) computed using the Lanczos algorithm with a stochastic approximation to the high mode component. We perform full time, spin, color and flavor dilution such that the random numbers are required only to induce a delta function in the spatial coordinate. We use a single random hit per configuration such that the delta function is generated under the ensemble average.
where Θ(| x − y|) = exp(−| x − y|/r), with radius r = 2, is the smearing function, and p = p 1 + p 2 is the pion momentum, with p 1 and p 2 chosen from the appropriate set of allowed momenta. For each operator and pion momentum we choose p 1 and p 2 as those that have the smallest magnitude subject to the constraints; the chosen momenta are given in Table VI. Here in each case the O − π form contains the lower magnitude quark momenta, and is therefore the natural primary choice of operator. (Note that for the parity partners of these four momenta, the O + π operator would be the natural choice.) We fit the correlation function to the form (153) C ± (t) = N ± π e −E ± π t + e −E ± π (L T −t) .  Table VII. We observe that the pion energies are all consistent as expected, but for both operators the amplitudes fall into two distinct groups: those of the three momenta perpendicular to the (1, 1, 1) diagonal are consistent, and that of the momentum parallel to the diagonal is ∼16% larger (about 9σ ) for the O − π operator and ∼4.5% smaller (about 3.5σ ) for the O + π operator, respectively. Thus we have clear and unambiguous evidence of the rotational symmetry breaking. The three momenta that have the same amplitude can be recognized as being related by the discrete group of 120 • rotations around the diagonal of the spatial box, which can readily be seen as a residual symmetry of the action from the discussion in Section IV G.
As part of this study we also measured the correlation function of the averaged operator The amplitudes and energies of the corresponding fits are given in Table VIII. We observe that in performing this average we significantly reduce the amount of rotational symmetry breaking, and in fact the amplitudes all now agree within statistics.
In order to test this further we also generated pion two-point functions of higher momentum  Table IX. To better discern the pattern of the symmetry breaking for this noisier data, we perform a simultaneous fit over all 12 correlation functions with a common energy. The fit is performed to  Table X. We observe that the amplitudes for each operator clearly fall into three discrete groups: two groups of three whose momenta are related by 120 • rotations about the (1, 1, 1) diagonal, and one group of six whose momenta are related by 60 • rotations about the same axis. These groups are recognizable as the sets of points that lie on discrete planes whose normals are parallel to this diagonal. In Table XI we show the results for the averaged operator O avg π where we again find that the symmetry breaking is heavily suppressed.
We make use of these averaged pion operators in Ref. [5] when constructing ππ operators that are intended to transform under particular representations of the cubic symmetry group. Details of the operators and further tests will be presented in a forthcoming paper.

D. Residual mass
The residual mass m res is a measure of the explicit chiral symmetry breaking in the domain wall formalism due to the finite extent of the fifth dimension. It acts as an additive renormalization of the bare quark masses,m = m + m res , and it is the limitm → 0 that defines the chiral limit. We make use of m res below when computing the pseudoscalar decay constants, and when testing the Ward Takahashi identity, hence it is important to understand the effects of the G-parity boundary conditions on this quantity.
The residual mass is measured by fitting the ratio to a constant in time, and extrapolating the results to the massless limit. Here J a 5 is the local where Ψ are the 5d fields, q the 4d surface fields and τ a are the generators of the flavor symmetry (SU(2) in our case). J a 5q is also a local pseudoscalar density, but defined from the midpoint fields in the fifth direction: (157) J a 5q (x) = −Ψ(x, L s /2 − 1)P L τ a Ψ(x, L s /2) +Ψ(x, L s /2)P L τ a Ψ(x, L s /2 − 1) .
Up to a constant we can compute m ′ res for a neutral pion (τ a = σ 3 /2) on the GP0 ensemble via the following (wall source) correlation function and the equivalent for J 5q in which the propagator sink is evaluated on the appropriate s-slices.
Here we have used the degeneracy of the quarks to relate the up and down quark propagators, and we have not shown the Coulomb gauge fixing matrices on the source timeslice.
We would, of course, obtain the same expression as the above for any of the three pion states due to the exact isospin symmetry, but the neutral pion can be created by a particularly convenient operator in our G-parity formalism, giving the following expression for C J 5 with GPBC: (159) where ψ are the usual two-flavor G-parity fields and we are using the same symbol G l for the one-flavor and two-flavor light quark propagators in Eq. (158) and (159). Note that we have not applied the 1+σ 2 flavor projection operator at the source here; the result is that we potentially have a poorer projection onto the desired moving pion state. Nevertheless, the ratio m ′ res is common to all pseudoscalar states, and therefore can be obtained just as well from a linear combination of pseudoscalar states as it can be from just the ground state.
We compute m ′ res using just the light-quark (m l = 0.01) propagators with antiperiodic temporal boundary conditions. For simplicity we use a uniform fit range of t = 4-28 and we perform uncorrelated fits. The values that we obtain are listed in Table XII. We observe very good agreement between the GP0 and GP2 ensembles, but the value on the GP1 ensemble is slightly lower by ∼2σ .
Examining this in more detail, we plot m ′ res as a function of time for the GP0 and GP1 ensembles in Figure 9. Here we see no evidence of any systematic deviation between the ensembles, suggesting the discrepancy is merely due to statistical fluctuations. In principle we would not expect m res to have any dependence on the boundary, it being simply a ratio of amplitudes separated in the fifth dimension, and our results are in line with this expectation.  Given that we measure with only a single valence quark mass, we cannot extrapolate to the massless limit. However, the result measured on the GP0 ensemble agrees very well with the value m ′ res = 0.003102 (25) quoted in Ref. [19]. Recall these measurements were performed with a 25% higher sea strange mass; this suggests the sea strange mass dependence is very weak, and therefore that we can use the value in the chiral limit of m res = 0.00308(4) given in that paper as our final value for the residual mass on all three ensembles.
where J 5 and J 5q were defined in the previous section. Here {λ a } are N f × N f flavor matrices, and a and b are flavor indices (of the conventional sort) .
is the discretized backwards derivative, withμ a unit vector in the µ-direction. The partially-conserved current A µ is composed of the 5d domain wall fields Ψ, and has the "point-split" form, i.e. it is defined along the link between x and x +μ rather than on the site x. This current can be related to the local axial current A a µ =qγ µ γ 5 λ a q, comprised of the domain wall surface fields q as follows where the relation becomes exact in the continuum limit. In the above, Z A is the renormalization factor relating the local current to the continuum-normalized Symanzik currents, and Z A is the same for the partially-conserved current (for further information we refer the reader to Ref. [28]).
Combining Eqs. (160) and (162) we obtain where we have used Eq. (155) to remove the J 5q terms. Inserting a complete set of states and taking the large time limit, the lhs of Eq. (163) becomes where we have allowed the ground state pion to have non-zero three-momentum p and corresponding Euclidean four-momentum p µ = (−iE p , p). Inserting the above into Eq. (163) and taking the spatial Fourier transform, we obtain Noticing that, in the large time limit, the source matrix element π( p)|J b 5 |0 appears on both sides of the equation, and that the time-dependence e −E p t is identical, we can transform the above into an expression that acts upon the ground-state amplitudes of general point-sink two-point functions with a pseudoscalar source: where S represents a general source smearing. We can therefore utilize the wall-point amplitudes determined in Section VIII A and m res determined above to obtain the relative normalization of the local and partially-conserved domain wall axial currents, Z A /Z A .
Notice that the rhs of Eq. (166) is pure-real, whereas the lhs contains a factor of i for the spatial components of the sum. This implies that the spatial amplitudes must be pure-imaginary, which we previously found to be the case. The expression therefore reduces to We combine the amplitudes in Table IV with m res obtained in the previous section to obtain Z A Z A from the above expression. We use a = b = 3, i.e. the third isospin component of the pion triplet.
Note that N LS J 5 P = −N LS PP where the minus sign occurs because in the latter we compute The values that we obtain are listed in Table XIII. The value measured on the GP0 ensemble agrees very well with Z A /Z A (m l = 0.01) = 0.71807 (14) given in Ref. [19], but those obtained from the G-parity ensembles are ∼2% higher. This discrepancy may arise due to the boundary conditions but also due to differing discretization errors as the G-parity Green's functions are evaluated at nonzero momentum. Nevertheless the effect is small. For use below, we take use Z A /Z A = 0.7162 (2) obtained in the massless limit in Ref. [19] for all ensembles.

F. Pion decay constant
The decay constants for a pseudoscalar meson P are defined on the lattice as [29,30] where the states are normalized as f P can be obtained from the following ratio of the pseudoscalar correlator amplitudes computed in Section VIII A: Here Z A is the renormalization factor relating the local domain wall current to the continuum rather than the ratio Z A /Z A that we obtained in the previous section. However, as can obtain a good approximation to f P using this ratio in place of Z A in the above.
For the pion we label the decay constant f ll . In Table XIV we list our measured values for this quantity, computed on all ensembles using the amplitudes and energies given in Table IV. We measure using the temporal (µ = 4) component of the axial-vector operator at the sink, and also using the non-zero-momentum spatial components on the G-parity ensembles (for which q µ = π/L). As expected we observe excellent agreement between the results on all three ensembles and between spatial and temporal determinations. and sink momentum directions ±2 p on the GP1 ensemble. The first column gives the sign of the source flavor projection, and the second the sign of the sink momentum.

G. Impact of the flavor projection
It is interesting to observe the effect of the the 1 ± σ 2 flavor projection for non-local sources.
To do so we consider the pseudoscalar point-sink correlation functions where we vary both the sign in the source projection operator and the sign of the sink momentum projection. We fit each independently to obtain the amplitude and ground-state energy. For simplicity we consider only the GP1 ensemble (with p = (π/2L, 0, 0)) and use a common fit range of 5-25.
We list the fit results in Table XV. The first line of the table gives the values with the 'correct' flavor and sink momentum projections, i.e. for which the sink momentum projection equals the total source momentum, and we use the correct 1 + σ 2 projector to make the source translationally covariant. The amplitude and energy for this correlation function agree within statistics with those obtained in our simultaneous fit and listed in Table IV. If we attempt to project onto the opposite sink momentum we find an amplitude (second line of the table) statistically consistent with zero as expected due to momentum conservation.
For the 'wrong' source flavor projection we observe from the third line of Table XV that this source also couples strongly to the forwards-moving pion state with momentum π/L. However we find (fourth line) that it also projects onto the backwards-moving pion with momentum −π/L, despite our explicit source momentum projection onto the forwards-moving state, seemingly violating momentum conservation. To see how this comes about, recall that the 1 − σ 2 operator projects out the ψ − field, for which the allowed momenta are q ∈ π 2L {. . ., −9, −5, −1, 3, 7, . . .} = −π(1 + 4n)/2L in each G-parity direction, where n is an integer.
In the following discussion we suppress the coordinates in the y and z and t directions. The phases exp(iπ[1 + 4n]x/2L) associated with the allowed momenta for ψ − , form an orthonormal basis: We can therefore define a Fourier transform and its inverse into a momentum space indexed by n: The Fourier transform of the quark field ψ − (x) with the non-allowed momentum π/2L can then be written as we observe that the coefficient is generally complex and is non-zero for all n. The same will be true for the ψ + field operator, which also has the same set of allowed momenta.
The source operator can then be written as which has non-zero projection onto all states of momentum −[2 + 4n + 4n ′ ]π/2L. This includes the pion state with momentum −π/L, for which n = −n ′ . For example, with n = n ′ = 0 we have which has a non-zero imaginary component. We therefore predict that the correlation function with source projection 1 − σ 2 and sink momentum sign −1 should also have an imaginary component.
Performing the fit to our data we indeed find a statistically-resolvable imaginary component with amplitude which is pure-real. We therefore predict that we will not observe an imaginary component for the correlation function with source projection 1 − σ 2 and sink momentum sign +1. In practise we found that we could not fit to a cosh-like form to the imaginary component (the fitter did not converge) unless we performed a simultaneous fit including the real-part with a shared energy; there we found an amplitude for the imaginary component of 1.7(3.9) × 10 2 , consistent with zero as expected. We can fit the imaginary part alone to a constant, for which we obtain a value 0.68(85), again consistent with zero.
We conclude that it is vital to perform the flavor projection if one is concerned with the specific direction of the resulting pion's momentum. This is important, for example, when constructing a two-pion state with zero total momentum.

H. Kaon mass and decay constant
We obtain the mass and correlator amplitudes for the neutral kaon on the GP0 ensemble and those of the mixed state |K 0 + = 1 √ 2 |K 0 + |K ′ 0 , discussed in Section VI A above, on the Gparity ensembles.
On the GP0 ensemble we compute wall-point correlators of the form using our Coulomb gauge fixed wall source propagators. Likewise, on the G-parity ensembles we compute (179) C LW Recall that the mixed combination |K 0 was introduced in Section VI A as the eigenstate of the QCD Hamiltonian which is an energy and momentum eigenstate with the energy of the kaon and zero spatial momentum. In order to obtain a zero-momentum state under the constraint that the quark momenta are odd-integer multiples of π/2L, we assign momentum + p to the anti-light quark and − p to the strange quark as indicated above, where p i = π 2L for each G-parity direction and zero otherwise. Here we must be careful with the momentum projection required to create a translationally covariant source operator: The momentum − p (n p = −1) of the strange quark requires the projection ψ h → 1 2 (1 − σ 2 )ψ h , and the momentum + p of the anti-light quark requires ψ l →ψ l 1 2 (1 − σ 2 ) (cf. Eq. (96)). If we swap the momentum assignment we must also swap the sign of the projection.
For the kaon wall-sink two-point function we measure the following: Computing the decay constant for the physical kaon requires a little more thought. The continuum temporal axial-vector operator that annihilates the K 0 is A 4 = −isγ 4 γ 5 d, where we have chosen a phase convention such that 0|A 4 |K 0 is real and positive for the operators specified in Section VI A. From the G-parity fields ψ h and ψ l we can construct such a bilinear operator (here for a generic spin-matrix Γ) as (181) sΓd = ψ h F 11 ΓF 11 ψ l where F 11 = 1 2 (1 + σ 3 ), such that A 4 = −iψ h F 11 γ 4 γ 5 ψ l . The axial-sink pseudoscalar source correlator is therefore where the coefficient of √ 2 compensates for the fact that only the physical component of the incoming state couples to the operator (up to exponential corrections) as explained in Section VI B.
(Note, as discussed in that section we have limited the observable operator appearing in Eq. 182 to one containing only the s-quark.) We can also compute the matrix element with a similar axial operator that connects to the unphysical K 0 ′ state, . The operators A 0 and A ′ 0 interchange under the G-parity operation, hence the G-parity symmetry of the action implies A 4 |K 0 = A ′ 4 |K 0 . This is easily verified numerically. In order to improve statistics we can therefore take the average of the two.
The data are fit to the functional form given in Eq. 147. The fitted masses and amplitudes as well as the corresponding values of the decay constant f K are given in Table XVI, alongside the fit ranges chosen by eye based on effective mass plots (used uniformly for all correlators on a given ensemble) and the uncorrelated χ 2 /dof. In Figure 10 we show effective mass plots for the PP point-sink channel. The errors shown in the table are only statistical but appear to be sufficient to explain any discrepancies between the results found for m K and f K among the three ensembles.
We observe good agreement between the values of both the masses and decay constants of all three ensembles, suggesting that masses of the degenerate K and K ′ particles are not significantly altered when they are allowed to mix by the boundary and that we are able to extract the decay constants of the physical kaon despite this mixing.
Here we have made use of Eqs. (67) and (43) which are identical in form up to the presence of the additional flavor matrix F 11 , further demonstrating the utility of our notation.
The expressions given in Eqs. (185) and (186) assume that the kaon state is created by a point source at (x 1 ,t 1 ) and absorbed by a point sink located at (x 2 ,t 2 ). In a modern calculation of neutral kaon mixing, improved source and sink wave functions are used which better project onto the kaon state and which fix the momentum of the kaon to be zero in order to reduce systematic errors arising from the contributions of excited or moving kaon states. Specifically we might use Coulomb gauge fixed wall sources and sinks for the light and strange quarks introducing explicit position-dependent phase factors and the associated projection operator (1 ± σ 2 ) to give these quarks the momenta required by the boundary conditions. If we do not show the Coulomb gauge fixing matrices, aK 0 interpolating operator at the time t would be written as where as before the momentum p has components +π/2L for directions in which G-parity boundary conditions are imposed and zero otherwise.
The bag parameter B K is constructed from the Green's functions containing the O VV+AA operator as follows: where the denominator serves to divide out the normalizations of the source and sink operators.
On the G-parity ensembles the K 0 operators in the above are replaced byK 0 operators and a factor of 2 is required in the numerator, as described above.
Following Ref. [31] we set t 2 = T and t 1 = 0 and utilize the forwards (F) propagators (obtained by using quark propagators which are the sum of propagators obeying periodic and anti-periodic boundary conditions in the time) to form the kaon between t 1 and t such that it falls off exponen- Similarly we form the kaon between t and t 2 from the backwards (B) propagators (obtained by using quark propagators which are the difference of propagators obeying periodic and anti-periodic boundary conditions in the time) such that it falls off exponentially in . Similarly the A 4 (t)K 0 (t 1 ) term in the denominator of Eq. (188) is constructed with the F propagators and the A 4 (t)K 0 (t 2 ) from the B propagators. The exponential time dependence then cancels exactly between the numerator and denominator such that B lat K (t) is constant up to excited state contamination. We will next present the results from such a calculation of B K . We begin by explaining that this calculation was carried at the start of our study of G-parity boundary conditions and the flavor projection matrix 1 2 (1 − σ 2 ) shown in Eq. (187) was not included. This omission might lead to the presence of additional kaon states with non-zero momenta in both the numerator and denominator of Eq. (188), resulting in systematic errors. However, for our particular operators and kinematics, we can argue that such effects are at or below the 1% level.
First we observe that the operatorK 0 (t) defined in Eq. (187) and its conjugateK 0 (t) have positive G-parity. Therefore, they obey periodic boundary conditions in each spatial direction and can only create states whose momentum components are integral multiples of 2π/L. Second, both the weak mixing operator O VV+AA and the axial current operators which appear in the numerator and denominator of Eq. (188) are summed over the full spatial volume. Since these operators do not contain a translationally invariant combination of our four quark flavors, this sum over space is insufficient to ensure that these operators cannot create or destroy momentum. However, for the case of momentum components that are multiples of 2π/L, (i.e. p = 2π n/L where n is a vector of integers) these spatial sums will give zero unless n = 0. Thus, for the factors in the denominator, only zero-momentum kaons can contribute. Similarly the spatial sum in the numerator guarantees that the momentum of the initial and final kaon must agree. Third, such an unwanted kaon with non-zero momentum traveling from the source to the sink will be suppressed relative to the kaon with zero momentum by the factor exp{−( M 2 K + (2π/L) 2 − M K )T } = 0.0027 for our T = 32, supporting our assertion that the resulting errors are at or below 1%.   In Table XVII we list the hand-chosen fit ranges, the fitted values of B K and the associated χ 2 /dof. Plots of B lat K (t) are given in Figure 11. We observe excellent consistency between all three ensembles, demonstrating our ability to compute physical kaon matrix elements using the G-parity mixed kaon states.

IX. CONCLUSIONS
This document contains a thorough investigation of the use of G-parity spatial boundary conditions (GPBC) [6][7][8] in zero temperature lattice simulations. The G-parity operation is a combination of charge conjugation and an isospin rotation by π radians around the y-axis under which the charged and neutral pions are all odd eigenstates. GPBC are therefore equivalent to antiperiodic boundary conditions for the pions and the allowed discretized pion momenta become odd-integer multiples of π/L, where L is the lattice spatial extent. This results in the removal of the stationary pion state from the spectrum, making this technique ideal for studying interactions of moving pions as it is no longer necessary to battle to isolate the weaker excited-state contribution involving moving pions from the dominant ground-state contribution. This is similar to the technique of applying antiperiodic boundary conditions to just the down or up quark [7], which also results in antiperiodic charged pions and has been used successfully to compute the ∆I = 3/2 K → ππ decay amplitude [4,14]. GPBC are superior however as they also give rise to moving neutral pions, and preserve the isospin symmetry which is explicitly broken if different boundary conditions are applied to the up and down quarks.
We demonstrated this technique numerically in Section VIII by performing measurements of a variety of quantities including the pion and kaon masses and decay constants, the bag parameter B K , as well as the domain wall residual mass and axial current renormalization, using two customgenerated 16 3 × 32 × 16, 2+1 flavor domain wall ensembles with the Iwasaki gauge action at β = 2.13 (a −1 = 1.73(3)) and G-parity boundary conditions in one and two spatial directions.
These were compared with the same quantities computed on an identical lattice but with periodic boundary conditions in all spatial directions. The generation of custom ensembles is necessary here as gauge invariance requires the gauge fields to obey charge conjugate (complex conjugate) boundary conditions (cf. Section III). Using these ensembles we show, for example, that the pion ground-state energies indeed obey the continuum dispersion relation as the number of G-parity directions increases, and that we can extract decay constants from both the spatial and temporal axial current operators that agree in all cases.
The implementation of GPBC is complicated by the flavor mixing under the isospin rotation at the boundary. In Section III C we demonstrated that GPBC in a single direction is equivalent to antiperiodic boundary conditions on a lattice of doubled spatial extent in the G-parity direction where we force the gauge links on the second half to be complex conjugates of those on the first.
While this provides a convenient implementation, extending this technique to multiple G-parity directions requires either further (unnecessary) doubling, increasing the computational cost, or applying unusual non-local boundary conditions in the directions perpendicular to the first doubling.
We concluded that a general implementation of the lattice action with an explicitly two-flavor fermion field that mixes at the lattice boundary is more practical, and it is this approach that we use for our later measurements and in Ref. [5].
In Section IV we investigated the symmetries of the lattice action. The mixing of quark flavor at the boundary explicitly breaks baryon number conservation as, for example, it transforms a proton (uud) into an anti-neutron (ddū). There is also a more subtle symmetry breaking of the flavor non-singlet axial current, which was originally recognized by Wiese [6]. This potentially gives rise to a change in the value of the chiral condensate and a corresponding shift in the energy spectrum of the pion states, but in Section VIII we argue that this is exponentially suppressed in m π L and indeed found the effect to be small on our 16 3 × 32 ensembles, for which m π L ∼ 4.
We also observed that, in contrast to anti-periodic boundary conditions, the imposition of Gparity boundary conditions in all three direction for a cubic box does not result in a system which is symmetric under cubic rotations. For example the quark free-field eigenstates of the system are not symmetric under cubic rotations but instead, for our conventions, must have the three components of their momentum equal modulo 2π/L. This implies, for example, that the momentum ( π 2L , π 2L , 0) (for GPBC in two directions) is allowed, but ( π 2L , − π 2L , 0) is not. This symmetry breaking does not extend to the pions, which obey standard antiperiodic boundary conditions in the G-parity directions as we demonstrated explicitly in Section VIII by showing that the energies of pions moving in orthogonal directions are in good agreement. However, the fact that we are restricted in our allowed choices of quark momentum components forces us to use operators that are inequivalent under cubic rotations in order to create pions moving in different orthogonal directions. This results in values of the two-point function amplitude that depend somewhat strongly on the direction of the pion momentum. We found that by choosing operators that are averages over multiple choices of quark momenta with fixed total momentum we drastically reduce the variation in the amplitude. This ultimately enabled us to to construct an interpolating operator that was rotationally symmetric within our measured accuracy to absorb the ππ state in our K → ππ calculation [5].
A further difficulty of the G-parity approach involves the treatment of the strange quark. In Section VI we describe how charge-conjugation boundary conditions, while consistent with the gauge field boundary conditions, are not suitable for constructing the stationary kaon state required for a K → ππ calculation. Such a state can be constructed if we instead impose G-parity boundary conditions between the strange quark and a fictional degenerate partner, s ′ , and comprises an admixture of the continuum kaon and an unphysical, degenerate partner state. We discuss how operators that interact only with the physical kaon retain their continuum values up to a known numerical factor and additional exponentially-suppressed finite-volume corrections that can be neglected. The effects of the unphysical strange sea quark introduced here are removed by using the square root of the resulting 2-quark fermion determinant in order to revert to a three-flavor simulation. This can be accomplished, for example by using the RHMC algorithm. We show that the non-local effects of applying this square root are exponentially suppressed in the spatial lattice size. In Section VIII we study the strange quark states numerically and observe that both B K and the kaon masses and decay constants remain in excellent agreement as we impose GPBC in a successively larger number of directions.
We conclude that G-parity boundary conditions provide a novel and useful means with which to study mesonic interactions and decays involving pions carrying above-threshold momenta. The RBC & UKQCD collaborations have already made use of these boundary conditions in order to compute the I = 0 K → ππ decay amplitude with physical kinematics [5], which requires moving pions. Here the statistical noise is very large due to disconnected contributions, which makes isolating an excited state contribution in such a setup difficult. Further investigations of the momentum dependence of the ππ scattering phase shift are among the future plans for this technique.  [15] for many of the highperformance optimized kernels. Measurements and ensemble generation were performed using the IBM Bluegene/Q computers [16] formerly at Brookhaven National Laboratory (BNL). We thank the RIKEN BNL Research Center and the Brookhaven National Laboratory for the use of these machines.

Appendix A: Domain wall fermions with G-parity boundary conditions
In Section III we formulated the lattice action for G-parity boundary conditions in terms of two four-dimensional fields ψ 1 (x) = d(x) and ψ 2 (x) = Cū T (x) that transform asT ψ 1 (L − 1)T −1 = ψ 2 (0) andT ψ 2 (L−1)T −1 = −ψ 1 (0) under translations across the boundary in a G-parity direction (suppressing orthogonal coordinates). In this Appendix we show how to extend this formulation to the five-dimensional quark fields in the domain wall fermion framework.
For domain wall fermions the fundamental fields Φ(x, s) are functions of the fifth dimensional coordinate s, and the four-dimensional fields are constructed as surface fields as follows: where φ and Φ here are generic 4d and 5d fields respectively, and we have used φ /Φ rather than the conventional ψ/Ψ to avoid confusion with our G-parity fields. Here and below we capitalize the field variables to indicate that they are five dimensional, and we explicitly display the coordinate in only one of the four space-time directions, which is assumed to have GPBC.
Starting with the known action of the charge conjugation operation on the 4d field, we can induce its action on those in five dimensions. This can be accomplished by expressing the chargeconjugated 4d field,Ĉφ (x)Ĉ −1 in two ways: Cφ (x)Ĉ −1 = P RĈ Φ(x, 0)Ĉ −1 + P LĈ Φ(x, L s − 1)Ĉ −1 (A2) and Cφ (x)Ĉ −1 = Cφ T (x) = P L CΦ T (x, 0) + P R CΦ T (x, L s − 1) (A3) and then equating the P L and P R terms found in these two equations. We conclude that charge conjugation acting on the 5d domain wall fields involves a reflection in the fifth dimension: where the second equation can be obtained by applying the same analysis to the conjugate fields φ andΦ.