How many vector charmonium(-like) states sit in the energy range from 4 . 2 to 4 . 35 GeV?

In recent years many vector charmonium(-like) states were reported by different electron-positron collider experiments above 4 . 2 GeV. However, so far, there not only exists sizable tension in the parameters of those states, but there is also no consensus on the number of the vector states in this energy range. To some extend, this might be caused by the fact that the experimental data were typically analyzed in single channel analyses employing overlapping Breit-Wigner functions, in particular ignoring the effect of opening thresholds. In this study, we focus on the mass range between 4 . 2 and 4 . 35 GeV, conducting a comprehensive analysis of eight different final states in e + e − annihilation. Our findings demonstrate that, within this mass range, a single vector charmonium-like state, exhibiting properties consistent with a D 1 ¯ D molecular structure and characterized by a pole location (cid:113) s Y (4230) pole = (cid:0) 4227 ± 3 − i 2 (50 +10 − 5 ) (cid:1) MeV, can effectively describe all the collected data. This is made possible by allowing for an interference with the well-established vector chamonium ψ (4160) along with the inclusion of the D 1 ¯ D threshold effect. Moreover, in contrast to experimental analyses, our study reveals that the highly asymmetric total cross sections for e + e − → J/ψππ and e + e − → J/ψK ¯ K around 4230 MeV stem from the same physics, rooted in the approximate SU(3) flavor symmetry of QCD.

In recent years many vector charmonium(-like) states were reported by different electronpositron collider experiments above 4.2 GeV.However, so far, there not only exists sizable tension in the parameters of those states, but there is also no consensus on the number of the vector states in this energy range.To some extend, this might be caused by the fact that the experimental data were typically analyzed in single channel analyses employing overlapping Breit-Wigner functions, in particular ignoring the effect of opening thresholds.In this study, we focus on the mass range between 4.2 and 4. 35 GeV, conducting a comprehensive analysis of eight different final states in e + e − annihilation.Our findings demonstrate that, within this mass range, a single vector charmonium-like state, exhibiting properties consistent with a D 1 D molecular structure and characterized by a pole location s Y (4230) pole = 4227±3− i 2 (50 +10 −5 ) MeV, can effectively describe all the collected data.This is made possible by allowing for an interference with the well-established vector chamonium ψ(4160) along with the inclusion of the D 1 D threshold effect.Moreover, in contrast to experimental analyses, our study reveals that the highly asymmetric total cross sections for e + e − → J/ψππ and e + e − → J/ψK K around 4230 MeV stem from the same physics, rooted in the approximate SU (3) flavor symmetry of QCD.

I. INTRODUCTION
Since the discovery of the first exotic state, i.e. χ c1 (3872) also known as X(3872), in the cc-sector in 2003, a large number of states was discovered in the charmonium and bottomonium mass range that show properties incompatible with expectations from the original quark model.For recent reviews see, e.g., Refs.[1][2][3][4][5][6].The amount of available data is especially rich in the J P C = 1 −− channel, since here states containing cc can be generated directly in e + e − -collisions and can therefore straightforwardly be studied at experiments like BaBar, Belle and BESIII.In this work, we focus on vector states in the mass range from 4.2 to 4. 35 GeV.This energy range hosts most prominently the ψ(4230) also known as Y (4230) and potentially one additional state located at 4.32 GeV.The latter was introduced in the analyses of the BESIII collaboration for the reaction e + e − → J/ψπ + π − to account for the highly asymmetric line shape seen in the experiments reported in Refs.[7,8].In particular, the most recent analysis [8] revealed for the Breit-Wigner mass and width of the Y (4230) in this channel M Y (4230) = 4221.4± 1.5 ± 2.0 MeV Γ Y (4230) = 41.8 ± 2.9 ± 2.7 MeV , and for the Y (4320) where the first and second uncertainty is statistical and systematic, respectively.The Y (4320) is also needed for analyzing the J/ψπ 0 π 0 channel [9], and the parameters above are consistent with the data in this channel.On the other hand, the Y (4230) is seen in the eight additional channels shown in Fig. 1, admittedly with largely inconsistent parameters, while the state dubbed Y (4320) shows up in none of them, at least within the mass range consistent with Eq. ( 2), not even in e + e − → J/ψK K which is connected to e + e − → J/ψππ by the approximate SU(3) flavor symmetry of QCD.In experiments by BaBar and Belle a state named Y (4360), with a mass of about 4345 MeV, was discovered in the ψ(2S)π + π − [10,11] final state.However, the recent BESIII measurement of the same channel revealed that the Y (4360) emerges due to a subtle interference of the Y (4230) and a state at 4390 MeV with a width of 140 MeV [12], which is thus in a mass range close to the ψ(4415), however, twice as wide.
A signal at 4390 MeV with the consistent parameters was also observed by BESIII in the h c ππ [13] and J/Ψη [14] final states.Since this state is outside the mass range in focus here, we do not discuss it any further.
In the mass range from 4.2 to 4. 35 MeV the findings just described raise the following questions: 1. Why does the observed width of the Y (4230) deduced from the J/ψππ channel, differ so significantly from that deduced from the D * Dπ channel, where the measured width is twice as large [15]?
2. What can we learn from the cross section differences for Y (4230) in its various decay channels?Note that the cross section in the D D * π channel is about one orderof-magnitude larger than those of hidden charm decays.
3. Why is the Y (4230) observed in final states with both cc spin 1 (i.e.J/ψππ and ψ ′ ππ channels) and cc spin 0 (i.e.h c ππ channel) at a similar rate, despite being produced via a photon, which leads to cc in spin 1 only?Can we understand this seemingly large violation of heavy quark spin symmetry?
4. Why is the Y (4320) seen only in a single channel?
5. Can the apparent asymmetry of the J/ψπ + π − line shape be generated by the opening of the D 1 (2420) D channel just below the nominal mass of the Y (4320)?
Here the D 1 (2420) is the narrow axial vector with a width of about 30 MeV, which decays to the πD * channel predominantly in D-wave; the nearby D 1 (2430) has a width of about 300 MeV and decays to the πD * channel predominantly in S-wave.Thus, the broad D 1 is not capable of producing structures as narrow as those discussed here, although its mixing with the narrow D 1 , emerging from spin symmetry violation, is relevant for the detailed description of the data.
In this work we address the mentioned issues starting from the assumption that the Y (4230) is a D 1 (2420) D hadronic molecule, proposed originally in Ref. [16], and refined in Ref. [17] by taking into account the the D 1 D cut properly and in particular the triangle singularity mechanism, which is crucial for the production of the Z c (3900), accounts for the lineshape of the J/ψππ and in this way leads to a pole position around 4. 23 GeV.Before we discuss in detail the observable consequences of this assumption we present the other structure assumptions put forward for this state in the literature, namely studies that do not need the pole of the Y (4230) at all as well as the three types that do call for a state in this mass range, namely the hybrid, the hadrocharmonium and the compact tetraquark interpretation -details are given in the following paragraphs.This allows us to demonstrate that the implications of a molecular structure for the Y (4230) are very specific and significant.
In Refs.[18,19] no pole for the Y (4230) needs to be introduced.In the former reference the structure near √ s = 4230 MeV is generated from the ψ(4160) coupling to the D s Ds channel.
While this provides a reasonable description of the J/ψππ final state, it is unlikely that the same scenario also allows a description of all the other final states, especially the D D * π channel.The same comment applies to Ref. [19], where the Y (4230) is generated from an interference of the neighbouring charmonium states.We therefore do not consider these mechanisms any further.
In the hadrocharmonium picture an exotic hidden charm state appears as a compact cc core surrounded by some typically excited light quark cloud [24] 1 .While this explains naturally that the Y (4230) decays into J/ψππ and not D ( * ) D( * ) as would be expected for a cc quark-model state, it appears at odds with the fact that the Y (4230) is observed also in the h c ππ final state, since heavy quark symmetry calls for a conservation of heavy quark spin.To overcome this problem it was proposed in Ref. [26] that the Y (4230) and the next higher state are in fact emerging from a mixing of two states, one with a spin 0 cc core and one with a spin 1 core.Thus, this scenario calls for a second nearby vector state -a currently good candidate being the above mentioned Y (4320).Moreover, this mixing scenario implies the existence of four spin symmetry partners [27].For example, there should be two exotic η c states, one in between the two vector states, one significantly lighter than the Y (4230).
Very early after its discovery, the Y (4230) was proposed to be a hybrid state based either on phenomenological calculations [28][29][30] or heavy quark effective field theory [31].In the hybrid picture, both quarks and gluons contribute as valence degrees of freedom.A study of the decays employing heavy quark effective field theory disfavors a pure hybrid interpretation of the Y (4230) [32].In any case, also the hybrid picture calls for a mixing of two close-by vector states with different spin of the cc component of the wave functions and thus for the existence of both Y (4230) and Y (4320) to accommodate the decays into final states with both spin 0 and spin 1 for the outgoing charmonium.In addition, for a hybrid vector the decays into J/ψππ and J/ψK K are connected by SU(3) flavor symmetry.The rate in the latter channel deduced in this way is however larger than what one finds experimentally.
In the compact tetraquark picture the states are typically made of heavy-light diquarks and anti-diquarks.This approach calls for four nonstrange vector states with masses in the range 4220 MeV and 4660 MeV [33,34], since the diquarks can have either spin one or spin zero allowing for the following spin couplings with positive C parity, [0, 0] 0 , [1,0] with the spins of diquark and antidiquark in the brackets and their total spin as subindex outside -note that a state that contains two spin 1 substructures coupled to total spin 1 has negative C parity.To get the negative parity needed for a vector state, an angular momentum of 1 needs to be introduced between the diquark and antidiquark that in addition flips the C parity to the needed -1.For example, the currently preferred fit of Ref. [33] includes both Y (4220) as well as Y (4320).An alternative approach to compact tetraquarks, similar in spirit, but different in the realization, is outlined in Ref. [35].Thus, we see that three of the non-molecular scenarios prefer the presence of both Y (4230) and Y (4320) while the remaining are challenged by the decay properties of the Y (4230).
In this study, we investigate the feasibility of a combined analysis involving eight different final states excited in e + e − annihilation, namely molecule.The main message of this work is that the data available in this mass range is consistent with the presence of a single exotic state predominantly of molecular nature, since such a state necessarily has a large coupling to the D 1 D channel.The molecular scenario for the Y (4230) was already advocated in Refs.[37,38] based on an analysis of older data in the J/ψππ and h c ππ channels.It is crucial to emphasize that, while certain properties of the data emerge naturally in the current analysis, there are cases where fine-tuned parameters are necessary.It turns out that in order to obtain a coherent picture, it is unavoidable to include the interference with an additional vector state whose properties we fix to those of the well known charmonium state ψ(4160).This is illustrated in Fig. 2, where we show the fit results with and without the ψ(4160): The narrow structure seen in the J/ψππ channel is by itself incompatible with the broad structure seen in D D * π, as well as some other channels discussed below.Ref. [39] puts forward the hypothesis that ψ(4160) and Y (4230) could actually be the same state.However, the data shown in Fig. 2 indicate that this conjecture is not compatible with the data.We regard this work as some exploratory study -accordingly some disclaimers need to be given, which will be overcome in subsequent publications: • We treat the effect of the interference of the ψ(4160) with the Y (4230) perturbatively.While this simplifies the fitting, it violates unitarity since only terms linear in the vector propagators are included in the evaluation of the hadronic cross sections.
• This is a phenomenological study.In particular, we cannot estimate uncertainties from a truncation error in some systematic expansion.This is appropriate, however, since we only aim at demonstrating what is possible with a single exotic particle in the mass range of interest.Accordingly, uncertainties of e.g. the pole parameters of the Y (4230) were only roughly estimated at this stage.
• We focus on the effect of the D 1 D intermediate state in the decays of the Y (4230), basically ignoring that heavy quark spin symmetry also calls for the coupled channels D 1 D * and the D 2 D * -this is the main limiting factor when considering the energy range.
• The channels with two pions or two kaons in the final states necessitate the proper inclusion of ππ/K K final-state interactions, as discussed in previous works [40][41][42][43][44].In this study we simplify the treatment of these effects.While our approximation shows qualitatively very reasonable results, the data for e + e − → ψ(2S)ππ, exhibit a very unusual energy dependence in the subsystem invariant mass distributions at √ s = 4230 and 4260 MeV, which seem to require a more refined treatment.Consequently, data from e + e − → ψ(2S)ππ are not included in the current fits.
• The data currently available do not show apparent peak structures of Y (4230) in D ( * ) D( * ) channels, which must appear in odd partial waves to reach J P C = 1 −− .This suggests that the couplings of Y (4230) to the two-body open charm channels are much smaller than those of the vector charmonium states.In Ref. [45] it was demonstrated that the dips seen in the data of e Note that with respect to exploiting the implications of the heavy quark spin symmetry there are more advanced studies than this one already published [46,47].However, both those works focus solely on the pole locations that emerge from solving the scattering equations for the members of the spin multiplet {D 1 , D 2 } scattering off those of { D, D * }.No attempt is made to investigate the resulting line shapes in the various decay channels.Contrary to those works, we here study the energy dependence of the cross sections in the various decay channels.This allows us to demonstrate that the inclusion of the ψ(4160) together with the strong coupling of the Y (4230) to the D 1 D channel that is a consequence of its assumed molecular nature, is sufficient to describe all data sets studied here without the need for an additional exotic state in the mass range of interest.This conclusion is in line with the preliminary results of this study announced in Ref. [48].
While preparing this manuscript we got aware of Ref. [49], where the degrees of freedom included in this work are in as well.The central finding of this work relevant for us is that in total three poles are needed in the energy range studied here.While one of them might well represent the ψ(4160), only with somewhat shifted parameters, and another one the Y (4230), there is still a resonance needed close to 4320 MeV, absent in our analysis.A more detailed comparison with that work will only be possible, once more details are published.
An alternative analysis to ours that includes both ψ(4160) as well as Y (4230) in the energy range of interest here but does not call for a state located at 4320 MeV is Ref. [50].In this work the asymmetric shape observed in the total cross sections of ππJ/ψ and D * Dπ can be reproduced as an interference effect between ψ(4160) and the higher energy state ψ(4415)-the latter state is beyond the energy range considered in our analysis-combined with a non-resonant background.Then the inclusion of Y (4230) is needed for fine-tuning the agreement with data near 4.2 GeV.Another striking difference between that work and ours is their omission of any threshold effects.As we argue below, the significance of the D 1 D threshold in the Y line shapes is a direct hint towards its molecular nature -accordingly Ref. [50] argues that their analysis is consistent with a cc structure of the Y (4230).Thus studying observable differences between the results of that work and ours is important to pin down the nature of the Y (4230).We come back to this when discussing the results.It should also be mentioned that the relatively large cross section seen in e + e − → h c ππ -where the final state contains a cc pair in spin zero, contrary to the production of a cc pair with spin one -suggests a considerable amount of heavy quark spin symmetry (HQSS) violation.This phenomenon is naturally explained in a molecular scenario for the Y , whereas it seems highly unnatural for a cc structure.
The paper is structured as follows: We start with some general considerations about the diagrams to be included in the molecular approach.Then, in Sec.III, we describe in some detail the formalism employed.Sec.IV contains the fitting results as well as their discussion.We close with a summary and outlook in Sec.V. Additional technical details of the calculations are delegated to Appendices.

