Pseudoscalar mesons in a finite cubic volume with twisted boundary conditions

We study the effects of a finite cubic volume with twisted boundary conditions on pseudoscalar mesons. We first apply chiral perturbation theory in the p-regime and calculate the corrections for masses, decay constants, pseudoscalar coupling constants and form factors at next-to-leading order. We show that the Feynman-Hellmann theorem and the relevant Ward-Takahashi identity are satisfied. We then derive asymptotic formulae a la Luscher for twisted boundary conditions. We show that chiral Ward identities for masses and decay constants are satisfied by the asymptotic formulae in finite volume as a consequence of infinite-volume Ward identities. Applying asymptotic formulae in combination with chiral perturbation theory we estimate corrections beyond next-to-leading order for twisted boundary conditions.


Introduction
The study of finite volume effects, besides a purely theoretical interest, is also motivated by the need to correct results from lattice simulations. These are necessarily performed in a volume of finite extent on which some form of boundary conditions are imposed. If one chooses periodic boundary conditions, momenta are discretized and can not take continuous values. To overcome such limitation twisted boundary conditions were proposed [1][2][3][4]. In this paper we study the effects of a finite cubic volume with twisted boundary conditions on observables related to pseudoscalar mesons by applying chiral perturbation theory (ChPT). As is well known ChPT is the low-energy effective field theory of QCD and can be formulated in finite volume thereby providing a systematic tool to study finite-size effects on observables calculated in lattice QCD. The first analytical study of finite-volume effects with twisted boundary conditions was published soon after the proposal to use these was made [4]. This relied on one-loop ChPT. Further analytical calculations have been made also by other groups [5][6][7][8], and more recent ones have appeared in the last few years [9,10]. All these studies rely on ChPT at one loop. In the case of periodic boundary conditions it has been shown that the combined used of asymptotic formulae [11,12] and ChPT is the most efficient way to estimate higher orders in ChPT [13][14][15], with only tiny deviations from the results of a full two-loop calculation [16].
Asymptotic formulae for twisted boundary conditions are not available in the literature yet, so that all estimates of finite-volume effects have to rely on one-loop ChPT. The latter estimates are known to suffer from large two-loop corrections [13,14] (even though the absolute size of finite-volume effects remains small). This paper fills this hole: its main aim is the derivation and application of asymptotic formulae for finite-volume effects with twisted boundary conditions for meson masses and decay constants. We also apply and extend a suggestion by Häfeli [17] to use the Feynman-Hellmann theorem [18,19] to derive an asymptotic formula for the scalar form factor of the pion at zero momentum transfer.
While for mesons most quantities have already been calculated at one loop in ChPT, we provide here for completeness also expressions for finite-volume effects for masses, decay constants and form factors. This gives us also the chance to discuss issues like the very definition of masses and decay constants in the presence of twisted boundary conditions -since this is subject to a certain degree of arbitrariness -and the role of chiral Ward identities in finite volume. We also clarify the meaning of the Ward-Takahashi identity for the electromagnetic current in finite volume and discuss its violation due to the breaking of Lorentz invariance.
With the help of these asymptotic formulae we make a numerical analysis of finitevolume effects which goes beyond one-loop ChPT. As in the case of periodic boundary conditions, also here two-loop effects can be very sizable, whereas the effect of the twist tends to be, for small twisting angles, which is the most relevant case, small. We trust that the results presented here will allow a more reliable correction of finite-volume effects in all lattice calculations using twisted boundary conditions. This work is structured as follows. In section 2 we present ChPT in finite volume: we give an overview of periodic boundary conditions and introduce the twisted ones. Section 3 focuses on finite-volume corrections at next-to-leading order. The corrections for masses, decay constants, pseudoscalar coupling constants are recalculated and new results for pion form factors are presented. We show that the corrections of the pion scalar form factor satisfy the Feynman-Hellman theorem [18,19]. In section 4 we derive asymptotic formulae à la Lüscher for twisted boundary conditions. We sketch the steps necessary to generalize the original derivation of Lüscher [11] and present asymptotic formulae for masses, decay constants and pseudoscalar coupling constants. These are related via chiral Ward identities as we will show below. We also derive asymptotic formulae for the pion scalar form factor at zero momentum transfer relying on the Feynman-Hellmann theorem. In section 5 we apply the asymptotic formulae in combination with ChPT. We use the chiral representation at one loop to express the amplitudes entering the formulae and present results beyond nextto-leading order. Section 6 contains the numerical analysis. Appendices give further details on analytical aspects. Appendix A provides a list of useful results for the evaluation of loop diagrams in finite volume. Appendix B is devoted to the (electromagnetic) gauge symmetry in finite volume. Therein, we construct an effective theory for charged pions which is invariant under gauge transformations and which reproduces results of section 2. As the gauge symmetry is preserved, we show that the Ward-Takahashi identity [20][21][22] holds in finite volume if the momentum transfer is discrete. Finally appendix C collects some long expressions related to results presented in section 5.

Chiral perturbation theory
QCD is the fundamental theory of the strong interaction [23,24]. It describes the dynamics of the strong interaction in terms of gluons and quarks. The Lagrangian can be written as where G µν a is the strength field tensor of gluon fields, q = q(x) represents the quark fields arranged as a vector of flavor space and M is the matrix of the quark masses. If quark masses are zero, L QCD exhibits a global chiral symmetry. It is well known that chiral symmetry spontaneously breaks down and gives rise to Goldstone bosons which can be identified with the lightest pseudoscalar mesons. For N f = 3 quark flavors, the fields of Goldstone bosons can be parametrized by a 3 × 3 unitary matrix, (2. 2) The parameter F 0 is the decay constant in the limit of zero masses. The effective chiral Lagrangian can be ordered as a series in powers of momenta and quark masses L eff = L 2 + L 4 + . . . , (2.3) where each quark mass counts as a momentum square, i.e. O(m f ) ∼ O(p 2 ). At leading order the effective Lagrangian consists of terms O(p 2 ) and can be written as The angular brackets . denote the trace in flavor space and with B 0 a parameter of the effective Lagrangian. We have introduced external fields v µ , a µ , s, p as sources for the chiral Noether currents and quark scalar and pseudoscalar bilinears. They also allow one to include electromagnetism, semileptonic weak interactions, as well as an explicit breaking of chiral symmetry via quark masses. In particular, we work in the isospin limit (wherem := m u = m d ) and include masses by setting s = M = diag(m,m, m s ). At next-to-leading order (NLO) the effective Lagrangian consists of terms O(p 4 ) and can be compactly written as (2.6) The coupling constants L j contain the so-called low-energy constants (LECs) and the monomials P j are contructed from U , v µ , a µ , s, p. Their explicit expressions is well known and can be found, e.g. in [25].

