Scrutinizing the η-η ′ mixing , masses and pseudoscalar decay constants in the framework of U ( 3 ) chiral effective field theory

We study the η-η′ mixing up to next-to-next-to-leading-order in U(3) chiral perturbation theory in the light of recent lattice simulations and phenomenological inputs. A general treatment for the η-η′ mixing at higher orders, with the higher-derivative, kinematic and mass mixing terms, is addressed. The connections between the four mixing parameters in the two-mixing-angle scheme and the low energy constants in the U(3) chiral effective theory are provided both for the singlet-octet and the quark-flavor bases. The axial-vector decay constants of pion and kaon are studied in the same order and confronted with the lattice simulation data as well. The quark-mass dependences of mη, mη′ and mK are found to be well described at next-to-leading order. Nonetheless, in order to simultaneously describe the lattice data and phenomenological determinations for the properties of light pseudoscalars π,K, η and η′, the next-to-next-to-leading order study is essential. Furthermore, the lattice and phenomenological inputs are well reproduced for reasonable values of low the energy constants, compatible with previous bibliography.

dependencies in terms of m π and study the observables as functions of m π . The η and η ′ lattice simulations have not been thoroughly analyzed in the chiral framework yet and it is the central goal of the present work. However, the numerical uncertainties resulting from our analyses in this work must be taken with a grain of salt as correlations between the different lattice data points and other systematic errors are not considered here.
In addition to lattice QCD, there are also phenomenological studies of the η and η ′ mixing, which has been extensively investigated in radiative decays of light-flavor vector resonances ρ, ω, φ and J/ψ → V P, P γ processes [14,15,[24][25][26][27][28][29][30]. In these works, the modern two-mixing-angle scheme for the η and η ′ mesons, which was first advocated in Refs. [4,5], was employed to fit various experimental data. The common methodology in these works is that the two-mixing-angle pattern for the η and η ′ is simply adopted to perform the phenomenological discussion and the mixing parameters are then directly determined from data. This is a bottom-up approach to address the η-η ′ mixing problem and it is quite useful for the phenomenological analysis. Contrary to the bottom-up method, it is also very interesting to study the η-η ′ mixing from a top-down approach in which one first constructs the relevant χPT Lagrangian and then calculates the η-η ′ mixing pattern and parameters in terms of the LECs. In this case, one can predict the η-η ′ mixing parameters once the values of the unknown LECs are given. The present work belongs to the latter category of top-down approaches.
Though the singlet η 0 meson, which is the main component of the physical η ′ state, is not a pNGB due to the strong U (1) A anomaly, it can be formally introduced into χPT from the large-N C point of view. The argument is that the quark loop induced U (1) A anomaly, which is responsible for the large mass of the singlet η 0 , is 1/N C suppressed and hence the η 0 becomes the ninth pNGB in the large N C limit [31]. Based on this argument, the leading-order (LO) effective Lagrangian for U (3) χPT, which simultaneously includes the pNGB octet π, K, η 8 and the singlet η 0 as dynamical fields, was formulated in Ref. [3]. Later on, a full O(p 4 ) U (3) chiral Lagrangian was constructed in Ref. [7] and the discussion on the O(p 6 ) unitary group chiral Lagrangian has been very recently completed in Ref. [32]. Subtle problems about the choice of suitable variables for the higher order U (3) χPT Lagrangian in the large N C framework were analyzed in Ref. [6].
The standard power counting employed in SU (2) and SU (3) χPT in powers of the external momenta and quark masses [1], is not valid any more in U (3) χPT, due to the appearance of the large η 0 mass. However, since the singlet η 0 mass squared behaves like 1/N C in large N C limit, the η 0 mass can be harmonized with the other two expansion parameters if one assigns the same counting to 1/N C , the squared momenta p 2 and the light quark masses m q . As a result of this, in order to have a systematic power counting, the combined expansions on momentum, light quark masses and 1/N C are mandatory in U (3) χPT [6,7]. We will work in this combined expansion in our study and denote it as δ expansion throughout the paper, where O(δ) ∼ O(p 2 ) ∼ O(m q ) ∼ O(1/N C ). This counting rule is different from the one proposed in Ref. [33], where the η 0 mass is counted as O(1) and the infrared regularization method is employed to handle the chiral loops.
Some recent works in Refs. [8,[34][35][36][37] have addressed the η-η ′ mixing in the chiral framework up to next-to-leading order (NLO). As an improvement, we will perform the systematic study of the η-η ′ mixing in the δ-expansion scheme up to next-to-next-to-leading order (NNLO) and take into account the very recent lattice simulation data, which are not considered in the previous works [8,[34][35][36][37]. In addition, we also simultaneously analyze the m π dependences of other physical quantities from lattice simulations, such as the axial π, K decay constants and the mass ratio of the strange and up/down quarks, in order to further constrain the χPT LECs.
This article is organized as follows. In Sect. II, we introduce the theoretical framework and calculate the relevant physical quantities. In Sect. III, the phenomenological discussions will be presented. Conclusions will be given in Sect. IV. Further details about the calculations up to NNLO are relegated to App. A.

A. Relevant chiral Lagrangian
At leading order in the δ expansion, i.e. O(δ 0 ), the U (3) χPT Lagrangian consists of three operators where the chiral building blocks are defined as [1,[6][7][8] appearing in Eqs. (1) and (2) corresponds to the pNGB axial decay constant in the large N C and chiral limits. The light quark masses are introduced by setting (s + ip) =diag{m,m, m s }, beingm the averaged up and down quark masses and m s that of the strange quark.
Notice the structure of the LO Lagrangian in Eq. (1): the first operator is of O(N C , p 2 ) type, the second one corresponds to the type of O(N C , m q ) and the last one stems from the QCD U (1) A anomaly and is of O(N 0 C , p 0 ) type, where U is counted as O(1),F 2 ∼ O(N C ) and M 2 0 ∼ O(N −1 C ) in the classification O(N j C , p k , m ℓ q ) of the EFT Lagrangian operators in Eq. (1). In the following, we will denote the chiral expansions in powers of squared momenta p 2 and quark masses m q simply as a generic expansion in p 2 .
The NLO U (3) chiral Lagrangian, i.e., O(δ), contains O(N C , p 4 ) and O(N 0 C , p 2 ) operators. The relevant ones in our work read [6] L (δ) = L 5 u µ u µ χ + + L 8 2 . Their explicit forms read [7,38] where the first line corresponds to the O(N −1 C , p 2 ) type, the second line is of the O(N 0 C , p 4 ) type and the last two lines are of the O(N C , p 6 ) type. The LECs carry the scalings v  (4) and (5) that are pertinent to our present study, not aiming at giving the complete sets of operators. The conventions to label the LO, NLO and NNLO operators in Eqs. (1), (4) and (5) follow closely the notations in Refs. [6,7,38]. Unless it is explicitly stated, the LECs will correspond to U (3) χPT and must not be confused with those in SU (3) χPT. The matching between these two EFTs can be found in Ref. [6]. The terms L j are denoted as β j in Ref. [7,8].
Comparing the U (3) and SU (3) theories one can observe that some terms have been reshuffled in the δ expansion of the U (3) Lagrangian. For example, the L i=4,5,6,7,8 terms are NLO in SU (3) χPT, but they are now split into NLO and NNLO in the δ expansion (see Eqs. (4) and (5)). We have several additional new operators, namely the last one in Eq. (1), the Λ i=1,2 in Eq. (4) and the v (2) 2 , L 18 , L 25 terms in Eq. (5), that are absent in the SU (3) χPT case. Finally, the chiral loops start contributing at NNLO in the δ expansion, while they appear at NLO in the conventional SU (3) case.

B.
The η-η ′ mixing at NNLO in δ expansion Next we calculate the η-η ′ mixing order by order in the δ expansion. In literature, there are two bases to address the η-η ′ mixing, namely the singlet-octet basis with η 0 and η 8 , and the quark-flavor basis with η q and η s . The relations between fields in these two bases are In the large-N C limit where the U (1) A anomaly is absent, η q and η s are the mass eigenstates and they are generated by the axial-vector currents with the quark flavors qq = (uū + dd)/ √ 2 and ss, respectively. The two bases are related to each other through an orthogonal transformation and provide an equivalent description for the η-η ′ mixing.
As noticed in Refs. [37,39], when doing the loop calculations with η and η ′ , it is rather cumbersome to work with the η 0 and η 8 states. The reason is that at leading order the Lagrangian in Eq. (1) gives the mixing between η 0 and η 8 , and the mixing strength is proportional to m 2 K − m 2 π , which in the δ expansion is formally counted as the same order as the diagonal terms in the mass matrix for η 0 and η 8 . As a result, the insertion of the η 0 -η 8 mixing in the chiral loops will not increase the δ order of the loop diagrams. This makes the loop calculation technically much more complicated, as one needs to consider the arbitrary insertions of the η 0 -η 8 mixing in the chiral loop diagrams. Nevertheless, Refs. [37,39] provide a simple recipe to handle this problem by expressing the Lagrangian in terms of the η and η ′ states which result from the diagonalization of η 0 and η 8 at leading order in δ. The main difference is that the mixing between η and η ′ is now at least a NLO effect in δ, while the η 0 -η 8 mixing was appearing at LO. The relation between the LO mass eigenstatesη andη ′ and the singlet-octet basis is given by the mixing angle θ: with c θ = cos θ and s θ = sin θ. The LO mixing angle θ and masses of η and η ′ are given by the leading order Lagrangian L (δ 0 ) in Eq. (1) (see e.g. Ref. [37]): Here m K and m π denote the LO kaon and pion masses, respectively. When higher order corrections are taken into account, the LO diagonalized η and η ′ will get mixed again. Up to the NNLO, a general parametrization of the bilinear terms involving the η and η ′ states can be written as where the δ ′ i s contain the NLO and NNLO corrections. Here these operators must be understood as the terms of the effective action that provide the pseudoscalar meson self-energies. The higher-derivative terms δ j=1,2,3 in the first line of Eq. (11) are exclusively contributed by the O(p 6 ) operator C 12 in Eq. (5), which belongs to the NNLO Lagrangian. The remaining δ ′ i s receive contributions from the NLO operators in Eq. (4), the NNLO ones in Eq. (5) and the one-loop diagrams, which contribute at NNLO. Their explicit expressions can be found in App. A.
At leading order, there is only the mass mixing term from Eq. (1) whereas at NLO and NNLO one has to deal in addition with the kinematic mixing terms in Eq. (11), apart from the mass mixing. The physical states of η and η ′ can be obtained from the perturbative-expansion (δ-expansion) in three steps: as a first step, we eliminate the higher-derivative terms through the field redefinitions of η and η ′ ; then we transform and rescale the fields resulting from the first step in order to write the kinematic terms in the canonical form; after the preceding two steps, there is only the mass mixing term left, which is straightforward to handle.
In the first step, we make the following field redefinitions for the η and η ′ states with the d'Alembert operator ≡ ∂ µ ∂ µ . After some algebra manipulations, it is straightforward to obtain so that the three higher-derivative terms in Eq. (11) will be eliminated. Notice that the α 1,2,3 are NNLO, i.e., O(δ 2 ). Substituting the field redefinitions from Eq. (12) into the general mixing structure in Eq. (11) and keeping the terms up to NNLO, the resulting bilinear Lagrangian reads In the second step, we need to eliminate the kinematic mixing term in Eq. (14), and then to rescale the fields to have them in the canonical forms. This can be done perturbatively. In the final step, we take care of the mass mixing term. The last two steps can be achieved through the following field transformations with η, η ′ the physical states and where δ η,NLO , δ η ′ ,NLO , δ k,NLO stand for the NLO parts of the three quantities respectively. We point out that δ η , δ η ′ , δ k receive both NLO and NNLO contributions, while δ 1 , δ 2 , δ 3 are only contributed by the NNLO effect, which is the C 12 operator in Eq. (5). Comparing with the NLO results in Eq. (15) from our previous paper [37], we have generalized the expression to the NNLO case in the present Eq. (15). Another way to treat the mixing of pseudoscalar mesons in χPT was also previously studied in Ref. [40] and applied to the π 0 -η case up to the two-loop level.
In the practical calculation, it is more often to use the inverse of the relations in Eq. (15), where the perturbative expansion leads to with The θ δ appearing in Eqs. (15) and (17) is determined through with where δ i,NLO stand for the NLO parts of δ i . In the phenomenological discussions, the popular two-mixing-angle parametrization in the singlet-octet basis [4,5] takes the form Combining Eqs. (7) and (15), it is straightforward to derive the relations between the four parameters in the two-mixing-angle scheme in Eq. (21) and the χPT LECs: where the χPT LECs are implicitly included in θ, θ δ , δ A , δ B and δ C . Since θ δ , δ A , δ B , δ C ∼ O(δ) or O(δ 2 ), at LO one has F 8 = F 0 = F and one mixing-angle θ 8 = θ 0 = θ. The relations between the physical η, η ′ states and the quark-flavor basis is commonly parametrized as Combining Eqs. (6), (7) and (15), it is straightforward to obtain the parameters in Eq. (23): where at LO in the δ-expansion one has F q = F s = F and φ q = φ s = θ id −θ, with the ideal mixing θ id = − arcsin 2/3.

C. Insights into previous studies of the η-η ′ mixing
In the previous subsection we have performed the full computation of the mixing up to NNLO in the δ expansion. It is interesting to make a brief summary of the assumptions made in previous works, where plenty of mixing formalisms have been proposed to address the η-η ′ system [8, 24, 34-37, 41, 42]. In Ref. [41], only the lowest order in the quark masses and 1/N C , i.e. the LO contributions in the δ expansion, were taken into account. Even though it provided a reasonable first approximation, it failed to give an accurate description of the experimentally observed mass ratio m 2 η /m 2 η ′ . The O(p 2 ) contributions were studied up to NLO in 1/N C in Ref. [42] (including the terms in Eq. (1) and Λ 1 and Λ 2 in Eq. (4)), perfectly explaining the experimental value of m 2 η /m 2 η ′ . However, it turned out to be inadequate to give a proper value for the η-η ′ mixing angle. On the other hand, the authors in Refs. [34][35][36] went up to NLO in the p 2 expansion but keeping just the LO in 1/N C (including the terms in Eq. (1) and L 5 and L 8 in Eq. (4)). Both the η-η ′ mixing angle and the ratio F K /F π were qualitatively reproduced in this case. The full set of NLO contributions in the δ-expansion (i.e., the effects up to NLO both in 1/N C and p 2 ) was analyzed in Ref. [8], together with the mixing angle and the π, K, η and η ′ axial-vector decay constants. In Ref. [37], the contributions from the tree-level resonance exchanges and partial NNLO effects, e.g. the loop diagrams, were considered for the masses of η and η ′ . In this work, we generalize the discussions up to the full NNLO study in the δ-expansion and confront our theoretical expressions with the very recent lattice simulation data and the phenomenological inputs from the two-mixing-angle scheme.
Reference [24] introduced a quark-model inspired approach to the η-η ′ mixing, which is commonly referred as the FKS formalism and used in many phenomenological analyses [43]. The essence of the FKS formalism is the assumption that the axial decay constants in the quark-flavor basis takes the same mixing pattern as the states where the decay constants are defined as the matrix elements of the axial currents From another point of view, the pattern of Eq. (25) employed in the FKS formalism relies on the assumption that there is no mixing between the decay constants of the flavor states η q and η s . In the χPT framework, the physical masses and decay constants can be obtained from the bilinear parts of nonet fields in the effective action with the correlation function of two axial currents. Since the correlation function is the second derivative with respect to the axial-vector external source a µ , and a µ always appears in the Lagrangian together with the partial derivative ∂ µ as shown in Eq. (2), the absence of the mixing for the η q and η s decay constants in Eq. (25) implies that there are no kinematic mixing terms for the quark-flavor states η q and η s in the FKS formalism. In fact, the assumption in Ref. [34] is in accord with the FKS formalism. This can be simply demonstrated by expanding the chiral operators considered in Refs. [34,35], i.e. those in Eq. (1) and L 5 , L 8 in Eq. (4), up to quadratic terms in η q and η s . 1 No kinematic mixing terms for the η q and η s fields result from these chiral operators. This also confirms the finding in Ref. [36] that only when the NLO of 1/N C operator is excluded the FKS formalism is recovered with their chiral Lagrangian calculations.
Since general terms up to NNLO in δ expansion are kept in our discussion, unlike in the previous works [8, 34-37, 41, 42] where different assumptions, such as the preference of the higher order p 2 and 1/N C effects, are made, it is important and interesting for us to justify these assumptions in later discussions.

D. Masses and decay constants of pion and kaon up to NNLO in δ expansion
The NLO expression of the pion decay constant in the δ expansion reads or, up to the precision considered, one can also use the physical F π in the expression inside brackets, The differences between Eqs. (27) and (28) are NNLO effects. We mention that at a given order there is always ambiguity in choosing the renormalized quantities in the higher order expressions. In contrast, there is formally no ambiguity in the expressions in terms of the quantity F , which is the pNGB axial decay constant in the chiral and large N C limits. For example, if we limit our analysis up to NLO, formally, it is equally good to use F π or F K in the denominators of the NLO part in Eq. (28), since the difference is beyond the NLO precision. A typical solution in the chiral study is to express the quantities, such as m π , F π , m K , F K , in terms of the renormalized F π in the higher order corrections, as done in the two-loop calculations in SU (3) χPT [44]. We follow this rule throughout the current work to estimate the uncertainty due to the truncation of the δ expansion when one works at a given order in perturbation theory. We mention that the notation of m 2 π in the above equations stands for the renormalized pion mass squared and the leading order mass squared is denoted by m 2 π . Notice the LO pion mass squared m 2 π is the one that is linear in the quark masses. The expressions relating m 2 π and m 2 π will be discussed below. Similarly up to NNLO, we can either use F or F π in the NLO and NNLO expressions for other quantities such as F K and the δ i 's in Eq. (11). In the NNLO expressions, the difference between using F or F π in the denominators is a next-to-next-to-next-to-leading order effect (N 3 LO). Since in this work we study lattice simulation data up to pion mass of 500 MeV, the convergence of the chiral series is expected to be much slower than that in the physical case with m π = 135 MeV. Therefore it is a priori not trivial to judge whether the two approaches-using 1/F 2 and 1/F 2 πare numerically equivalent or the lattice data prefer one of them. Indeed in Ref. [45], it is already noticed that to use F or F π could cause some noticeable effects. We will use the difference between both approaches as an estimate of the truncation error at a given order in δ.
We take the pion decay constant as an example to illustrate the differences of using F and F π in the higher order expressions. Using F in the higher order corrections, its expression reads The one-point loop function A 0 (m 2 ) is calculated in dimensional regularization within the M S − 1 scheme [1] and it reads corresponds to our Λ 2 operator in Eq. (4). The Λ term, though introduced from the beginning in these references, is dropped in their later discussions, since it is 1/N C suppressed.
with the renormalization scale µ fixed at 770 MeV throughout. Using Eq. (27) to replace F by F π in the NLO and NNLO corrections, the resulting form is In the δ expansion, the expressions for a physical quantity with F or F π in the higher order chiral corrections differ only for the L 5 L j=5,8 and L 5 Λ j=1,2 terms, since the differences by replacing F by F π are originated from the NLO expressions of F π in Eq. (28) and we only retain terms up to NNLO in this work. It is clear that the difference between Eqs. (29) and (31) is the L 2 5 term. Notice that in the δ expansion scheme, the terms like L 5 L 4 are N 3 LO and will be dropped throughout the article.
The corresponding expression for the kaon decay constant when one uses F to express the NLO and NNLO corrections reads On the other hand, expressing the NLO and NNLO contributions in terms of F π yields The expanded expression for the ratio of F K /F π in terms of F up to NNLO in δ expansion, takes the form When expressing the previous result in terms of F π , it reads which differs from Eq. (34) in the L 2 5 term. The pion squared mass up to NNLO is given by with When expressing the renormalized m π in terms of F π , the only differences are the L 5 L 8 and L 2 5 terms in Eq. (39) and the other parts are the same as in Eq. (36) with the explicit replacement of F by F π in Eq. (27). Therefore we only give the different parts for simplicity when expressing in terms of F π and they read The mass squared for kaon up to NNLO is provided by When expressing Eq. (41) in terms of F π , the differences are the L 5 L 8 and L 2 5 terms in Eq. (44) and the new expressions are Notice that the masses of pion and kaon appearing in NLO and NNLO parts in the above equations correspond to the renormalized quantities, instead of their LO expressions. In addition, this gives the quark mass ratio relation m s /m = 2m 2 K /m 2 π − 1. When performing the chiral extrapolation of the lattice data, instead of the renormalized m 2 K as in the previous equations, it is convenient to use the LO kaon mass squared in the higher order corrections. In this way, we do not need to iteratively solve Eq. (41) in order to give the value of m K for a given m π . The result in terms of m K in the NLO and NNLO expressions becomes with When expressing Eq. (46) in terms of F π , the differences are the L 5 L 8 and L 2 5 terms in Eq. (48) and the new expressions are When confronting with the lattice data, we only consider the simulated points with the physical strange-quark mass, i.e. the lattice ensembles that when extrapolating to the physical pion masses lead simultaneously to physical kaon masses. In this case, we can express the LO kaon mass squared as where m 2,Phy K and m 2,Phy π can be obtained through Eqs. (36) and (41) by substituting the physical masses of π, K, η, η ′ in the NLO and NNLO expressions. For m 2 π , which varies in the lattice simulation, we can extract its value by using Eq. (36). In this case, m π in Eq. (36) takes the value from lattice simulation, and m K , m η , m η ′ , which only appear in the NNLO part, can be approximated by their LO expressions.
In the above discussions, we have distinguished the situations of using 1/F 2 and 1/F 2 π in the higher order corrections for various observables. Similarly, we can also generalize the discussions by replacing the renormalized masses (m π and m K ) with the LO ones (m π and m K ) in the higher order corrections. We take the observables F π and F K as examples to illustrate the differences. The renormalized m π and m K have been used in Eqs. (29) and (32) for F π and F K with 1/F 2 in the higher order terms, respectively. After replacing m π and m K in Eqs. (29) and (32) with their expressions in terms of the LO masses m π and m K through Eqs. (36) and (41) respectively, the corresponding expressions are found to be As in the discussion of 1/F 2 versus 1/F 2 π up to the NNLO precision, the expressions for a specific observable by using the renormalized masses m π , m K and the LO m π , m K only differ in the terms like L i L j , being L i and L j the NLO LECs in Eq. (4). This can be clearly seen when comparing Eqs. (29) and (51). E.g. the differences caused by using the renormalized masses and the LO ones are the L 2 5 and L 5 L 8 terms, apart from the explicit replacement of m π and m K by m π and m K respectively. Similar rules are also applied to Eqs. (32) and (52).
To replace m π , m K by m π and m K in the NLO and NNLO corrections in Eq. (36), the only changes happen for the L i L j terms and the corresponding new expressions read In principle, we should also present the results expressed with the LO masses m π , m K and the renormalized decay constants F π , F K , which can be straightforwardly obtained by substituting the relations in Eqs. (36) and (41) into the corresponding observables. We consider the expressions given in terms of the renormalized masses and 1/F 2 as our preferred ones in this work. The reason to choose the renormalized masses is for practical purpose, since in lattice simulations the different observables are typically given as functions of the renormalized m 2 π . Also most of the chiral studies choose to express the quantities with the renormalized masses in the higher order corrections, such as in Refs. [44,45]. Following this rule we consider the results with the renormalized masses and 1/F 2 π as an estimate of systematic errors due to the truncation of the δ expansion when one works at a given order in perturbation theory. While for the case with the LO masses, we shall also comment the results in the following numerical discussions.

III. PHENOMENOLOGICAL DISCUSSIONS
The big challenge in the present general discussions on the η-η ′ mixing is the determination of the unknown LECs in Eqs. (1), (4) and (5). The recent lattice simulations on the light pseudoscalar mesons provide us valuable sources to constrain these free parameters. The considered lattice simulations include the m π dependences of the masses of η, η ′ [9-13] and kaon [16,17], and the π, K decay constants [16,17] and their ratios [18]. Moreover, relevant phenomenological results and experimental data will be also included to constrain the LECs.
Since we do not consider the isospin violating effects, we will take the values for the physical pion and kaon masses in the isospin limit from Ref. [46], where the corrections from the electromagnetic contributions are removed, These values will be used in later chiral extrapolations, while for the physical masses of η and η ′ and the decay constants of pion and kaon, we will take their world-average values from Ref. [47].
In order to show the results step by step, we present the discussions in the following sections split in three parts: we consider fits performed at leading order, next-to-leading order and next-to-next-to-leading order.

A. Leading-order analyses
At leading order, the η-η ′ mixing is described by one free parameter, namely the singlet η 0 mass M 0 in Eq. (1) and the explicit expressions for the masses and mixing angle are given in Eqs. (8), (9) and (10). At this order, the π, K decay constants are degenerate and given by their chiral and large N C limits, i.e. F π = F K = F . Therefore we shall not take the lattice simulations of the decay constants into account for the LO discussion as they clearly show the need of higher order corrections for a suitable description. Also at leading order, F will not enter the masses and mixing angle, as shown in Eqs. (8), (9) and (10). As a result of this, we do not need to distinguish the two situations with F or F π discussed previously in the expressions of different observables. Apart from the lattice simulation data, we also fit the physical values of the η and η ′ masses. Nonetheless, fitting the physical masses with the experimental precision at the level of several hundred-thousandth is too ambitious. Since the ultimate goal of the present work is the NNLO study, the ballpark estimate of our theoretical uncertainty, starting from the N 3 LO part, should be around 3%. This value is obtained from the general rule that each higher order correction in δ expansion, either the SU (3)-flavor breaking or the 1/N C effect, is around 30%. In fact, the estimated three-percent uncertainty is also similar to the typical error bars reported in many lattice simulations, in the range from 3 − 10% [9][10][11][12][13]. Consistently, we assign a 1% uncertainty to the physical values of m η and m η ′ in the fits.
The value of the singlet mass M 0 from the LO fit is The physical masses for the η, η ′ and their LO mixing angle θ from the fit are found to be The resulting plots can be seen in Fig. 1. We verify that if the physical masses are excluded in the fit, M 0 = 813 ± 11 MeV results. If we only include the physical masses and exclude the lattice simulation data in the fit, M 0 = 859 ± 11 MeV is obtained. These determinations of M 0 lie within the broad range summarized in Ref. [43] and are quite close with the commonly used values of M 0 = 850 MeV [43]. Taking into account the large uncertainties of the lattice simulation data, specially for m η ′ , and the concise formalism of the LO mixing, it is impressive that the lattice simulation data can already be qualitatively described with the LO analysis, as shown in Fig. 1. This also indicates that the higher order mixing effects can only give moderate corrections to the masses of η and η ′ .
Nevertheless, in order to describe the lattice data more accurately, specially the η masses, the chiral corrections beyond the leading order are needed. For the physical masses, it has also been shown that the LO description fails to explain the mass ratio of η and η ′ accurately enough [41]. Therefore it is essential to generalize the discussions to NLO and NNLO in order to achieve a precise description both for lattice simulations and physical data.  [12,13] (ETMC), [11] (UKQCD), [10] (RBC/UKQCD), [9] (HSC), where we only take into account the simulation points with mπ < 500 MeV. The shade area surrounding each curve stands for the statistical uncertainty from the fit.

B.
Next-to-leading order analyses At next-to-leading order, in addition to the parameter M 0 at leading order, there are five additional free parameters: the decay constant F at chiral and large N C limits, and the four NLO LECs L 5 , L 8 , Λ 1 and Λ 2 in Eq. (4). At this order, as well as at next-to-next-to-leading order, one can rewrite the chiral expansion of the observables in various equivalent ways up to the perturbative order in δ under consideration. In the following discussion, we will perform two types of fits: one using F in the theoretical NLO and NNLO expressions and the other employing F π , as discussed in Sect. II D. Since the differences of the theoretical expressions used in the two types of fits are beyond the considered precision, the variances of the outputs from the two fits can be considered as systematic errors from the theoretical models by neglecting higher order contributions. In the following, we will explicitly present the fit results by using F in the theoretical expressions, which is the most straightforward option, as discussed in Sect. II D. The outputs of the fits with the theoretical formulas expressed in terms of F π will be used to estimate the systematic errors: the difference between the central values of the two types of fits will be used to estimate the truncation uncertainty due to working just up to a given order in the δ-expansion, providing the second error for each quantity in the following tables.
In Refs. [34][35][36], it is argued that at each chiral order, the leading N C effects are dominant, or in other words that the Λ 1 and Λ 2 terms are assumed to be much less irrelevant than the L 5 and L 8 terms in the NLO δ expansion. This assumption has been more or less confirmed when focusing on the masses of η, η ′ and the LO mixing angle at the physical points [34][35][36]. In Ref. [37], the local higher order LECs were estimated by the tree-level resonance exchanges and it was found that with those LECs Λ 2 seems to be more important than Λ 1 when focusing on the physical masses for η and η ′ . It is interesting to check how these assumptions work when including the lattice simulations and the phenomenological results of the two-mixing-angle parameters, which are not considered in Refs. [34][35][36][37]. Different sets of fits to the lattice data and phenomenological inputs from the two-mixing-angle scheme are performed either by fixing Λ i=1,2 to zero or releasing their values, in order to reexamine the assumptions. Interestingly we do not find qualitative changes between the fits with fixed Λ i=1,2 = 0 and the ones with free values for these parameters. This tells us that indeed the Λ 1 and Λ 2 terms do not significantly improve the fit results, even after taking into account the lattice simulations. Nevertheless, we find that these two terms are quite important to reproduce the phenomenological mixing angles θ 0 and θ 8 in the fits where M 0 is fixed at its LO value. If M 0 is released in the fits we find that including Λ 1 and Λ 2 improves the descriptions of m η ′ from lattice simulations. Therefore, we will not further discuss fits with Λ 1 and Λ 2 set to zero in the following. Instead, we focus on the results given in Table I with all the four NLO LECs in the fits, namely L 5 , L 8 , Λ 1 and Λ 2 in Eq. (4).
For the parameter M 0 , we take two strategies to estimate its value in NLO analysis. In one of them we fix M 0 = 835.7 MeV from its LO determination (NLOFit-A) and in the other case we free its value for the NLO fit (NLOFit-B). These two NLO fits are given in Table I. The first error bar for each fitted parameter corresponds to the statistical one from the fits and the second error bar is estimated from the variation of the fits between those using F and F π in the NLO (and later also NNLO) theoretical expressions. From the two fits shown in Table I, one can see that releasing M 0 in the fits barely changes the fit quality with respect to the cases when its value is fixed, although there are slight variations in the determinations of M 0 and Λ 2 .
Concerning the results of the LECs in Table I, the resulting values for F from the two fits are quite compatible and close to the physical pion decay constant. For Λ 1 and Λ 2 , their values are poorly known in literature and it is helpful to compare our values with the following estimate for their ranges: we take the LO determination M 0 = 835.7 MeV, and we then separately include the Λ 1 and Λ 2 terms in the η-η ′ mixing and vary their values to obtain new results for m η and m η ′ with the physical m π . Since the Λ 1 and Λ 2 terms are NLO 1/N C effects, it is reasonable to assume that their corrections to m 2 η or m 2 η ′ should be at most around 30% of the LO results. In this way we can set up conservative and rough estimates for the ranges of Λ 1 and Λ 2 , which are found to be The resulting magnitudes of Λ 1 in our fits are tiny and consistent with zero, as shown in Table I. For the parameter Λ 2 , our determinations lie within the ranges estimated in Eq. (57). Its value, specially the one from NLOFit-A, is close to the one used in Ref. [48], where the mixing was discussed at next-to-leading order. However the determinations for Λ 2 in Table I become much more precise than those given in Refs. [37,39], where the lattice simulations for m η and m η ′ are not included, indicating the usefulness of incorporating the lattice data in the U (3) χPT study. Our determinations of L 5 and L 8 are in good agreement with the leading N C predictions from resonance chiral theory [49], the SU (3) one-loop results in Ref. [1] and the one-loop resonance chiral theory determination for L 8 [50]. But the values here are clearly larger than those from the recent two-loop determinations [20,21], the results from Kπ scattering in the scalar channels [51], and the one-loop resonance chiral theory estimates for L 5 [23]. The discrepancies of L 5 and L 8 , comparing with the recent two-loop determinations [20,21], can be eliminated once the O(p 6 ) LECs are taken into account, as we will show in the NNLO discussion.
The values of the parameters in the two-mixing-angle scheme and the mass ratio of strange and up/down quarks resulting from the fits are given in Table II. Similarly, the first error bar for each quantity is the statistical error and the second one corresponds to the systematic error, which is obtained in the same way as the one in Table I. Notice that these inputs have already been satisfactorily reproduced in NLO analyses.
The other quantities in the fits are presented in Figs. 2, 3, 4 and 5, together with the lattice simulation data and the experimental inputs. We find that the final outputs from NLOFit-A and NLOFit-B are quite similar, so only the plots from NLOFit-B are given explicitly. The shaded area surrounding each curve corresponds to the statistical error band for each quantity. In Fig. 2, we show the resulting figures from NLOFit-B for the masses of η and η ′ . In Figs. 3, 4 and 5, we show the corresponding plots for m 2 K , F π,K and F K /F π as functions of m 2 π , respectively.    (24). The phenomenological values for the mixing parameters are taken from Ref. [15] and we triple the error bands here in order to make a conservative estimate. The input of ms/ m is taken from the FLAG working group in Ref. [46] and we assign the 10% error bar as done in Ref. [20]. For the error bars of each quantity, the first one corresponds to the statistic error and the second one is for the systematic error, which are explained in detail in the text.

C. NLO fits focusing on the masses
In this section, we present another kind of NLO fits by focusing on the masses of η, η ′ , K and excluding the decay constants F π , F K and their ratio. This kind of discussion is well motivated, since it is known that the NNLO corrections in δ counting, such as the LEC L 4 , are important to simultaneously describe F π and F K [20][21][22]. But this LEC is absent in NLO study. We have also provided another independent confirmation on this finding in Fig. 4, where one can see that the decay constants of pion and kaon are poorly reproduced at next-to-leading order in δ expansion. When only focusing on the η, η ′ and K masses and the ratio m s /m at next-to-leading order the parameter F can not be resolved, because it always appears in the form L 5 /F 2 or L 8 /F 2 . We will fix its value to F = 90 MeV, close to the values given in Table I. For the mixing parameters we consider the mixing angles of θ 0 and θ 8 , but exclude the constants F 0 and F 8 . This is because F 0 and F 8 are dependent on the parameter F and should be determined together with F π and F K . For simplicity in later discussion, we call the fits performed in this section as the mass-focusing type throughout.
As in the previous section, we present the fits with F in the denominators of the theoretical expressions (e.g. Eq. (27)) and use the fits with F π to estimate the systematic errors. For each case, we perform the fits either by fixing M 0 at its LO determination (NLOFit-C) or by freeing its value (NLOFit-D). The fitted parameters are given in Table III and the m s /m ratio and mixing angles are given in Table IV. The resulting figures from NLOFit-C and NLOFit-D are quite similar and we explicitly show one set of them, e.g. NLOFit-D in Figs. 2 and 3 for the η ( ′ ) and kaon masses, respectively. A significant difference between the results in Table I and the mass-focusing fits in Table III is that much larger statistical error bars are obtained in the latter case, especially for the LECs L 5 and L 8 , as they are constrained by fewer data. Likewise, there are large systematic errors for the values of L 5 and L 8 in Table III, indicating a larger truncation uncertainty due to higher orders. We do not see a significant improvement when freeing the value of M 0 in the fits.     Table II for the phenomenological inputs. The first error for each quantity corresponds to the statistical one and the second error denotes the systematic uncertainty. See the text for details.

D.
Next-to-next-to-leading order analyses From the NLO discussions in the previous two sections, we observe that the phenomenological results and the lattice simulations on η and η ′ states can be reasonably reproduced. This is an important improvement comparing with the LO study, since at this order we only have the conventional one-mixing-angle formalism. The two-mixingangle formalism only shows up beyond LO. However, observing m K , F π , F K and their ratio in Figs. 3, 4 and 5, it is clear that the NLO analysis is still inadequate. We need to include higher order contributions beyond NLO in order to further improve the descriptions. Moreover, the chiral logarithms predicted by χPT at one loop start at NNLO in the δ expansion. Due to their importance in other observables, we consider it is relevant to discuss the impact of these chiral logs.
As in the NLO case, we perform two types of fits, using the NLO and NNLO theoretical expressions given in terms of F and F π for various observables. We explicitly present the fit results with F in the theoretical expressions and use the alternative fits expressed in terms of F π to estimate the systematic errors, due to working up to NNLO in δ and neglecting higher orders. According to the Lagrangian in Eq. (5), eleven additional unknown LECs appear at NNLO and there will be seventeen parameters in total for the NNLO study. At the present precision of the lattice simulations and phenomenological inputs, it is impossible to obtain sensible and stable fits if we free all of the seventeen parameters. Therefore, we need to take other independent determinations for some of the LECs in order to proceed the NNLO study.
We mention that the state-of-art determinations of the O(p 4 ) LECs in SU (3) χPT suffer uncertainties from the many poorly known O(p 6 ) LECs [20,21]. Because of the large number of barely known O(p 6 ) LECs, it is rather difficult to get conclusive results in the present two-loop SU (3) χPT studies [20,21]. In the present work, there are five O(p 6 ) LECs, i.e. C 12 , C 14 , C 17 , C 19 , C 31 , in Eq. (5) and we cannot make precise determinations of these C i parameters here. Maybe when taking into account the scattering data, one can make more stringent constraints on the C i LECs in U (3) χPT. But this is beyond the scope of current work. Instead we take the C i values from the Dyson-Schwinger-like approach given in Ref. [52], where all of the O(p 6 ) C i at leading N C are predicted. In order to show the dependences of the final results on the C i values, we also perform other fits by using their updated determinations [53]. Like in Ref. [20], we multiply the O(p 6 ) C i from Refs. [52,53] by a global factor α and consider α as a free parameter in the fits. In this way, we partially compensate the large uncertainties of the C i parameters.
For the operators proportional to v 2 , L 18 and L 25 in Eq. (5), they are not present in SU (3) χPT and purely contribute to the η-η ′ mixing, being irrelevant to the pion and kaon observables. Since the η-η ′ mixing parameters have already been satisfactorily described in the NLO fits, we do not further include v (2) 2 , L 18 and L 25 at NNLO study. 2 Their inclusion in the present analysis tend to make the fit unstable. Clearly studying more η ( ′ ) related observables it would be possible to extract these parameters but this is out of the reach of the present analysis. A global fit is too unconstrained, being unstable and producing values of the latter couplings compatible with zero within uncertainties. Then we are left with three O(N 0 C , p 4 ) operators: L 4 , L 6 and L 7 , which have corresponding parts in SU (3) χPT. Since U (3) and SU (3) χPT contain different dynamical degrees of freedom, the corresponding LECs from the two theories can be different. A typical example is the L 7 parameter in SU (3) χPT, which is demonstrated to be dominated by the singlet η 0 state [1]. Since in U (3) χPT the singlet η 0 has been explicitly introduced, the value of L 7 in this theory 2 Indeed, in this work M 0 and v 2 (2m 2 K + m 2 π ), which is the parameter we are actually extracting. The contribution v in general can not be recovered in the present numerical fits, due to the presence of far too many parameters in the problem and the large uncertainties of the lattice simulation data, specially for the determinations of m η ′ .
can be totally different from L SU(3) 7 in SU (3) case. While for other O(p 4 ) LEcs, such as L i=4,5,6,8 , the differences between U (3) and SU (3) χPT are not expected to be as large as the L 7 case, since they do not receive the tree-level contributions from the η 0 state.
Another subtlety to take into account is that m η and m η ′ appear in the chiral loops and, at the same time, the final expressions of m η and m η ′ depend on the these loops as well. In order to avoid making the complicated iterative procedure to obtain the η-η ′ mixing parameters, we use the LO formulas for m η and m η ′ in the chiral loops. The differences caused by this simple treatment and the strict iterative procedure are beyond the NNLO precision in δ expansion, since the chiral loops themselves are already NNLO. Our simple solution is also justified by the fact that the LO description of m η and m η ′ is in qualitative agreement with the lattice simulation data, as shown in Sect. III A. Since the qualitative agreement between the LO formulas and the lattice simulation data requires the value of M 0 to be around 835.7 MeV, as given in Eq. (55), we fix M 0 = 835.7 MeV in the following discussions. This also helps to stabilize the NNLO fits, with its many free parameters. Other useful criteria to discriminate reasonable fits are the a priori ranges estimated in Eq. (57), since the fits with large magnitudes of Λ 1 and Λ 2 imply unphysically large corrections to the η-η ′ mixing parameters and the breakdown of the δ expansion. In the following we only present the fit results that are consistent with Eq. (57).
With all of the above setups, the values of parameters from the NNLO fits are summarized in Table V in the same units as before.
It is clear that the parameters resulting from fits with different C i inputs slightly differ from one another. We remind that the first error bar for each parameter in Table V corresponds to the statistical one directly from the fits and the second error bar stands for the systematic one, which is estimated, as usual, from the variation of the parameter from the NNLO fits with the theoretical expressions in terms of F and those expressed as functions of F π .
At NNLO, one has the contributions from the chiral loops and the O(p 6 ) LECs, which make our determinations in Table V closer to the recent two-loop results of the SU (3) χPT LECs, comparing with the NLO determinations in Table I. Some typical trends of the values of parameters from the NLO study in Table I to the NNLO one in Table V are summarized now. The axial-vector decay constant F at leading N C and chiral limit is reduced at NNLO, which is mainly due to the inclusion of L 4 . Our conclusion is based on the fact that strong correlations between F and L 4 always appear, which has been confirmed in previous study [22,23]. For L 5 and L 8 , we find that their values are obviously reduced compared to the NLO determination and become closer to the two-loop results in Ref. [21]. As mentioned in the former reference, the discussions in the two-loop SU (3) χPT are sensitive to the value of the 1/N C suppressed LEC L 4 . The present study provides an independent determination for this parameter and for the 1/N C suppressed LEC L 6 as well. We mention that our determinations of L 4 have opposite signs with respect to that in Ref. [21], which may be the source of the smaller F obtained in that reference. Notice that the present values of L 4 , L 5 , L 6 , L 8 are rather compatible with the combinations of 2L 8 − L 5 and 2L 6 − L 4 given in Ref. [54]. Fit solutions with larger Λ 1 and Λ 2 than those in Eq. (57) (out of the a priori range (57)) are discarded: they are not considered as reasonable physical solutions and will not be discussed any further. According to the values of α in the two fits, it seems that our study somewhat prefers smaller magnitudes of the O(p 6 ) LECs than those from the Dyson-Schwinger approach given in Refs. [52,53] and also prefers a global change of sign with respect to Eqs. (58) and (59). We have investigated the impact of fitting α but releasing one of the O(p 6 ) LECs as an independent parameter (e.g., C 14 ), but no definitive conclusion could be extracted. These puzzles cannot be resolved here and it is definitely interesting and necessary to further investigate the values of the O(p 6 ) LECs in the future.
The various plots from the NNLO fits are shown in Fig. 2 for m η and m η ′ , Fig. 3 for m K , Fig. 4 for F π and F K , and Fig. 5 for the ratio F K /F π , together with the NLO results and the lattice simulation data and experimental inputs. The shaded area surrounding each curve represents the statistical error band. The figures from NNLOFit-B are compatible with those from NNLOFit-A within the uncertainties, so we only show the results for the former in Figs. 2, 3, 4 and 5.
In addition, to demonstrate the effects by using the LO masses in the higher order corrections, instead of the renormalized ones, we explicitly show the results for m η and m η ′ expressed in terms of the LO masses m π , m K and 1/F 2 in Fig. 2, with the lines labeled as NNLOFit-B-m π,K . The values of the LECs when plotting these lines are exactly the same as those from the NNLOFit-B column in Table V. In this way, one can directly see the differences due to the N 3 LO truncation uncertainty caused by using the renormalized masses and the LO ones at the NNLO level. According to Fig. 2, we conclude that the differences for m η and m η ′ caused by using different types of masses in the higher order corrections are rather within the statistical uncertainties from the fits and therefore the differences should be perfectly compatible within the total uncertainties after taking into account the systematic ones in Table V. We verify that similar conclusions are obtained for other cases. In order not to overload the plots in other figures, we shall not explicitly show the results given in terms of m π and m K .
From Figs. 2, 3, 4 and 5, we observe, when compared with the curves of the NLO study, slight improvements in the reproduction of the masses for η, η ′ and significant ones for m K , F π , F K and the ratios of F K /F π . Moreover the χ 2 for the NNLO fits are greatly reduced compared with χ 2 for the NLO ones, indicating that the NNLO corrections are important at the present level of precision and essential to simultaneously describe the lattice simulation data and experimental inputs of the light pseudoscalar mesons π, K, η and η ′ .    Table II for the explanation of the phenomenological inputs. The first error for each quantity corresponds to the statistical one and the second error denotes the systematic one. See the text for details.

IV. CONCLUSIONS
In this article we have performed a thorough study on the η-η ′ mixing, and axial-vector decay constants for the pion and kaon, up to next-to-next-to-leading order in δ expansion within U (3) chiral perturbation theory. We have carried on a detailed scrutiny and discussions of our results, which have been carefully compared to other works in literature for the η-η ′ mixing. A general mixing formalism, including the higher-derivative terms and kinematic mixing cases, has been addressed in detail. The connections between the mixing parameters from the popular two-mixing-angle scheme and the low energy constants from chiral perturbation theory have been established, both for the singlet-octet basis and the quark-flavor basis.
The considered quantities, including the masses of η, η ′ and K, the quark mass ratio of m s / m, the parameters in the two-mixing-angle scheme and the π, K decay constants have been confronted with recent lattice simulations and phenomenological inputs. We find that the next-to-leading-order fits yield satisfactory descriptions for the masses of the three pseudoscalar mesons as functions of m 2 π and the four mixing parameters (F 0 , F 8 , θ 0 , θ 8 ), producing in addition reasonable values of low energy constants. Nonetheless, when the π and K decay constants are included together with the masses and mixing parameters in the fits, the next-to-leading-order analyses are inadequate and it is necessary to step into the next-to-next-to-leading-order study. Using the O(p 6 ) LECs determinations from a Dyson-Schwinger-like approach [52,53] multiplied by a global factor, we are able to achieve a reasonable description for all of the physical quantities considered above and the resulting values for the leading N C O(p 4 ) low energy constants L 5 and L 8 turn to be compatible with the very recent two-loop determinations in Ref. [21]. Therefore we conclude that the large N C U (3) chiral perturbation theory offers a concise theoretical framework that is able to simultaneously reproduce accurately the general η-η ′ mixing and to provide sophisticated enough expressions to describe the chiral extrapolations of the π and K decay constants and masses.
Our results are also useful for future phenomenological studies of different processes involving η and η ′ . Combining Eq. (21) or Eq. (23) with Table VI, one can directly find the relations between the physical states η, η ′ and the octet-singlet bases η 8 , η 0 or the quark-flavor bases η q , η s . These relations are consistent with the requirements from the recent lattice simulations and phenomenology.
Finally, it is worthy to remark that some of the parameters in our best analysis (NNLOFit-B) in Table V have been determined with relatively small errors. For instance, the NLO parameters Λ 1,2 , which are fitted up to O(N −2 C ) in the NNLO analysis, become Λ 1 = −0.04 ± 0.06 ± 0.13 , Λ 2 = 0.14 ± 0.10 ± 0.40 .
For completeness, we also give the results in terms of the LO masses m π and m K and 1/F 2 . Only the terms with L i L j , being L i and L j the NLO LECs in Eq. (4), will be different, comparing with the expressions in terms of m π and m K and the other parts remain the same, apart from the obvious replacement of the renormalized masses by the LO ones. Therefore, we only give the parts that are different from those in terms of m π , m K and 1/F 2 and it turns out that in this case all of the L i L j terms for δ η , δ η ′ , δ k , δ m 2 η , δ m 2 η ′ , δ m 2 vanish.