II. GENERAL CONSIDERATIONS
In the molecular scenario, the coupling of some physical Y states to the nearby continuum channel, h 1 h 2 , that forms the molecule is maximal [51,52], see also [4] for a review.In fact, this large coupling is what encodes the molecular nature of a given state.Accordingly, the transition of Y to the channel h 1 h 2 dominates over the others and the diagrams containing this coupling appear always at leading order.In effec-tive field theories (EFT), the relative importance of the diagrams is controlled by power counting rules, as presented for similar systems, e.g., in Refs.[53,54].However, our case involves several additional complexities such as the presence of the unstable particle D 1 in the transition, exploration of a relatively wide energy range, and the analysis of three-body final states.To illustrate the second point, we note that for the energies near the Y (4230) peak, triangle diagram (c) of Fig. 3 is potentially important, as long as we look at D * D invariant masses close to the mass of the Z c (3900).Indeed, not only the D * D intermediate state is in this case nearly on shell, but also the nearby triangle singularity significantly enhances the contribution of this diagram [55].However, this diagram is suppressed over a large fraction of the Dalitz plot apart from this range.Therefore, in what follows, we employ a more pragmatic strategy to consider the most natural and phenomenologically motivated production mechanisms and investigate their relative importance in different energy regimes.Due to this, we postpone an estimation of the theoretical uncertainties to a later publication.
A. e + e − → D 0 D * − π + The most direct access to a molecular state is provided by its imprint on the near threshold cross section of the channel that forms the molecular state, since as outlined above the coupling of the molecule to its constituents is large.The same reason is also the origin of the unnaturally large nucleon-nucleon scattering lengths [4,51,56].This phenomenon arises from the existence of nearby molecular states in both the spin-1 channel (with the deuteron as a true bound state) and the spin-0 channel (with a closely located virtual state).The propagator -after all we are near the pole of a narrow state.In our study we also treat the Z c (3900) as a hadronic molecule, in line with Ref. [17].Accordingly the coupling of the Z c to the D * D channel is also large.Moreover, the triangle diagram, which is part of diagram (c), is enhanced by a very close by triangle singularity [55].Thus we expect diagrams (a) and (c) to contribute significantly to the observables.
The pole of the Y (4230) is located about 60 MeV below the nominal D 1 D threshold, which is about twice the width of the D 1 (2420).This width still allows for a resonance signal even at √ s = 4230 MeV [57,58], however, with a significant kinematic suppression -the large coupling of the Y (4230) to the D 1 D channel with the D 1 → D * π decay in D-wave develops its effect mostly above the D 1 D threshold.However, although violating HQSS, the narrow D 1 (2420) is expected to mix with the much broader D 1 (2430).In this way the narrow D 1 (2420) also gets an S-wave decay [59], which does not so strongly suffer from the above mentioned kinematic suppression allowing it to contribute significantly to the Y (4230) peak in the πD * D channel.In particular, the decays of the D 1 (2420) to D * π both in S-and D-waves are therefore included in diagrams a) and c).In this paper, whenever referring to D 1 without mass number, we talk about the narrow D 1 (2420), to simplify notation.
In addition, because of the large width of the D 1 (2430), its residual effect acts effectively like a very short ranged contribution.We thus do not calculate loop contributions involving this broad state explicitly but parametrise it by a point coupling of the Y (4230) to πD * D with S-waves in all subsystems.Since the D * D in S-wave also undergoes final state interactions, the just men- tioned point coupling cannot occur in isolated form, but needs to get dressed by the Z c (3900) propagator that parametrises the D * D S-wave interaction.This results in an expression that is represented by diagram 3(b) -details for the expressions employed are given in Sec.III as well as in the Appendix.This construction is automatically consistent with the Watson theorem [60].
Finally, in the experimental data for the πD * D channel the width of the structure around 4230 MeV is notably broader than that observed, e.g., in the J/ψππ channel -see Fig. 2. Because of this, the parameters extracted for the Y (4230) in the two channels by the BESIII collaboration are inconsistent with each otherc.f.Fig. 1.A possible mechanism that allows for a combined fit of the various channels is that the ψ(4160) also has some small coupling to the D * Dπ.The experimental signal observed could then be interpreted as the result of an interference of the signatures from the two resonances.Also for the ψ(4160) we assume that the coupling is in S-wave with respect to all subsystems and, as before, also here the direct transition ψ(4160) → D * Dπ gets dressed by the D * D final state interaction parametrised by the Z c (3900) propagator.The corresponding diagram is shown in Fig. 3(d).We are aware that, if the ψ(4160) were (predominantly) a Dwave charmonium, there should also be angular momenta in the final state as a consequence of HQSS.However, the data do not call for an additional coupling structure and we thus omit it from our study.
B. e + e − → J/ψ(ππ/ KK) Next we turn to the discovery channel of the Y (4230), e + e − → J/ψππ, where the highly asymmetric lineshape lead to the claim for the existence of an additional state called Y (4320) [7,8].Again, driven by the assumed molecular nature of the Y (4230), contributions that run through the D 1 (2420) D intermediate state are sizable and need to be considered.Then, to reach the J/ψππ final state, possible topologies are either box diagrams (see Fig. 4,  (a) and (b) as well as Fig. 5 for a complete set of box diagrams) or a triangle followed by a Z c (3900) propagator (Fig. 4(c) and (d)).As before we need to allow for additional processes and also here include a diagram for the contact transition of the Y (4230) to the J/ψππ final state (Fig. 4(e)), as before dressed by the final state interaction that leads to the occurrence of the Z c (3900) -see Sec.III B for a detailed discussion.Furthermore also in this channel we allow for a contribution of the ψ(4160), shown in Fig. 4(f ).To come to the full amplitudes, the ππ final state interaction needs to be taken into account as well.Since the initial photon generates a ccpair, which is isoscalar, and the final cc pair is isoscalar as well, the pion pair must be isoscalar with even angular momentum (the latter also follows from parity conservation).In the vincinity of the Y (4230) pole, the ππ system is probed in the energy range from its threshold up to about 1.1 GeV.Since the scalar-isoscalar ππ interaction has a strong coupling to the KK system, the final state interaction is included by employing a formalism that explicitly treats the coupled channels.Since the full treatment of the system is technically very demanding [44] (see Refs. [40][41][42] for related studies) because of the intricate singularity structure of the pertinent integrals, in this exploratory study we employ an approximate treatment that still allows for a sensible description also of the ππ spectra -details are given in the next section and in Appendix B.
FIG. 5. Decomposition for the box topology of e + e − → J/ψπ + π − strangeness sources are also included in the calculation of the J/ψππ final state.
C. e + e − → h c ππ The diagrams contributing here are in principle analogous to those for the J/ψππ channel, shown in Fig. 4.However, in contrast to that channel, we exclude diagrams containing a Z c (3900).This is based on the observation that Z c (3900) does not show a significant contribution to the h c π invariant mass distribution.Additionally, we point out that the Z c (4020) is not included in this work, since this would require a complete treatment of the coupled channels, and of the {D D * , D * D * } subsystems, which is postponed to future work.Moreover, the contact terms that drive the contributions shown in diagrams (e) and (f ) of Fig. 4 in the J/ψππ channel are omitted here since they violate spin symmetry.This symmetry violation is overcome by the loop diagrams as a result of the spin symmetry violation that enters through the mass differences of D and D * as well as D 1 and D 2 -the former one being included explicitly in the calculation, the lat-ter one by choosing an energy range where the D 2 contribution should be negligible.For a detailed discussion on how the spin symmetry gets restored in the heavy quark limit even in the presence of hadronic molecules, see Ref. [61].In summary, for the h c ππ channel we only include the box topologies shown in Fig. 7, expecting some deviations from experiment as a result of the omission of the Z c (4020).On the other hand, it is not expected that the Z c (4020) will generate significant structures in the total cross section of h c ππ, which is the focus of the current work, since in this case, the narrow peak from Z c (4020) in the πh c subsystem is smeared.The same effect is demonstrated explicitly in this work, where the narrow structures of the Z c (3900) seen in the J/ψπ subsystem do not modify visibly the energy dependence of the total cross section for J/ψππ.
D. e + e − → X(3872)γ If Y (4230) is a D 1 D hadronic molecule and both Z c (3900) and X(3872) are D * D hadronic molecules with I(J P C ) = 1(1 +− ) and I(J P C ) = 0(1 ++ ), respectively, the production mechanism of the latter pair in Y (4230) decays must be analogous [62].Only that the particle radiated off in the course of the Y (4230) decay must have positive C parity for the transition to the Z c and negative C parity for the transition to the X(3872).Thus, all that needs to be done to get from the diagram that generates the Z c in Y (4230) → πZ c to the one that generates the X(3872), is to replace the pion in the final state by a photon.The resulting diagrams are shown in Fig. 8.
For each reaction discussed so far the electromagnetic production mechanism and the strong decay were entangled in a special way.What makes the e + e − → µ + µ − especially interesting is, that here we may isolate production from decay, since the total cross section is by far dominated by the real valued tree-level diagram (first diagram in Fig. 9) and the hadronic cross sections only contribute significantly through their interference with the mentioned dominating one.
Moreover, the decays of Y (4230) and ψ(4160) into the same hadronic channels induce some mixing of these in the γ * → γ * transition amplitudes.The diagrams contributing to the process are shown in Fig. 9.The mentioned mixing of the two vector resonances is depicted here as the hatched blob.The imaginary part of this mixing amplitude is given by the respective interference terms that contribute also to the various exclusive hadronic channels discussed above.It is dominated by the transitions Y (4230) → D D * π → ψ(4160), since the D D * π cross section is by far the largest hadronic cross section.The details of the calculations can be found in Sec.III C 8. Therefore, the simultaneous study of the hadronic channels and the e + e − → µ + µ − channel provides a sanity check for the size of the induced mixing of the vector states, which turn out to be significant.