Periodic boundary conditions
Numerical simulations of lattice QCD are by necessity performed in a volume of finite extent. The volume is usually a spatial cubic box of the side length L on which boundary conditions are imposed. Mostly employed are periodic boundary conditions (PBC) which require the periodicity of fields within the cubic box, q(x + Lê j ) = q(x), j = 1, 2, 3. (2.7) Here,ê µ j are unit Lorentz vectors pointing in the j-th spatial direction. In momentum space, the periodicity of the fields corresponds to a discretization of the momenta. The spatial components of momenta are discrete and read p = 2π m/L, m ∈ Z 3 . As a consequence, Lorentz invariance is broken. Still, the subgroup of spatial rotations of 90 • (the so-called cubic invariance) remains intact.
The momentum discretization also introduces a new scale in the theory: 1/L. To apply ChPT one must consider momenta smaller than Λ χ = 4πF π -where F π is the pion decay constant -and this provides the following quantitative condition [26]: In refs. [27][28][29] Gasser and Leutwyler showed how to apply ChPT in finite volume. They proved that for large enough volume the effective chiral Lagrangian and the values of LECs remain the same as in infinite volume. The most relevant change concerns the counting scheme to be applied to the effective Lagrangian and the propagators. The counting scheme has to take into account also 1/L together with momenta and quark masses, moreover propagators are modified by the discrete momenta. For what concern the counting, there are two possible ones corresponding to different regimes: Since here we will work in the p-regime we do not discuss the -regime any further and refer the reader to [28,[30][31][32] for more details about it. In the p-regime, M π is larger than 1/L: a pion fits well inside the box and behaves almost as if it were in infinite volume. Here, the counting scheme can be applied with the additional rule 1/L ∼ O(p), see ref. [29]. The expressions for propagators are similar to the infinite-volume ones, but integrals over spatial components must be replaced by sums over discrete values. This makes propagators periodic and dependent on L. Physical observables can be then calculated as in infinite volume: tree graphs produce exactly the same contributions and just the loop diagrams generate a finite-volume dependence.

Twisted boundary conditions
A serious limitation of PBC is the momentum discretization which makes it difficult to access very small, finite momenta without using huge volumes. Twisted boundary conditions (TBC), see refs. [1][2][3][4], have been introduced to overcome this difficulty. They require that fields are periodic up to a global symmetry transformation, Here, the subscript T indicates that fields satisfy TBC. The transformation U j has to be a symmetry of the action and as such depends on the form of the Lagrangian [4]. For QCD with 3 light flavors one can consider where v µ ϑ = diag(ϑ µ u , ϑ µ d , ϑ µ s ) is a traceless matrix commuting with M. The Lorentz vectors, are called twisting angles and their spatial components can be arbitrarily chosen. It is convenient to redefine quark fields as periodic ones by means of The periodicity of q(x) follows from the condition (2.10). After this field redefinition the twisting angles appear in the Lagrangian as a constant vector field, (2.14) The momentum of the flavor f is shifted by the corresponding twisting angle ϑ µ f , which is a free parameter and can therefore be varied continuously. Note that one can either impose the condition (2.10) and work with the original form of the Lagrangian or redefine quark fields as periodic and introduce the twisting angles through the constant vector field v µ ϑ in L QCD . The two approaches are equivalent.
Since they point in specific directions, twisting angles break various symmetries. In particular cubic invariance in momentum space is broken. More generally, all symmetries whose generators do not commute with v µ ϑ are broken. For three different twisting angles these are: the vector symmetry SU(3) V and the isospin symmetry. Note that the third isospin component I 3 , the strangeness S and the electric charge Q e are still conserved quantities in this case. In addition, the symmetry U j induces a new one: at each vertex the sum of incoming and outgoing twisting angles is conserved and equal to zero, see ref. [4].
In the effective theory the condition (2.10) implies that the unitary matrix parametrizing the fields of pseudoscalar mesons satisfies Here, the repetition of j does not imply any sum and T specifies that fields satisfy TBC. Through a field redefinition the unitary matrix can be made periodic, The twisting angles enter the effective Lagrangian as a constant vector field v µ ϑ where each derivative is replaced by ∂ µ . → ∂ µ . − i v µ ϑ , . . At leading order the Lagrangian reads , and U (x) now satisfying periodic boundary conditions. The commutator [v µ ϑ , U ] acts in different ways on the fields of pseudoscalar mesons. Pseudoscalar mesons sitting in the diagonal of U commute with v µ ϑ and their momenta are unshifted. Pseudoscalar mesons off the diagonal do not commute with v µ ϑ and their momenta are shifted by the twisting angles, Note that twisting angles reflect the flavor content of the particles. A pseudoscalar meson with the flavor content q fqf has the twisting angle ϑ µ q fqf = (ϑ q f − ϑ q f ) µ whereas its antiparticle has a twisting angle of opposite sign.
As pointed out in ref. [4] twisting angles enter the expressions of external states and modify internal propagators. As an example we consider charged pions (kaons have similar expressions). The field redefinition (2.16) implies that the propagators read The propagators are still periodic but obey modified Klein-Gordon equations: 1 We observe that twisting angles shift the poles in the denominator of propagators. Moreover, the substitution k 0 → −k 0 and k → − k reverses the propagation direction and the sign of twisting angles. Since antiparticles have twisting angles of opposite sign, we conclude that the propagation of a positive pion with ϑ µ π + in the forward direction of space-time is equivalent to a propagation of a negative pion with ϑ µ π − in the backward direction. To close this section, we remark that TBC are a generalization of PBC. Setting all twisting angles to zero, the condition (2.10) reduces to eq. (2.7), which means that by taking this limit in a TBC calculation and comparing the results with those for PBC, one gets a non-trivial -albeit partial -check.

Masses, decay constants and pseudoscalar coupling constants
In general, we define the corrections of an observable X as the exponential function of the numerator. In that case, propagators are periodic up to a phase [e.g. ∆ π ± ,L (x + Lêj) = e −iLê j ϑ π ± ∆ π ± ,L (x) for j = 1, 2, 3] and obey the usual Klein-Gordon equations.
where ∆X := X(L) − X is the difference among the observable evaluated in finite and in infinite volume. The corrections of masses and decay constants of pseudoscalar mesons were first calculated in ChPT with TBC in ref. [4]. For what concerns the masses, before giving explicit expressions for the finite-volume corrections, we need to define what we mean by "mass" in the presence of TBC. Having introduced twisting angles and with them a breaking of Lorentz symmetry, when we calculate the self-energy of a particle we will get Lorentz non-invariant contributions proportional to the twisting angles. In general, the self-energy of a meson P will have the form: The equation which determines the mass is: Since the solutions of the equation do not lie on a constant p 2 surface the mass would seem to depend on the direction of the momentum p µ . A possible solution for this non-invariance of the pole position can be obtained by completing the square and interpreting ∆ϑ µ Σ P as a renormalization of the twisting angle [7]: where ∆θ µ Σ P := ∆ϑ µ Σ P /(1 + B P ). The pole position is then given by: which is the mass definition we adopt, in agreement with refs. [4,7]. In contrast, the authors of ref. [10] adopt a mass definition which is momentum-dependent and treat the terms proportional to ∆ϑ µ Σ P as part of the finite-volume correction to the mass. We stress that our choice of the mass definition is consistent with the idea to treat the twisting angle as part of the momentum -which is the basic point of TBC. In general the mass is defined as the energy of a particle at zero spatial momentum, which is what is measured on the lattice. It appears therefore natural to take as definition of the mass the energy of a particle at zero total momentum (kinetic + twisting angle). The interpretation of ∆ϑ µ Σ P is of course somewhat arbitrary. As we will see below, the fact that at NLO exactly the same contribution also appears as an additive correction to F P (p + ϑ P ) µ in the matrix element of the axial current strongly suggests its interpretation as a renormalization of the twist. The definition of mass given in eq. (3.5) then naturally follows.
Indeed a similar situation occurs for the decay constants: the matrix element of the axial current is not just proportional to momentum, but there is an extra shift, defined as follows: Here too we consider ∆θ µ A P as a twisting-angle renormalization and not as part of the finite-volume correction to the decay constant (again in agreement with refs. [4,7] and in contrast to ref. [10]).
The corrections decay exponentially in λ P and depend on twisting angles through phase factors. Note that twisting angles can change the overall sign. This is a consequence of the breaking of the vector symmetry SU(3) V . For instance, mass corrections can turn negative whereas corrections of decay constants can turn posistive depending on ϑ µ π + , ϑ µ K + , ϑ µ K 0 . With an appropriate choice (or averaging over randomly chosen twisting angles) one can even suppress the corrections as discussed e.g. for nucleons in ref. [9].
In addition, the evaluation of the corrections involves the terms which remormalize the twisting angles. Such terms are not present for PBC and are generated by the breaking of the cubic invariance, see ref. [7]. In our notation, they read The function f µ 1 (λ P , ϑ) is defined in eq. (A. 4). Note that extra terms emerge only in the evaluation of corrections of charged pions and kaons. They are non-vanishing in the directions where twisting angles are non-vanishing and disappear for ϑ µ π + = ϑ µ K + = ϑ µ K 0 = 0. This intimately relates ∆ϑ µ Σ π ± , ∆ϑ µ Σ K ± , ∆ϑ µ Σ K 0 to twisting angles. Expressions for ∆ϑ µ A π ± , ∆ϑ µ A K ± , ∆ϑ µ A K 0 would also be needed to complete the formulae at NLO, but at this order they exactly coincide with the extra terms of eq. (3.9).
The corrections of the pseudoscalar coupling constants were first calculated in ref. [10]. At NLO the results read , (3.10) Note that δG π 0 (resp. δG η ) correspond to ∆ V G π 0 3 /G π (resp. ∆ V G η8 /G η ) of ref. [10]. At this order, no extra terms like eq. (3.9) appear. We show that these results can be obtained with the mass definition [4,7] relying on chiral Ward identities. We just illustrate the case of charged pions. The relevant chiral Ward identities read Here, the subscript L indicates that the matrix elements are evaluated in finite volume. The operators are linear combinations of the pseudoscalar densities resp. axialvector currents: The matrix elements on the left-(resp. right-) hand side of the chiral Ward identities are proportional to the pseudoscalar coupling (resp. decay) constants. Working out both sides and retaining terms up The mass definition [4,7] implies that charged pions lie on the following mass shells, Hence, the momentum squares on the right-hand side of eq. (3.12) value  and -when multiplied with F π ± (L) -produce contributions that exactly cancel 2F π (p + ϑ π ± ) µ ∆ϑ µ A π ± as at NLO, ∆ϑ µ A π ± coincide with ∆ϑ µ Σ π ± . Dividing bymG π = M 2 π F π we get δG π ± = δM 2 π ± + δF π ± + O(ξ 2 π ).

Pion form factors
In infinite volume the pion form factors are defined by the matrix elements, (3.16) They depend on the square of the momentum transfer, q 2 = (p − p) 2 and their expressions are known in ChPT at NLO [33] and at NNLO [34]. At vanishing momentum transfer, they satisfy the relations which follow from the Feynman-Hellmann theorem [18,19] and the Ward identity [35]. In finite volume the matrix elements of pion form factors receive additional corrections that can be defined as (3.18) Here, L (resp. q 2 = 0) indicates that matrix elements are evaluated in finite volume (resp. in infinite volume at vanishing momentum transfer). These corrections still depend on the momentum transfer. The twisting angles shift the momenta of charged pions but not necessarily induces a continuous momentum transfer. If the incoming and outgoing pions are the same, the twisting angles cancel out from the spatial components of the momentum transfer and We use this fact to work out the matrix elements in finite volume and evaluate the corrections. However, keep in mind that -depending on the kinematics chosen -the zeroth component q 0 may contain twisting angles of external pions and hence, vary continuously.
To study the pion form factors in finite volume we consider only N f = 2 light flavors. The corrections of masses, decay constants and pseudoscalar coupling constants can be obtained from eqs. (3.7, 3.8, 3.10) discarding the contributions of the virtual eta meson and kaons. Note that in this case, the extra terms read We first discuss the corrections of the matrix elements of the scalar form factor. At NLO the corrections can be evaluated from the loop diagrams of figure 1. The tadpole diagram generates corrections similar to those encountered before. The fish diagram generates additional corrections which can be calculated with the Feynman parametrization (A.5). Altogether, we find with λ z = λ π 1 + z(z − 1)q 2 /M 2 π . The functions g 2 (λ z , q, ϑ), g 2 (λ z , q) := g 2 (λ z , q, 0) originate from the fish diagram and can be evaluated by means of the Poisson resummation formula (A.2). Note that g 2 (λ z , q, ϑ) is even in the second and third argument. This is a consequence of the fact that the spatial components of the momentum transfer are discrete. The last term in eq. (3.21) consists of the product among (3.22) The Lorentz vectors ∆Θ µ π ± have non-vanishing components in the directions where both ϑ µ π ± and q µ are non-vanishing. They disappear for ϑ µ π ± = 0. The function f µ 2 (λ z , q, ϑ) is defined in eq. (A.7) and originates from the fish diagram. Note that as q is discrete, the function f µ 2 (λ z , q, ϑ) is even in the second argument and odd in the third one. The corrections of the matrix elements of the scalar form factor decay exponentially in λ π = M π L and disappear for L → ∞. As a check we set ϑ µ π ± = 0 and find the result for PBC [17,36]. In that case, the corrections are negative. For small twisting angles, the corrections stay also negative as the dependence on twisting angles is roughly a phase factor. They may turn positive for large twisting angles. Note that P µ ∆Θ µ π ± depend linearly on ϑ µ π ± . This dependence increases δΓ π ± S at large twisting angles. Thus, in order to keep the corrections under control, it is important to employ small twisting angles, e.g. | ϑ π ± | < π/L. At vanishing momentum transfer, the corrections reduce to (3.23) The functions g 2 (λ π , ϑ π + ) := g 2 (λ π , 0, ϑ π + ), g 2 (λ π ) := g 2 (λ π , 0, 0) and f µ 2 (λ π , ϑ π + ) := f µ 2 (λ π , 0, ϑ π + ) are defined in eq. (A.4). At vanishing momentum transfer the Feynman-Hellman theorem [18,19] relates the scalar form factor with the derivative of the pion mass, see eq. (3.17). This relation can be extended to finite volume. However, one must make some specification. As pointed out in ref. [33] the Feynman-Hellmann theorem states that the expectation value π b | S 0 |π a q 2 =0 is related to the derivative of the energy level describing the pion eigenstate. In finite volume the energy levels are additionally shifted by twisting angles and by the corrections of self-energies: Here, ∆Σ π ± := Σ π ± (L) − Σ π ± just contain corrections in finite volume. Taking the derivative ∂m = (∂M 2 π /∂m) ∂ M 2 π , and retaining only contributions in finite volume, we obtain This relation extends the statement of the Feynman-Hellmann theorem in finite volume: at q 2 = 0 the corrections of the matrix elements of the scalar form factor are related with the derivative of the self-energies with the respect to the mass. One can show that deriving the expressions of the self-energies at NLO one obtains the expressions (3.23). The corrections of the matrix elements of the vector form factor can be similarly evaluated from the loop diagrams of figure 1. In this case, some attention must be paid as the evaluation involves tensors in finite volume. At NLO we find (3.26) Here, Q e = ±1 represents the electric charge of π ± in elementary units and the functions g 1 (λ z , q, ϑ), h µν 2 (λ z , q, ϑ) are defined in eq. (A.7). At this order, ∆ϑ µ Γ π ± exactly coincide with the extra terms of eq. (3.20).
The corrections of the matrix elements of the vector form factor decay exponentially in λ π = M π L and disappear for L → ∞. As a check we set ϑ µ π ± = 0 and find the result for PBC obtained by Häfeli [17]. This result differs from the expression published in ref. [37] by a term proportional to q µ . We stress that such term contributes for PBC and only disappears when the momentum transfer is zero, as well. In this sense, we disagree with ref. [7] where it is claimed that the contribution G rot F V disappears for vanishing twisting angles. Actually, such contribution disappears when both the momentum transfer and the twisting angles are zero. In ref. [10] the corrections were calculated in SU(3) ChPT with TBC. Their results coincide with eq. (3.26) if contributions of the virtual eta meson and kaons are discarded. In general, the corrections roughly depend on twisting angles through a phase factor. For small twisting angles the corrections stay negative (resp. positive) for positive (resp. negative) pions. Note that in (3.26) there are terms linear in ϑ µ π ± . This increases the absolute value of (∆Γ π ± V ) µ at large twisting angles. At vanishing momentum transfer the corrections reduce to 27) and disappear only for L → ∞. Here, h µν 2 (λ π , ϑ π + ) := h µν 2 (λ π , 0, ϑ π + ), see eq. (A.4). Setting ϑ µ π ± = 0 the corrections (3.27) reduce to the result obtained for PBC [38]. In general, the fact that (∆Γ π ± V ) µ q 2 =0 differ from zero indicates that the (electromagnetic) gauge symmetry or Lorentz invariance are broken. In infinite volume, the vector form factor equals unity at vanishing momentum transfer, see eq. (3.17). This result follows from the Ward identity [35] which relies on both the gauge symmetry as well as Lorentz invariance and can be derived from the Ward-Takahashi identity [20][21][22]. The derivation implies a continuous limit on the momentum transfer which in finite volume can not be taken due to the discretization of spatial components. This invalidates the Ward identity in finite volume. It turns out that the corrections (3.27) are consequences of the breaking of Lorentz invariance. In ref. [38] these considerations were presented for PBC. The authors demonstrated that at vanishing momentum transfer the corrections respect the gauge symmetry and that the Ward-Takahashi identity holds for PBC. In appendix B we generalize their derivations for TBC. In particular, we construct an effective theory invariant under gauge transformations which reproduces eq. (3.27) and we show that the Ward-Takahashi identity holds for TBC as long as the spatial components of the transfer momentum are discrete.

