Twisted Boundary Conditions in Lattice Simulations

By imposing twisted boundary conditions on quark fields it is possible to access components of momenta other than integer multiples of 2pi/L on a lattice with spatial volume L^3. We use Chiral Perturbation Theory to study finite-volume effects with twisted boundary conditions for quantities without final-state interactions, such as meson masses, decay constants and semileptonic form factors, and confirm that they remain exponentially small with the volume. We show that this is also the case for"partially twisted"boundary conditions, in which (some of) the valence quarks satisfy twisted boundary conditions but the sea quarks satisfy periodic boundary conditions. This observation implies that it is not necessary to generate new gluon configurations for every choice of the twist angle, making the method much more practicable. For K->pipi decays we show that the breaking of isospin symmetry by the twisted boundary conditions implies that the amplitudes cannot be determined in general (on this point we disagree with a recent claim).


Introduction
In lattice simulations of QCD on a cubic volume (V = L 3 ) with periodic boundary conditions imposed on the fields, the hadronic momenta (p) are quantized according to p i = 2π/L × n i , where i = 1, 2, 3 and the n i are integers. For currently available lattices this means that the lowest non-zero momentum is large (typically about 500 MeV or so) and there are big gaps between neighbouring momenta. This limits the phenomenological reach of the simulations. In ref. [1] Bedaque proposed the use of non-periodic boundary conditions which would allow hadrons with arbitrarily small momenta to be simulated (see also the references cited in [1] for earlier related ideas). We refer to these boundary conditions as twisted boundary conditions 1 . This technique has subsequently been used in a quenched study of the energy-momentum dispersion relations of pseudoscalar mesons [3] and the finite-volume corrections for two-particle states with twisted boundary conditions have been calculated [4].
In this letter we use chiral perturbation theory (χPT) to analyse some of the properties of twisted boundary conditions and show that: 1. For physical quantities without final state interactions, such as masses or matrix elements of local operators between states consisting of the vacuum or a single hadron, the flavour symmetry breaking induced by the twist only affects the finite-volume corrections, which nevertheless remain exponentially small.
2. For amplitudes which involve final-state interactions, such as those for K → ππ decays, in general it is not possible to extract the physical matrix elements using twisted boundary conditions (see sec. 4). On this point we disagree with ref. [4]. 3. For Partially Twisted boundary conditions, in which (some of) the valence quarks satisfy twisted boundary conditions but the sea quarks satisfy periodic ones, one also obtains the physical quantities described in item 1 with exponential precision in the volume. This implies that in unquenched simulations it is not necessary to generate new gluon configurations for every choice of boundary condition, thus making the method much more practicable.
In ref. [5] Kim and Christ propose H-and G-parity boundary conditions in which the minimum non-zero hadronic momenta are reduced from 2π/L → π/L (see also ref. [6]). These authors impose H-parity boundary conditions for K → ππ decays in which the two-pions are in an I = 2 state. This is a particular case of twisted boundary conditions, corresponding to the specific choice of π for the twisting angle (as stated in item 2 above and explained in sec. 4 below, it is not possible to study K → ππ decays with a general choice of twisting angle). Kim and Christ also show that G-parity boundary conditions, which exploit the discrete charge conjugation transformations, can be used for an I = 0 two-pion state (in unquenched simulations), but the formalism will have to be extended to incorporate the kaon. Although we do comment below on H-and G-parity boundary conditions in order to illustrate our discussion, the main focus of this letter is on boundary conditions based on continuous symmetries.
When considering K → ππ decays, throughout this letter we restrict our discussion to the centre-ofmass frame for the two pions. For this case and with periodic boundary conditions, the finite volume corrections which decrease as powers of the volume have been derived for the two-pion spectrum [7] and matrix elements [8,9]. At present the theory of finite-volume corrections in a moving frame has not been developed for matrix elements (but for a discussion of finite volume corrections to the two-pion spectrum in a moving frame see ref. [10]). We therefore do not generalise our discussion to the moving frame at this stage.
The plan of the remainder of this letter is as follows. In the following section we define twisted boundary conditions in QCD and briefly review their properties. In sec. 3 we impose twisted boundary conditions on the chiral lagrangian and demonstrate that their effect is to shift the momenta of internal propagators and external mesons by amounts corresponding to the twists. Sec. 4 contains a discussion of finite-volume effects when twisted boundary conditions have been imposed. We discuss partially twisted boundary conditions in sec. 5 and present our conclusions in sec. 6. There are two appendices in which we derive the finite-volume corrections with twisted boundary conditions at one-loop order in χPT (appendix A) and present the corresponding results for meson masses and decay constants (appendix B).