F. Further channels
As shown in Fig. 1, in addition to the channels discussed in detail above, the Y (4230) is  seen also in the final states ωχ c0 , ηJ/ψ and ψ(2S)ππ.In this work we do not study this last decay channel as the ψ(2S)π invariant mass distributions vary so dramatically when the total energy is changed mildly form 4.226 to 4.258 MeV [63] that there must be some highly nontrivial interplay of different mechanisms at work that to our understanding are not yet under-stood microscopically (while in Ref. [43] a description of the invariant mass distributions is provided, no attempt is made to understand the energy dependence of the total cross section).
For the first two channels, both triangle diagrams as well as direct transitions contribute as shown in Figs. 10 and 11, respectively.Below we discuss the results for these channels as well.

III. FORMALISM
In this section the formalism underlying the calculations is presented in some detail with additional material provided in the appendix.Those readers most interested in the results and their physical interpretation may want to jump to Sec.IV immediately.

A. The Y(4230) as a D 1 D state
We can write the D 1 D wavefunction as a negative C-eigenstate As the Y (4230) is predominantly produced by cc it must be an iso-singlet.Following the convention ) the iso-singlet wavefunction is given by where the couplings g Y 0 and g 1 include the heavy quark mass normalization of the fields.Typically a proper field redefinition allows one to absorb the effect of non-perturbative hadronhadron scattering into a pole term.This is not possible only if there is more than one pole on the physical sheet in the mass range of interest [64].Since this is not the case here we can safely set the parameter g 1 to zero2 .Thus we get for the D 1 D scattering potential where the bare Y propagator reads with ω Y for the on-shell energy of the Y (4230) from the field normalization and E = √ s.Here we dropped the spin indices although G 0 and various other propagators below refer to the propagation of a spin one particle.The reason is that in our non-relativistic treatment the spin structure simply refers to a δ ij -the spin simply runs through unchanged.The relation of the bare propagator G 0 (E) to the full propagator G Y (E) is given by the Dyson equation From this one finds for the D 1 D scattering amplitude with Note that the last term in the denominator was added to account for the contribution to the width of the Y (4230) from the various inelastic channels.The self-energy Σ for a resonance R can be derived from the standard, scalar one-loop diagram, which reads in dimensional regularization for the intermediate two-body state a, up to terms irrelevant in what follows .
Here s = E 2 .The masses in the expression refer to the masses of the two particles propagating in channel a.To come from this to the expression for the self-energies employed in the propagators, we use With this subtraction, the real part of the inverse Y propagator vanishes at E = m 0 and it reduces significantly the correlations between couplings and bare masses [65].
Using the D 1 D scattering amplitude and the Y (4230) propagator G Y , one is in the position to derive the pointlike production operator M Y via the Y (4230) to D 1 D (see Fig. 12 for the graphical illustration) where c is the direct coupling of the photon to D 1 D in the quantum numbers J P C = 1 −− , which vanishes in the HQSS limit, and α is the source term coupling of the photon to the bare Y state.Eq. ( 13) gives the impression as if it had a pole at the bare mass m 0 , however, from Eq. ( 8) one gets that which allows us to rewrite Eq. ( 13) as Here, in the purely one channel D 1 D problem, unitarity requires the parameters to be real.However, allowing for additional complex phases at the photon-resonance couplings enables us to effectively include other effects, such as interference between ψ(4040) and ψ(4160), as will be discussed below.