Asymptotic formulae for TBC
Asymptotic formulae represent another method to estimate finite-volume corrections. They relate the corrections of a given physical quantity to an integral of a specific amplitude, evaluated in infinite volume. The method was introduced by Lüscher [39] and it has been widely applied in combination with ChPT as it allows one to get a chiral order almost for free at the price of neglecting exponentially suppressed contributions. Currently, there are asymptotic formulae for pseudoscalar mesons [12-14, 40, 41], nucleons [15,[42][43][44][45][46] and heavy mesons [15]. These formulae are valid in the p-regime and for PBC. Here, we generalize the method to TBC and derive asymptotic formulae for small twisting angles estimating the corrections for masses, decay constants, pseudoscalar coupling constants and scalar form factors of pseudoscalar mesons.

Generalization of Lüscher's derivation
The derivation of the asymptotic formulae for TBC can be led back to the original derivation of Lüscher [11]. In the following we outline the necessary steps to generalize the Lüscher's derivation to TBC and refer the reader to his paper for details.
In the first part of the proof, Lüscher showed by means of abstract graph theory that, for a generic loop diagram contributing to the self-energy, the dominant corrections are obtained if one takes all propagators in infinite volume but the ones of the lightest particles (in this case, pions) are taken in finite volume. 3 For TBC the propagators depend also on twisting angles. In that case, the twisting angles flowing along internal lines of a generic loop diagram add up to match the total external twisting angles entering the diagram. The situation can be illustrated by the graph in figure 2. This is a consequence of the conservation law pointed out in ref. [4]: at each vertex the sum of twisting angles is conserved and equal zero. In position space the contribution J(D, L) of the diagram D to the self-energy takes then the general form, J(D, n, L), Here, x µ (v) is the space-time coordinates of the vertex v ∈ V := V \ { b } and ∆x µ is the difference among the final and initial vertex of the line . The quantity V is a product of differential operators and generates the vertex functions of the diagram D. The pion propagators G π ∆x + Ln( ) are in infinite volume and can be expressed, e.g., with the heat kernel representation, see ref. [11]. For every line we have assigned an integer Lorentz vector n µ ( ) = 0 n( ) . Note that the summation over all possible sets of integers [n] for all internal propagators is a well-defined operation also in the case of TBC. The term [n] = [0] corresponds to the contribution in infinite volume with external momenta shifted by ϑ µ P . Discarding this term one finds that |J(D, L)| is exponentially bound so that the self-energy decays as O(e − √ 3 λπ/2 ) at asymptotically large L. The dominant corrections of I P Figure 3. Skeleton diagram contributing to asymptotic formulae for masses. Solid lines stand for a generic pseudoscalar meson P and dashed lines for virtual pions. The spline indicates that the pion propagator is in finite volume and accounts for integer vectors n ∈ Z 3 with | n| = 1. The blob corresponds to the vertex functions defined by the one-particle irreducible part of the amputated four-point function, see ref. [11].
the self-energy is then given by the contribution where for only one propagator | n( * )| = 1 and all others | n( = * )| = 0. This is represented by the skeleton diagram of figure 3.
The second part of the derivation consists in showing that by modifying the integration countour in the complex plane, the dominant corrections can be written as an integral of the forward P π-scattering amplitude evaluated in Minkowski space and analytically continued to complex values of its arguments. For TBC the pion propagators as well as the vertex functions depend on twisting angles. Through an integration shift one can express the dependence on twisting angles of virtual particles as a phase factor multiplying the vertex functions. The dependence on external twisting angles may be worked out expanding the vertex functions around small twisting angles. This is hardly a limitation since the main goal of the introduction of TBC is precisely to be able to access small momenta, which requires the use of small twisting angles. The first term of the expansion contributes to the dominant corrections of masses whereas a part of the second term provides the dominant contribution to the extra terms of the self-energy. The results are asymptotic formulae which, in the case of the neutral pion and the eta meson are valid for arbitrary twisting angles, and in the case of charged pions and kaons are valid for small external twisting angles only. Note that similar argumentations can be extended to the derivation of asymptotic formulae for decay constants and pseudoscalar coupling constants.

Analytical results
From the generalization of the Lüscher's derivation we obtain the following asymptotic formulae for the masses of pseudoscalar mesons, , Here, we display the resummed version of the formulae for whichλ =M L withM = ( [14,26]. Each symbol on the right-hand side refers to a quantity in infinite volume. The amplitudes are all defined in similar way. For instance, represent the isospin components of the Kπ-scattering in the t-channel with zero isospin [14]: In general, the asymptotic formulae depend on twisting angles through the phase factor exp(iL n ϑ π + ) in the amplitudes. The formulae for δM π ± , δM K ± , δM K 0 additionally depend on external twisting angles through We stress that the latter formulae are only valid for small external twisting angles. At tree level the chiral representation of F π ± (ν, ϑ π + ) = −M 2 π /F 2 π , F K ± (K 0 ) (ν, ϑ π + ) = 0 does not depend on twisting angles and if inserted in the asymptotic formulae provides the results (3.7) obtained with ChPT at NLO. The formulae for δM π 0 , δM η are valid for arbitrary twisting angles. Inserting the chiral representation at tree level of F π 0 (ν, ϑ π + ) = M 2 π 1 − 2 e iL n ϑ π + /F 2 π and F η (ν, ϑ π + ) = M 2 π 1 + 2e iL n ϑ π + / 3F 2 π , one recovers the results obtained with ChPT at NLO. While the asymptotic formulae are in principle valid up to terms O(e −λ ), at this chiral order they give the full result. Note that setting all twisting angles to zero, the asymptotic formulae reduce to the formulae valid for PBC [14].
Along with the formulae for mass corrections, we also derive asymptotic formulae for the extra terms breaking the Lorentz invariance of the self-energies. These formulae read . In that case, the amplitudes are given by differences of isospin components. For instance, and similarly for other pseudoscalar mesons. These asymptotic formulae are valid for small external twisting angles and provide the results (3.9) obtained with ChPT at NLO if one inserts the tree-level chiral representation of the amplitudes, namely G π ± (ν, Note that for ϑ µ π + = 0 the sums in eq. (4.6) are odd in n and hence, vanish -as they should for PBC. The asymptotic formulae for decay constants are similar to those for masses, , where D π ± , D K ± (K 0 ) are given in eq. (4.5). The amplitudes are defined from the matrix elements of the axialvector decay after a pole subtraction described in refs. [12,14]. For instance, the amplitudes of charged kaons read are the isospin components of the matrix elements describing the kaon decay into a two-pion state with zero isospin [14]: Inserting the chiral representation at tree level of the amplitudes, the asymptotic formulae provide the results (3.8) obtained with ChPT at NLO and if all twisting angles are set to zero, they reduce to the formulae valid for PBC [14].
We also derive asymptotic formulae for the extra terms arising in the matrix elements of the axial vector current. These read .
The amplitudes are given by differences of isospin components, e.g.
For ϑ µ π + = 0 the sums in eq. (4.11) are odd in n and hence vanish, as expected for PBC. The asymptotic formulae for the pseudoscalar coupling constants are , , The amplitudes are all defined in similar way for the various pseudoscalar mesons. For charged kaons, they are where the isospin componentsC K ± π 0 (s, t − u),C K ± π + (s, t − u), can be determined from the matrix elements, , after the pole subtraction, are the isospin components of the Kπ-scattering in the s-channel, see eq. (4.4). Note that the isospin component can be determined in a similar way from C K ± π + by exchanging p µ 1 ↔ p µ 2 . Furthermore, we derive asymptotic formulae for extra terms of the matrix elements These formulae are valid for small external twisting angles and read The amplitudes are given by the difference of the isospin components, e.g.

Chiral Ward identities
In ref. [12] it was pointed out that the asymptotic formulae for masses, decay constants and pseudoscalar coupling constants are related by means of chiral Ward identities. We are now going to show that the relation can be generalized to TBC. For convenience, we only illustrate the case of charged pions. We start from the relevant chiral Ward identities which in momentum space read where (4.21) The matrix elements (4.21) are in infinite volume though momenta are shifted by a twisting angle. Note that π c (p 1 )π c (p 2 )| is a two-pion state with zero isospin. We can leave out the twisting angles of that state as they will appear as phase factors after a shift of the loop momentum, see eq. (4.27). This is why the twisting angle just appears in the initial states. According to ref. [47] the matrix elements A π ± πc (ϑ π ± ) µ have a pole that does not enter the amplitudes in the asymptotic formulae. We subtract that pole expanding the matrix elements around Here, T π ± πc s, t(ϑ π ± ) − u(ϑ π ± ) correspond to the isospin components of the ππ-scattering in the s-channel with zero isospin and are Mandelstam variables shifted by the twisting angles of p 3 andQ. Note that the first variable does not depend on ϑ µ π ± as twisting angles exactly cancel out. The bar of Ā π ± πc (ϑ π ± ) µ indicates that the pole has been subtracted from the matrix elements.
Similarly, the pole in the matrix elements C π ± πc (ϑ π ± ) should also be removed, see [12]. We do this by expanding the matrix elements around We insert (4.22, 4.24) in the identities (4.20) and divide by M π . Using the relation (4.25) The last term can be brought on the left-hand side: we can set p µ 1 = −p µ 2 = k µ and rewrite the momenta where ν ± = u θ π ± −t θ π ± /(4M π ) = ν +k µθ µ π ± /M π . Now, we expand for small external twisting angles (i.e. aroundθ µ π ± = 0 or ν ± = ν) and multiply both sides by The integration over k can be performed to an accuracy of O(e −λ ) and what remains is the integral over y = −ik 0 /M π appearing in the asymptotic formulae. This provides us with where δM π ± , δF π ± , δG π ± resp. ∆ ϑ Σ π ± , ∆ ϑ A π ± , ∆ ϑ G π ± are given in terms of eqs. (4.2, 4.8, 4.13) resp. eqs. (4.6, 4.11, 4.18). These relations hold if the amplitudes entering the asymptotic formulae satisfŷ which they actually do in general. As a check one can insert the chiral representation of these amplitudes provided below and explicitly verify that both relations hold.
The relations for other pseudoscalar mesons can be proved in an analogous way and give conditions on the amplitudes similar to eq. (4.29), see ref. [48]. From these conditions it is possible to determine unknown amplitudes starting from the explicit representation of the related amplitudes. For instance, one can determine N π 0 from the chiral representation of F π 0 , C π 0 (see refs. [33,49] [50,51]. Later, we will use this fact to work out the asymptotic formulae for the decay constant of the neutral pion and for the pseudoscalar coupling constants of charged kaons.

Pion form factors
As proposed by Häfeli [17] one can rely on the Feynman-Hellmann theorem to derive asymptotic formulae for the matrix elements of the scalar form factor. In finite volume the Feynman-Hellmann theorem relates the corrections of the matrix elements of the scalar form factor with the derivative of the self-energies, see eq. (3.25). Starting from such relation we can derive asymptotic formulae valid at vanishing momentum transfer via .
Taking the derivative of the asymptotic formulae δM π 0 , δM π ± , ∆ ϑ Σ π ± one finds These asymptotic formulae depend on the amplitudes entering the expressions for δM π 0 , δM π ± , ∆ ϑ Σ π ± . Here, the dependence on the twist is threefold. The formulae depend on the twisting angle of the virtual positive pion through the phase factor exp(iL n ϑ π + ) in the amplitudes. Furthermore, they depend on the external twisting angles through the parameter D π ± and through the product n ϑ π ± in the last line of eq. (4.31). Note that the formulae for charged pions are only valid for small external twisting angles. The formula for the neutral pion is valid for arbitrary twisting angles.
We have checked the above asymptotic formulae in two ways. First, we have set all twisting angles to zero and we have recovered the formula valid for PBC, originally proposed by Häfeli [17]. Second, we have inserted the tree-level chiral representation of F π 0 , F π ± , G π ± and we have obtained the results (3.23) found at NLO with SU(2) ChPT.
The general derivation of asymptotic formulae for form factors is complicated by the presence of a non-zero momentum transfer. We anticipate that in this case we were not able to derive asymptotic formulae. The first part of the derivation is similar to the one outlined in section 4.1.1. In this case, there is an additional external line due to the insertion of the pseudoscalar densities (resp. vector currents). Taking appropriate modifications one can still show that the argumentation holds and the matrix elements of form factors decay as O(e − √ 3 λπ/2 ) at asymptotically large L. In the second part of the derivation, complications arise due to the injected non-zero momentum transfer, which implies that the integration over the loop momentum can not be performed as for the self-energy. Note that these complications should not arise in matrix elements with two different external states, as e.g. K + | S 4+i5 π 0 . In such matrix elements, both external particles can be taken in the rest frame and the momentum transfer remains non-zero. The vertex functions could be then expanded around small external angles as outlined in section 4.1.1 and the integration performed in a similar way as for the self-energy [52].

Amplitudes in ChPT
We apply the asymptotic formulae in combination with ChPT and estimate finite-volume corrections beyond NLO. We use the chiral representation of the amplitudes at one loop Asymptotic formula Amplitude Process ref. Table 1. We summarize the quantities for which asymptotic formulae have been derived in this work, the corresponding infinite-volume amplitudes needed and the processes from which the chiral representation can be determined. We also provide references where the results for the amplitudes at one loop can be found. Here, P a −→ π a (π c π c ) denotes the pseudoscalar decay into three pions in which two pions form a state with zero isospin.
which is known for most of the amplitudes derived in this work. In table 1 we summarize all quantities for which we have derived an asymptotic formula and the infinite-volume amplitudes needed therein. If available we give the reference where the chiral representation at one loop can be found. We also list the process from which the chiral representation can be determined. Note that currently, N π 0 , C K ± , K K ± are unknown in ChPT but they can be determined from F π 0 , C π 0 , F K ± , N K ± and G K ± , H K ± relying on chiral Ward identities.

Chiral representation at one loop
Pions The amplitudes F π 0 , F π ± , G π ± are defined similarly to eqs. (4.3, 4.7). Their chiral representation can be determined from the ππ-scattering, π(p 1 ) + π(p 2 ) −→ π(p 3 ) + π(p 4 ), for forward kinematics: All isospin components of the ππ-scattering can be given in terms of the invariant amplitude A stu := A(s, t, u) which is known up to NNLO in ChPT [49]. In terms of A stu the isospin components of F π 0 , F π ± , G π ± read Inserting the expression of A stu at NLO [33] and evaluating the isospin components in the kinematics (5.1) one obtains the chiral representation of F π 0 , F π ± , G π ± at one loop.
The amplitudes N π ± , H π ± have similar definitions as eqs. (4.9, 4.12). Their chiral representation can be determined from the matrix elements of the τ -decay, τ (p τ ) −→ π(p 1 )+ π(p 2 ) + π(p 3 ) + ν τ (p τ − Q), for forward kinematics: In ChPT the matrix elements of this τ -decay are known at NLO [47] and can be written as According to [47] these matrix elements can be decomposed as 4 (5.6) As described in ref. [12] the matrix elements (5.5) have a pole that does not contribute to the amplitudes of asymptotic formulae and must be subtracted. The subtraction occurs expanding the matrix elements around Q 2 = (p 1 + p 2 + p 3 ) 2 = M 2 π and isolating the pole appearing in F 123 , G 123 resp. F where T π − π 0 (s 3 , s 1 −s 2 ), T π − π − (s 3 , s 1 −s 2 ) are isospin components of the ππ-scattering amplitude in the s 3 -channel, see eq. (5.2). The bar indicates that the pole has been subtracted from the matrix elements. The amplitudes can be obtained contracting with p µ 3 /M π and setting the momenta to p µ The isospin componentĀ π − π + (s 3 , s 1 − s 2 ) can be obtained from A π − π − µ exchanging p µ 1 ↔ p µ 2 . Note that isospin symmetry relates the amplitudes of the negative pion to those of the positive pion via and H π + (ν, ϑ π + ) = −H π − (ν, ϑ π + ). (5.9) The amplitudes C π 0 , C π ± , K π ± are defined similarly to eqs. (4.14, 4.19). Their representation can be determined from the point function containing four pseudoscalar densities, where three of them serve as interpolating fields for pions. In ChPT such point function can be calculated from eq. (16.2) of ref. [33]: one must take the off-shell amplitude A(s, t, u; p 2 1 , p 2 2 , p 2 3 , p 2 4 ) and set the momenta to p 2 1 = p 2 2 = p 2 3 = M 2 π resp. p 2 4 =Q 2 with Q µ = (p 3 − p 1 − p 2 ) µ . Defining C stu := A(s, t, u; M 2 π , M 2 π , M 2 π ,Q 2 ), the isospin components of C π 0 , C π ± , K π ± can be expressed as after having subtracted the pole atQ 2 = M 2 π . In section 4.1.2 we describe how to subtract this pole and how to define the isospin components, see eq. (4. 16). From that definition, one obtains the chiral representation evaluating the isospin components in the forward kinematics: In section 4.2 we have derived asymptotic formulae for the matrix elements of the scalar form factor of pions at vanishing momentum transfer. The amplitudes entering those formulae are F π 0 , F π ± , G π ± . The chiral representation is given by the isospin components (5.2) evaluated in the forward kinematics. In this case, the asymptotic formulae also contain ∂ M 2 π F π 0 , ∂ M 2 π F π ± , ∂ M 2 π G π ± , for which attention must be paid during their calculation. Kaons The amplitudes F K ± (K 0 ) , G K ± (K 0 ) can be determined from the Kπ-scattering, π + (p 1 ) + K + (p 2 ) −→ π + (p 3 ) + K + (p 4 ), in the case of forward kinematics: In ChPT this scattering is given by the amplitude T 3/2 stu := T 3/2 (s, t, u) which is known at NLO from ref. [50]. 5 In terms of T The chiral representation of F K ± (K 0 ) , G K ± (K 0 ) can be then obtained evaluating these expressions in the kinematics (5.12). Note that from eq. (5.13) it turns out that the isospin 5 As noted in [14] there are two misprints in eq. (3.16) of [50]. The prefactor of [M r πK (u) + M r Kη (u)] should read (M 2 K − M 2 π ) 2 and the factor of 3 components of the negative kaon are equal to those of the neutral and hence, (5.14) The amplitudes N K ± , H K ± are defined in eqs. (4.9, 4.12). Their chiral representation can be determined from the matrix elements of the K 4 -decay, K + (p 3 ) −→ π(p 1 ) + π(p 2 ) + + (p ) +ν (Q − p ), for forward kinematics: In ChPT the matrix elements of the K 4 -decay are known at NLO [51] and read According to [51] these matrix elements can be decomposed as 6 As described in ref. [14] the matrix elements (5.17) have a pole that must be subtracted. Expanding aroundQ 2 = M 2 K the pole appears in R + stu resp. R stu and can be expressed in terms of the isospin components of the Kπ-scattering amplitude. One gets The amplitudes of asymptotic formulae can be obtained contracting with p µ 3 /M K and setting the momenta to p µ 2 = −p µ 1 resp. p µ 3 = M K 0 . Thus, the isospin components of N K + , H K + readĀ andĀ K + π − (s, t−u) can be obtained from Ā K + π + µ exchanging p µ 1 ↔ p µ 2 . Note that isospin symmetry relates the amplitudes of the positive kaon to those of the negative one as and H K − (ν, ϑ π + ) = −H K + (ν, ϑ π + ). (5.21) 6 The factors Fstu, Gstu, Rstu are 1/ √ 2 times smaller than those of ref. [51].

Eta meson
The amplitude F η is defined similarly to eq. (4.3). Its chiral representation can be determined from the ηπ-scattering, π(p 1 ) + η(p 2 ) −→ π(p 3 ) + η(p 4 ), for forward kinematics: In ChPT the ηπ-scattering is given by the invariant amplitude T πη (s, t, u). In terms of T πη (s, t, u) the isospin components of F η are all the same, Inserting the expression of T πη (s, t, u) at NLO [53] and evaluating the isospin components for forward kinematics one obtains the chiral representation of F η at one loop.

Chiral expansion
We now apply the asymptotic formulae of section 4. The results are presented in sections 5.2-5.4 with long expressions relegated to appendix C. To better organize these results we follow ref. [14] and make use of the chiral expansion. We first keep the discussion general and consider a pseudoscalar meson P with the mass M P and the twisting angle ϑ µ P . The amplitudes of table 1 can be expressed in the generic forms, where ν = (s − u)/(4M P ) andν = ν/M π . The functions Z P π 0 (t, u − s), Z P π + (t, u − s), Z P π − (t, u − s) are isospin components in infinite volume. Let us assume that X P (ν, ϑ π + ) enters the asymptotic formula for the observable X P and so, provides an estimate for the corrections δX P . The asymptotic formula has then the form where λ P = M P L and D P = M 2 P + | ϑ P | 2 − M P . Analogously, assume that the amplitude Y P (ν, ϑ π + ) enters the asymptotic formula for the extra term ∆ϑ µ X P which has the form For convenience, we rewrite the amplitudes as X P (ν, ϑ π + ) = X P (X P , π 0 ) + X P (X P , π ± ) e iL n ϑ π + , Y P (ν, ϑ π + ) = Y P (ϑ X P ) e iL n ϑ π + , where we collect the isospin components in X P (X P , π 0 ) := Z P π 0 (0, −4M P ν), X P (X P , π ± ) := Z P π + (0, −4M P ν) + Z P π − (0, −4M P ν), Y P (ϑ X P ) := Z P π + (0, −4M P ν) − Z P π − (0, −4M P ν).
Each contribution is rescaled by X π /X P where X π denotes the observable X P in infinite volume for P = π. The integrals I (j) (X P , π 0 ), . . . , I D (X P , π ± ) can be determined from the terms X (j) P (X P , π 0 ), X (j) P (X P , π ± ) of eq. (5.30). Note that the contributions R D (X P , π 0 ), R D (X P , π ± ) are proportional to the parameter D P and are not present if the particle P has no external twisting angle (as e.g. for the neutral pion and the eta meson).
The asymptotic formula (5.26) can be expanded in a similar way as The integrals I (j) (ϑ X P ) can be determined from the term Y (j) P (ϑ X P ) of eq. (5.30). Note that R(ϑ X P ) is not present if the particle P has no twisting angle (as e.g. for π 0 and η) and disappears for ϑ µ π + = 0.

Masses
We start with the asymptotic formulae for pion masses. In the generic form of eq. (5.32) the formula of the neutral pion exhibits two contributions: R(M π 0 ) = R(M π 0 , π 0 )+R(M π 0 , π ± ). At one loop, R(M π 0 , π 0 ) and R(M π 0 , π ± ) are given by the first two expressions of eq. (5.33) multiplied by (−1/2) and replacing X P = X π = M P = M π . The integrals I (j) (M π 0 , π 0 ) resp. I (j) (M π 0 , π ± ) can be determined from the chiral representation of F π 0 and read The functions B 2k = B 2k (λ π | n|) were defined in ref. [13] and can be evaluated analytically where K r (x) are modified Bessel functions of the second kind. The constants¯ j were originally introduced in ref. [33] and depend logarithmically on the pion mass, Here, M phys π = 0.140 GeV and¯ phys j are listed in table 2. The terms S (4) (M π 0 , π 0 ), S (4) (M π 0 , π ± ) contain integrals that can not be evaluated analytically but just numerically. Their explicit expressions are given in eq. (C.1).

Decay constants
The asymptotic formula for the decay constant of the neutral pion is given in eq. (4.8).

(5.46)
These relations follow from eq. (5.10) once the pole has been subtracted. The integrals I (j) D (G π ± , π 0 ) resp. I (j) D (G π ± , π ± ) can be determined from the derivative of the chiral representation ∂ y C π ± and are related to those of masses and of decay constants through The asymptotic formulae for the extra terms ∆ϑ µ G π ± can be expressed in the form of eq. (5.34) replacing ϑ X P = ϑ G π ± , X P = X π = G π and M P = M π . In this case, the equation must be divided by M 2 π . The integrals of R(ϑ G π ± ) can be evaluated from the chiral representation of K π ± and are related to those of eqs. (5.40, 5.44) through

Scalar form factors at vanishing momentum transfer
In section 4.2 we have presented asymptotic formulae for the matrix elements of the scalar form factor valid at vanishing momentum transfer. The formula for the neutral pion is given in eq. (4.31). If we insert the chiral representation of F π 0 and expand according to eq. (5.31a), the formula exhibits two contributions, R(Γ π 0 S ) = R(Γ π 0 S , π 0 ) + R(Γ π 0 S , π ± ). At one loop, R(Γ π 0 S , π 0 ) and R(Γ π 0 S , π ± ) are given by the first two expressions of eq. (5.33) multiplied by (−1/2) and replacing X P = X π = F S (0) as well as M P = M π . The integrals I (j) (Γ π 0 S , π 0 ) resp. I (j) (Γ π 0 S , π ± ) can be evaluated from the chiral representation of F π 0 . The evaluation involves terms with 1 + y 2 F π 0 and we make use of R dy y 2k λ π | n| 1 + y 2 e −λπ| n| where B 2k are defined in eq. (5.36). The derivative ∂ M 2 π F π 0 must be evaluated with care. The operator ∂ M 2 π acts on all quantities depending on the pion mass. In particular, it acts on the decay constant F π and on the constants¯ j of eq. (5.37). This leads to supplementary terms which must be integrated and added according to their chiral order to the corresponding integral. Altogether, we find where S (4) (Γ π 0 S , π 0 ), S (4) (Γ π 0 S , π ± ) are explicitly given in eq. (C.7). The asymptotic formulae for the matrix elements of the scalar form factor of charged pions are given in eq. (4.31). If we insert the chiral representation of F π ± , G π ± and expand according to eq. (5.31a) the formulae exhibit five contributions which can be written as . At one loop, the first four contributions are given by the expressions of eq. (5.33) multiplied by (−1/2) where X P = X π = F S (0), M P = M π and D P = D π ± . The integrals I (j) (Γ π ± S , π 0 ) resp. I (j) (Γ π ± S , π ± ) can be evaluated from the first group of terms in the square brackets of eq. (4.31). Using the chiral representation of F π ± we find where S (4) (Γ π ± S , π 0 ), S (4) (Γ π ± S , π ± ) are given in eq. (C.8). The integrals I (j) S , π ± ) can be determined from the second group of terms in the square brackets of eq. (4.31) after having taken the derivative of the chiral representation ∂ y F π ± . We find where C π ± = M π /(M π + D π ± ) with D π ± defined by eq. (4.5). The terms S (4) S , π ± ) are explicitly given in eq. (C.8). The fifth contribution R(Θ π ± ) can be expressed in the form of eq. (5.34) replacing ϑ X P = Θ π ± , X P = X π = F S (0) and M P = M π . In this case, the equation must be divided by 2M 2 π . The integrals of R(Θ π ± ) can be evaluated from the last group of terms in the square brackets of eq. (4.31). Using the chiral representation of G π ± we find where S (4) (Θ π + ) can be found in eq. (C.8).

Masses
The asymptotic formulae for the kaon masses are given in eq. (4.2). They differ in terms of the parameters D K ± , D K 0 and in terms of the amplitudes F K ± , F K 0 . In section 5.1.1 we have seen that the chiral representation of F K 0 is equal to that of F K − , see eq. (5.14). Hence, we can just consider the asymptotic formulae of charged kaons as the results of the neutral kaon can be obtained replacing D K ± with D K 0 .
We insert the chiral representation of F K ± and expand according to eq. (5.31a). The asymptotic formulae of charged kaons exhibit four contributions which can be written as . At one loop, these contributions are given by eq. (5.33) multiplied by (−1/2) where X P = X π = M P = M K and D P = D K ± . The integrals I (j) (M K ± , π 0 ) resp. I (j) (M K ± , π ± ) can be evaluated from the chiral representation of F K ± and are related to the integrals I (j) M K of ref. [14] through The integrals I (j) can be determined from the derivative of the chiral representation ∂ y F K ± and read Here, N = (4π) 2 , P = 2 log(M P /µ), x P Q = M 2 P /M 2 Q for P, Q = π, K, η and L r j are the renormalized LECs, see [25]. The expression of S (4) D (M K ± , π 0 ) is given in eq. (C.12). The asymptotic formulae for the extra terms ∆ϑ µ Σ K ± can be expressed in the form of eq. (5.34) replacing ϑ X P = ϑ Σ K ± and X P = X π = M P = M K . The integrals of R(ϑ Σ K ± ) can be evaluated from the chiral representation of G K ± and read where S (4) (ϑ Σ K + ) is given in eq. (C.13). Note that the asymptotic formula for ∆ϑ µ Σ K 0 differs from that for ∆ϑ µ Σ K − only in terms of G K 0 , G K − . As the chiral representation of G K 0 coincides with that of G K − , the asymptotic formulae are equal in this case, and we can use R(ϑ Σ K − ) to estimate ∆ϑ µ Σ K 0 .

Decay constants
The asymptotic formulae for the decay constants of charged kaons exhibit four contributions: . Their expressions at one loop are given by eq. (5.33) where X P = F K , X π = F π , M P = M K and D P = D K ± . The integrals I (j) (F K ± , π 0 ) resp. I (j) (F K ± , π ± ) can be evaluated from the chiral representation of N K ± and are related to the integrals I (j) F K of ref. [14] by virtue of 7 (5.57) 7 In ref. [14] there are two missprints: in eq. (57) the term 2 π (xπη − 9 4 ) should read 2( η − π The integrals I (j) can be determined from the derivative of the chiral representation ∂ y N K ± and read is given in eq. (C.14). The asymptotic formulae for the extra terms ∆ϑ µ A K ± can be expressed in the form of eq. (5.34) replacing ϑ X P = ϑ A K ± , X P = F K , X π = F π and M P = M K . The integrals of R(ϑ A K ± ) can be evaluated from the chiral representation of H K ± and read where S (4) (ϑ A K + ) is given in eq. (C.15).

Pseudoscalar coupling constants
The asymptotic formulae for the pseudoscalar coupling constants of charged kaons are given in eq. (4.13). In section 4.1.3 we have mentioned that C K ± is related to F K ± , N K ± by virtue of chiral Ward identities. We use that relation to determine the chiral representation of C K ± at one loop. Expanding the representation as in eq. (5.31a) the asymptotic formula exhibits four contributions: R(G K ± ) = R(G K ± , π 0 ) + R(G K ± , π ± ) + R D (G K ± , π 0 ) + R D (G K ± , π ± ). Their expressions at one loop are given by eq. (5.33) where X P = G K , X π = G π , M P = M K and D P = D K ± . The integrals can be evaluated from those of eqs. (5.54, 5.57) by means of where j = 2, 4 and Here, similar relations hold also for I (j) . In the evaluation, some attention must be paid: because of the prefactors (i.e.

•
x πK x −1 πK and F K /F π ) the integrals with j = 2 on the right-hand side generate terms that contribute to the integrals with j = 4 on the left-hand side.
The asymptotic formulae for the extra terms ∆ϑ µ G K ± are presented in eq. (4.18). As mentioned in section 4.1.3 the chiral representation of K K ± is related to that of G K ± , j¯ phys j ref.   Table 2. Values of LECs used in the numerical analysis. Here, the renormalization scale is µ = 0.770 GeV. The central value of L r 7 has been determined evaluating the expression of M η at the physical point (see text) whereas its uncertainty is taken from table 5 (column "All") of ref. [56].
H K ± by means of chiral Ward identities. We use that relation to determine the chiral representation of K K ± at one loop. The asymptotic formulae for ∆ϑ µ G K ± can be then expressed in the form of eq. (5.34) replacing ϑ X P = ϑ G K ± , X P = G K , X π = G π and M P = M K . In this case, the equation must be divided by M 2 K . The integrals of R(ϑ G K ± ) can be evaluated from those of eqs. (5.56, 5.59) by means of Note that because of the prefactors (viz.

Numerical set-up
We adopt the numerical set-up of refs. [13,14] and express the quantities in infinite volume appearing in the formulae (i.e. F π , F K , M K , M η ) as functions of M π . For F π we use the expression at NNLO obtained with SU(2) ChPT while for F K , M K , M η we use expressions at NLO of SU(3) ChPT. In any cases, the values of the relevant LECs are summarized in table 2. If available, these values are taken from results of lattice simulations with N f = 2 + 1 dynamical flavors. For LECs of SU(2) ChPT we take the averages of the FLAG working group [55] which are obtained from refs. [58][59][60][61]. For LECs of SU(3) ChPT we take the results of ref. [57] as recommended by FLAG [55]. The remaining values come from phenomenology [54,56] or have been determined evaluating the expressions of F π , F K , M K , M η at the physical point (see later).

Pion mass dependence in infinite volume
In infinite volume, the expression of the pion decay constant at NNLO reads This expression is obtained with SU (2) ChPT, see ref. [34]. Here, ξ = M 2 π /(4πF ) 2 and π = 2 log(M π /µ). The constants¯ j are given in eq. (5.37) and the values of¯ phys j are in table 2. The parameter r F (µ) is a combination of LECs at NNLO. In ref. [13] it was estimated as r F (µ) = 0 ± 3.
We can determine F numerically, evaluating eq. (6.1) at the physical point. Taking M π = M phys π ± and inverting the expression, we find F = (86.6 ± 0.4) MeV. This value agrees with the result of ref. [13] and provides the ratio F π /F = (1.065 ± 0.006) which in turn agrees with the FLAG average obtained from simulations with N f = 2 + 1 dynamical flavors, see refs. [55,57,58,60].
In figure 4 we represent the pion mass dependence of decay constants. We observe that the dependence of F π is rather mild. As in ref. [13], we conclude that this dependence is too mild to violate the condition of eq. (2.8). Thus, we can rely on ChPT and apply the formulae of sections 2 and 5 to estimate finite-volume corrections.  At NLO the expressions of M 2 K , M 2 η , F K obtained with SU(3) ChPT read [14]: • η (6.2b) In these expressions, k 1 = 2L r 8 − L r 5 , k 2 = 2L r 6 − L r 4 , k 3 = 3L r 7 + L r 8 and • P = 2 log( • M P /µ) resp. P = 2 log(M P /µ) for P = π, K, η. The decay constant F π is expressed as in eq. (6.1). Here, the circled masses are not just at leading order but are of hybrid nature, The part containing M 2 π is at NLO while the part containing B 0 m s is at leading order. This is unavoidable if we want to study the dependence on the pion mass of M K , M η , F K in ChPT. In practice, we use the first two expressions (6.2a, 6.2b) to determine B 0 m s , L r 7 and the third one (6.2c) to check the numerical results.
Taking the LECs from table 2 we evaluate eq. (6.2a) at M π = M phys π ± . Requiring M K = M phys K ± we find B 0 m s = 0.241 GeV 2 . The uncertainty on this value is of order 10 −5 GeV 2 and may be neglected. Then, we evaluate eq. (6.2b) at M π = M phys π ± and require M η = M phys η . We find L r 7 = −0.16 · 10 −3 . This value agrees with the result presented in table 5 (column "All") of ref. [56]. As uncertainty estimated we take the one given in [56], yielding L r 7 = (−0.16 ± 0.15) · 10 −3 . To check our numerical results on B 0 m s , L r 7 we insert the expressions of M 2 K , M 2 η in eq. (6.2c) and evaluate F K at M π = M phys π ± . We find F K = (0.106 ± 0.004) GeV which agrees with the PDG result [62] and with the FLAG average obtained from simulations with N f = 2 + 1 dynamical flavors, see refs. [55,59,63,64].
We stress that F π is expressed here with SU(2) ChPT even in the SU(3) expressions of M 2 K , M 2 η , F K . This choice was already made in ref. [14] and for m s = m phys s it exactly reproduces what one would get in the SU(3) framework. As lattice simulations are usually performed at m s ≈ m phys s we expect that such choice remains a valid approximation also in our numerical analysis.
In figure 5 we show the pion mass dependence of M K , M η . We observe that the dependence on M π is mild. The same holds for F K as one sees from figure 4. Note that for M π ≈ 0.500 GeV we have M K ≈ 0.610 GeV and M η ≈ 0.640 GeV. In that case, the values of M π , M K , M η are all similar. If we consider ξ P = M 2 P /(4πF π ) 2 for P = π, K, η we expect that the expansion parameter stays small for all pion masses in [0.1 GeV, 0.5 GeV]. This is confirmed by figure 5 where the pion mass dependence of ξ P is represented graphically. In the numerical analysis we will use the expansion parameter as in figure 5 and consider ξ P as exact, ignoring its uncertainty.

Twisting angles and multiplicity
To perform the numerical analysis we use the following configuration of twisting angles In principle, the angle θ can be arbitrarily chosen. Since we rely on ChPT, we require θ/L 4πF π . If we consider L 1 fm -as requested by eq. (2.8) -such condition is certainly satisfied for θ ∈ [0, 2π].
The configuration (6.4) allows us to simplify a lot of calculations. As twisting angles are aligned to one specific axis, we may use the formula, to rewrite the three sums over integer vectors as a nested sum over n ∈ N. In general, this speeds up the numerical evaluation by a factor 15. The notation . indicates the floor function and m(n, n 1 ) is the multiplicity (i.e. the number of possibilities to construct a vector n ∈ Z 3 with n = | n| 2 , having previously fixed the value of the first component to n 1 ∈ Z). As an illustration, we list in table 3 the values of the multiplicity for n ≤ 10.
In the following we present our numerical results. We plot the dependences of the corrections on M π and θ. The pion mass dependence is plotted for different values of L and θ. Lines of different colors refer to different values of L whereas lines of different hatchings refer to different values of θ. The dark (resp. light) yellow areas refer to the region M π L < 2 m(n, n 1 ) n 1 0 ±1 ±2 ±3  Table 3. Multiplicity of vectors n ∈ Z 3 with n := | n| 2 ≤ 10.
for θ = 0 (resp. θ = π/3). The dependence on the angle θ is plotted for different values of M π , at fixed L = 2 fm. Lines of different hatchings refer to different values of M π . For L = 2 fm, the region M π L < 2 begins for masses smaller than M π = 0.197 GeV. We remind the reader that in that region formulae obtained in the p-regime are not reliable any more.

Finite-volume corrections at NLO
The corrections of masses, decay constants and vector form factors were numerically evaluated at NLO in refs. [7,10]. As we have reproduced their plots we refrain from presenting results for these quantities at this order. We just focus on pseudoscalar coupling constants and scalar form factors for which no numerical analysis was published, yet. We start with the pseudoscalar coupling constants. At NLO the corrections of the pseudoscalar coupling constants exhibit all a dependence on twisting angles. In section 3.1 we have seen that the corrections of the pseudoscalar coupling constants are given by the sum of the corrections of masses and decay constants, see eq. (3.15). As δM π ± , δM K ± and δF π ± , δF K ± were numerically evaluated in ref. [10] one can use those results to determine δG π ± , δG K ± . Moreover, δG K 0 can be determined from δG K ± substituting ϑ K 0 ↔ ϑ K + , cfr. eq. (3.10). Thus, we just study the dependence of δG π 0 , δG η on the pion mass.
In figure 6a (resp. 6b) we represent the pion mass dependence of −δG π 0 (resp. −δG η ). The logarithmic graphs illustrate the exponential decay O(e −λπ ) of the corrections. The line slopes depend on L while the y-intercepts on θ. In figure 6a the lines are so close that they overlap in the graph: δG π 0 is practically insensitive to θ. On the contrary, δG η is noticeably sensitive to θ. In general, the corrections are negative and for θ ∈ { 0, π/8, π/4, π/3 } their absolute values decrease with the angle. Note that δG π 0 , δG η reach the percentage level before entering yellow areas and are thus comparable with the statistical precision of lattice simulations.
In figure 7a and 7b we represent the corrections of the matrix elements of the pion scalar   form factor at vanishing momentum transfer. The pion mass dependence is represented in linear graphs where yellow areas refer to the region M π L < 2. The solid lines (θ = 0) reach that region when -starting from the right-hand side of the figure -they first touch the dark yellow area. The dotted lines (θ = π/3) reach this region when they enter in the light yellow area. We observe that the corrections decay exponentially as O(e −λπ ) and are mainly negative. They may turn positive depending on the pion mass. In the region where the p-regime is guaranteed, the corrections are less than the percentage level and thus negligible.
To estimate the corrections at a non-zero momentum transfer we consider the incoming pion at rest (i.e. p = 0) and the outgoing pion moving along the first axis carrying the first non-zero momentum (i.e. | p | = 2π/L). This kinematics provides the momentum transfer, The zeroth component corresponds to the energy transfer among external pions, For the configuration (6.4) the pion energies take the following forms, E π + (L) = M 2 π ± (L) + (2π + θ) 2 /L 2 + 2 (2π + θ)∆ϑ/L + O(ξ 2 π ). (6.8) Here, M π 0 (L), M π ± (L) are the pion masses in finite volume and for SU(2) ChPT, their corrections are given from eq. (3.7) discarding the contributions of kaons and the eta meson. The quantity, m(n, n 1 ) n 1 e in 1 θ , (6.9) corresponds to the first component of the extra term (3.20) in the configuration (6.4). We refrain from presenting numerical results of the square radius as according to ref. [65] it is more effective to correct the matrix elements of form factors and from those, extract the square radii. In figure 8a (resp. 8b) we represent the pion mass dependence of −δΓ π 0 S (resp. −δΓ π + S ) at q 2 = q 2 min . The logarithmic graphs illustrate the exponential decay O(e −λπ ) of the corrections. The corrections are mainly negative. For θ ∈ { 0, π/8, π/4, π/3 } the absolute value of δΓ π 0 S increases (resp. that of δΓ π + S decreases) with the angle. Note that the corrections reach the percentage level before entering yellow areas and they should be subtracted when the scalar form factor is extracted from lattice data.   Table 4. Comparison of corrections evaluated at NLO with ChPT and estimated beyond NLO with asymptotic formulae. Here, we use L = 2.83 fm so that M π L = 2 for M π = M phys π ± . The uncertainties of R(M π 0 ) originate from the errors of LECs contained in the integrals I (4) (M π 0 , π 0 ), I (4) (M π 0 , π ± ). Note that results with M π < 0.140 GeV should be taken with a grain of salt as they are in the region where the p-regime is no more guaranteed.

Masses and extra terms of self-energies
In figure 9 and 10 we represent the mass corrections estimated with R(M π 0 ), R(M π ± ). The pion mass dependences of R(M π 0 ), R(M π ± ) are represented in logarithmic plots whereas the angle dependences in linear ones. In these graphs the bands of the uncertainty are displayed for solid lines (i.e. for θ = 0 in figure 9 resp. for M π = 0.140 GeV in figure 10). Bands of different colors refer to different values of L. The uncertainty bands are calculated with the usual formula of the error propagation. As unique source of error, we have taken the uncertainties on the LECs contained in the integrals I (4) (M π 0 , π 0 ), . . . , I   figure 9b the lines are exactly straight and are so close that they overlap in the graph. This indicates that for small angles, R(M π 0 ) is more sensitive to θ. The corrections are significantly bigger than those at NLO. As illustration, we list in table 4a (resp. table 4b) numerical values for δM π 0 , R(M π 0 ) at L = 2.83 fm for θ = 0 (resp. for θ = π/3). In some cases, R(M π 0 ) is of order 50% with respect to δM π 0 . Such significant subleading effects were already observed in finite volume with PBC, see ref. [40]. However, the comparison of numerical results obtained from asymptotic formulae with amplitudes at two loops [13,14] has showed that subsubleading effects are small and that the expansion does have a good converging behaviour for M π L > 2. We are confident that this is also true for TBC and that our numerical estimates are reliable, at least for small angles (i.e. θ < π).
In figure 10 we display the angle dependences of R(M π 0 ), R(M π ± ). We observe that R(M π 0 ) depends on θ as a cosine function. The corrections oscillate with a period of 2π and have maxima (resp. minima) at even (resp. odd) integer multiples of π. If we consider M π = 0.197 GeV, the difference among maxima and minima is 5.4% at L = 2 fm. This is a sizable effect and should be taken into account when physical observables are extrapolated from lattice data. In figure 10b we observe that R(M π ± ) depends on θ as (a + cos θ) √ 1 + θ 2 with a > 0. This dependence originates from contributions O(ξ 2 π ) and provides large corrections at large angles. However, one should here retain only results in the interval θ ∈ [0, π]. The reason is, R(M π ± ) is derived by means of an expansion which is valid for small external twisting angles. Considering M π = 0.197 GeV, the difference among the local maximum and the local minimum in θ ∈ [0, π] is 0.2% at L = 2 fm. This is a negligible effect. Note that the corrections estimated with R(M K ± ), R(M η ) are less than a percent for small angles and can be neglected.
In figure 11 we plot the extra term ∆ϑ µ Σ π + estimated by means of R µ (ϑ Σ π + ). We represent the first spatial component as it is the only one which is non-zero for configuration (6.4). The values on the y-axis are in GeV since R µ (ϑ Σ π + ) is a dimensionful quantity. Uncertainty bands are displayed for θ = π/8 in figure 11a and for M π = 0.140 GeV in figure 11b. From the logarithmic graph of figure 11a we observe that the extra term decays exponentially as O(e −λπ ). For θ ∈ { π/8, π/4, π/3 } its absolute value increases with the angle. This can also be seen in figure 11b where R µ (ϑ Σ π + ) is represented as a function of θ. We observe that the extra term depends on θ almost exactly as a (negative) sine function. The zeros correspond to integer multiples of π and extrema are close to half-integer multiples of π. In this graph, one should retain only results for θ ∈ [0, π] as by derivation, R µ (ϑ Σ π + ) is valid for small external angles. Note that R µ (ϑ Σ K + ) has similar dependences on M π resp. θ and its absolute value is in general smaller than that of R µ (ϑ Σ π + ).

Decay constants and extra terms in axialvector matrix elements
In figure 12 we represent the corrections of decay constants estimated with R(F π ± ), R(F K ± ). The logarithmic graphs illustrate the exponential decay O(e −λπ ) of the corrections. The corrections are negative and for θ ∈ { 0, π/8, π/4, π/3 } their absolute values decrease with the angle. Note that the corrections may reach more than 10% before entering yellow areas and they should be subtracted before decay constants are extracted from lattice data.     Figure 13. Pion mass dependence of −R µ (ϑ A π + ) for µ = 1. Note that the dark yellow area refers to the region M π L < 2 for θ = π/8.
In figure 13 we represent the extra term ∆ϑ µ A π + estimated with R µ (ϑ A π + ). We represent the first component as it is the only one which is non-zero for configuration (6.4). Also this extra term decays exponentially with λ π . Comparing figure 13 with figure 11a we observe that for a fixed angle −R µ (ϑ A π + ) is smaller than −R µ (ϑ Σ π + ). This difference is O(ξ 2 π ) and is proportional to the extra term R µ (ϑ G π + ). A similar observation can be made for charged kaons where the difference among R µ (ϑ A K + ) and R µ (ϑ Σ K + ) is O(ξ 2 π ) and is proportional to R µ (ϑ G K + ).

Pseudoscalar coupling constants
In figure 14 we represent the pion mass dependence of −R(G π 0 ) which is as well exponential. In general, the corrections are negative and for θ ∈ { 0, π/8, π/4, π/3 } their absolute value increases with the angle. In this case, the corrections are smaller than those at NLO. This can be explained if we look at the contributions O(ξ 2 π ), see eq. (5.45). For θ ∈ { 0, π/8, π/4, π/3 } the contribution originating from integral I (4) (G π 0 , π 0 ) is negative but that from I (4) (G π 0 , π ± ) is positive. As the negative contribution is smaller than the positive one, the corrections estimated with −R(G π 0 ) are smaller than those evaluated with −δG π 0 at NLO.

Pion form factors
We consider pions at rest and estimate the corrections of the matrix elements of the scalar form factor with R(Γ π 0 S ), R(Γ π + S ). In figure 15a [resp. 15b] we represent the pion mass  dependence of R(Γ π 0 S ) [resp. R(Γ π + S )] at q 2 = 0. From the graphs we observe that the corrections decay exponentially as O(e −λπ ). In general, they are negative but they may turn positive depending on the pion mass. Note that the corrections reach percentage level before entering yellow areas and they should be subtracted when the scalar form factor is extracted from lattice data.

Conclusions
In this work we studied the effects of a finite cubic volume with twisted boundary conditions on pseudoscalar mesons. We applied chiral perturbation theory (ChPT) in the p-regime and introduced twisting angles by means of a constant vector field, see ref. [4]. The corrections for masses, decay constants, pseudoscalar coupling constants were recalculated at next-to-leading order (NLO) and new results for pion form factors were presented. In the calculations we adopted the mass definition of ref. [4,7] which treats new extra terms as renormalization terms of twisting angles, and argued in some detail about the reasons behind this choice. These extra terms originate from the breaking of the cubic invariance and can be reabsorbed in the on-shell conditions modifying the mass definition in finite volume. We found that the Feynman-Hellmann theorem [18,19] as well as the Ward-Takahashi identity [20][21][22] are satisfied. To prove the Ward-Takahashi identity we constructed an effective field theory for charged pions invariant under gauge transformations which reproduces results obtained with ChPT.  We generalized the derivation of Lüscher [11] and derived asymptotic formulae for twisted boundary conditions. We showed that the asymptotic formulae for masses, decay constants, pseudoscalar coupling constants are related by means of chiral Ward identities where extra terms satisfy such relations in an independent way. Applying asymptotic formulae in combination with ChPT, we estimated corrections beyond NLO and found that, as in the case of PBC, NNLO corrections can be very significant, indeed almost as large as NLO corrections. This underlines the importance of using asymptotic formulae combined with NLO chiral calculations of the relevant infinite-volume amplitude to reliably estimate finite-volume corrections. From our numerical analysis we see that the corrections can be comparable (or even larger) than the statistical precision reached in simulations of lattice QCD and hence, should be taken into account.

A Sums in finite volume
We list some results which are useful in the evaluation of loop diagrams in finite volume. For convenience, we define where g is a generic function in momentum space and L is the side length of the finite cubic box. The right-hand side of the equation represents the difference among contributions in finite and infinite volume. For loop diagrams encountered in this work, this difference is finite and can be calculated by means of the Poisson resummation formula [30]: The first group of results is g µν 2 g r−1 (λ P , ϑ) + h µν r (λ P , ϑ) , for r ∈ N. Here, Γ(r) is the gamma function, g µν is the metric of Minkowski space-time, ϑ µ = 0 ϑ is a twisting angle and λ P = M P L. The functions on the right-hand side can be expressed in terms of modified Bessel functions of the second kind, K r (x). They read with n µ = 0 n . Corrections of masses and decay constants were calculated using the results (A.3) for r = 1.
To evaluate loop diagrams with two different propagators one can use the Feynman parametrization, Here, we consider A = M 2 P − (k + ϑ) 2 and B = M 2 P − (k + ϑ + q) 2 where q µ is an external momentum. The second group of results we present, is with λ z = M P L 1 + z(z − 1)q 2 /M 2 P . The functions on the right-hand side can be ex-pressed as For q 2 = 0 the functions g r (λ z , q, ϑ), f µ r (λ z , q, ϑ), h µν r (λ z , q, ϑ) reduce to the expressions of eq. (A.4). Note that if q = 2π L l with l ∈ Z 3 , the results (A.6) can be simplified by means of substitutions z → (1 − z) and n → − n. This leads to the results of section 3.2 and appendix B.2.

B Gauge symmetry in finite volume
To explain the results of section 3.2 we construct an an effective theory for charged pions which is invariant under electromagnetic gauge transformations. The theory reproduces the expression obtained at vanishing momentum transfer and indicates that the gauge symmetry is preserved in this case. Relying on this observation, we show that the Ward-Takahashi identity [20][21][22] holds in finite volume as long as the momentum transfer is discrete. Only the differential form of the identity -the Ward identity [35] -is violated due to the discretization of the spatial components. These considerations were presented for PBC in ref. [38] and are here generalized to TBC.

B.1 Construction of a gauge invariant effective theory
We consider a finite cubic box of side length L on which we impose TBC. In presence of two light flavors, we can introduce the electromagnetic gauge field through the external vector field, v µ = −eA µ (x) Q.
Here, e is the elementary electric charge of the positron and Q = diag(2/3, −1/3). As long as Q is diagonal we may redefine the fields so that they are periodic and introduce the twisting angles by means of a constant vector field. Since A µ (x) as well as other fields are periodic, we can proceed in a similar way as in ref. [38] to construct the effective theory. The only difference is that the effective Lagrangian contains additional terms due to the constant vector field proportional to the twisting angles. At low energies, the relevant degrees of freedom are pions and for simplicity, we just consider the charged ones in the following. In absence of the electromagnetic interaction the Lagrangian of the effective theory reads and M π ± (L) is the mass of charged pions in finite volume. The kinetic term contains the The constant vector field w µ ϑ is proportional to the third Pauli matrix, τ 3 = diag(1, −1) and introduces the twisting angle ϑ µ π + as well as the extra term ∆ϑ µ Γ π + . Here, ϑ µ π + , ∆ϑ µ Γ π + break Lorentz invariance. For ϑ µ π + → 0 the field w µ ϑ disappears and the cubic invariance is restored: in this case the theory respects PBC. Note that M π ± (L), ∆ϑ µ Γ π + implicitly depend on parameters of the effective theory (like the LECs).
To add the electromagnetic interaction we must include all possible operators which are invariant under gauge transformations. This can be achieved using Wilson loops, see ref. [38]. We limit ourselves to include operators containing the zero mode of the gauge field A µ (x) as these allow us to study the electromagnetic form factor at vanishing momentum transfer. Proceeding in a similar way as ref. [38] we obtain the following effective Lagrangian in presence of the electromagnetic interaction, is constructed from Wilson loops, see ref. [38]. The expression (B.4) needs some explanations. The dots at the end indicate that we have just written down the relevant terms of the effective Lagrangian. The most general effective Lagrangian contains terms with arbitrary many insertions of W µ − , which we are not writing explicitly. The expansion of W µ − starts with a term linear in the zero mode which allows us to study the electromagnetic form factor at vanishing momentum transfer. The tensor Q(L) µν breaks the Lorentz as well as the cubic invariances and must be determined by matching. For ϑ µ π + = 0 we expect that Q(L) µν reproduces the result for PBC [38] and that it disappears for L → ∞.
The effective theory of eq. (B.4) reproduces the expression of the vector form factor at vanishing momentum transfer. The presence of Wilson loops ensures that the theory is invariant under gauge transformations. As long as A µ (x) is periodic this invariance is preserved. Starting from this observation, we show in appendix B.2 that the Ward-Takahashi identity holds for TBC and that the corrections to the vector form factor are related to inverse propagators.

B.2 Ward-Takahashi identity
In infinite volume gauge symmetry implies that the electromagnetic vertex function Γ µ satisfies the Ward-Takahashi identity [20][21][22]: (B.7) Here, q µ = (p − p) µ is the momentum transfer, ∆(p ) resp. ∆(p) are the propagators of outgoing and incoming particles and Q e = Q/e is the electric charge of external particles in units of the positron charge. In the limit q µ → 0, the identity tends to a differential form, known as Ward identity [35]: For external charged pions, we can calculate the electromagnetic vertex function from the matrix elements, iΓ µ = π ± (p )| V µ 3 |π ± (p) . In ref. [33] these matrix elements are evaluated in ChPT at NLO and amount to whereJ(q 2 ) is the finite part of the loop integral (C.11). Here, we display all terms, even those that disappear as external momenta are on-shell. For on-shell momenta only the term proportional to (p + p) µ contributes and provides the vector form factor, see ref. [33]. One can show that the vertex function (B.9) satisfies the Ward-Takahashi identity by contracting with q µ and arranging the surviving terms in inverse propagators. Taking the limit q µ → 0, the same vertex function satisfies the Ward identity, indicating that the electromagnetic current as well as the electric charge are conserved.

(B.13)
if q µ is non-vanishing and if q = 2π L l with l ∈ Z 3 \ { 0 }. These relations can be obtained by partial integration and by using properties of the derivatives of the modified Bessel functions of second kind.

C Terms S (4)
We list some explicit expressions of the terms S (4) introduced in section 5 indicating the equation where they appear. Other terms S (4) can be found in appendix A of ref. [14].

C.2 Kaons
We list the terms S (4) appearing in the asymptotic formulae for kaons. Note that the functions S kl P Q are defined in eq. (C.16). In the next expressions, we denote the ratio of mass squares as x P Q = M 2 P /M 2 Q for P, Q = π, K, η. The functions S k,l P Q entering the above expressions are defined as S k,l P Q = S k,l P Q (λ π | n|) (C.16)
Here, N = (4π) 2 and g (1) P Q (x) =J P Q (x), g P Q (x) = M 2 KJ P Q (x), g P Q (x) = M 4 KJ P Q (x), g P Q (x) = K P Q (x), g P Q (x) = M 2 K K P Q (x), g P Q (x) = M 4 K K P Q (x), g P Q (x) =M P Q (x), g P Q (x) = M 2 KM P Q (x), g P Q (x) = M 4 KM P Q (x). (C.17) The explicit forms of g (l) P Q (x) were presented in ref. [51]. They can be expressed in terms of the loop-integral functionJ P Q (q 2 ) = J P Q (q 2 ) − J P Q (0) evaluated in d = 4 dimensions, We conclude with a remark on the use of the loop-integral functions in the asymptotic formulae. The loop functions (C.19) need to be evaluated for complex values of their arguments. ForJ P Q (M 2 P +M 2 Q +2iM P M Q y) there is an ambiguity due to the negative value of ρ = −4M 2 P M 2 Q (1 + y 2 ), which (C. 19) does not resolve explicitly. An explicit analytic continuation was provided in ref. [14] but unfortunately was not correct. The correct prescription is as follows: take the positive value of the square root √ ρ = 2iM P M Q 1 + y 2 for which the logarithm in (C. 19) becomes (t = M 2 P + M 2 Q + 2iM P M Q y) (C.20)