Twisted Boundary Conditions in QCD
In this section we will define the twisted boundary conditions and derive some of the constraints they have to satisfy. Since the choice of boundary conditions is a non-local effect, we can present the discussion, without any loss of generality, within the framework of continuum quantum field theory. It should be noted however, that the discussion also applies to every lattice discretization. Although local discretization artefacts may affect the constant pre-factors, they do not affect the functional behaviour with the volume. For definiteness we present the discussion in Euclidean space with an infinite time dimension and a finite cubic spatial volume of size L 3 .
When formulating quantum field theory in a finite cubic volume, in order to avoid boundary terms, periodic boundary conditions are frequently imposed on the fields. This is equivalent to defining the theory on a torus and the periodicity of the fields ensures that the fields are single valued. However requiring that the fields be single valued is not necessary; it is sufficient instead to require that the observables be single valued, which is equivalent to the condition that the action be single valued on the torus. This means that the generic field φ has to respect the following boundary conditions: where the index i is not summed and U i is a symmetry of the action. Imposing the condition in eq. (1) is sufficient to cancel the boundary terms.
Consider now the Dirac term in the (Euclidean) QCD Lagrangian, where for our discussion it will be convenient to consider the quark field q(x) to be a vector in flavour space and the quark mass matrix M to be a diagonal matrix. The possible boundary conditions depend on the symmetries of the action, and in particular on the form of M , i.e. on whether there is any degeneracy assumed for the quark masses. Here we will consider the most general continuous symmetry, i.e. the flavour symmetry group U (N ) V and its subgroups, and will not discuss the possible use of discrete symmetries (and charge conjugation in particular [5]). An advantage of the use of continuous symmetries is that the minimum momentum can take any value less that 2π/L, whereas with discrete symmetries such as G-parity the lowest momentum is π/L. Eq. (1) then implies that U i has to commute with the Dirac operator, and in particular with the quark mass matrix. For general values of the quark masses this implies that U i is a diagonal matrix. In the isospin limit one could in principle, take U i to be non-diagonal in the u -d sector, however this choice breaks the conservation of electric charge, and will not be considered explicitly here. We can therefore write the boundary condition for the quark fields in the form: where the T a 's are the generators in the Cartan subalgebra of the flavour U (N ) V group commuting with the quark mass matrix. It is convenient to redefine the quark fields by: The fieldsq(x) satisfy periodic boundary conditions, and the Lagrangian (2) is given in terms of these fields by: whereD This is the Lagrangian of QCD with quark fields satisfying periodic boundary conditions interacting with a constant background vector field which couples to quarks with charges determined by the phases in the twisted boundary conditions. The external field, in addition to breaking the cubic group symmetry, breaks also all the symmetries which do not commute with it. For generic diagonal B i the broken symmetries are flavour SU(3) and I 2 , but not I z , strangeness and the electric charge.
To illustrate some of the above points and the kinematic nature of the symmetry breaking induced by the twisted boundary conditions, we end this section by exhibiting the propagator of a free quark using both the q andq definitions. For compactness of notation we drop the flavour index and take B = θ/L to be the twist corresponding to the flavour represented by q with mass M . The propagators are then where in both cases the sum is over momenta k = (2π/L) n and n is a vector of integers. S(x) satisfies the twisted boundary condition (S( x,0 whereasS(x) satisfies periodic boundary conditions (S(x i + L) =S(x)) and x,0 . The momentum in the denominators is shifted (or boosted) by θ/L.

The Effective Chiral Lagrangian
In this section we derive the low-energy effective Lagrangian for QCD in the presence of twisted boundary conditions and study its properties. The derivation of the chiral Lagrangian could be performed directly by coupling the Gasser-Leutwyler Lagrangian to the external vector field B µ introduced above. Here instead, we choose to follow the steps of sec. 2 in order to show the equivalence with QCD explicitly. Imposing the boundary conditions of eq. (1) on the fields, implies that: where Σ is the coset representative of SU . Note that this relation is completely fixed once the quark boundary conditions are imposed, and the results below are implied unambiguously by this relation. Following the presentation in sec. 2, we redefine the fields by: so thatΣ satisfies periodic boundary conditions. Eq. (11) corresponds to a local symmetry transformation so that only the derivative terms are affected: In terms ofΣ the chiral Lagrangian reads and represents the trace. The Lagrangian in eq. (15) is the standard χP T Lagrangian with periodic fields coupled to the vector external field B µ . Note that the long-distant nature of the boundary conditions allows χP T to take their effects completely into account through the simple modification in eqs. (15) and (16). The low energy constants are not affected by the twist (analogously to the arguments in [11]).
The effects of the twist on the mesons can be obtained directly from eq. (16). The neutral-meson fields commute with B µ (recall that Θ is diagonal) and do not receive any boost. The charged-meson fields, on the other hand, are boosted by an amount proportional to the difference of the twists of the two valence quarks (v 1 and v 2 ): and the spectrum of allowed meson momenta is shifted accordingly, both in external states and in internal propagators. Figure 1: Auxiliary one-loop diagram with n external mesons, used in the demonstration that the effect of twisting is to shift the internal and external momenta accordingly. The unprimed and primed variables represent the external and internal lines respectively.
From the chiral Lagrangian in eq. (15), and its extensions to higher order in the momentum expansion, it follows that the only effect of the twisted boundary conditions is to shift all the momenta consistently in order to recover the correct boost corresponding to the flavour of each propagator and external line. We illustrate this by considering the loop contribution represented in fig. 1, which may be a one-loop contribution to an n-body process or an insertion in a higher-order diagram. The contribution from this diagram is of the form: where the sum is over momenta k i = (2π/L) n i and the {n i } are integers 2 . The factor in the numerator represents the derivative terms at vertices in the chiral Lagrangian and B l ′ refers to the momentum shift due to the external field on the meson in the l ′ propagator of the loop. We now perform the trivial change of variables k → k ′ = k + B n ′ to rewrite the sum in eq. (18) as where now the sum is over momenta k ′ i = (2π/L) n i + B n ′ i with integer n i . B i is the twist corresponding to the flavour of the external line i and we have used the fact that at each vertex the sum of the twists is zero (e.g. . This condition is a consequence of the invariance of the action under the twist transformations. So far we have considered the theory on a single volume, where there are finite-volume artefacts, and in sec. 4 below we investigate the size of these corrections (which do depend on the boundary conditions which have been imposed). In phenomenological applications we generally wish to eliminate FV corrections by taking, in principle at least, the infinite-volume limit, so that the sum in eq. (18) goes over into the integral .
(20) As required, this is precisely the expression for the contribution from this diagram with external momenta p i + B i . For fixed external momenta (p i + B i in the notation of eq. (20)), the integral is independent of the boundary conditions which have been used in the finite-volume calculations.
We conclude this section with some brief comments about the way that the infinite volume limit might be taken in principle. We start of course by studying a physical quantity in a finite-volume. For illustration imagine that the process depends on a component of momentum p i which is smaller than 2π/L for a particular lattice simulation and so we introduce a twist θ i for the corresponding flavour in direction i. Now we can envisage taking the infinite volume limit in a number of ways. For example we may keep θ i fixed so that B i → 0 as the volume is increased. p i is kept fixed in physical units, and since p i = (2π/L) n i + θ i /L for some integer n i , as we increase the volume we take higher excitation levels n i . The effect of the twist decreases as the volume increases, and the results approach those obtained with periodic boundary conditions. This is also true for momentum sums such as that in eq. (18), which are dominated by momenta of order of some physical scales and hence the relevant n i increase as L increases. Thus again we see that the effects of the twist decrease as the volume is increased. This feature is generally true as long as the infinite-volume limit is taken keeping the physics fixed.

Finite Volume effects with Twisted Boundary Conditions
Finite-volume corrections in general, and those due to the choice of boundary conditions in particular, are long-distance effects which can be studied using χPT (for sufficiently light pseudo-Goldstone mesons and large volumes). We start by considering processes without any final-state interactions, such as particle masses or matrix elements of local operators with external states which consist of either the vacuum or a single hadron. For these quantities finite-volume corrections are known to be exponentially suppressed with the volume, due to the fact that in the absence of branch cuts (which is the case for these quantities), the Poisson formula allows us to replace the sums over the discrete momenta in finite volume by infinitevolume integrals. Differences between the two are exponentially small in the volume and this remains true in the presence of twisted boundary conditions. As shown in appendix A, the finite-volume correction can be calculated in terms of elliptic-ϑ functions, and decrease exponentially at large volumes.
We now report the asymptotic finite-volume corrections (in the limit L → ∞) for pion masses and decay constants with twisted boundary conditions; the results demonstrate explicitly the isospin breaking at finite volume. For each physical quantity we present the results in the form where X(L) and X(∞) are the results in finite and infinite volume respectively. The full expressions for the finite-volume corrections (at NLO in χPT) can be found in appendix B, and their asymptotic behaviour as L → ∞ is as follows: For matrix elements involving two or more mesons in the final state the situation is more complicated: there are both exponential and power finite-volume corrections. The latter are parametrically larger (and in most cases numerically larger). Two-particle energy shifts due to the finite volume, Lellouch-Lüscher factors relating finite-volume matrix elements to physical amplitudes and finite-volume corrections to the two-particle interpolating operators at the sink which contain final state interactions, arise as power corrections in the volume. If twisted boundary conditions affect these terms, then they inevitably affect unitarity with obvious consequences for the extraction of the matrix elements.
Consider, for instance, the case of K → ππ decays, and imagine that only the u-quark is twisted, so that the charged pions are boosted but not the neutral ones. In such a situation I 2 is no longer a good quantum number, so that the energy eigenstates are no longer states with a definite isospin I; in particular they are no longer states with I=0 or I=2 as is the case when isospin is a symmetry. This can be seen even in the free-theory. Since we require that the two-pion state is at rest, at tree-level the lowest energy state is |π 0 π 0 with both pions at rest and the first excited state is |π + π − with the momenta of each of the two pions p π ± = ± θ/L, where θ is the twist on the up-quark. In the interacting theory the presence of π + π − ↔ π 0 π 0 transitions complicates the analysis very significantly and, as explained in the next paragraph, it is not possible to determine physical K → ππ amplitudes from simulations on finite volumes using twisted boundary conditions with the power corrections in the volume kept under control. These issues were not considered in ref. [4] and we therefore do not agree that the formulae for finite-volume corrections presented in that paper can be applied to K → ππ decays.
We now briefly demonstrate the difficulties in studying quantities involving two-pion states using boundary conditions which break isospin invariance. Consider, for example, the correlation functions 0|π 0 (t)π 0 (t) σ(0)|0 and 0|π + (t)π − (t) σ(0)|0 , where σ is some operator which can create two pions from the vacuum and π i is an interpolating operator for a pion with charge i. σ is placed at the origin and we have taken a Fourier transform at zero momentum of each of the π i fields so that only their time dependence is exhibited (of course the boundary conditions induce a momentum of O( Θ/L) for charged pions). Fitting the measured behaviour to two exponentials we would find: where, at tree-level in chiral perturbation theory, E 0 = 2m π and E 1 = 2 m 2 π + p 2 π± . The ellipses represent terms with higher energies and we assume here these can be neglected. By fitting the correlation functions above, the constants A 00 , A +− , B 00 , B +− can, in principle at least, be determined numerically and we would then know which combinations of the two-pion operators have no overlap with states with energies E 0 and E 1 respectively (we call these states |s 0 and |s 1 and denote by Π 2 0 and Π 2 1 the operators with no overlap with |s 1 and |s 0 respectively). Note that in order to extract B 00 and B +− it is necessary to include the non-leading exponential in the fit, which eliminates a major potential advantage of using twisted boundary conditions for K → ππ decays. Combining the results from these fits, together with those of four-pion correlation functions of the form 0|π(t)π(t) π(0)π(0)|0 , we can determine the matrix elements 0|Π 2 0 |s 0 , 0|Π 2 1 |s 1 , 0|σ|s 0 and 0|σ|s 1 . Unfortunately, even if we are able to carry out the procedure described above with reasonable accuracy, it is still not clear how to relate the finite-volume eigenstates |s 0 and |s 1 (which have different energies) to the infinite-volume eigenstates |(ππ) I=0 and |(ππ) I=2 since the known procedures for doing this [7,8,9] rely on isospin symmetry. In some respects this problem resembles the one of extending the discussion of refs. [7,8,9] above the KK threshold. We conclude that new ideas would be necessary before K → ππ matrix elements could be determined with twisted boundary conditions. In order to overcome the above difficulties one should introduce twisted boundary conditions which preserve isospin symmetry and this is not possible in general. Christ and Kim [5] however, have pointed out that one can make some progress for (ππ) I=2 states if one restricts the twist angle to π (they call this case H-parity). The K + → π + π 0 matrix element can be related by the Wigner-Eckart theorem to a matrix element of a ∆I z = 3/2 operator into a π + π + final state. By choosing θ = (π, 0, 0) for the down quark and 0 for the up quark and performing the Fourier transforms over the positions of the two pions with weights 1 and exp{i(−2π/L)x} respectively, we obtain a π + π + final state with the two pions having momenta π/L and −π/L in the x direction (hence remaining in the centre-of-mass frame). With this procedure we are restricted to θ x = π but the need for extracting terms corresponding to non-leading exponentials is avoided. Note that this procedure is possible, because the required matrix element can be related to one in which the final state only contains π + mesons. For I = 0 final states this is not possible, although in ref. [5] it is also shown that by introducing discrete (G-parity) boundary conditions one can treat I = 0 two-pion states with the two pions having momenta ±π/L (but additional ideas will have to be introduced to incorporate kaon states at rest into the formalism).

Partially Twisted Boundary Conditions
Until now we have assumed that the twisted boundary conditions are applied consistently to both the valence and sea quarks. In lattice simulations this implies that a new set of gauge configurations must be generated for each choice of the twist. In addition, if different twists are imposed on the u-and d-quark fields then one must use formulations of lattice fermions for which the light quark determinant is positive definite for each flavour. It would clearly be very welcome if one could avoid new simulations for every value of Θ and in this section we analyse the consequences of introducing different boundary conditions for sea and valence quarks. In particular we consider the case in which the valence quarks satisfy twisted boundary conditions and the sea quarks satisfy periodic boundary conditions. In this case the QCD Lagrangian can be conveniently written as: where the subscripts v, g, s stand for valence, ghost and sea and q g are commuting spinors. Moreover in order to have a cancellation of valence loops we require that D g = D v and M v = M g . Eq. (27) can be rewritten in the form and now both B µ and M take values in the graded algebra of The derivation of the Feynman rules for both QCD and χP T is standard and we do not present it here. Having different twists for valence and sea quarks breaks the valence-sea symmetry. This is clearly a finite-volume effect, but the relevant question is whether the corrections induced by this asymmetry decrease like powers of the volume or exponentially 4 . We find that the situation is analogous to the violation of unitarity in partially quenched QCD [14] and that for many physical quantities (including those with at most a single hadron in the initial and final states) the use of partially twisted boundary conditions induces errors which are exponentially small. This is because sea quarks appear in loops and the sums over the loop-momenta can be approximated by integrals with exponential precision.
We start by considering processes with at most one hadron in the external states. As long as the shift does not induce cuts in the correlation function, the correction is still exponentially suppressed in the volume. For example, if we consider an unquenched simulation with three flavours, in the asymptotic limit, we find (in the notation defined in eq. (21) ): for the three cases in which the d and s quarks satisfy periodic boundary conditions but the up quark is (a) untwisted, (b) fully twisted (both valence and sea u-quarks satisfy twisted boundary condition) and (c) partially twisted (only the valence u quark is twisted). This shows that, in general, finite volume corrections could be different for the three cases but they are always exponentially small.
In [15] Golterman and Pallante demonstrated that (partial) quenching can induce "unphysical" mixing among weak operators because of their different transformation properties under the flavour group and its graded extension. These effects are proportional to the difference M s − M v and have to disappear in full QCD. One could also ask whether imposing different boundary conditions for sea and valence quarks could lead to similar effects. Such effects are proportional to θ v − θ s and again appear as exponentially small finite-volume corrections.
Not surprisingly the case of amplitudes with multiparticle external states is much more complicated. We have seen in sec. 4 that it is not possible to isolate ππ states with a given isospin using periodic boundary conditions. We therefore restrict our consideration here to the H-parity and G-parity cases for I = 2 and I = 0 two-pion states respectively. Since now the twist angle is fixed to be π, the practical advantage of using partial twisting to avoid generating new gluon configurations for every value of the angle is much less compelling, but it is interesting nevertheless to consider the theoretical issues. The effects of different boundary conditions for sea and valence quarks are analogous to those discussed in [14] for the partially quenched theory. For I = 2 the sea quarks do not enter the FSI directly and the difference between imposing H-parity boundary conditions fully or partially is exponentially small. (The power corrections in the volume with H-parity boundary conditions are different of course from those with periodic ones, but they are calculable [4,10].) For I = 0 and partial G-parity boundary conditions on the other hand, the mesons in intermediate states in correlation functions necessarily include both sea and valence quarks, whereas the external states are made of valence quarks only. The lack of degeneracy between external and internal states implies a breakdown of unitarity and Watson's theorem and we are then unable to extract the physical matrix elements.
In conclusion we have found that for a large class of processes (those without final state interactions) it is possible to neglect the twist of the determinant avoiding the need to generate new gauge configurations for each twist. This is not true however, for all processes. In particular, for K → ππ matrix elements, with the two pions in an I = 0 state, if G-parity boundary conditions are used they must be implemented for both the valence and sea quarks.

Conclusions
In this letter we have used χPT to study the finite-volume corrections with twisted boundary conditions. For quantities without final-state interactions, such as meson masses, decay constants or semileptonic and other form-factors, we confirm that these corrections remain exponentially small in the volume. This remains true with partially twisted boundary conditions for which only the valence quarks are twisted, thus eliminating the need to generate a new set of gluon configurations for each choice of twisting angle and makes the technique much more useful.
We have also demonstrated that twisted boundary conditions cannot be applied in general to processes with final-state interactions, such as K → ππ decays. This is disappointing since twisted boundary conditions would have been particularly useful for lattice studies of these decays, extending very significantly the kinematic range accessible in a simulation. In spite of this particular disappointment, we look forward to the implementation of twisted boundary conditions to the wide range of processes for which they are applicable and phenomenologically useful.

A Finite-Volume Corrections in Chiral Perturbation Theory
In this appendix we derive the finite-volume corrections with twisted boundary conditions at one-loop order in χPT. The generic expression for tadpole diagrams in finite volume is given by the left-hand side of 1 The first term on the right-hand side of eq. (30) is the corresponding infinite-volume integral and ξ θ s (L, M ) contains the finite-volume corrections. We now generalise the procedure of ref. [16] to twisted boundary conditions and demonstrate that these corrections are exponentially small in the volume.
If cos(θ i ) = 0, then the leading finite-volume corrections are given by the m = ±2 terms in eq. (36) and hence decrease with a larger exponent. For the generic case in which cos(θ i ) = 0 for i = 1, 2, 3 the behaviour of the finite-volume corrections as L → ∞ is given by where ξ 0 s are the finite-volume corrections with (θ 1 , θ 2 , θ 3 ) = 0. Eq. (38) demonstrates that finite-volume corrections are exponentially small with twisted boundary conditions.
The second diagram which appears at one loop level contains two propagators and in finite volume gives a contribution proportional to: where the numerator N is a function of momenta and masses and q is the injected momentum. For illustration we consider the simplest case for which N = 1 (terms in N containing the loop momentum can be reduced to the tadpole integrals of the form in eq. (30) ) and q = 0 (the finite-volume effects in integrals with non-zero q can readily be obtained from the expressions below by the substitution A 2 → A 2 + q ).
After introducing the Feynman parameter x and performing the k 4 integration eq. (40) reduces to: where where in the last line we have made the replacement q 4 → iE and E is the physical (Minkowski) injected energy. Eq. (41) has the same form as the left-hand side of eq. (30), so we can proceed just as for the tadpole integral to obtain the expression for the corresponding ξ θ -function: which is exponentially small in the volume as long as M 2 (x) > 0, i.e. as long as no branch cuts appear.
In partially quenched chiral perturbation theory there are also contributions with double poles in which one or both propagators in eq. (40) are squared. These can be written in terms of derivatives of (40) w.r.t. the masses and the finite-volume corrections therefore remain exponentially small.

B Masses and Decay Constants
In this appendix we present the full finite-volume corrections for meson masses and decay constants at NLO in χPT with twisted boundary conditions (using the notation of eq. (21) ):