B. Production in the presence of coupled channel final state interactions
Through hadronic final state interactions, unitarity links contact terms to resonance propagators -a special example of this was already demonstrated above: Eq. ( 13) contains both a contact term to the final state as well as the resonance contributions collecting the interactions in that final state.As demonstrated there, employing unitarity makes the tree level production term vanish and the final amplitude, Eq. ( 15), is proportional to the dressed resonance propagator.
Analogously one cannot discuss other treelevel operators or contact terms involved in the decay transitions without the inclusion of the non-perturbative final state interactions in the relevant subsystem parametrized via the pertinent resonance propagators.In analogy to the Y propagator provided in Eq. ( 10) we find for the propagator of the Z c (3900) from solving the related Dyson equation where the sum in the denominator runs over all relevant channels, which for the Z c (3900) are D * D and J/ψπ [66] (denoted as channels 1 and 2 respectively) and E is the energy in these subsystems.Furthermore, g i stands for the couplings of the Z c (3900) with the channel i, and Σ i refers to the self energy in the corresponding channel.As before the trivial spin structure of the propagator is not shown.As the energy range studied in this work is far above the threshold of J/ψπ, the contribution of this channel to the self energy is well approximated by a constant whose real part can be absorbed into the bare mass m 0 .
The production amplitude F 2 for channel 2, shown in Fig. 13, can now be expressed as with M i denoting the production operator for channel i.The expression for F 1 can be easily obtained from Eq. ( 17) by interchanging 1 and 2. Defining M 1 = g 1 M1 and M 2 = g 2 M2 , the form factor F can be therefore expressed by which takes the form already expressed in the Feynman diagrams in Sec.II, namely that it is given by some vertex structure for the source term, times the Z c propagator times the respective channel coupling.As indicated in the lower line of Fig. 14, the effective coupling of the Z c to the J/ψπ channel, here abbreviated as g 2 , contains besides a contact term also a triangle topology.The same is true for M 2 , as shown by the upper line of Fig. 14.The triangles for M 2 and g 2 in this figure are essentially identical, except for the couplings of Y and Z c to D-mesons, which are evidently different.In particular, they incorporate the D ( * ) D( * ) J/ψ vertex in P-wave, causing the principal value part of these triangles to depend on a regulator that must be renormalized by a contact term, consistent in both cases.In the picture advocated here, where the decay of Y (4230) is predominantly governed by diagrams involving the D 1 D intermediate state rather than those depicted in Fig. 14, it is reasonable to assume that the overall coefficient M2 connecting M 2 and g 2 is real-valued.While formally present in the transition amplitude, we observed that the fits to the experimental data do not need the term proportional to ( M2 − M1 ), since it was consistently found to be zero.We thus omit the corresponding terms from the start and employ for the production amplitude 1 (α being free parameters to be determined in the fit.These form factors appear in both the Y and the ψ decays.The corresponding strength parameters of the latter resonance are denoted as β  3a, to provide the most significant contribution followed by the triangle loop and contact interactions.As argued in Appendix A, see the discussion below Eq. (A16), the D 1 (2420) can decay into D * π in both S-and D-wave, such that the spin structure of the Y → D 0 D * − π + amplitude can be written as where we introduced as short hand notation s and h ′ are fixed from the D 1 decay properties -details are given in App. A. To respect the Goldstone theorem, stating that the pion amplitude has to vanish in the chiral limit for p π → 0, the S-wave vertex and other amplitudes below scale with the on-shell pion energy ω π = m formalism described in section III B is used, 2 +E DD * )δ kj (21) This expression agrees to Eq. ( 19), only that the generic name g 1 used there was now adapted to the notation employed in this section.
The relative factor 2 between the tree-level and Z c contributions comes from the isospin coefficient of the Y (4230) wave function shown in Eq. ( 4 Formally, the Y (4230) should be treated as emerging from a D D * π three-body system, which can be most conveniently handled using time-ordered-perturbation-theory (TOPT) [67].In preparation for this more complete treatment, that we will attack in a subsequent publication, we evaluate also the loop integrals in this work using the same formalism.Thus, the scalar triangle with pion emission is given by where the D 1 energy is given by with m D 1 and Γ D 1 for the mass and width of the D 1 , respectively.The other particle energies are defined analogously, however, with their widths neglected.The width of the D 1 can here be treated as constant, since the D 1 pole is sufficiently high above the πD * threshold [58].We checked that the energy dependencies of the various loop diagrams included in this study agree to the analogous loops evaluated covariantly.As argued above, a simultaneous treatment of Y → J/ψππ and Y → D D * π is possible only if also the interference with the ψ(4160) is included.The contribution of ψ(4160) → D 0 D * − π + is parameterized as again in line with Eq. (19).The free parameters that appear in the equations above were fixed in a fit to data -the resulting values are listed in Tab.II.
Here, a comment is in order.As our focus lies in examining the interference effect between Y (4230) and ψ(4160) on the line shapes in various channels-specifically, for e + e − → D 0 D * − π, J/ψπ + π − , J/ψK + K − , h c π + π − , χ c0 (1P )ω, J/ψη, and X(3872)γ-we derive the corresponding observables by multiplying the amplitudes M Y and M ψ , discussed in this and subsequent sections, by the same complex couplings g γR = exp(iδ Rγ )em 2 R /f R of the photon with the resonance R as defined in Eq. (A23), where R represents both Y (4230) and ψ(4160).Clearly, for all channels listed above, only the relative phase of the two resonances plays a role.This is, however, not the case for e + e − → µ + µ − , where the two phases enter individually, see Sec.III C 7 for details.
2. e + e − → J/ψπ + π − The Feynman diagrams for e + e − → J/ψπ + π − are shown in figure 4. The dominant contributions corresponding to the molecular nature of the Y (4230) are the box, below denoted as M □ , and triangle, M △ , topologies, since those contain the D 1 D intermediate state.
As the second triangle in figure 4 c) is divergent, due to the internal P -wave vertex that is connected to a J/ψ coupling to a pair of D ( * ) mesons, a counter term M △ CT is also introduced 2 where g Zc J/ψπ is the coupling of Z c → J/ψπ, given by the triangle transition shown in figure 14 g To reduce the run-time of the numerical evaluation of the loop integrals only two out of four contributions of the Y (4230) wave function shown in Eq. ( 4) with p π 1 ↔ p π 2 are considered, as the differences due to isospin breaking are negligible small.For example, for the box topologies shown in figure 5 only the particle content spelled out first at each line in the boxes is evaluated explicitly, while those in brackets are included via a multiplication with a factor 2. The Y (4230) contact term M J/ψππ Y CT has two contributions, one in the ππ invariant mass from the subtraction polynomial of the ππ final state interaction and the other in J/ψπ from the chiral contact term and intermediate Z c (3900) with c △ CT denoting the free parameter of the triangle counterterm.The amplitudes of the loop diagrams are given below, where the notation and numerical implementation are discussed in the appendix C where the q I , q II , q III , q ′ I , q ′ II , q ′ III denote the relative momenta at the J/ψD ( * ) D ( * ) vertex for the different box and triangle topologies.Additional free parameters come from the production polynomials of the Y (4230) and Ψ(4160) contact terms, namely α  In general, one expects the diagrams for h c π + π − to be analogous to J/ψππ, apart from the fact that the Y (4230) contact term is omitted as it violates HQSS.In addition, since the h c π subsystem does not show any prominant signal of the Z c (3900), no triangle operators are included in this study.Meanwhile, the h c π subsystem shown in Ref. [68] shows a strong peak from the Z c (4020), which would, however, require to include the coupling of Y (4230) → D 1 D * , as the Z c (4020) couples strongly to D * D * .On the other hand, the Z c (4020) is not anticipated to generate significant structures in the total cross section of h c ππ, which is part of the current analysis.The inclusion of this state will be postponed for the upcoming full coupled channel analysis, such that for now we only consider the box topologies, where the free parameters are completely fixed by D 0 D * − π + , J/ψπ + π − and the two-body final states.The amplitude there-fore reads with M hcππ □ given by where B hcππ I and B hcππ II are given in the appendix C.
4. e + e − → J/ψK + K − With J/ψππ included in the study, we can also easily access J/ψK K, as the main contribution is expected to go via the ππ → K K final state interaction in the S-wave, where no new parameters need to be introduced.Here the contributions of the triangle topologies are negligible, as the partial wave projection on the ππ system contains a tiny S-wave piece due to the presence of the near on-shell Z c (3900) in the J/ψπ subsystem.The amplitude is given by where we collected the loop diagrams in the amplitude and (32) Furthermore, M jl Y →J/ψKK is a strange source shown in Fig. 6.We postpone the inclusion of strange triangles, including the Z cs (4000), to a later, more complete analysis.In this sense we regard this channel in this analysis as a consistency check.On the other hand, the Z cs (4000) can only appear in conjunction with an additional kaon within the triangular mechanism.Consequently, this state is expected to contribute significantly only in the energy range around 4470 MeV, well exceeding the energy range of interest in this study, even when accounting for the Z cs width.
5. e + e − → χ c0 ω The Feynman diagrams are shown in figure 10.The main contribution is expected from the triangle, which scales like the scalar traingle as both the D 1 → Dω and DD → χ c0 are S-wave at leading order.Additionally, there are two Swave contact terms for the Y (4230) and ψ( 4160) respectively where c △ χ c0 ω , c Y χ c0 ω and c ψ χ c0 ω are free parameters.The width of the ω is included by convolving the cross section for a fixed ω mass with the ω spectral function -see, e.g., Ref. [69].
6. e + e − → J/ψη For J/ψη, the couplings of the triangle shown in figure 11 are fixed.The vector-vector-axial vector vertex of the contact terms must couple via ϵ µνρσ which reduces to a three dimensional ϵ mjl in the rest frame of the incoming vector particles: where q denotes the relative momentum at the J/ψ vertex and c Y J/ψη and c ψ J/ψη are free parameters.We do not consider the mixing of the sin-glet η 1 and octet η 8 to the physical η and η ′ states, but just match η 8 = η, as the mixing effects are small.The diagrams for Y (4230) → X(3872)γ are shown in Fig. 8, and are analogous to Y (4230) → Z c π as well as J/ψη.However, the quality of data for X(3872)γ does not allow one to distinguish between the triangle and contact transition of Y (4230) → X(3872)γ, such that we omit the latter from the start3 .The vector-vector-axialvector coupling of D 1 → D * γ scales with ϵ kjl , such that the amplitude is given by with c Y Xγ and c ψ Xγ being free parameters to be determined in a fit.
8. e + e − → µ + µ − As already explained in section II E we consider three main contributions for e + e − → µ + µ − , namely for the tree level amplitude and we introduced and where g γR = exp(iδ Rγ )em 2 R /f R defined in Eq. (A23) with δ Rγ denoting a phase factor discussed in Sec.IV.The individual terms in Eq. ( 36) represent the different diagrams shown in Fig. 9.The imaginary part of M RR ′ mix is fixed by unitarity and can be reconstructed from the optical theorem where f = DD * π, χ c0 ω, J/ψη and X(3872)γ are all final states with significant contributions from both ψ(4160) and Y (4230) studied in this work.Note that the sum runs over all allowed final states with the given particle contentaccordingly f = [DD * π] should be understood as ) Since all those channels are connected via isospin symmetry, they can be included via a proper multiplicity factor -clearly for that we need to neglect, e.g., the mass differences between the different channels.For example, for DD * π we denote decay amplitudes for the transition of Y (4230) and ψ(4160) to the experimentally measured channel D 0 D * − π + as A and B, respectively, where in accordance to Eq. (39) A and B do not contain the resonance propagators, but only the decay vertices.Summing over all channels one therefore obtains where f ∈ [DD * π] was defined in Eq. ( 41).The factor 4 in Eq. ( 43) arises from the four different decay modes of the Y (4230) wave function given in Eq. ( 4).For each mode the subsequent D 1 decay can produce a charged or a neutral pion, e.g.D 0 1 can decay into D * 0 π 0 and D * + π − , where the amplitudes scale as 1 and 1/ √ 2, respectively, due to the isospin factors.The additional factors arising in the other channels are 3/2 for J/ψππ and 1 for χ c0 ω, J/ψη and X(3872)γ.The real part of M RR ′ mix can in principle also be constructed dispersively, however, there is still freedom in the subtraction constant.So for now we just approximate it via a real constant

