On the transition form factors of the axial-vector resonance $f_1(1285)$ and its decay into $e^+e^-$

Estimating the contribution from axial-vector intermediate states to hadronic light-by-light scattering requires input on their transition form factors (TFFs). Due to the Landau-Yang theorem, any experiment sensitive to these TFFs needs to involve at least one virtual photon, which complicates their measurement. Phenomenologically, the situation is best for the $f_1(1285)$ resonance, for which information is available from $e^+e^-\to e^+e^- f_1$, $f_1\to 4\pi$, $f_1\to \rho \gamma$, $f_1\to \phi \gamma$, and $f_1\to e^+e^-$. We provide a comprehensive analysis of the $f_1$ TFFs in the framework of vector meson dominance, including short-distance constraints, to determine to which extent the three independent TFFs can be constrained from the available experimental input, a prerequisite for improved calculations of the axial-vector contribution to hadronic light-by-light scattering. In particular, we focus on the process $f_1\to e^+e^-$, evidence for which has been reported recently by SND for the first time, and discuss the impact that future improved measurements will have on the determination of the $f_1$ TFFs.


Introduction
The interaction of an axial-vector resonance A with two electromagnetic currents is subject to the venerable Landau-Yang theorem [1,2], which states that a spin-1 particle cannot decay into two on-shell photons. Accordingly, the decay A → γγ is forbidden, and the simplest process from which information on the general A → γ * γ * matrix element can be extracted is the singly-virtual process. Such measurements are available from the (spacelike) reaction e + e − → e + e − A for A = f 1 (1285) and A = f 1 (1420) [3][4][5][6][7][8], providing results for the so-called equivalent two-photon decay widthΓ γγ as well as some constraints on the momentum dependence of the process. Assuming U (3) symmetry then allows some inference for A = a 1 (1260), but other direct phenomenological input is scarce.
As a first step to improving this situation, a systematic analysis of the axial-vector TFFs has been presented recently in Ref. [80], including the decomposition into Lorentz structures that guarantee the absence of kinematic singularities in the TFFs, following the recipe of Bardeen, Tung, and Tarrach (BTT) [81,82], and the derivation of short-distance constraints in analogy to the light-cone expansion of Brodsky and Lepage (BL) [83][84][85]. Here, we provide a comprehensive analysis of the TFFs of the f 1 (1285), for which the most phenomenological input is available. In addition to e + e − → e + e − f 1 [5][6][7], there are data for f 1 → 4π [86], f 1 → ργ [86,87], f 1 → φγ [86,88], and, most recently, f 1 → e + e − [89], all of which probe different aspects of the TFFs, as we will study in detail in this paper.
Given that there are three independent TFFs, in contrast to just one in the case of pseudoscalar mesons, a full dispersive reconstruction as in Refs. [26,27,[90][91][92][93][94] for the π 0 or in progress for η, η [95][96][97][98][99] appears not feasible given the available data. Accordingly, we will study the simplest vector-meson-dominance (VMD) ansatz, to elucidate which parameters can presently be determined from experiment. In contrast to previous work [100,101], our parameterization ensures the absence of kinematic singularities, includes short-distance constraints, and accounts for the spectral function of the isovector resonances. In particular, we critically examine which of the processes listed above do allow for an unambiguous extraction of TFF properties. We focus on the f 1 → e + e − decay, evidence for which has been observed only recently by the SND collaboration [89], with future improvements possible in the context of the ongoing program to measure e + e − → hadrons cross sections. Further, since this process involves a loop integration that depends on all three TFFs, it should provide some sensitivity also to the doubly-virtual TFFs, which are particularly difficult to measure otherwise.
The outline of this article is as follows: in Sec. 2, we review the BTT decomposition of the A → γ * γ * matrix element as well as the asymptotic constraints. In Sec. 3, we then construct a minimal VMD ansatz, an extended version, and study their asymptotic behavior. The tree-level processes e + e − → e + e − f 1 , f 1 → 4π, and f 1 → V γ (V = ρ, φ, ω) used to constrain the parameters are discussed in Sec. 4, followed by the f 1 → e + e − decay in Sec. 5. The full phenomenological analysis is provided in Sec. 6, before we summarize our findings in Sec. 7. Further details are provided in the appendices.

.2 Asymptotic constraints
In analogy to the asymptotic limits of the pseudoscalar TFF derived in Refs. [83][84][85], one can use a light-cone expansion to obtain the asymptotic behavior of the axial-vector TFFs.
Using the distribution amplitudes from Refs. [103,104], the asymptotic behavior is given by [80] F 1 (q 2 1 , q 2 2 ) = O(1/q 6 i ), where we generically denoted powers of asymptotic momenta by q i = q 1 , q 2 and the wave function φ(u) = 6u(1 − u) is the asymptotic form that already contributes to the pseudoscalar case. In writing Eq. (2.12), we furthermore defined an effective decay constant where the decay constants F a A are defined via (2.14) The Gell-Mann matrices λ a and the conveniently normalized unit matrix λ 0 = 2/3 1 determine the flavor decomposition, with the flavor weights C a in the effective decay constant given by C a = 1/2 Tr(Q 2 λ a ), i.e., C 0 = 2/(3 √ 6), C 3 = 1/6, and C 8 = 1/(6 √ 3). In Eq. (2.12) we retained the leading mass effects in the denominator, but stress that this does not suffice for a consistent treatment of such corrections. We will thus mostly set m A = 0 in the denominators when implementing the short-distance constraints, but address the treatment of the leading mass effects in App. A. Rewriting the results in terms of the average photon virtuality Q 2 and the asymmetry parameter w,  Figure 1: Asymmetry functions f 2 (w) and f 3 (w), Eq. (2.17), with values for the limiting cases w ∈ {−1, 0, 1} of f 2 (w), corresponding to q 2 1 = 0, q 2 1 = q 2 2 , and q 2 2 = 0, respectively. The analogous limits for f 3 The asymmetry functions f 2/3 (w) are shown in Fig. 1, where we also illustrate the values of the function f 2 (w) for the limiting cases w = −1 (q 2 1 = 0), w = 0 (q 2 1 = q 2 2 ), and w = 1 (q 2 2 = 0); since f 2 (−w) = −f 3 (w), the analogous limits for f 3 (w) follow accordingly. More specifically, the symmetric doubly-virtual and singly-virtual asymptotic limits of the TFFs-the latter often being referred to as the BL limit-become while the expressions for F 2 (0, q 2 ) and F 3 (q 2 , 0) diverge. Given that the derivation of Eq. (2.12) can only be justified from the operator product expansion for |w| < 1/2 [105,106], the singly-virtual limits need to be treated with care. 3 However, physical helicity amplitudes only depend on the well-defined limits in Eq. (2.18), in such a way that the problematic limits F 2 (0, q 2 ) and F 3 (q 2 , 0) do not contribute to observables. We will return to this point in the context of the f 1 → e + e − loop integral.

Vector meson dominance
Given the scarcity of data for axial-vector resonances, we will perform our phenomenological analysis in the context of a VMD description, which has proven to provide successful approximations for a host of low-energy hadron-photon processes [111][112][113][114][115][116]. Most notably, the underlying assumption that the interaction is dominated by the exchange of vector mesons predicts the charge radius of the pion at the level of 10%. Even though the ensuing model dependence is hard to estimate a priori, this approach allows us to analyze all experimental constraints simultaneously in a common framework, which could be refined as soon as improved data become available.
To construct VMD representations of the TFFs as defined in Sec. 2, it is convenient to recast them in terms of their symmetric (s) and antisymmetric (a) combinations with the indicated symmetry properties under the exchange of momenta, q 2 1 ↔ q 2 2 . Consequently, the basis of structures transforms according to where these functions fulfill the same symmetry properties under photon crossing. Given this alternative basis, the equivalent two-photon decay width, Eq. (2.11), becomes and the tensor matrix element of Eq. (2.5) takes the form

Quantum numbers and mixing effects
Since by far the best phenomenological information is available for the f 1 ≡ f 1 (1285), we will focus on this resonance in the remainder of this work, but remark that information on the f 1 ≡ f 1 (1420) and the a 1 (1260) can be derived when assuming U (3) flavor symmetry.
As a first step towards constructing our VMD ansatz for the TFFs, 4 we review the relevant quantum numbers and mixing patterns. From the G-parity G = + of the f 1 , it is immediately clear that both photons have to be either in their isoscalar or isovector state when neglecting isospin-breaking effects. Hence, the VMD coupling can only proceed via ρρ-like or via some combination of an ωand φ-like vector meson, each of which will be discussed in turn in Sec. 3.2 and Sec. 3.3, respectively. As we will show in the following, it is the isovector channel that dominates, with isoscalar corrections typically at the level of 5%.
To this end, we have to take into account mixing effects between the (physical) mesons of the corresponding J P C = 1 ++ axial-vector nonet, i.e., the mixing pattern [86] where f 0 and f 8 denote the isoscalar singlet and octet states of the J P C = 1 ++ nonet and θ A is the corresponding mixing angle. Pure octet/singlet mixing is reproduced for θ A = π/2, whereas ideal mixing is obtained for θ A = arctan(1/ √ 2). Including only the two resonances f 1 and f 1 , the U (3) parameterization of the J P C = 1 ++ axial vectors reads and when splitting the charge matrix into isovector and isoscalar components according to Q = Q 3 + Q 8 , one finds Using the mixing angle θ A = 62(5) • as determined by the L3 collaboration [7,8], see Sec. 4.1, one thus finds that the ratio R S/V of isoscalar to isovector contributions for the f 1 γγ coupling is given by

Isovector contributions
For the isovector contributions to the TFFs in Eq. (3.1) we include the ρ ≡ ρ(770) and the ρ ≡ ρ(1450), since this is the minimal particle content that produces a non-vanishing contribution for the antisymmetric TFFs. We propose the minimal parameterizations , (3.10) where Γ ρ (q 2 ) and Γ ρ (q 2 ) are yet to be specified energy-dependent widths. 5 Moreover, ρρ and ρ ρ terms will be added to F s (q 2 1 , q 2 2 ) below, to help incorporate the asymptotic constraints from Sec. 2.2. We adopt the dispersion-theoretical point of view to model the singularities of the TFFs based on vector-meson poles, and refrain from constructing these using effective Lagrangians in order to facilitate the implementation of high-energy constraints.
Concerning the energy-dependent width Γ ρ (q 2 ), the decay ρ → ππ is described by where γ ρ→ππ (q 2 ) is constructed to be in accord with the behavior of the decay width for variable M 2 ρ = q 2 , see Eq. (B.9), and Γ ρ is the total width of the ρ meson. For the energydependent width Γ ρ (q 2 ) on the other hand, we will consider two different parameterizations. First, we assume the decay channel ρ → 4π to be dominant and thus adopt the nearthreshold behavior of the four-pion phase space [117,118]. Second, we construct a spectral shape from the decay channels ρ → ωπ (ω → 3π) and ρ → ππ, neglecting, however, another significant contribution from ρ → a 1 π (a 1 → 3π) [86]. These parameterizations read where γ ρ →4π (q 2 ) is taken from Refs. [117,118] and Γ ρ is the total decay width of the ρ meson, and where (3.14) Estimates for the branching fractions required to evaluate these expressions are provided in App. B. Finally, the standard form of the ρ → ππ spectral function in Eq. (3.11) proves disadvantageous for the evaluation of superconvergence relations in Sec. 3.4 due to its highenergy behavior. We thus follow Refs. [119,120] and introduce barrier factors according to where concurrent adjustments to the ρ → ππ channel of Γ (ωπ,ππ) ρ (q 2 ), Eq. (3.13), are implied. In the end, the numerical impact of the choice of the ρ spectral function is subdominant, and our results will be shown for Γ (2) ρ (q 2 ) (both for the ρ and the 2π component of Γ (ωπ,ππ) ρ (q 2 )), which is identified as the best phenomenological description for the ρ meson in Ref. [119].
For the one-loop process f 1 → e + e − discussed in Sec. 5 we will use dispersively improved variants of the isovector form factors to ensure the correct analyticity properties when inserting the TFFs into the loop integral. The corresponding spectral representations are constructed from the energy-dependent widths, i.e., where the dispersive ρ and ρ propagators are given by The spectral functions are and the threshold s thr ∈ {16M 2 π , 4M 2 π } depends on the choice of Γ ρ (q 2 ), Eq. (3.12) or Eq. (3.13). The normalization constants N a and N s are introduced in order to retain the form factor normalizations C a 1/2 and C s from Eq. (3.10),  i.e., to ensure that the constants C a 1/2 and C s carry the same meaning in the original and the dispersively improved VMD parameterizations, see Table 1. With these conventions, we will drop the distinction between F i (q 2 1 , q 2 2 ) and F i (q 2 1 , q 2 2 ) in the following, the understanding being that f 1 → e + e − is evaluated with the dispersively improved variants.
Given that excited ρ mesons need to be introduced for the antisymmetric TFFs, it is natural to consider an extended VMD parameterization of the symmetric form factor including ρρ and ρ ρ terms, , (3.20) which is normalized in such a way that F I=1 s (0, 0) = C s = F I=1 s (0, 0). Here, 1 and 2 could be treated as additional free parameters, but instead we will use this freedom to match to the asymptotic constraints in Sec. 3.4. Similarly to Eq. (3.16), the spectral representation for F I=1 s (q 2 1 , q 2 2 ) is given by with normalization see Table 1.

Isoscalar contributions
In the following, we estimate the isoscalar contributions to the TFFs of Eq. (3.1) under the assumption of U (3) flavor symmetry, where we will include the resonances ω ≡ ω(782) and φ ≡ φ(1020) as well as their excited equivalents ω ≡ ω(1420) and φ ≡ φ(1680) into our parameterization. Mixing effects between the (physical) mesons of the corresponding J P C = 1 −− vector-meson nonets are taken into account via the pattern [86] ω ( ) where ω 0( ) and ω 8( ) denote the isoscalar singlet and octet states of the respective vectormeson nonet with mixing angle θ V ( ) . For our considerations, we assume both nonets to be ideally mixed, i.e., θ V = arctan(1/ √ 2) = θ V . Finally, we need the U (3) parameterization of the J P C = 1 −− vector mesons, which reads when including only the aforementioned resonances. Since the U (3) couplings f 1 ωφ, f 1 ω φ, and f 1 ωφ vanish for ideally mixed vector mesons, we propose the minimal parameterizations .

(3.25)
The resonances ω and φ should be well described by a narrow-resonance approximationwith M 2 V → M 2 V − i for time-like applications-while for a realistic description of the excited-state isoscalar resonances their widths would need to be taken into account. Due to the expected smallness of the isoscalar contributions, see Eq. (3.9), we refrain from giving an extended VMD parameterization analogous to Eq. (3.20).
With the U (3) parameterization of the axial-vector mesons, Φ A µ , and the charge matrix Q from Sec. 3.1, the ratios of isoscalar to isovector couplings are found to be 6 C ωω 6 The notation is to be understood in such a way that for each term the prefactor of the fields indicated as a subscript is taken, with the U (3) parameterizations from Eq. (3.6), Eq. (3.7), and Eq. (3.24). In the ratios only the traces are relevant, as the common Lagrangian parameters cancel.
which, using the mixing angle θ A = 62(5) • as determined by the L3 collaboration [7,8], see Sec. 4.1, implies (34). (3.27) The additional suppression in Eq. (3.9) then results from a cancellation between ω and φ contributions In practice, we will restrict the analysis of isoscalar contributions to the symmetric TFF. First, F s (q 2 1 , q 2 2 ) gives the dominant contribution to the observables, so that the most important isoscalar correction is expected from there. In addition, for the antisymmetric TFFs we would need to include the excited ω and φ states, incurring significant uncertainties from their spectral functions and, especially for the f 1 → e + e − application, the asymptotic matching due to their large masses. Alternatively, isoscalar antisymmetric TFFs could be produced via deviations from ideal φ-ω mixing, but again the uncertainties would be difficult to control. For these reasons we conclude that the isoscalar contributions to the antisymmetric TFFs should be irrelevant at present, with potential future refinements once better data become available.

Asymptotics
The VMD representations for the TFFs should comply with the asymptotic constraints reviewed in Sec. 2.2, mainly to ensure that the f 1 → e + e − loop integral does not receive unphysical contributions in the high-energy region. We will focus on the isovector amplitudes, given the strong suppression of the isoscalar contributions. Translated to the basis of (anti-)symmetric TFFs, we have see Fig. 2. The symmetrical doubly-virtual limits become (λ ≈ 1) but upon symmetrization all singly-virtual limits of F a 2 /s (q 2 1 , q 2 2 ) diverge. For this reason, the asymptotic limits for F a 2 /s (q 2 1 , q 2 2 ) cannot be considered in isolation, but need to be implemented in such a way as to reproduce the physical behavior of F 2/3 (q 2 1 , q 2 2 ). Figure 2: Asymmetry functions f a 2 (w) and f s (w), Eq. (3.29), with values for the limiting cases w ∈ {−1, 0, 1}, corresponding to q 2 1 = 0, q 2 1 = q 2 2 , and q 2 2 = 0 respectively.
We first consider the asymptotic behavior of the minimal VMD parameterization, Eq. (3.10), In this case, the scaling is correct in the doubly-virtual direction of F I=1 a 1 /s (q 2 1 , q 2 2 ), while F I=1 a 2 (q 2 1 , q 2 2 ) drops too fast and the singly-virtual limits too slowly, see Table 2. Phenomenologically, the symmetric TFF gives the dominant contribution to f 1 → e + e − , see Sec. 5, so that here also the coefficient deserves some attention. Comparing the asymptotic limit of Eq. (3.10) with Eq. (3.30), the VMD ansatz for F s (q 2 1 , q 2 2 ) implies the following estimate for the effective decay constant defined in Eq. (2.13): where we already used the L3 result for C s including the isoscalar contribution, see Eq. (4.7) below. Within uncertainties, this value agrees with the result from light-cone sum rules (LCSRs) [80,104] F eff so that even the minimal VMD ansatz should display a reasonable asymptotic behavior. To go beyond this minimal implementation, we now turn to the extended VMD ansatz for F s (q 2 1 , q 2 2 ). We follow the strategy from Refs. [26,27] and add an explicit asymptotic term that incorporates the correct doubly-virtual behavior, obtained by rewriting Eq. (2.12) in terms of a dispersion relation; see also Ref. [121]. Accordingly, we need to ensure that the isovector VMD contribution to F s (q 2 , q 2 ) behaves ∝ 1/q 6 , resulting in This leaves the freedom to choose 1 , which we use to implement the physical singly-virtual Further, the coefficient of 1/q 4 in the resulting F I=1 2 (q 2 , 0) only depends on C s , and matching to Eq. (2.18) implies F eff reasonably close to the LCSR estimate of Eq. (3.33). In general, the choice for 1 in Eq. (3.35) enforces the expected singly-virtual behavior at the expense of a large coefficient, e.g., for C a 2 = 0 one has 1 = −1.08, so that a better low-energy phenomenology might be achieved when considering 1 a free parameter instead. We will continue to use Eq. (3.35) as a benchmark scenario in comparison to the minimal VMD ansatz, keeping this caveat regarding 1 in mind. In choosing the above 1/2 , we did not take the spectral representations of Eq. (3.16) and Eq. (3.21) into account, which would lead to a set of superconvergence relations that need to be fulfilled, but instead made an approximate choice in terms of Eq. (3.20) and  Table 3: Numerical values of P 0 ρ and P 0 ρ , Eq. (3.38), as obtained with the parameterizations Γ Solving this for 2 and 1 , we find in accordance with Eq. (3.34) and Eq. (3.35) upon the replacements Numerical values for P 0 ρ and P 0 ρ are collected in Table 3. These results show that most correction factors are close to unity, in which case the only potentially significant correction arises from the different normalizations N a and N s for 1 , see Table 1. However, our central results will employ Γ (ωπ,ππ) ρ (q 2 ), and given the abovementioned caveats in the choice of 1 , we conclude that at the current level of accuracy the naive VMD expressions Eq. (3.34) and Eq. (3.35) are sufficient.
The doubly-virtual behavior is implemented as follows [26,27]: first, we rewrite the asymptotic form factors F 2 (q 2 1 , q 2 2 ) and F 3 (q 2 1 , q 2 2 ) from Eq. (2.12) into a double-spectral representation, which allows us to isolate the different energy regions, in particular those that give rise to the correct asymptotic limits. Setting m A = 0 in the respective integrands of Eq. (2.12), we observe that take exactly the same form as for the pseudoscalar case, except for the partial derivatives with respect to q 2 i . Accordingly, the same arguments as in Refs. [26,27,121] apply, and the integral over the wave function can be formally expressed by a double-spectral representation , The asymptotic form arises from the high-energy part of these integrals, so that, to avoid overlap with the VMD contribution at low energies, we impose a lower cutoff s m , which, in the language of LCSRs, could be identified with the continuum threshold. Evaluating the partial derivatives and dropping surface terms in the evaluation of the δ distribution [26,27], we find By construction, the asymptotic contributions in this form saturate the doubly-virtual limits of Eq. (3.30), while not affecting the singly-virtual contributions F 2 (q 2 , 0), F 3 (0, q 2 ) already taken into account via the extended VMD representation. The opposite-unphysical-cases F 2 (0, q 2 ), F 3 (q 2 , 0), which do not contribute to helicity amplitudes, are equally suppressed in the f 1 → e + e − loop integral, see Sec. 5. Given that m A > 1 GeV, it is also worthwhile to consider the potential impact of mass corrections to the asymptotic constraints. A formulation in terms of a generalized double-spectral density is given in App. A.
In conclusion, the extended VMD ansatz together with the asymptotic contribution of Eq. (3.44) complies with the short-distance constraints of Eq. (2.12), apart from the singly-virtual behavior of F a 1 (q 2 1 , q 2 2 ) and small violations due to the isoscalar contributions of the form factors, see Eq. (3.9). As we will demonstrate below that F a 1 (q 2 1 , q 2 2 ) gives the smallest contribution to the f 1 → e + e − loop integral, see Eq. (5.5), the resulting VMD representation should provide a decent approximation to its high-energy part. In particular, the sensitivity to the high-energy assumptions can be monitored by comparing the two VMD variants constructed in this section.

Tree-level processes
The VMD parameterizations constructed in the previous section involve the free parameters C a 1 , C a 2 , and C s (and, for the extended variant, the onset of the asymptotic contributions s m ). In the following, we collect the available data that can, in principle, determine these parameters, starting with the processes in which the TFFs appear at tree level: 1. e + e − → e + e − f 1 , which mainly determines the equivalent two-photon decay width 3. f 1 → ργ, whose branching fraction and helicity components encode information on the TFFs, see Sec. 4.3.
In a more rigorous, dispersive, reconstruction of the TFFs, the (partially) hadronic final states would serve as input to a determination of their discontinuities. The strategy to investigate the impact of these reactions on a determination of the various TFFs has already been followed in Refs. [100,101], albeit with rather different form factor parameterizations. Moreover, we investigate the following tree-level decays: 4. f 1 → φγ and f 1 → ωγ, where the measured branching fraction of the former allows for a consistency check of our U (3) assumption for the isoscalar TFFs and the latter predicts a branching ratio that can be confronted with potential future measurements, see Sec. 4.4.
In contrast to (pseudo-)scalar or tensor resonances, axial-vector resonances are only visible in e + e − collisions, see Fig. 3, as long as at least one of the photons is off shell, a direct consequence of the Landau-Yang theorem [1,2]. The required challenging measurements have been performed for the f 1 and f 1 , by the MARK II [3,4], the TPC/Two-Gamma [5,6], and, more recently, by the L3 [7,8] collaborations. With both measurements required to constrain the mixing angle θ A from the data, we will restrict our analysis to the L3 data, given that they are more accurate than the results from the preceding experiments. The L3 analyses are based on the model of Ref. [122], which assumes F 1 (q 2 1 , q 2 2 ) = 0 for the first form factor from Eq. (2.5) and uses a dipole ansatz for Under the assumption B(f 1 → KKπ) = 1-which appears justified in light of the smallness of the other available channels [86]-the measured parameters are where the quoted uncertainties are statistical and systematic, respectively. Employing the two-photon decay widths of the f 1 and f 1 , the mixing angle of the J P C = 1 ++ axial-vector nonet as defined in Eq. (3.5) can be extracted as follows: one calculates the coupling of the axial-vector mesons f 1 and f 1 to two photons in analogy to Eq. (3.8), yielding so that using the formula for the equivalent two-photon decay width Γ γγ , Eq. (2.11), one finds where θ 0 = arcsin(1/3). Solving for θ A and inserting the above values for Γ f 1 γγ and Γ f 1 γγ , one finds the result of Refs. [7,8], where the statistical and systematic uncertainties have been added in quadrature. Next, the measurement of Γ f 1 γγ determines the normalization of the symmetric TFF, |C s | = |F I=1 s (0, 0)| when neglecting the isoscalar contributions, according to Eq. (3.3), |C s | = 0.89 (10). (4.6) Taking into account the isoscalar contributions and, in particular, the ratios R ω and R φ of isoscalar to isovector couplings, Eq. (3.27), the normalization of the symmetric TFF becomes |F I=1 which is slightly larger than Eq. (4.6), as expected from the negative ratio found in the estimate of Eq. (3.9). In the following, we will use Eq. (4.7) for the normalization of the symmetric TFF.
In addition, Eq. (4.2) determines the slope of F 2 (q 2 , 0), based on the assumption of a dipole form. The asymptotic behavior matches onto Eq. (2.18) with [80] F eff below both the LCSR estimate, Eq. (3.33), and the effective decay constant implied by VMD, Eq. (3.32), and close to the scale derived from the singly-virtual behavior of the extended VMD representation, Eq. (3.36). 7 The uncertainty in Eq. (4.8) is mainly driven by the dipole parameter Λ D . In fact, most of the data points measured by the L3 collaboration lie well below the obtained dipole scale, in such a way that the data should be similarly well described by a monopole ansatz, when adjusting the slopes of the parameterizations to coincide at q 2 = 0. The corresponding monopole scale becomes thus providing strong motivation for the VMD representation constructed in Sec. 3.
To constrain the singly-virtual VMD limits further, we need to match the L3 parameterization onto the full description of the e + e − → e + e − f 1 cross section, which depends on the combination [ (4.11) The normalization agrees by construction, while matching the slopes at q 2 = 0 leads to for the minimal VMD representation, and for the extended one. The factor N ωφ = 1+R ω +R φ arises from accounting for the isoscalar terms in the normalization, see Eq. (4.7).

f 1 → 4π
In addition to e + e − → e + e − f 1 , the normalization of the symmetric TFF would be accessible in the process f 1 → ρρ → 4π if the ρ intermediate states largely saturated the decay within 7 Matching the effective decay constant in the doubly-virtual direction to the quark model of Ref. [122] instead, one would obtain F eff (42) MeV, closer to Eq. (3.32) and Eq. (3.33). This reflects the factor 3/2 by which the relative coefficients of the singly-and doubly-virtual limits differ between the quark model and the BL prediction [80]. Figure 4: Feynman diagrams for f 1 → π + π − π + π − via two ρ mesons. Since the two π + and π − are respectively indistinguishable, there exist two contributions (left and right).
regions of the phase space reasonably close to their mass shell. In fact, up to corrections due to the two-pion channel ρ → π + π − , such an identification appears natural within the VMD approach. In constructing an amplitude M(f 1 → π + π − π + π − ), which can be obtained by means of M(f 1 → ρ 0 * ρ 0 * ) and the ρππ coupling dictated by Eq. (B.8), only the symmetric form factor F I=1 s (q 2 1 , q 2 2 ) and the symmetric Lorentz structure T µνα s (q 1 , q 2 ) are relevant under the above assumptions and when restricting to the minimal VMD parameterization. More specifically, we use the amplitude M(f 1 → γ * γ * ), in the decomposition of Eq. (3.4), and remove the external photons by dropping the relevant ρ-meson propagator poles and the factors of e, at the same time dividing by the ργ coupling g ργ , Eq. (B.7), for each cut photon. In doing so, we arrive at where we defined . Observing that there exist two diagrams for f 1 → π + π − π + π − due to the indistinguishability of the two π + and π − -see Fig. 4-we use the ρππ coupling as prescribed by Eq. (B.8) to deduce Here, the momenta are defined as in Fig. 4 and the pions are on shell, p 2 1/2 = M 2 π = k 2 1/2 . Given this amplitude, one can calculate the decay width and thus branching ratio via the four-body phase-space integration of We use the differential four-body phase space dΦ 4 (P, p 1 , p 2 , k 1 , k 2 ) in the form [86] dΦ 4 (P, p 1 , where dΦ 2 (P ; q 1 , q 2 ), dΦ 2 (q 1 ; p 1 , p 2 ), and dΦ 2 (q 2 ; k 1 , k 2 ) are the respective two-body phase spaces of the subsystems Since the integration volumes of the phase spaces are Lorentz invariant, each two-body phase space can be evaluated in the corresponding center-of-mass frame and we have to perform an explicit Lorentz transformation from the center-of-mass frames of {π + (p 1 )π − (p 2 )} and {π + (k 1 )π − (k 2 )} into the one of {ρ(q 1 )ρ(q 2 )} in order to evaluate scalar products of the kind (p i · k j ), i, j ∈ {1, 2}, appearing in |M(f 1 → π + π − π + π − )| 2 -see, e.g., Ref. [123] for more details. 8 We perform the phase space integration numerically with the Cuhre algorithm from the Cuba library [124], where the energy-dependent width Γ ρ (q 2 ) is as specified in Eq. (3.15), and obtain [125] Γ(f 1 → π + π − π + π − ) = |C s | 2 |g ργ | 4 |g ρππ | 4 × 0.63 × 10 −10 GeV. Eq. (B.11), we find the branching ratio to be given by The comparison with the experimental ratio B(f 1 → π + π − π + π − ) = 10.9(6)% [86] yields 20) in serious disagreement with Eq. (4.7). Including ρ contributions within the minimal VMD representation, there are four additional diagrams as compared to Fig. 4 and the corresponding master formula takes the form where see Table 4 for the numerical values of the decay rates. The numerical pattern shows that even though the coupling κ itself is O(1), ρ contributions are significantly suppressed, both due to the propagators in Eq. (4.15) and because the ρ can never be on shell in the available phase space. For the solutions of the global phenomenological analysis in Sec. 6, we find that the interference effects tend to even slightly reduce the branching ratio in the minimal VMD case, while the large values of (1 − 1 − 2 ) in the extended VMD fits can increase B(f 1 → π + π − π + π − ) to the level of 1%, still far below the experimental value. The reason for this incompatibility can be understood as follows. The available phase space prohibits the two ρ mesons from being simultaneously on-shell, and the corresponding loss of resonance enhancement for two intermediate ρ mesons implies that other decay mechanisms become more important. A candidate for such a mechanism is given by the decay f 1 → a 1 π → ρππ → 4π, see App. D for an estimate of this decay channel. From this analysis, we indeed infer that the intermediate state a 1 π likely saturates the decay width to a large extent, so that we have to conclude that the decay f 1 → 4π does not allow one to extract further information on the f 1 TFFs. We will thus disregard this input entirely and adopt Eq. (4.7) for the symmetric normalization. With C a 1 , C a 2 , and C s all real couplings, we will further fix the global sign by demanding that C s be positive, C s = 0.93 (11).

f 1 → ργ
The construction of the amplitude for f 1 → ργ proceeds along the same lines as for f 1 → 4π, via M(f 1 → γ * γ * ), either by using the minimal or the extended VMD parameterization. By definition, this decay channel only probes the isovector contribution, up to negligible isospin-breaking effects.
For the amplitude M(f 1 → ργ), we then proceed as stated above, starting with the minimal VMD ansatz, and consider the ρ meson and photon on shell, q 2 1 = M 2 ρ , q 2 2 = 0, and * (q 1 ) · q 1 = 0 = * (q 2 ) · q 2 , which also implies Γ ρ (q 2 2 = 0) = 0 = Γ ρ (q 2 2 = 0) according to Eq. (3.11)-Eq. (3.15). The corresponding diagram is depicted in Fig. 5 and we find where we introduced C f ργ = eM 2 ρ /( g ργ m 2 f 1 ). The branching ratio of the decay is given by where-as throughout this work-the coupling constants are assumed to be purely real and we defined the coefficients As depicted in Fig. 6, the solution of Eq. (4.26) in terms of the unknown couplings C a 1 and C a 2 represents an ellipse, where we used the central values of C s = 0.93 (11) and B(f 1 → ργ) = 4.2(1.0)%, see Sec. 6, to illustrate the cut surfaces. Although it is straightforward to actually solve Eq. (4.26) for such an equation, we refrain from doing so here since there is no unique solution anyway without further input. The equivalent amplitude in the extended VMD representation reads the only difference compared to the minimal VMD parameterization being that C s → C s = (1 − 1 /2 − 2 )C s . Hence, the branching ratio given in Eq. (4.26) becomes which, when inserting 1 and 2 from Sec. 3.4, simplifies to where we defined the coefficients In this variant, the dependence on C a 2 thus disappears from the branching fraction, which is a subtle consequence of the correlation between C a 2 and C s imposed via the singly-virtual high-energy behavior, see Eq. (3.35).
Another measured quantity of interest with regard to f 1 → ργ is the ratio of the ρmeson's helicity amplitudes in its rest frame, which is accessible through the subsequent decay ρ → π + π − . In a similar manner to how we obtained the f 1 → ργ amplitudes in Eq. (4.25) and Eq. (4.28), we can construct an amplitude for f 1 → ργ → π + π − γ, where we indeed consider the subsequent decay of an on-shell ρ meson and furthermore use the ρππ coupling given by Eq. (B.8); the process is depicted in Fig. 7.
In the extended VMD case, one again needs to replace C s → C s = (1 − 1 /2 − 2 )C s , which then further simplifies to when inserting 1 and 2 from Sec. 3.4. The coupling C a 2 therefore does not contribute to either f 1 → ργ observable in the extended VMD ansatz. The solution of Eq. (4.34) in terms of the unknown couplings C a 1 and C a 2 is given by four unconnected straight lines, as apparent from Fig. 8, where we used the central values 4.4 f 1 → φγ and f 1 → ωγ The branching ratio of f 1 → φγ has been measured experimentally, B(f 1 → φγ) = 0.74(26) × 10 −3 [86,88], and thus allows for another consistency check of our VMD representations, in particular, the U (3) assumptions for the isoscalar TFFs. Similarly, we can predict the branching fraction for f 1 → ωγ once all the parameters are determined, which could be confronted with potential future measurements. In complete analogy to Sec. 4.3, we construct amplitudes for f 1 → V γ, V = φ, ω, i.e., where we defined C f V γ = eM 2 V /( g V γ m 2 f 1 ). In terms of the ratio R V = R φ , R ω of isoscalar to isovector couplings, Eq. (3.27), the branching ratio of the decay is given by The generalization to the extended VMD representation would be straightforward, once applied to the isoscalar sector.
As the discussion in Sec. 4 shows, in general the constraints from e + e − → e + e − f 1 , f 1 → 4π, and f 1 → ργ do not suffice to reliably determine all three free VMD parameters, with the branching fraction of f 1 → 4π not able to provide any additional input at all due to significant contamination from decay channels not related to the TFFs. In this way, the evidence for the decay f 1 → e + e − reported by the SND collaboration [89] is extremely interesting as future improved measurements of the decay have the potential to overconstrain the system of C a 1 , C a 2 , and C s , as we will demonstrate in Sec. 6. In this section, we provide the required formalism to extract information on the f 1 TFFs from its decay into e + e − ; cf. also Ref. [100]. The Feynman diagram for the one-loop process is depicted in Fig. 9. The general form of the amplitude is which implies for the spin-averaged amplitude squared and a decay width of Here and in the following, the arguments of the reduced amplitude A 1 will be suppressed and we will work in the limit m e = 0. To extract A 1 from the full amplitude, we first consider the amplitude M(f 1 → γ * γ * ) and recast it into the more convenient form Inserting this amplitude into the QED loop, the full amplitude can be written as where we have used the on-shell condition for the fermions, neglected their masses, and written the loop integration in the most symmetric way. In particular, rewriting the TFF combinations as (q 2 2 − q 2 1 )F a 2 (q 2 1 , q 2 2 ) − (q 2 2 + q 2 1 )F s (q 2 1 , q 2 2 ) = −2q 2 1 F 2 (q 2 1 , q 2 2 ) + 2q 2 2 F 3 (q 2 1 , q 2 2 ), 2q 2 1 q 2 2 + k 2 (q 2 1 + q 2 2 ) F s (q 2 1 , q 2 2 ) + k 2 q 2 1 − q 2 2 F a 2 (q 2 1 , q 2 2 ) = 2(k 2 + q 2 2 )q 2 1 F 2 (q 2 1 , q 2 2 ) − 2(k 2 + q 2 1 )q 2 2 F 3 (q 2 1 , q 2 2 ) (5.6) shows that the BL limits that are not well-defined-see Eq. (2.18) and the subsequent comment-always appear suppressed by the respective on-shell virtuality, as expected from the form of the physical helicity amplitudes. We conclude that these integration regions will therefore be of minor importance. Moreover, all remaining integrals are ultraviolet and infrared convergent by inspection of the parameterization of the form factors in Eq. (3.10) and Eq. (3.20). However, inserting the (isovector) VMD expressions directly into the loop integral would produce unphysical imaginary parts, which can be avoided by using the spectral representations of Eq. (3.16) and Eq. (3.21) instead, to ensure the correct analytic properties.
For the numerical analysis we further write the coefficients in Eq. (5.7) according to where the prefactor for D asym is motivated from Eq. (3.32) to ensure that the resulting dimensionless coefficients can be compared in a meaningful way. Our numerical results are  shown in Table 5, including the uncertainties from the variation in Γ ρ . Even after taking the change in the normalizations into account, see Table 1, these results show that the uncertainties due to the spectral shape and the width itself can lead to comparable effects.
To be able to better compare the various contributions, we also show the coefficients including their normalizations, see Table 6, where we used the value of Eq. (3.33) for the asymptotic contribution. These numbers show that the symmetric contribution still produces the largest coefficient, but not by much. Accordingly, the f 1 → e + e − decay proves sensitive to the antisymmetric TFFs, about which not much is known at present. For the extended VMD ansatz, this observations implies an important caveat regarding the numbers shown in the table, which have been produced under the assumption that C a 2 = 0. In this case, one observes distinct differences between the two VMD versions, which can be traced back to the different weight given to the ρρ contribution. Finally, the real part of the isoscalar coefficient comes out larger than expected from Eq. (3.9). This is due to the fact that the loop integral is effectively regularized by the vector-meson mass, and the masses of ω and φ differ by a sufficient amount that the cancellation in Eq. (3.28) between the two contributions becomes less effective. The imaginary part of the loop integral is finite also in the infinite-mass limit, so that its size complies better with the expected isoscalar suppression.
Since the coupling constants are real, we use the decay width from Eq. (5.3) to obtain a branching ratio of where we defined and the terms involving D asym are only included in the extended VMD representation. Similarly to Eq. (4.26), the solution of Eq. (5.13) in terms of the unknown couplings C a 1 and C a 2 represents an ellipse in the minimal VMD case, which, however, changes for the extended VMD representation, see Fig. 10. Here, we used the central value of C s = 0.93 (11), Eq. (4.24), to remove one unknown and set √ s m = 1.3 GeV for the asymptotic contribution [26,27]. In fact, the results in Table 5 and Table 6 show that D asym remains small for a wide range of matching scales s m , so that the details of the matching do not play Reference  a role in view of the present experimental uncertainties. For definiteness, we will continue to use √ s m = 1.3 GeV in the following, with the understanding that the matching can be refined once improved data become available, along the lines described in App. A. In order to solve for all couplings, we need to consider a combined analysis of all constraints, see Sec. 6. However, given that the biggest contribution tends to come from the symmetric term, see Table 5, it is instructive to study the case C a 1 = C a 2 = 0 and consider the f 1 → e + e − decay as an independent determination of C s . For the minimal VMD ansatz we find where the isoscalar contribution implies an increase by about 0.3 (1). The extended variant gives 9 C s = 1.9 +0.8 −0.6 , (5. 16) where the uncertainties from the dependence on the ρ spectral function, its width, and the asymptotic contribution, ∆C s 0.03, are negligible compared to both the experimental error and the uncertainty from the isoscalar contribution. Both values are larger than the L3 result given in Eq. (4.24), indicating that indeed a significant contribution from the antisymmetric TFFs should be expected, which in view of the results from Table 6 is well possible with plausible values of C a 1/2 . Finally, the difference between Eq. (5.15) and Eq. (5.16) gives a first estimate of the sensitivity to the chosen VMD ansatz.

Combined phenomenological analysis
In this section, we perform a global analysis of the experimental constraints from e + e − → e + e − f 1 , f 1 → ργ, and f 1 → e + e − . We will also consider f 1 → φγ due to its relation via U (3) symmetry, but not include f 1 → 4π for the reasons stated in Sec. 4.2 and App. D. Most of the input quantities follow in a straightforward way from the experimental references and the compilation in Ref. [86], see Table 7, except for the branching fraction of the ργ   Ref. [146] tends to further reduce the average a little, which together with a slightly increased scale factor when including Refs. [87,146] leads us to quote as our final average, which we will use in the subsequent analysis, see Table 7. While our main argument in favor of this procedure is the avoidance of a fit bias towards the larger ργ branching fractions, one may also compare to theoretical expectations. The models considered in Refs. [140,[147][148][149][150] in general do prefer smaller ργ branching fractions, but the spread among the models is too large to make that comparison conclusive. The results of the global analysis are shown in Table 8 and Table 9, restricted to the parameterization Γ (ωπ,ππ) ρ (y) due to the dominant experimental uncertainties. The latter are propagated as given in Table 7, except for B(f 1 → φγ), for which we use B(f 1 → φγ)/(R φ ) 2 = 3.0(1.6)% as data point in the minimization, including the uncertainty on R φ from Eq. (3.27). As a side result, Table 8 and Table 9 also contain predictions for   Table 9: Same as Table 8, but for the extended VMD case, including the resulting parameter 1 .
the branching fraction of the yet unmeasured decay f 1 → ωγ. The outcome in the four cases considered-minimal and extended VMD representations each with and without the constraint from B(f 1 → φγ)-is illustrated in Fig. 11 and Fig. 12. In all cases the parameter C s is by far best constrained, its value hardly changes compared to the L3 reference point given in Eq. (4.24), with a slight preference for a small upward shift. The main distinctions concern the couplings C a 1 and C a 2 , with qualitative differences between the two VMD scenarios. In each case, however, we find two sets of solutions, corresponding to a small negative value of C a 1 (Solution 1) or a sizable positive one (Solution 2), respectively, both of which are shown in the tables and figures. In most cases, Solution 1 is strongly preferred, the exception being the minimal VMD fit including B(f 1 → φγ), in which case Solution 2 displays a slightly better fit quality.
In the minimal VMD representation, all constraints are sensitive to C a 2 , but especially once including B(f 1 → φγ) there is significant tension among the different bands. In Solution 2, the region preferred by all constraints but B(f 1 → e + e − ), which thus dominate the fit, would imply a much smaller value of B(f 1 → e + e − ) than reported by SND [89], while Solution 1 is better in line with the SND result. An improved measurement of B(f 1 → e + e − ) could therefore differentiate between these scenarios. In addition, we compare the resulting relevant form factor combination to the L3 dipole fit-see Sec. 4.1-in Fig. 13. While some tension is expected due to the singly-virtual asymptotic behavior of F a 1 (q 2 1 , q 2 2 ), see Table 2, the resulting curves for Solution 2 start to depart from the L3 band already around Q = 0.5 GeV, which further disfavors this set of solutions.   Figure 11: Contours in the C a 1 -C a 2 plane for the best-fit value of C s in the minimal VMD representation: for Solution 1 (left) and Solution 2 (right), without (upper ) and including (lower ) the constraint from B(f 1 → φγ). The best-fit region is indicated by the ∆χ 2 = 1 ellipse (inflated by the scale factor).
In the extended VMD representation, the dependence on C a 2 disappears in all observables apart from B(f 1 → e + e − ) and, potentially, B(f 1 → φγ). Accordingly, in the fit without the latter, the value of C a 2 is solely determined by B(f 1 → e + e − ), and the best-fit value of this branching fraction thus coincides with the input. There is good consistency among the other constraints, as reflected by a reduced χ 2 around unity. In this case, an improved measurement of B(f 1 → e + e − ) could thus be interpreted as a determination of C a 2 . Once B(f 1 → φγ) is included, one obtains an additional constraint on C a 2 , which, however, needs to be treated with care. First, the uncertainties on R φ have been included in the fit, but in addition there are U (3) uncertainties that are difficult to quantify. Moreover, the isoscalar contributions have been treated in their minimal variant throughout, but if excited ω and φ states were included, the dependence on C a 2 would again change, even disappear in a scenario similar to the extended VMD representation for the isovector contributions. Since the fit including B(f 1 → φγ) favors a value of B(f 1 → e + e − ) smaller than SND (for Solution 1 similar in size to the ones for Solution 1 in the minimal VMD  Figure 12: Contours in the C a 1 -C a 2 plane for the best-fit value of C s in the extended VMD representation: for Solution 1 (left) and Solution 2 (right), without (upper ) and including (lower ) the constraint from B(f 1 → φγ). The best-fit region is indicated by the ∆χ 2 = 1 ellipse (inflated by the scale factor). We do not consider equivalent solutions with very large negative C a 2 , as arise without the B(f 1 → φγ) constraint. Further local minima when including B(f 1 → φγ) mirror the indicated Solutions 1 and 2 on the lower branch of the ellipse, but display a worse χ 2 /dof and are thus discarded. case), an improved measurement of B(f 1 → e + e − ) would also allow one to differentiate between these scenarios. In addition to the worse χ 2 , Solution 2 is again disfavored by the comparison to L3, see Fig. 13.
In contrast, for Solution 1 of both the minimal and the extended VMD fit departures from the L3 dipole only arise around Q = 1 GeV, which implies agreement with all but the last data point of Ref. [7] (centered around Q = 1.8 GeV, where the curves still agree within uncertainties). In fact, a large part of the pull is a result of the slightly increased value of C s from the global fit, while the impact of the asymptotic behavior of F a 1 (−Q 2 , 0) remains small. Finally, we observe that most extended VMD fits require a substantial ρ contribution, as reflected by the large values of 1 shown in Table 9. In fact, for the fit without B(f 1 → φγ) it even exceeds the coefficient of the ρ contribution, which could With f 1 → φγ (Solution 1) With f 1 → φγ (Solution 2) Figure 13: Comparison of the fit solutions for the effective form factor probed in e + e − → e + e − f 1 to the L3 measurement [7], according to Eq. (4.11), for the minimal VMD representation (left) and the extended one (right). The L3 dipole band includes the uncertainties on |F D (0, 0)| and Λ D as given by Eq. (4.2), added in quadrature; ours propagate the uncertainties from Table 8 and Table 9, respectively.
be considered an indication that smaller values of B(f 1 → e + e − ) are preferred. We also implemented a variant of the extended VMD fit in which 1 was allowed to float freely, but this did not improve the fit quality, with a resulting 1 consistent with the ones imposed via Eq. (3.35).

Summary and outlook
In this paper, we performed a comprehensive analysis of the TFFs of the axial-vector resonance f 1 (1285), motivated by its contribution to HLbL scattering in the anomalous magnetic moment of the muon. Our study is based on all available constraints from e + e − → e + e − f 1 , f 1 → 4π, f 1 → ργ, f 1 → φγ, and f 1 → e + e − , all of which are sensitive to different aspects of the f 1 → γ * γ * transition. Since the amount of data is limited, a completely model-independent determination of all three TFFs is not feasible at present, leading us to consider parameterizations motivated by vector meson dominance. To assess the sensitivity to the chosen parameterization, we constructed two variants, a minimal one that produces non-vanishing results for all TFFs, and an extension that improves the asymptotic behavior by matching to short-distance constraints. In each case this leaves three coupling constants as free parameters, C s , C a 1 , and C a 2 , for the symmetric and the two antisymmetric TFFs, in terms of which the analysis is set up. As a first step, we derived master formulae for all processes in terms of these couplings and performed cross checks when analyzing each process in terms of the dominant coupling C s . This reveals that the decay f 1 → 4π does not provide further information on the TFFs, as the mechanism f 1 → a 1 π → ρππ → 4π likely dominates with respect to f 1 → ρρ → 4π, and only the latter can be related to the f 1 TFFs. The process is thus discarded in the subsequent analysis. For the remaining observables we performed detailed uncertainty estimates, including the subleading isoscalar contributions, the properties of the ρ meson and its spectral function, and the matching to short-distance constraints. In all cases we conclude that the dominant uncertainties are currently of experimental origin. Combining all constraints in a global fit, we found that the symmetric coupling C s is by far best determined, with substantial differences in C a 1 and C a 2 among the different scenarios, see Table 8, Table 9, Fig. 11, and Fig. 12 for our central results. Out of two sets of solutions-Solution 1 with a small negative value of C a 1 , Solution 2 with a sizable positive one-the former is in general preferred by the fit, with Solution 2 further disfavored by the comparison to space-like e + e − → e + e − f 1 data, see Fig. 13. In the case of the minimal VMD representation, we observed some tension between B(f 1 → e + e − ) and the remaining constraints especially when including B(f 1 → φγ) in the fit, leading to a preference for a branching fraction below the value recently reported by the SND collaboration. In the extended parameterization, the dependence on C a 2 drops out in all observables but B(f 1 → e + e − ) and, potentially, B(f 1 → φγ), but limited information about the isoscalar sector together with necessary U (3) assumptions render the latter constraint less reliable. While the f 1 → φγ branching fraction seems to prefer a smaller value of B(f 1 → e + e − ) (similar to the minimal VMD fit), we conclude that the parameter that controls its size, C a 2 , is largely unconstrained at the moment, and would thus profit most from an improved measurement of B(f 1 → e + e − ).
In general, new measurements of B(f 1 → e + e − )-as possible in the context of e + e − → hadrons energy scans at SND and CMD-3-would be highly beneficial to further constrain the f 1 TFFs, given that the resulting constraints are complementary to other observables, in particular, providing sensitivity to doubly-virtual kinematics and the antisymmetric TFFs. Apart from a more reliable determination of C a 2 , one could also validate and, if necessary, refine the underlying VMD assumptions. Furthermore, improved measurements of e + e − → e + e − f 1 would be valuable to further constrain the singly-virtual TFFs-in particular, the asymptotic behavior of F a 1 (q 2 , 0)-ideally adding new data points above 1 GeV and being analyzed using the full momentum dependence given in Eq. (4.11), to avoid the corresponding limitation in the interpretation of the L3 data. Such analyses are possible at BESIII [151] and Belle II [152]. To go beyond VMD parameterizations, the energy dependence in the (dispersively improved) Breit-Wigner propagators would need to be constrained by data, which would require differential information on f 1 decays. At the moment, our analysis summarizes the combined information on the f 1 TFFs that can be extracted from the available data in terms of simple parameterizations, which we expect to become valuable for forthcoming estimates of the axial-vector contributions to HLbL scattering.

Acknowledgments
We thank Gilberto Colangelo, Stefan Leupold, and Peter Stoffer for valuable discussions and comments on the manuscript. Financial support by the DFG through the funds provided to the Sino-German Collaborative Research Center TRR110 "Symmetries and the Emergence of Structure in QCD" (DFG Project-ID 196253076 -TRR 110) and the SNSF (Project No. PCEFP2_181117) is gratefully acknowledged.

A Asymptotic behavior including mass effects
In this appendix, we generalize the considerations of Refs. [26,27] regarding a doublespectral representation of BL scaling to include mass effects that arise from the kinematic variables in the denominator. Starting from see Eq. (2.12), we see that the asymptotic behavior of the axial-vector TFFs can still be deduced from the simpler pseudoscalar case, which leads us to study the generic integral which, in the case q 2 1 = q 2 2 = −Q 2 , evaluates tõ Given the large masses of the axial-vector mesons, m = m A , such corrections in ξ may become relevant and Eq. (A.2) defines a convenient test case to study their impact. As a first step, we observe that Eq. (A.2) can still be formulated as a single dispersion relation [121] via the transformation x = − 1−u u q 2 2 − m 2 u , where y = q 2 2 has been assumed to be space-like. Analytic continuation in q 2 2 then allows one to rewrite the imaginary part in Eq. (A.4) in terms of another dispersion relation, leading toĨ Restricting the integration in x, y should then allow one to isolate the asymptotic contributions while keeping the leading mass corrections. In the subtracted version, the singlyvirtual limits become explicit since Further, to make connection with the massless limit of Eq. (3.44), which amounts tõ see Eq. (3.42), we first note that this variant had been constructed in such a way that the singly-virtual contributions are removed, suggesting a matching in the limit q 2 (A.9) To evaluate Eq. (A.5) in the same limit, we symmetrize the integration to v = x + y, w = x − y and introduce a step function θ(v − v m ). In these variables, the w integration extends between w ± = ± √ 2m 2 v − m 4 , which shows that in the massless limit the doublespectral density indeed collapses to a δ function, see Eq. (3.43). For q 2 1 = q 2 2 = −Q 2 , the w integration can be performed analytically, leading tõ The first three terms in the expansion thus match upon the identification v m = 4s m .

B Phenomenological Lagrangians
In this appendix, we define the Lagrangians used for the ργ, ρππ, and ρωπ couplings and discuss the information that can be extracted for their ρ analogs. In particular, we derive estimates for the branching ratios B(ρ → ππ) and B(ρ → ωπ), which are necessary inputs for the construction of the energy-dependent width Γ (ωπ,ππ) ρ (q 2 ) in Eq. (3.13). For the coupling of photons to the vector mesons {ρ, ω, φ, ρ , . . .}, we use the effective interaction Lagrangian [115] µ }, {ω µν , ω µ }, and {φ µν , φ µ } are the respective vector meson equivalents, and the ellipsis refers to excited isoscalar vector mesons that we omit from the following discussion for simplicity. The couplings of the three ground-state vector mesons are linked via SU (3) symmetry according to g ργ : g ωγ : g φγ = 1 : 3 : 3/ √ 2 [115], with the sign of g φγ adjusted according to Eq. (3.24). In the following, we neglect complex phases associated with actual pole residues (which are known to be tiny [153]), and work with the phase convention sgn g ργ = +1. From the Lagrangian, the partial decay width of the vector mesons into e + e − follows as For the ρ meson, one can solve for the coupling and insert the (experimental) value Γ(ρ → e + e − ) = 7.04 keV [86] to find g ργ = 4.96.

(B.3)
This value agrees well with the residue |g ργ | = 4.9(1) extracted from the pion vector form factor [153], and is also close to the expectation from SU (3) symmetry, g For the VMD application considered in this work, we also need a formulation in which the coupling of photons to vector mesons is momentum independent, with the respective vector meson considered on shell. Such a coupling can be formally defined via the Lagrangian L Vγ = eA µ g ργ ρ µ + g ωγ ω µ + g φγ φ µ + g ρ γ ρ µ + . . . , where matching the amplitudes resulting from Eq. (B.1) and Eq. (B.6) for on-shell mesons determines In particular, we carry over the sign convention for the coupling constants g Vγ from g Vγ above.
Analyses of the pion vector form factor using improved variants of Eq. (B.10) suggest a ρ contribution relative to the dominant ρ therein of an approximate strength [91,155,156] g ρ ππ /g ρ γ g ρππ /g ργ ≈ − 1 10 . (B.14) On the other hand, the ω → πγ * TFF [91,157] can be approximated in a VMD picture according to .

(B.15)
The asymptotic behavior f ωπ (q 2 ) = O(q −4 ) [83,84,158,159] implies a superconvergence sum-rule constraint on the couplings of Eq. (B.15) according to which is consistent with the experimental analysis of Ref. [157]. From the experimental width Γ(ω → πγ) = 0.71 MeV [86] and the corresponding formula [91] we furthermore obtain the normalization |f ωπ (0)| = 2. so that under the assumption Γ ρ ≈ Γ(ρ → ππ) + Γ(ρ → ωπ)-neglecting another significant contribution from ρ → a 1 π (a 1 → 3π) 13  The branching ratios then become and, for completeness, the ρ γ coupling is estimated as The estimate Eq. (B.21) agrees with the expectation that the ρ should be largely inelastic, and the resulting spectral function in Eq. (3.13) thus essentially defines an estimate of the 4π channel dominated by ωπ. We stress that these considerations should only be considered rough estimates, the main point being to define another plausible variant that allows us to assess the sensitivity of our results to the assumptions made for the ρ spectral function. Finally, for our analysis of f 1 → 4π including effects of the ρ , we require the ratio of coupling constants g ρ ππ × g ρ γ g ρππ × g ργ ≈ −0.7. (B.23)

C Comparison to the literature
In this appendix, we briefly compare the basis of Lorentz structures and TFFs as well as the parameterization of the latter for the f 1 used in this work to the previous analysis of Refs. [100,101]. Since the TFFs are not (anti-)symmetrized in Ref. [100], we use the basis introduced in Sec. 2 for our comparison, that is, in particular, the structures from Eq. (2.6). When using Eq. (2.2) to translate the amplitude M(f 1 → γ * γ * ) from Ref. [100] to the tensor matrix element given in Eq. (2.5), we find the structures to be related by and the TFFs to be linked via While the structures are thus identical to ours except for two global signs and a permutation, the additional factor of 4π in the TFFs appears due to the fact that the fine-structure constant α is used in the definition of their matrix element instead of the factor e 2 . The symmetry properties of the TFFs in their basis are given by F The strategy that is used in Ref. [100] to determine the explicit parameterization of the TFFs in accord with a VMD model is, in fact, quite different from our approachthe model does not correspond to a strict VMD ansatz. Instead of proposing a VMD-like parameterization for the form factors F [100] i (q 2 1 , q 2 2 ) as we did in Eq. (3.10), three form factors h i (q 2 1 , q 2 2 ) are introduced, based on which an amplitude M(f 1 → ρ 0 * ρ 0 * ) is constructed by replacing F ; analogously, two complex coupling constants g 1 and g 2 are introduced to construct an amplitude M(f 1 → ργ). We disagree that such complex couplings are allowed since the resulting imaginary parts need to reflect the actual analytic structure of the amplitude. Moreover, the explicit form of the h i (q 2 1 , q 2 2 ), to account for an off-shell dependence of the ρ mesons, introduces unphysical kinematic singularities.
By employing a ργ coupling similar to the one we introduced by means of Eq. (B.1), the form factors F [100] i (q 2 1 , q 2 2 ) and h i (q 2 1 , q 2 2 ) are then related to each other in Ref. [100], where the latter can further be linked to the coupling constants g 1 and g 2 . Using the ργ coupling in the convention of the present work, the form factors are found to be , the width Γ ρ being the (energy-independent) total width of the ρ meson, as opposed to our energy-dependent parameterization of Eq. (3.11) and Eq. (3.15). Moreover, the magnitude of the couplings g 1 and g 2 is determined in Ref. [100] by making use of experimental data on f 1 → ργ, see Sec. 3; the relative phase between these coupling constants remains undetermined, despite using, in addition, input from f 1 → 4π. By rewriting Eq. (C.4) as one observes that F (q 2 1 , q 2 2 ) does not correspond to a VMD ansatz in the strict sense, but rather arises from two diagrams, each being composed of one direct photon coupling and one VMD-like ρ coupling. As we argued in Sec. 3, an actual VMD representation of the antisymmetric TFFs requires the introduction of a second multiplet. Further, Eq. (C.4) shows that the second and third TFFs are parameterized symmetrically, i.e., the antisymmetric part is neglected. In either case, we believe that the f 1 → 4π decay does not allow one to extract information on the f 1 TFFs, for the reasons described in Sec. 4.2 and App. D.
Finally, we would like so stress that, in addition to using complex couplings, energyindependent widths are problematic when inserted into the f 1 → e + e − loop integral, leading to imaginary parts below the respective thresholds and thus distorting the analytic structure. Given, in addition, the appearance of kinematic singularities and different highenergy behavior, it is difficult to compare our phenomenological results to the ones of Refs. [100,101].

E Constants and parameters
In this appendix, we collect the particle masses and decay widths used throughout this work, see Table 10. Isospin-breaking effects can be safely neglected, in particular, the pion mass is identified with the mass of the charged pion. Some comments are in order, however, regarding the treatment of broad resonances, most notably the ρ(1450) and, to  Table 10: Selected masses and decay widths from Ref. [86], in comparison to the ρ(770) and ρ(1450) parameters from Refs. [155,164].
a lesser extent, the ρ(770). Especially for the former, the quoted masses and widths are strongly reaction dependent, as referring to Breit-Wigner parameters, not to the modelindependent pole parameters. We thus need to make sure that we use determinations that apply to the channels that we consider here. Since the main application concerns the description of multi-pion decay channels in the VMD propagators, both for the ρ(770) and the ρ(1450), it appears most natural to consider reactions that provide access to both resonances, which points towards τ → ππν τ from Ref. [155] and e + e − → ππ from Ref. [164].
In particular, this allows us to see if there are relevant systematic differences between the charged and neutral channel. For the ρ(770), the mass parameter agrees well between all channels, but while there is also good agreement between Refs. [155,164] for the width, the compilation from Ref. [86] quotes a significantly lower value for the neutral channel.