IV. FIT STRATEGY, RESULTS AND DISCUSSION
In 2022 and 2023 BESIII published new XYZ data sets for J/ψπ + π − [8] and D 0 D * − π [36] with very impressive statistics.Those data clearly highlight the asymmetric lineshapes of the total cross sections in these two channels.It turns out that from those channels most of the parameters specific for the Y (4230) are fixed.The Z c (3900) shows up prominently only in the D D * [70] and J/ψπ ± [71] subsystems of those channels.To get a better constraint on the light quark SU(3) singlet and octet components (for details see Appendix B) we also include J/ψK + K − in the first fit.This may overestimate the contributions of the contact term in J/ψK + K − to some extend as it needs to compensate for a possible contribution from the missing Z cs (4000) triangle, but allows us to reduce the correlation of the parameters.We do not include the data for the J/ψπ 0 π 0 channel in the fit, due to their reduced statistics in comparison to J/ψπ + π − .Since µ + µ − is the only channel showing a clear separation of the Y (4230) and ψ(4160) signals and their interference, it is also included in the first fit.This further allows us to properly separate photon and strong couplings, since in the hadronic channels they only appear as a product.With this in mind, our fit strategy is the following: 1.The resonance parameters of the Y (4230) and Z c (3900), as well as the channel dependent parameters of D 0 D * − π, J/ψπ + π − , J/ψK + K − and µ + µ − are fitted simultaneously to the D 0 D * − π, J/ψπ + π − , J/ψK + K − and µ + µ − total cross sections, the D D * , J/ψπ ± and π + π − invariant mass distributions, and the pion Jackson angle extracted from D 0 D * − π + .
2. With the resonance and channel dependent parameters of D 0 D * − π + , J/ψπ + π − and µ + µ − being fixed, the remaining parameters in the channels χ c0 ω, J/ψη and X(3872)γ are fitted to the corresponding cross sections data.
3. At last, the parameters obtained in the previous steps are used as initial parameters for a global fit to all observables.
If we were working with a complete formalism, with all relevant channels dynamical and unitarity imposed, all parameters would necessarily be real.However, here some ingredients are approximated.e.g. as shown in Ref. [72] the direct transition of a photon to the D 1 (2420) D intermediate state that predominantly couples to the Y (4230), if it is a hadronic molecule, is suppressed by heavy quark spin symmetry, since this narrow D 1 state has a light quark cloud with j = 3/2.On the other hand, there is no such suppression for the transition of the photon to D 1 (2430) D, where the broad D 1 (2430) has its light quark cloud with j = 1/2.The D 1 (2430) D intermediate state may thus act as a doorway state to feed the production of the molecule.This effect can be included effectively via a complex coupling of the Y (4230) to the photon.Moreover, the ψ(4160) production from a photon sits in the tail of the ψ(4040) [22] an effect which may also be included by allowing for a complex coupling.It is worth noting that, while all hadronic cross sections are sensitive to the difference of those two phases only, the leptonic cross section e + e − → µ + µ − probes the phases individually, as shall be discussed below.It turns out that all other parameters of the model can be chosen real valued.
The results of the fits are shown in Fig. 15 to Fig. 20, the parameters and statistical uncertainties that emerge from the fit are listed in Table II.We see that the interplay of the ψ(4160) and Y (4230) is important and shows a non-trivial impact in almost all final states.This naturally explains the large scatter of the resonance parameters of the Y (4230) in the single channel analyses of BESIII -c.f.Fig. 1.
With the central values of the parameters fixed in the fits, the pole parameters of the Y (4230) can be extracted from its propagator.We find MeV .
(45) To estimate the uncertainties we varied the parameters that control the pole location of the Y (4230), namely m 0 and Γ in in Eq. (10), and − (133.9 ± 4) g 1 − (14.9 ± 0.9) 10 −3 g 8 (24 ± 1) 10  I. Parameters of the model as determined in the fit.We find the value of f J/ψ to be strongly dependent on the fit range in D 0 D * − π + , such that we did not assign an uncertainty to this quantity.refitted all other parameters for either of them fixed.These two pole parameters turn out to be correlated very little.The ranges of the variation of the pole location given above, reflect the variations of the parameters allowed by still providing an acceptable fit to the data for the J/ψππ final state, which turned out to be the most constrain- ing for the Y pole parameters.While a more rigorous uncertainty estimate would require the implementation of Bayesian methods, clearly beyond the scope of this study, the uncertainties provided above should still provide some guidance on how well the pole can be pinned down in an analysis of the kind presented here.We now discuss the results for the various channels in some detail.The results for the D 0 D * − π channel are shown in Fig. 15.The apparent peak structure around 4. 22 GeV emerges in our study from the interplay of the ψ(4160) and Y (4230).Remarkably, this interplay manifests differently in the D 0 D * − π and J/ψππ channels -we refer to Fig. 2  = (3884 − i44/2) MeV.In comparison to Ref. [73] we find a slightly higher mass, however, double the width for the Z c (3900).It remains to be seen if this feature is caused by the incomplete ππ − K K final state interaction, used in this work.The data for the pion Jackson angle are also reproduced well.Contrary to Ref. [37], in this study the S-wave is more prominent due to the presence of the ψ(4160) as well as the S-wave decay of the D 1 (2420).
Naturally, a prominent contribution from the D 1 D intermediate state not only influences strongly the energy dependence of the total cross section but also the D * π invariant mass distributions.Our predictions for those at total energies near the Y (4230) pole location and near the nominal D 1 D threshold are shown in the left and right panel of Fig. 16, respectively.In both panels the peak from the D 1 is clearly visible at the upper end of the spectrum.While the data currently available do not allow us to provide an unambiguous determination of the vari-  ous parameters leaving some freedom in the actual height of the D 1 signal, the presence of such a peak is a model independent prediction of the molecular scenario.Any model that does not account for the D 1 D as a prominent component of the Y (4230) wave function will not show such a structure -as such this invariant mass distribution is a crucial observable to either support or disprove the molecular picture.
The results for the J/ψππ final state are shown in Fig. 17.A linear non-interfering background of 9 pb is added due to the presence of the J/ψπ continuum.The loop contributions, dominant in the molecular scenario, enhance the cross section at the D 1 D threshold, allowing for a description of the highly asymmetric lineshape with just a single pole -in the experimental analyses of Refs.[7,8] this asymmetry was generated by the additional state Y (4320) as described in the introduction.It should be noted that also in Ref. [50] an analysis is presented, where, in particular, the J/ψππ data is described with essentially the same resonance content as presented here, along with the addition of a ψ(4415) state, but without including threshold effects.In this case, the asymmetry of  17. Fit results for the J/ψπ + π − cross section and the J/ψπ ± and π + π − invariant mass distributions.J/ψπ + π − XYZ data from Ref. [8], J/ψπ ± and π + π − invariant mass distribution from Ref. [71].The data for the J/ψK K channel are taken from Ref. [23].
the peak in the total cross section is driven by interferences, predominantly involving ψ(4160) and ψ(4415).While we regard our explanation of the data as more natural, since the data indicate some structure right at the D 1 D threshold, at some point experiment will allow us to choose between the two explanations, not only since the energy dependencies in the mentioned energy range are different (but not sufficiently to be distinguished given the current quality of the data), but also since an analysis of the type presented in Ref. [50] will provide completely different D D * and D * π invariant mass distributions compared to the ones shown in Figs. 15 and 16, respectively.
The J/ψπ ± invariant mass distribution shows a prominent peak, generated by the Z c (3900) pole, the D * D threshold and the nearby triangle singularity, and its reflection.In principle, the J/ψπ ± and D 0 D * − lineshape can also be described by just including the triangles and introducing a contact interaction for D D * → D D * , where the cusp is then generated simply by kinematic effects of the D D * rescattering without any resonance structure.However, in agreement with Ref. [74], we find that the strength of the D D * → D D * transition potential, necessary for producing the pronounced structure in the D * D invariant mass distribution, becomes too large to justify a perturbative approach, as the next to leading order loop contribution becomes larger than leading order.In this case the full non-perturbative hadron-hadron interaction naturally generates a pole in the S-matrix, corresponding to the Z c (3900).
As pointed out in the introduction, the ψ(4160) needs to be included to get a consistent simultaneous description of the J/ψππ and D * Dπ final states.We allow this well established vector charmonium state to contribute to both of these channels (as well as to all other channels included in the analysis), however, the fit reveals that no significant coupling of the ψ(4160) to the J/ψππ is needed.Indeed, the fit puts the parameter β (2) 1 , introduced in Eq. ( 24), for the production of the ψ(4160) to a value consistent with zero.
The results for the µ + µ − final state are shown Fit results for the µ + µ − cross section.Upper panel: The measured born cross section; lower panel: The same, however, with the cross section from the tree level amplitude subtracted.The data are taken from Ref. [22], where the data points with an uncertainty smaller than 32 pb are shown in black to better highlight the structure in the data.
in the upper panel of Fig. 18.The cross section is completely dominated by the real tree level amplitude, shown as the first diagram in Fig. 9. Accordingly, following Eq.( 36), the signal of interest to us reads to very good approximation This quantity is shown in the lower panel of Fig. 18.As argued above, we allow for complex couplings at the resonance-photon vertices.Contrary to all observables studied so far, where only the relative phase of the Y (4230) and ψ(4160) contributions played a role, this leptonic cross section is sensitive to the individual phases of these resonance contributions.The phases of those couplings are thus fixed by the fit to the µ + µ − cross section.
In Ref. [22] a sum of Breit-Wigner terms with complex couplings is used to parameterize the data including the ψ(4040), ψ(4160), S(4220) ≡ Y (4230) (and the ψ(4415), which is, however outside our energy range of interest).We find that the complex phase δ ψγ called for by the fit in the production vertex of the ψ to the photon agrees within uncertainties to the one of the experimental paper.With this phase included we can reproduce the µ + µ − lineshape in the energy range studied.We see that the contribution of the ψ(4160) is dominant in comparison to that of the Y (4230) at least in the energy below 4.2 GeV, as expected in the molecular scenario, since the coupling of a photon to the D 1 D channel violates spin symmetry [26].One should, however, keep in mind that there should also be some suppression of the coupling of the ψ(4160) to the photon, if it indeed is predominantly a D-wave charmonium state.The peak in the data near 4230 MeV in our fits emerges from both the interference of the two resonances and the Y (4230) itself.The main contribution to the imaginary part of the pertinent mixing matrix element of the ψ(4160) and the Y (4230) is generated from the D D * π intermediate state -this part is fixed completely by the data for e + e − → D D * π.As outlined above, the corresponding real part is here taken as a free parameter that is adjusted in the fit.It is reassuring, however, that the real and the imaginary part of the mixing amplitude contribute with comparable strength, as shown in figure 19.Although the energy dependence emerging from the real and imaginary part of the mixing amplitude, A mix , resembles that of a single resonance structure, it emerges from an interplay of the different resonance propagators as well as the mixing amplitude, M RR ′ mix , as outlined in Eq. ( 36) and below.One may naively expect that the imaginary part of the mixing matrix element does not contribute to the total cross section significantly, as only interferences of the strong amplitudes with the real tree level amplitude matter quantitatively.However, the phases of the resonance propagators in Eq. ( 36) non-trivially mix real and imaginary parts of the  4160) in e + e − → µ + µ − , denoted by A mix in Eq. (39).The brown dash-dotted curves here and in Fig. 18 are identical.mixing amplitude, allowing both contributions to interfere with the tree level amplitude.
The fit result for J/ψK + K − is shown in the top right of Fig. 17.Note that the line shape emerging for this channel is closely linked to that of the J/ψππ channel -there are no new independent parameters entering for this hidden strangeness channel.In our fit the contact term is the dominant contribution.A possible reason is that it needs to absorb the effects of the Z cs (4000) and the corresponding triangle diagrams not included in this analysis, though their main effect is expected at the energies above those considered here.The boxes again show a very strong enhancement in the cross section at the D 1 D threshold explaining the apparent asymmetry in the data.We find the Y → J/ψπ + π − → J/ψK + K − contribution generated by the ππ/K K FSI to be by far dominant in the studied energy range, in comparison to the box with strange D-mesons, as shown in Fig. 6, where only the D 1 D cut goes on-shell.At higher energies, above the D * s DK threshold at about 4.47 GeV, the D * s DK intermediate state in this box will go on shell.Consequently, we expect a more pronounced contribution from the strange source in this mass range.Starting from this energy also the Z cs (4000) generated via the triangle mechanism should contribute considerably.The data are taken from Ref. [13] It should be stressed that in our analysis the unusual energy dependencies of the J/ψππ and J/ψK K cross sections emerge from the same physics, which is natural given the approximate SU(3) flavor symmetry of QCD, while in the experimental analyses the former is driven by an interference of the Y (4230) with the new resonance Y (4320) [8] and the latter predominantly by the shifted threshold for J/ψ KK vs. J/ψππ with some small distortion from an interference with another new resonance, called Y (4500) [23].
To complete the discussion of the final states with three hadrons, in Fig. 20 we show the cross sections with an h c in the final state.These rates are of particular interest, since the h c has its cc pair in the spin singlet state, while it is originally produced in a spin triplet via its coupling to the photon.Thus, in this transition the heavy quark spin changes, at odds with heavy quark spin symmetry.However, besides violations of that symmetry due to the relatively small charm quark mass, spin symmetry transitions can also be circumvented by the presence of hadronic molecules.In the molecular picture it is therefore natural to expect the h c ππ and the J/ψππ cross sections of similar order of magnitude as is confirmed by the data.Moreover, by using values for both the J/ψD ( * ) D( * ) and the h c D ( * ) D( * ) couplings available in the literature (details on how the various couplings were determined are given in App.A), we can describe the cross sections in both channels, providing addi- Ref. [20],and Ref. [21], respectively.tional support for the molecular picture.
We now turn to a discussion of remaining hadronic two-particle final states, also included in Fig. 21.As one can read from the figure, the energy dependences of the χ c0 ω, the J/ψη, and the X(3872)γ cross section are rather different: While the first one shows a very narrow structure, the structure in the second is already a lot broader and the one in the last is more than four times as broad as the first -this is also reflected in the resonance parameters extracted in the single channel analyses of the BESIII collaboration collected in Fig. 1.In contrast to this, our model allows us to describe all three cross sections with consistent resonance parameters as a result of an interplay of the two vector resonances ψ(4160)and Y (4230): the narrow peak in the χ c0 ω channel emerges from a destructive interference of the triangle diagram shown in Fig. 10a and the ψ(4160) contact term, shown in Fig. 10c, since the energy dependencies of the two contributions are quite different, as can be clearly seen in Fig. 21 -we included the width of the ω by a convolution of the cross section with the omega spectral function as explained above which is the origin of the not vanishing cross section below the nominal χ c0 ω cross section.The mechanism we propose here is different to that studied in Ref. [75], however, the energy dependence found there appears inconsistent with that of the newest data set for this channel measured at BESIII [20].Also for the J/ψη and X(3872)γ final states the interplay of the two resonances is crucial, but less dramatic.

V. SUMMARY AND OUTLOOK
In this work we simultaneously analysed the lineshapes of the cross sections for e + e − to seven hadronic channels in the mass range from 4.2 to 4.35 GeV as well as data for e + e − → µ + µ − .We show that a description of all those channels is possible with consistent resonance parameters for a single Y state.This was made possible mainly because of two features of our model: We included the interference of the exotic Y (4230) with the more conventional ψ(4160) and considered the effects of the D 1 D intermediate states.The prominence of the latter is natural in a scenario, where the Y (4230) is a hadronic molecule in this channel.The interference of the Y (4230) and ψ(4160) is especially important to get a consistent description of both channels D 0 D * − π + and J/ψπ + π − , making a substantial impact on the former.It is at the same time necessary to describe the µ + µ − cross section, where the mixing of the two vector states deduced from the strong decay channels is in fact consistent with what is needed for the leptonic final state.We interpret this as providing additional support for the mixing scenario advocated here.The ex-plicit inclusion of the D 1 D intermediate states reflects itself in a significant distortion of line shapes, which are especially prominent in final states with a J/ψ and two light mesons.In particular, contrary to the experimental analyses, in our study the energy dependencies of the total cross sections for e + e − → J/ψππ and e + e − → J/ψK K emerge from the same physics as expected from the approximate SU(3) flavor symmetry of QCD.
For the other final states, within our model especially the energy dependence of the χ c0 ω cross section is non-trivial, which emerges from the distinct energy dependencies of the triangle diagram, influenced by the D 1 D intermediate state, and the ψ(4160) contact term, which does not.Moreover, within our analysis we understand that the very different lineshapes of the χ c0 ω cross section and, e.g., the J/ψη cross section emerge naturally through the interference of the Y (4230) with the ψ(4160).
To summarize, we have demonstrated that the data for electron-positron annihilation in to various final states in the mass range from 4.2 to 4.35 GeV are consistent with the existence of just a single vector charmonium-like state in this mass range -the ψ(4230) also known as Y (4230), with properties consistent with its being a D 1 (2420) D hadronic molecule.Moreover, we show that a consistent description of all channels with the same resonance content is possible only, if we allow for an interference with the conventional ψ(4160) resonance.
The non-trivial insights of this work were possible only because we studied various final states simultaneously -to get access to reliable resonance parameters this appears to be unavoidable, while single channel analyses have the tendency to provide resonance parameters with a wild scatter as shown in Fig. 1. arXiv:1902.10957 [hep-ph].

Appendix A: Lagrangians and parameter determination
To construct the Lagrangian we define superfields representing the different light-quark spin multiples [76].The ground states of heavy mesons with light quark quantum numbers j P = 1 2 − will be denoted by H a .In the presence of one unit of angular momentum there are 2 spin multiplets, one with j P = 3 2 + and one with j P = 1 2 + .
In the following the former is of relevance, which is denoted by T a .The states in the latter have widths of the order of 300 MeV and are only included implicitly.All together we may thus write where a is the SU(3) flavor index.We have, e.g., for j P = 1 (A2) The superfields creating heavy mesons are given by (A3) The corresponding superfields containing an anti-heavy quark Q can be constructed using the charge conjugation operator C = iγ 2 γ 0 , where we are following the convention The heavy field operators contain a factor √ M H and therefore have dimension 3/2.Pseudoscalar mesons couple through the vector V µ and axialvector A µ current containing an even and odd number of boson fields, respectively, conserving chiral symmetry.Chiral symme-try violation is introduced via constructions of the kind with χ = 2BM, where M is the quark mass matrix and B is a scale parameter related to the chiral condenstate.Here the exponential parametrisation is employed for the light Goldstone boson fields: with f π = 92 MeV denoting the pion decay constant in the chiral limit.The Lagrangian is constructed by imposing invariance under heavyquark spin and chiral transformation [77][78][79].The kinetic terms are and the relevant terms for the interaction are given by The relation between the decay width and the effective coupling of a resonance R with total angular momentum J R decaying into the two-body final state a in the narrow width approximation is given by [66] Γ where m R denotes the resonance mass, the phase space factor is ρ a (m 2 R ) = 2p a (m R )/(16πm R ) and p a denotes the relative momentum of the decay particles in the rest frame of the resonance, with m a,i for the masses of the particles in channel a.The summation runs over the polarisations of the final and initial state, respectively, in necessary.The pertinent matrix elements can be read off the Lagrangian given in Eq. (A9) straightforwardly allowing one to determine the couplings from the experimentally measured decay widths.
The squared matrix element for the transition of D * a → D b ϕ ab , summed over the D * polarisations, is given by pol.

|M D
where the coefficient c ab can be read off from the Goldstone boson matrix provided in Eq. (A7): c +0 = 1 and c ++ = 1/ √ 2. Using eq.(A10) we extract for where the central values listed in the Review of Particle Physics by the Particle Data Group [66] were used Analogously from D * + → D + π 0 we find with Γ(D * + → D + π 0 ) = 25.6 keV.It is this value that we use in our calculations in line with Ref. [80].Since in this work we do not aim at a calculation with controlled uncertainties but more at demonstrating the consistency of the existing data with just a single molecular state in the mass range studied, we feel safe to not keep track with the individual uncertainties of the parameters employed.The interaction of the j P l = 3 2 + doublet {D 2 , D 1 } with {D * , D} and the Goldstone bosons Φ given in Eq. (A9) can be reexpressed as The decay of the narrow D 1 into D * π is predominately in a D-wave, since the S-wave is suppressed by heavy-quark spin symmetry, which calls for the conservation of the light quark total angular momentum in the decay.However, violations of this symmetry in the charm sector can be sizeable.To get an estimate for the S-wave strength in the D 1 decay, we can use the fact that the spin partner of the D 1 , the D 2 , can only decay into D * π and Dπ in a pure D-wave due to the total angular momentum conservation [59].Adding the partial widths, according to Eq. (A10), the total width of the D 2 is given by ) where the factor 3/2 in front of each term results from adding the partial widths of the D ( * )+ π 0 and D ( * )0 π + in line with what was done for the decay of the D * .Using Γ D 2 = 47.3MeV [66], one can extract h ′ = 0.82 GeV −1 .Our calculation is not sensitive to the sign of this coupling which we chose positive.
Allowing for a D * π S-wave, the expression for the total width of the D 1 reads From Eq. (A16) one finds pol.
A20) where again a factor 3/2 was included to account for the two possible final states.With h ′ fixed above, one finds Γ d-wave D 1 = 15 MeV, in agreement with Ref. [59].Since the total width of the D 1 (2420) is 31 MeV [66], the remainder must be generated by the S-wave decay.Using where ω π denotes the energy of the pion and again the factor 3/2 accounts for the two decay channels D + 1 → D * 0 π + and D * + π 0 .This leads to h ′ s = 0.57.Below we study a pion angular distribution, which is sensitive to the relative sign of h ′ and h ′ s .We here already account for the observation that the data call for equal signs of the two.Photons couple via the field-strength tensor where A µ denotes the photon field.In this way gauge invariance is preserved automatically.The production of a vector resonance from a photon is thus described by where V µν = ∂ µ V ν − ∂ ν V µ and V denotes either the field for the Y (4230) or the ψ(4160).The implications of heavy quark spin symmetry on charmonium production from photons are discussed later in the chapter, but as the production of the Y (4230) must go via the broad D 1 (2430), thus we may allow for an additional phase in case of a point-like production.For the decay of Y (4230) → X(3872)γ we can describe the E1 transition of D 1 going to D * γ with the following Lagrangian where E i denotes the electric component of the photon field.
We now come to the description of the ground state doubly heavy vector fields of relevance to this study.Heavy quark spin symmetry allows us to write Q Q superfields [81].The ℓ = 0 superfield R where the interaction with D/D * is given by the Lagrangian The resulting vertex factors are = denoting the relative residual momentum between the D mesons.At leading order the masses of the multiples are degenerate m D * = m D = m H and q simplifies to q The coupling is traditionally parametrised as which includes the leading spin-symmetry violating effects via the mass factors.The ℓ = 1 superfield P (Q Q) contains the spin triplet χ c0 , χ c1 , χ c2 and the singlet h c (A29) Due to the Proca-constraint, P µ v µ , the leading order Lagrangian for the interaction of the ℓ = 1 spin-multiplet with D and D * contains only a single term: From this the vertex factors evaluate to A31) where we fixed α = 0 such that p α ≈ m hc .The coupling is parametrised as where f J/ψ = 416 MeV was determined using vector-meson dominance (VMD) [54] and f χ c0 = 510±40 from numerical results of QCD sum rules [82].Those parameters carry a systematic uncertainty which is difficult to quantify.In the fits we allow f J/ψ to vary within 10% of its value, while we fix f χ c0 to its central value, since our fit is not sensitive to this quantity.where σ k = 1 − (4m 2 k /s)Θ(s − 4m 2 k ) and the subscript indices j, k, ... refer to the coupled channels, which in our case correspond to ππ and K K. Furthermore, T jk corresponds to the meson-meson coupled-channel amplitude.The solution is given by Muskhelishvili Omnés function encoding the ππ/K K final state interaction.Therefore, the full amplitude is given by the sum of the amplitudes M j and Γ, involving the lefthand and right-hand (unitarity) cuts only, respectively, which reads for the amplitude T and the Omnés matrix is taken from Refs.[83,84].We now have a closer look at the principal value integral.Using the short hand notation kl Ω −1 jk T kl σ l M l = f j (z) we get The integral in the last line can be evaluated analytically P.V. and that in the second needs to be done numerically.However, in the exploratory study we aim at here, we neglect this part which we expect to play a role of a correction.The quality of the fits we achieve might be taken as justification of this treatment a posteriori, although some shift in importance of different contributions in a complete analysis cannot be excluded.Thus, the modified ππ and KK s-wave is given by such that the full amplitude (B4) is approximated as This also enables us to project the FSI onto the K K-channel, allowing us to also determine Y (4230) → J/ψππ → J/ψ KK.
The subtraction polynomial (P n−1 ) j in Eqs.(B4) and (B5) is matched to the Y (4230) → J/ψϕϕ chiral contact term.In Ref. [85] it was found that both SU

(B8)
The Lagrangian L Y ψϕϕ at leading order in chiral expansion is given by [86] with J = (ψ/ √ 3)1.The resulting s-wave projected chiral contact terms relevant for the J/ψππ J/ψK K final states are given by with q 2 = λ 1/2 (m 2 Y , m 2 J/ψ , s)/(2m Y ), resulting in To summarize, the amplitudes, incorporating the ππ/K K FSI, used in our calculations, are provided by Eqs.(B6), (B5) and (B11). with (C2) In the center of mass frame of the Y (4230) one can choose the momenta in a way that only the D * has an angular dependency, where p π = p π (0, 0, 1) T , such that z = cos θ denotes the cosine of the polar angle of the loop momentum l.Due to the width of the D 1 only the last propagator in eq.(C1) has poles on the real axis, as in comparison the D * width is negligible small.We define such that f is regular in l.The integral can now be rewritten as where the trivial integration over the loop momentums azimuthal angle is performed.Do- The inverse factor of l is canceled by f (E, l), while p π is canceled by the phase space integration.
With the relation Finally we arrive at C13) where the remaining radial integration can be preformed numerically.
In case of the J/ψD ( * ) D( * ) vertex, which scales with the loopmomentum itself, the integrand is modified ac- In this work the second cut corresponds for most box topologies to DD * π, which can go onshell.Analogous to the triangle we isolate the singularity and define the remaining part in a function f (E, l, φ, z, p 1 , p 2 ) that is regular in l and z.However, different to the triangle, it is not possible to to perform the integration of the polar angle analytically as f is also dependent on z, such that we perform a numerical subtraction + f (E, l, φ, z 0 , p 1 , p 2 ) where z 0 is the pole of the propagator.The integration of the second summand con now be done analytically, resulting in (C17) The remaining φ and l integration are performed numerically using Gauss-Legendre integration.To decrease the number of sample points needed to get a stable result, it is useful to further split the l integration at the poles of the propagator, as the distribution of Gauss-Legendre sample points is more dense at the integration borders  where T stands for the sum over the different time orderings and G i denotes the propagators for the different cuts, e.g.G 1 = 1/(E − w D 1 − ω D ).

FIG. 1 .
FIG.1.Mass and width of the two Y -states discussed in the introduction as extracted from the experimental analyses of the individual channels shown by the labels.All data below a mass value below 4240 MeV is interpreted a Y (4230) and the one data point above refers to the Y (4320).The experimental values are taken from[8,9,[12][13][14][15][20][21][22][23].The red dot denotes the pole location for the Y (4230) as extracted in this work.
ψη and X(3872)γ, in the mass range from 4.2 to 4.35 GeV, under the assumption that the Y (4230) is a D 1 D

FIG. 2 .
FIG.2.Fit result for D 0 D * − π + and J/ψπ + π − channels, including the ψ(4160) (solid red line) and omitting it (dashed blue line).The data for the D 0 D * − π + channel are from Ref.[36] and those for the J/ψπ + π − channel from Ref.[8].The vertical red dotted lines indicate the positions of the nominal mass of the ψ(4160) and D 1 D threshold, respectively.The black dotted line at 4.33 GeV denotes the upper limit of the D 0 D * − π + fit range.
+ e − → D * D * and D * s D * s are consistent with an interference from the D 1 D molecular nature of the Y (4230).
D 1 (2420) is unstable with a width of about 30 MeV.It decays predominantly into D * π in D-wave, thus the final state with closest connection to a possible molecular nature of the Y (4230) is the D * Dπ channel.The corresponding diagrams are shown in Fig. 3.Both diagram (a) and diagram (c) scale directly with the large Y D 1 D coupling.Moreover, they are both enhanced by the near on-shell D 1

FIG. 3 .
FIG. 3. Diagram contributing to e + e − → DD * π. a) tree level, b) Y (4230) contact term, c) Triangle, d) ψ(4160) contact term, where for the last two the final state interactions in the doubly heavy subsystem are included.

FIG. 4 .
FIG. 4. Diagram contributing to e + e − → J/ψππ.The thin lines in the box and triangle denote D * or D mesons.a) and b) boxes, c) triangle, d) triangle counter term, e) Y (4230) contact term, f) ψ(4160) contact term, where for the last two the final state interactions in the doubly heavy subsystem are included.
)The effective Lagrangian for the coupling of D 1 D to Y (4230) and D 1 D self-interactions reads[38]

FIG. 12 .
FIG.12.The Y (4230) induced production of the D 1 D pairs from a pointlike source.The solid lines denote D 1 and D mesons as well as the bare propagator G 0 , double line stands for the dressed propagator of the Y(4230), and the wiggly line corresponds to the initial photon.

FIG. 13 .
FIG.13.Feynman diagram for production of channel 2. The scattering amplitudes T ij in the channels ij (i, j = 1, 2) are related to the Z c propagator as T ij = g i G Z g j .
1. e + e − → D 0 D * − π + With D 0 D * − π + being the channel with the most direct access to the molecular nature of the Y (4230), one expects the tree-level decay, shown in Fig.

2 π + p 2 πFIG. 14 .
FIG. 14. Upper line: Feynman diagrams for production operator for Y (4230) → J/ψππ.In the full amplitude the J/ψ and one of the pions undergo final state interactions driven by the Z c .Lower line: The corresponding transition Z c → J/ψπ.
), as the Z − c (3900)π + pair is produced via Y (4230) → D + 1 D − and Y (4230) → D D1 .The coefficients of the iso-triplet production of the Z c (3900) from D * D and the C = −1 eigenstate are absorbed into the coupling g Z0 of the Z c (3900) with D D * .
well as the triangle counterterm c △ CT .The inclusion of the ππ − KK final state interaction is discussed in Appendix B.

7 .
e + e − → X(3872)γ 017 ± 0.003) GeV µ + µ − c mix (0.6 ± 0.01)TABLE for an illustration.In addition to this, we find a strong enhancement at the D 1 D threshold in the cross section, mainly driven by the prominent D 1 decay in Dwave.The deviations of our results from the data, starting around 4.35 GeV, are expected, as the molecular scenario predicts an additional bound state in the D 2 D * channel [25, 46, 47] 4 , which will be included in a subsequent study.The peak at low D D * invariant masses is generated by the interplay of the tree-level two-step decay Y (4230) → D 1 D → D * π D, the contact mechanism and the triangle operator.The last two mechanisms involve the rescattering of D D * into the Z c (3900).The resonance parameters of the Z c (3900) are very poorly constrained.The fit seems to prefer masses slightly above the D D * threshold, however, for the whole mass range of approximately m Z ∈ [3.86, 3.9] GeV, the data are described with similar quality.In the current fit the pole closest to the real axis of the Z c (3900) appears at the +− sheet with respect to the J/ψπ and D D * channels, respectively (where +(−) denotes the sign of the imaginary part of the three-momentum in each channel), with s Zc(3900) pole
FIG.18.Fit results for the µ + µ − cross section.Upper panel: The measured born cross section; lower panel: The same, however, with the cross section from the tree level amplitude subtracted.The data are taken from Ref.[22], where the data points with an uncertainty smaller than 32 pb are shown in black to better highlight the structure in the data.

FIG. 19 .
FIG.19.Contributions to the cross section difference from the real and imaginary parts of the mixing of Y (4230) and Ψ(4160) in e + e − → µ + µ − , denoted by A mix in Eq.(39).The brown dash-dotted curves here and in Fig.18are identical.

Appendix B : √ 2 α 1 − 1 dz
Inclusion of ππ − KK final state interactionAn amplitude M corresponding to the given isospin I (the isospin index is omitted in what follows) can be projected to a partial wave M l with definite angular momentum lM l (s) = 1 2 M(s, z)P l (z) , (B1)where P l denotes the Legendre polynomial of degree l, z the scattering angle and α = 0, 1, 2 is a symmetry factor for identical particles in initial and final states (e.g.α = 2 for AA → BB, α = 1 for AA → C C and α = 0 for AB → AB) .The full amplitude can be reconstructed using the orthogonality relation of the Legendre polynomialsM(s, z) = √ 2 α l (2l + 1)M l (s)P l (z) .(B2)On the other hand, from the unitarity of the Smatrix the discontinuity of the production amplitude M l is given by discM l j (s) = 2i k T * jk (s)σ k (s)M l k (s) , (B3)

Mn π dz z n Ω − 1 Ω
full j (s) = M j + Γ j = M j + k Ω jk (P n−1 ) k + lm s kl (z)T lm (z)σ m (s)M 0 m (z) jk (P n−1 ) k + s n π P.V. [...] + iT jk σ k M 0 k ,(B4)where P n−1 is a polynomial of the order n − 1, which is discussed below, and Ω is the S-wave Omnés matrix.The amplitudes M j correspond tocdiagrams discussed in Secs.II B and II C, while the right-hand cut in Γ emerges from the ππ/K K FSI in an S wave.The necessary input

C d 3 l 2π 3 Num(l, p 1 , p 2 ) 16ω 1 ω 2 ω 3 ω 4 G 1 G 2 G 3 ,
where l P i corresponds the poles of the propagator in l, which can be calculated analytically.The general notation used in this work is B(C, Num)= T