Production of $X(3872)$ and a Photon in $e^+e^-$ Annihilation

If the $X(3872)$ is a weakly bound charm-meson molecule, it can be produced in $e^+ e^-$ annihilation by the creation of $D^{*0} \bar D^{*0}$ from a virtual photon followed by the rescattering of the P-wave charm-meson pair into the $X$ and a photon. A triangle singularity produces a narrow peak in the cross section for $e^+ e^- \to X \gamma$ 2.2 MeV above the $D^{*0} \bar{D}^{*0}$ threshold. We predict the normalized cross section in the region of the peak. We show that the absorptive contribution to the cross section for $e^+ e^- \to D^{*0} \bar D^{*0} \to X \gamma$, which was calculated previously by Dubynskiy and Voloshin, does not give a good approximation to the peak from the triangle singularity.


I. INTRODUCTION
Since early in this century, a large number of exotic hadrons whose constituents include a heavy quark and its antiquark have been discovered in high energy physics experiments [1][2][3][4][5][6][7][8][9][10]. The first of these exotic heavy hadrons to be discovered was the X(3872) meson. It was discovered in 2003 in exclusive decays of B ± mesons into K ± X by observing the decay of X into J/ψ π + π − [11]. The J P C quantum numbers of X were eventually determined to be 1 ++ [12]. Its mass is extremely close to the D * 0D0 threshold, with the difference being only 0.01 ± 0.18 MeV [13]. This suggests that X is a weakly bound S-wave charm-meson molecule with the flavor structure However, there are alternative models for the X [1][2][3][4][5][6][7][8][9][10]. The X has been observed in 7 different decay modes, many more than any of the other exotic heavy hadrons. Despite these many decay modes, a consensus on the nature of X has not been achieved.
There may be aspects of the production of X that are more effective at discriminating between models than the decays of X. If the X is a weakly bound charm-meson molecule, it can be produced by any reaction that can create its constituents D * 0D0 and D 0D * 0 . It can be produced by the creation of D * 0D0 and D 0D * 0 at short distances of order 1/m π , where m π is the pion mass, followed by the binding of the charm mesons into X at longer distances. The X can also be produced by the creation of D * D * at short distances followed by the rescattering of the charm-meson pair into X and a pion at longer distances [14,15].
One way in which the nature of a hadron can be revealed in its production is through triangle singularities. Triangle singularities are kinematic singularities that can arise if three virtual particles that form a triangle in a Feynman diagram can all be on their mass shells simultaneously. There have been several previous investigations of the effects of triangle singularities on the production of exotic heavy mesons [16][17][18][19]. Guo has recently pointed out that any high-energy process that can create D * 0D * 0 at short distances in an S-wave channel will produce Xγ with a narrow peak near the D * 0D * 0 threshold due to a charmmeson triangle singularity [20]. One such process is electron-positron annihilation, which can create an S-wave D * 0D * 0 pair recoiling against a π 0 . Guo suggested that the peak in the line shape for Xγ due to the triangle singularity could be used to determine the binding energy of X more accurately than a direct mass measurement. Because of the charm-meson triangle singularity, a high-energy process that can create an S-wave D * D * pair at short distances can also produce Xπ with a narrow peak near the D * D * threshold. We noted previously the existence of a such a narrow peak in the production of Xπ in hadron colliders [15] and in B meson decays into KXπ [14], but we did not recognize the connection to triangle singularities.
The quantum numbers 1 ++ of the X imply that Xγ can be produced by e + e − annihilation into a virtual photon. The virtual photon can create D * 0D * 0 at short distances in a P-wave channel, and the charm-meson pair can subequently rescatter into Xγ. The production of Xγ in e + e − annihilation near the D * 0D * 0 threshold was first discussed by Dubynskiy and Voloshin [21]. They calculated the absorptive contribution to the cross section from e + e − annihilation into on-shell charm mesons D * 0D * 0 followed by their rescattering into Xγ. They predicted that the cross section has a narrow peak only a few MeV above the D * 0D * 0 threshold. In retrospect, the narrow peak comes from a triangle singularity. In Ref. [22], we calculated the cross section for e + e − → Xγ in the energy region near the D * 0D * 0 threshold, including the dispersive contributions as well as the absorptive contributions. The dispersive contributions have a significant effect on the shape of the narrow peak. We also predicted the normalization of the cross section for e + e − → Xγ using a fit to the cross section for e + e − → D * + D * − by Uglov et al. [23].
The production of Xγ in e + e − annihilation has been studied by the BESIII collaboration [24,25]. The cross section was measured at the center-of-mass energies ranging from 4.008 GeV to 4.6 GeV. The cross section was not measured at energies near the D * 0D * 0 threshold at 4.014 GeV, which is where the narrow peak from the charm-meson triangle singularity is predicted to appear. The observation of this peak would provide strong evidence in support of the identification of the X as a charm-meson molecule.
In this paper, we describe in detail the calculation of the cross section for e + e − → Xγ. In Section II, we calculate the cross section for e + e − → D * 0D * 0 near the threshold. We determine the normalization of the cross section by using a previous fit to Belle data on e + e − → D * + D * − . In Section III, we calculate the cross section for e + e − → Xγ from the creation of a P-wave D * 0D * 0 pair by a virtual photon followed by the rescattering of the charm mesons to Xγ. We reduce the amplitude to a scalar loop integral. In Section IV, we calculate the loop amplitude analytically and show that there is a triangle singularity. We predict the normalized cross section for e + e − → Xγ at energies in the region near the peak from the triangle singularity. In Section V, we express the loop amplitude in terms of the Schrödinger wavefunction for the bound state. We point out the differences from the wavefunction used by Dubynskiy and Voloshin in Ref. [21]. In Section VI, we calculate the absorptive contribution to the cross section for e + e − → Xγ from intermediate charm mesons D * 0D * 0 that are on their mass shells. We show that it does not provide a good approximation to the peak in the cross section from the triangle singularity. Our results are summarized in Section VII. In an Appendix, we present a diagrammatic derivation of the Schrödinger wavefunction of the X in a frame where its momentum is nonzero. We also present a prescription for calculating the cross section for producing the X resonance feature in cases where the X is not a narrow bound state. A pair of spin-1 charm mesons D * 0D * 0 can be produced from the annihilation of e + e − into a virtual photon. The Feynman diagram for this process is shown in Fig. 1. We use nonrelativistic normalizations for the charm mesons in the final state. In the center-of-momentum (CM) frame, the matrix element has the form where √ s is the invariant mass,v and u are the spinors for the colliding e + and e − , and J is the matrix element of the electromagnetic current between the QCD vacuum and the D * 0D * 0 state. Near the threshold for producing D * 0D * 0 , the charm-meson pair is produced in a P-wave state with total spin 0 or 2. The matrix element of the electromagnetic current that creates D * 0 andD * 0 with momenta +k and −k and with polarization vectors ε andε can be expressed as J i = A ijkl k j ε * kε * l . The Cartesian tensor A ijkl is where A 0 and A 2 are amplitudes for creating D * 0D * 0 with total spin 0 and 2, respectively. The numerical prefactor of A 2 in Eq. (3) was chosen for later convenience. The differential cross section for producing D * 0D * 0 with scattering angle θ is where M * 0 is the mass of the D * 0 and k is the relative momentum of the D * 0D * 0 pair: The cross section for e + e − annihilation into D * 0D * 0 near the threshold is The absolute values of the two amplitudes A 0 and A 2 could in principle be determined experimentally from the value of the cross section and from the angular distribution at a single energy near the threshold. The Belle collaboration has measured exclusive cross sections for e + e − annihilation into several pairs of charm mesons, including D * + D * − [26,27]. Uglov et al. have analyzed the Belle data using a unitary approach based on a coupled channel model [23]. They included a spin-2 F-wave amplitude for e + e − → D * D * as well as spin-0 and spin-2 P-wave amplitudes. Their fit to the cross section for e + e − → D * + D * − as a function of the center-of-mass energy W relative to the D * + D * − threshold is shown in Fig. 2, along with the spin-0 and spin-2 P-wave contributions. 1 The fitted cross section increases to a local maximum of 3.8 nb at W = 36 MeV, and then decreases to a local minimum of 2.1 nb at W = 97 MeV. At the energy of the first local minimum, the spin-2 F-wave contribution has increased to 21% of the cross section. Thus the P-wave contributions alone give a good approximation to the cross section for W up to about 100 MeV.
Near the D * + D * − threshold at 4020.5 MeV, the spin-0 and spin-2 P-wave contributions to the cross sections have the k 3 behavior in Eq. (5). A fit to the two terms in the cross section in Eq. (5), with M * 0 replaced by the mass M * 1 of the D * + and k = [ 1 The curves in Figure 1(c) of Ref. [23] are mislabeled in the figure caption. We wish to relate the amplitudes A 0 and A 2 for e + e − → D * 0D * 0 to the corresponding amplitudes for e + e − → D * + D * − . An amplitude A i for D * + D * − is the sum of an isopin-0 amplitude and an isospin-1 amplitude, while the corresponding amplitude for D * 0D * 0 is the difference. Fig. 2 shows that the spin-2 P-wave contribution to the cross section for e + e − → D * + D * − has a strong peak near 4050 MeV, and that the spin-0 P-wave contribution also has a peak near that energy. This strongly suggests that the amplitudes A 2 and A 0 are dominated by the ψ(4040) charmonium resonance, which has mass 4039 MeV and width 80 MeV. Since the width of this isospin-0 resonance is much larger than the 6.8 MeV difference between the charm-meson pair thresholds, the amplitudes A 2 and A 0 for D * 0D * 0 must also be dominated by the ψ(4040) resonance. We assume the isospin-1 amplitudes are negligible compared to the resonant isospin-0 amplitudes. We therefore approximate the amplitudes A 0 and A 2 for e + e − → D * 0D * 0 by the corresponding amplitudes for e + e − → D * + D * − . We can predict the cross section for e + e − → D * 0D * 0 near its threshold at 4013.7 MeV by inserting the values of |A 0 | and |A 2 | in Eq. (6) into Eq. (5).
The Belle data on e + e − annihilation into charm-meson pairs in Refs. [26] and [27] has also been analyzed by Du, Meißner, and Wang using an approach that takes into account P-wave coupled channel effects by solving Lippmann-Schwinger equations with contact interactions between the charm mesons [28]. They presented their fit to the cross section for D * + D * − in the form of histograms with 20 MeV bins. It is therefore not possible to determine |A 0 | and |A 2 | from their results. A comparison of their results with those of Ref. [23] can however Cross section for e + e − → D * + D * − near threshold as a function of the center-of-mass energy W relative to the D * + D * − threshold. The dots are the spin-2 P-wave contribution (higher red dots) and the spin-0 P-wave contribution (lower blue dots) from Uglov et al. [23], which are shown as curves in Fig. 2. The curves are fits to the analog of Eq. (5), and they determine the coefficients in Eq. (6).
give some indication of the possible size of theoretical errors. In Ref. [28], the spin-0 and spin-2 P-wave contributions to the cross section were given separately only at the single invariant mass √ s = 4.040 GeV, which is about 17 MeV above the D * + D * − threshold. The cross sections are 1.52 nb and 1.23 nb, respectively, compared to 0.75 nb and 2.19 nb from the fit in Ref. [23]. The sum of the two cross sections in Ref. [28] is only about 6% smaller than their sum from the fit in Ref. [23]. However the ratio 0.81 of the spin-2 and spin-0 cross sections in Ref. [28] is significantly smaller than the ratio 2.92 from the fit in Ref. [23]. We conclude that |A 0 | 2 + |A 2 | 2 can be determined more accurately by fitting cross sections than |A 2 |/|A 0 |. The phases of the amplitudes A 0 and A 2 can be chosen so that a linear combination of the amplitudes has no interference. For this choice, the values in Eq. (6) can be expressed as We will take these to be our preferred values for the amplitudes. However we will consider all possible complex values for the ratio A 2 /A 0 that are consistent with the value of |A 0 | 2 +|A 2 | 2 in Eq. (7). The decay width of the D * 0 can be predicted from measurements of D * decays: Γ * 0 = (55.9 ± 1.6) keV [29]. The present value of the difference E X between the mass of the X and the energy of the D * 0D0 scattering threshold is [13] The central value corresponds to a charm-meson pair just above the scattering threshold. The value lower by 1σ corresponds to a bound state with binding energy |E X | = 0.17 MeV. The X can be produced in e + e − annihilation through the creation of D * 0D * 0 by a virtual photon followed by the rescattering of the charm-meson pair into Xγ. The Feynman diagrams for this process are shown in Fig. 4. The vertex for the virtual photon with vector index i to create D * 0 andD * 0 with momenta +k and −k and with vector indices m and n is eA ijmn k j , where the Cartesian tensor is given in Eq. (3). The vertex for the transition of D * 0 to D 0 γ with a photon of momentum k is −eν ijm k m , where i and j are the vector indices of D * 0 and γ. The transition magnetic moment eν can be determined from the radiative decay width of D * 0 : where the photon energy ω satisfies ω + ω 2 /2M 0 = δ. The radiative width of the D * 0 can be predicted from measurements of D * decays: [29]. This determines the transition magnetic moment: ν = 0.95 GeV −1 . The binding of D * 0D0 or D 0D * 0 into X can be described within an effective field theory called XEFT [30,31]. The vertices for the couplings of D * 0D0 to X and D 0D * 0 to X can be expressed as where γ X is a parameter with dimensions of momentum and k and l are the vector indices of the spin-1 charm meson and the X [32]. In the Appendix, this vertex is derived from the residue of the pole in the elastic scattering amplitude for the charm-meson pair. We will assume in the body of this paper that the X is a narrow bound state whose binding energy |E X | is significantly larger than the decay width Γ * 0 ≈ 56 keV of the D * 0 . The energy of the X can therefore be expressed in terms of the positive real binding momentum γ X as For E X = −0.17 MeV, the binding momentum is γ X = 18 MeV. In the simplest plausible model for the resonant scattering amplitude of the charm mesons, this binding momentum coincides with the parameter γ X in the vertex in Eq. (10). We will therefore use Eq. (11) to determine the real parameter γ X from the binding energy |E X |. The production of the X resonance feature in the more general case where X is not a narrow bound state is considered in the Appendix. The matrix element for e + e − → Xγ is the sum of the two diagrams in Fig. 4. We use nonrelativistic propagators for the charm mesons. The matrix element for producing X and γ with momenta q and −q and with polarization vectors ε X and ε γ can be expressed as where F (W ) is a function of the center-of-mass energy W = √ s − 2M * 0 relative to the D * 0D * 0 threshold. We will refer to F (W ) as the loop amplitude, because it can be expressed as an integral over the undetermined energy and momentum in the charm-meson loops in Fig. 4: To obtain the scalar loop integral in Eq. (13), we used rotational symmetry to replace a factor of k i inside the momentum integral by (q · k)q i /q 2 . The resulting expression for the current J i is In the production of Xγ invariant mass √ s = 2M * 0 + W , the center-of-mass energy W relative to the D * 0D * 0 threshold is determined by energy conservation: where q is the photon energy and δ = M * 0 − M 0 . We have used the Galilean-invariant approximation for the kinetic energy of the X, in which its kinetic mass is the sum of the masses of D * 0 andD 0 . Assuming W is less than or of order δ, we can solve Eq. (15) for the photon energy q as an expansion in powers of δ/(M * 0 +M 0 ) and E X /δ: The differential cross section for producing X with scattering angle θ is We have used nonrelativistic phase space for X and relativistic phase space for the photon. The cross section for e + e − annihilation into Xγ near the D * 0D * 0 threshold is The factor that depends on A 0 and A 2 differs from the value of |A 0 | 2 + |A 2 | 2 in Eq. (7) by a multiplicative factor that depends on the value of the ratio A 2 /A 0 . For our preferred values A 2 /A 0 = ±1.9 i in Eq. (7), the multiplicative factor is 0.73. If we allow all possible complex values of A 2 /A 0 consistent with the value of |A 0 | 2 + |A 2 | 2 in Eq. (7), the multiplicative factor can range from 0.34 to 1.31.

IV. PEAK FROM THE TRIANGLE SINGULARITY
The loop amplitude F (W ) in Eq. (13) has a kinematic singularity called a triangle singularity [16][17][18][19] in the limit where the binding energy of X(3872) is 0 and the decay width of the D * 0 is 0. The singularity arises from the integration region where the three charm mesons whose lines form triangles in the diagrams in Fig. 4 are all on their mass shells simultaneously. The two charm mesons that become constituents of the X are both on their mass shells in the limit where the binding energy is 0. There is a specific energy W where the spin-1 charm meson that emits the photon can also be on its mass shell. If γ X = 0 and Γ * 0 = 0, F (W ) has a logarithmic divergence at W . When either γ X or Γ * 0 is nonzero, F (W ) has a narrow peak near W .
If the integral over the loop energy ω in Eq. (13) is evaluated by closing the contour in the lower half-plane, the only pole in ω is from the propagator of the spin-1 charm meson that becomes a constituent of the X. The denominators of the propagator of the spin-1 charm meson that emits a photon and the propagator of the spin-0 charm meson that becomes a constituent of the X can be combined into a single denominator by introducing an integral over a Feynman parameter x. The integral over the loop 3-momentum k can be evaluated analytically. The resulting loop amplitude can be expressed as The coefficients of the polynomial inside the square root are where k 2 = M * 0 W , q is the real photon energy given by Eq. (16), and γ X is the real binding momentum. The dependence of F on the real energy W is through k 2 and the energy q of the photon. In Eq. (20b), we have used the conservation of energy in Eq. (15) to express b as a sum of four terms that are all order δ 2 or smaller when W is order δ 2 /M * 0 . Note that the sum of the three coefficients in Eqs. (20) does not depend on the energy: The integral over x in Eq. (19) can be evaluated analytically: The loop amplitude in Eq. (22) is particularly simple in the limit Γ * 0 → 0: where The triangle singularity arises from the logarithm in Eq. (22). The denominator of the argument of the logarithm vanishes at a complex energy that approaches the real axis in the limits Γ * 0 → 0 and γ X → 0. In these limits, the denominator is zero at the energy for which k = (µ/M 0 )q. The energy W at which the triangle singularity occurs can be obtained by solving the equation k = (µ/M 0 )q, where k = √ M * 0 W and q is the function of W in Eq. (16). The solution can be expanded in powers of δ/(M * 0 +M 0 ) and E X /δ: The prediction from the leading term is W = 2.7 MeV. The Argand diagram for the amplitude F (W ) at equally spaced values of W is shown in Fig. 5. As W increases towards 0 from below, F (W ) increases along the positive real axis. As W passes through 0, Im[F (W )] begins to increase. The amplitude F (W ) then follows a roughly circular path. The value of F (W ) is (1.20 + 0.86 i) |F (0)| at W = 2.2 MeV, where |F (W )| 2 has its maximum value. At large W , the decreasing amplitude F (W ) approaches the positive imaginary axis. Thus the amplitude F (W ) moves counterclockwise around a loop in one quadrant of the complex plane. This should be compared to the path followed by the amplitude A(E) for an ideal resonance as a function of the energy E. As E approaches the resonance, A(E) increases from 0 along the positive real axis. It moves counterclockwise around a circle in the upper half of the complex plane, crossing the positive imaginary axis as E passes through the maximum of |A(E)| 2 . As E increases further, A(E) decreases towards 0 along the negative real axis. The path of the triangle-singularity amplitude is qualitatively similar to that of an ideal resonance, except that F (W ) traces out a loop that remains entirely in one quadrant of the complex plane.
The cross section for e + e − → Xγ as a function of the invariant mass √ s is obtained by inserting the loop amplitude F (W ) in Eq. (22) into the expression in Eq. (18). The cross section near the D * 0D * 0 threshold is shown in Fig. 6 for three values of the binding energy: |E X | = 0.30 MeV, 0.17 MeV, and 0.10 MeV. The peaks of the line shapes are produced by the triangle singularity. The position of the peak is 2.2 MeV above the D * 0D * 0 threshold, and it is insensitive to the binding energy. The height of the peak is also insensitive to |E X | provided |E X | Γ * 0 . In Fig. 6, we have used the preferred values of the amplitudes A 0 and A 2 given by Eq. (7). For |E X | = 0.17 MeV, the cross section at the peak is 0.51 pb. If we allow all possible complex values of A 2 /A 0 consistent with the value of |A 0 | 2 +|A 2 | 2 in Eq. (7), the cross section can be larger by a factor of 1.80 or smaller by a factor of 0.47. Beyond the peak, the cross section decreases to a local minimum at an energy W near 40 MeV before increasing because of the k 3 dependence of the P-wave cross section for producing D * 0D * 0 . For |E X | = 0.17 MeV, the cross section at the local minimum is about 0.07 pb. The BESIII collaboration has measured cross sections for e + e − annihilation into Xγ at center-of-mass energies ranging from 4.008 GeV to 4.6 GeV by observing X in the final states J/ψ π + π − and J/ψ π + π − π 0 [24,25]. They did not measure the cross section at energies between 4.009 MeV and 4.178 MeV, which includes the predicted energy 4.016 GeV of the peak from the triangle singularity. The BESIII collaboration measured the cross sections in 10 MeV steps between 4.178 GeV and 4.278 GeV [25]. The largest measured value of the product σ Br of the cross section and the branching fraction Br of X into J/ψ π + π − was about 0.5 pb. We derived upper and lower bounds on the branching fraction Br for the X bound state in Ref. [33]. The BaBar collaboration has recently measured the inclusive branching fraction of B + into K + plus the X resonance feature [34]. It implies a branching fraction into J/ψ π + π − from the X resonance feature that provides the loose lower bound Br > 4% on the branching fraction from the X bound state. An upper bound Br < 33% can be derived from measurements of branching ratios of J/ψ π + π − over other short-distance decay modes of the X. Thus the height of the peak from the charm-meson triangle singularity could be a significant fraction of the cross section that has been measured in the higher energy region.  7). If A 2 /A 0 is allowed to vary with |A 0 | 2 + |A 2 | 2 fixed, the normalizations can change by a factor ranging from 0.47 to 1.80.

V. BOUND-STATE WAVEFUNCTION
In Ref. [21], Dubynskiy and Voloshin (DV) presented an approximation for the absorptive contribution to the cross section for e + e − → X(3872)γ that involved the Schrödinger wavefunction of the X. The momentum-space wavefunction ψ(k) in the rest frame of the bound state is a function of the relative momentum k of the constituents. The standard normalization for the wavefunction is The universal momentum-space wavefunction for a weakly bound S-wave molecule whose constituents have short-range interactions is [35] In the case of the X, this wavefunction should be a good approximation out to momenta k of order m π , beyond which it should decrease more rapidly with k.
An expression for the loop amplitude F (W ) involving a wavefunction can be obtained by closing the contour for the integral over ω in Eq. (13) in the lower half-plane. The resulting loop amplitude can be expressed as where the last factor in the numerator is We have simplified the denominator by using the conservation of energy in Eq. (15) to eliminate W . The function ψ(q − k, k) in Eq. (28) can be interpreted as the wavefunction for a bound state whose constituents have momenta q − k and k. If we set Γ * 0 = 0, this is just the universal wavefunction in Eq. (26) with the relative momentum k replaced by |k − (µ/M 0 )q|, which is a function of the velocity difference k/M * 0 − (q − k)/M 0 of the constituents. The wavefunction is expected to depend on the velocity difference in a Galilean-invariant effective field theory. The wavefunction in Eq. (28) is the appropriate wavefunction for the bound state in a frame where its momentum is q. A diagrammatic derivation of this wavefunction is presented in the Appendix. The wavefunction in Eq. (28) satisfies the normalization condition in Eq. (25) if γ X is real and Γ * 0 = 0. The wavefunction for the X in its rest frame used by DV in Ref. [21] was where Λ is an adjustable parameter. DV identified the parameter γ X in Eq. (29) with the binding momentum of X. This identification can be justified rigorously only in the limit Λ → ∞. The DV wavefunction in Eq. (29) is a regularized form of the universal wavefunction ψ X (k) in Eq. (26) that reduces to ψ X (k) in the limit Λ → ∞. The universal wavefunction decreases as 1/k 2 when k is much larger than γ X . The subtraction in the DV wavefunction in Eq. (29) makes it decrease as 1/k 4 when k is much larger than the momentum scale Λ. At k = 0, the DV wavefunction is larger than ψ X (0) by a factor of 1 + 3γ X /2Λ. In Ref. [21], DV illustrated their results using the values Λ = 200 MeV and Λ = 300 MeV. The DV wavefunction was used previously by Voloshin in a study of the decays of X into D 0D0 γ [36]. In Ref. [21], DV assumed that the wavefunction for X in a frame where it has momentum q could be obtained from Eq. (29) simply by replacing the relative momentum k by |k − 1 2 q|, which is a function of the momentum difference k−(q −k) of the constituents. DV used that wavefunction to calculate the absorptive contribution to the loop amplitude F (W ). Their prescription for the wavefunction can be used to calculate the full amplitude F DV (W ) by replacing ψ(q − k, k) in Eq. (27) by ψ DV (|k − 1 2 q|). The amplitude F DV (W ) can be obtained most easily by first calculating the amplitude F X (W ) obtained by replacing ψ(q − k, k) in Eq. (27) by ψ X (|k− 1 2 q|), where ψ X is the universal wavefunction in Eq. (26). The expression for F X (W ) as a Feynman-parameter integral is the same as Eq. (19) except that a, b, and c are where k 2 = M * 0 W and q is the photon energy given by Eq. and by setting µ/M 0 to 1 in the coefficient of M * 0 Γ * 0 in b. The analytic result for F X (W ) in the limit Γ * 0 → 0 is The loop amplitude F DV (W ) in the limit Γ * 0 → 0 can now be obtained by subtracting from F X (W ) the expression with γ X replaced by Λ except in the factor √ πγ X , and then multiplying by the first prefactor in Eq. (29): In Fig. 7, we show the cross sections for e + e − → Xγ in the triangle singularity region calculated using the approximation for the loop amplitude F DV (W ) in Eq. (32) with Λ = 200 MeV. We have set the binding energy to |E X | = 0.17 MeV. The cross section is compared to the cross sections calculated using the complete result for F (W ) in Eq. (22) and the result for F (W ) with Γ * 0 = 0 in Eq. (23). The approximation using the DV wavefunction is significantly lower for E < 0, it is a little higher at the peak, and it rapidly approaches the exact result as W increases.

VI. ABSORPTIVE CONTRIBUTION
The loop amplitude F (W ) in Eq. (13) has an absorptive contribution that corresponds to e + e − annihilation into on-shell charm mesons D * 0D * 0 followed by the rescattering of the charm-meson pair into X(3872)γ. In the limit Γ * 0 → 0, the absorptive contribution is simply the imaginary part of F (W ). The imaginary part can be obtained by cutting rules that replace the propagators of the spin-1 charm mesons in Eq. (13) by delta functions. The absorptive contribution to the loop amplitude in the limit Γ * 0 → 0 is where W = √ s − 2M * 0 . The integral over k can be evaluated analytically: (34) where k = √ M * 0 W . This expression, which is nonzero only for W > 0, agrees with the imaginary part of the expression for F (W ) in Eq. (23), which was obtained by taking the limit Γ * 0 → 0.
In Fig. 8, we compare the cross section for e + e − → Xγ in the triangle singularity region with the absorptive contribution obtained by replacing F (W ) in Eq. (18) by Im[F (W )] in Eq. (34). The absorptive contribution is zero below the D * 0D * 0 threshold. Unlike the complete cross section, the position of the peak of the absorptive contribution depends on the binding energy. For |E X | = 0.17 MeV, the position of the peak is about 1.3 MeV higher than that of the full cross section, and the height of the peak is about 58% of that of the full cross section. Thus the absorptive contribution is not a good approximation to the cross section for e + e − → Xγ in the triangle singularity region. At larger energies, the absorptive contribution quickly approaches the full cross section.
In Ref. [21], Dubynskiy and Voloshin derived an approximation for the absorptive contribution to the cross section for e + e − → Xγ using a Schrödinger wavefunction, as discussed in Section V. Their approximation for the imaginary part of the loop amplitude F (W ) in Eq. (13) can be expressed as where ψ DV is the DV wavefunction in Eq. (29). The integral over k in Eq. (35) can be evaluated analytically. If we replace ψ DV in Eq. (35) by the universal wavefunction ψ X in Eq. (26), the imaginary part of F (W ) becomes This agrees with the absorptive part in Eq. (34) if µ/M 0 = 0.518 is set to 1 2 . The corresponding integral with the DV wavefunction in Eq. (29) can be obtained by subtracting from the factor in square brackets the corresponding factor with γ X replaced by Λ and then multiplying by the first prefactor in Eq. (29): In Ref. [21], Dubynskiy and Voloshin estimated the peak in the absorptive cross section to be "numerically of the order of 1 pb". They used as input a measurement of the cross section for e + e − → D * 0D * 0 by the CLEO-c collaboration. They quoted the cross section as 0.15 nb at a center-of-mass energy 1.6 MeV above the D * 0D * 0 threshold. This measurement does not seem to appear in the conference proceeding they gave as a reference [37]. If we insert that cross section into Eq. (6), we get |A 0 | 2 + |A 2 | 2 = 410 GeV −2 . This differs only by a factor of about 1.5 from the value in Eq. (7) that we obtained from a fit to the cross section for e + e − → D * +D * − .

VII. SUMMARY
In this paper, we presented details of the calculation of the cross section for X(3872)γ from e + e − annihilation. A pair of P-wave spin-1 neutral charm mesons is created by the virtual photon from e + e − annihilation, and the charm mesons then rescatter into Xγ. The cross section for producing Xγ is given in Eq. (18), and the loop amplitude F (W ) is given in Eq. (22). We predicted the normalization of the cross section by using a previous fit to the cross sections for e + e − → D * D * from Belle data by Uglov et al. [23] to determine the amplitudes A 0 and A 2 in Eq. (7). The cross section has a narrow peak at an energy 2.2 MeV above the D * 0D * 0 threshold as shown in Fig. 6. We presented the cross section in this region in a previous paper [22]. The peak is caused by a charm-meson triangle singularity. The height of the peak is predicted to be between 0.2 pb and 0.9 pb if the amplitude ratio A 2 /A 0 is allowed to vary with |A 0 | 2 + |A 2 | 2 fixed at the value in Eq. (7). If the X is observed in the decay mode J/ψ π + π − , the cross section must be multiplied by the branching fraction Br for the bound state to decay into J/ψ π + π − . The loose lower bound Br > 4% and the upper bound Br < 33% on that branching fraction were derived in Ref. [33].
The Argand diagram for the loop amplitude from the triangle singularity is shown in Fig. 5. It is qualitatively similar to that of an ideal resonance, with the amplitude tracing out a counterclockwise loop in the complex plane as the energy W increases. Unlike the case of an ideal resonance, the loop is restricted to a singe quadrant of the complex plane.
The loop amplitude is expressed in terms of the Schrödinger wavefunction of the bound state with nonzero momentum in Eq. (27). Such an expression was previously used by Dubynskiy and Voloshin in their calculation of the production of Xγ [21]. Their wavefunction is a function of the relative momentum of the two constituents. In the Appendix, we derived the wavefunction for a bound state with nonzero momentum from the transition amplitude for the case of a near-threshold S-wave resonance. The wavefunction is a function of the relative velocity of the two constituents.
In the previous study of the production of Xγ from e + e − annihilation by Dubynskiy and Voloshin, they calculated only the absorptive contribution to the cross section. The comparison between the full cross section and the absorptive contribution is shown in Fig.  8. The absorptive contribution is not a good approximation near the triangle singularity region. In the absorptive contribution, the widths of the spin-1 charm mesons are set to 0. The narrow peak in the cross section comes from the triangle singularitiy that arises when all three charm mesons that form the triangle are on shell. The binding energy of the X prevents both of its constituents from being on shell simultaneously, and this was taken into account in DV. However the widths of the spin-1 charm mesons prevent them from being on shell, so the widths are also important. We calculated the cross section for e + e − → Xγ using the amplitudes from the charmmeson triangle diagrams in Fig. 4, which begin with the creation of a D * 0D * 0 pair. Dubynskiy and Voloshin have pointed out that there are also short-distance contributions to that amplitude that begin with the creation of DD, D * D , or DD * [21]. Those amplitudes will be essentially constant in the region of the peak from the triangle singularity. If the shortdistance amplitudes are larger than the amplitude from the triangle diagrams, the cross section near the D * + D * − threshold will have only a small peak or even a dip on top of a larger smooth background cross section. We have assumed the short-distance amplitudes are negligible compared to the amplitude from the triangle diagrams. Quantitative estimates of the short-distance amplitudes would be useful.
In our prediction for the cross section near the peak from the triangle singularity in Fig. 8, we assumed the X is a narrow bound state. A prescription for calculating the cross section in the case of a resonance with a different character is presented in the Appendix. As illustrated in Fig. 11, the peak from the triangle singularity in the case of a zero-energy resonance or a virtual state can be qualitatively similar to that in the case of a narrow bound state.
The BESIII collaboration has measured the cross section for Xγ from e + e − annihilation at center-of-mass energies ranging from 4.008 GeV to 4.6 GeV [24,25]. The cross section seems to have a broad peak near 4.2 GeV. The cross section has not been measured in the region near the D * 0D * 0 threshold at 4.014 MeV. The height of the narrow peak near 4.016 MeV from the charm-meson triangle singularity is predicted to be large enough that it could be observed by the BESIII detector. The observation of this peak would provide strong support for the identification of the X(3872) as a weakly bound charm-meson molecule.

ACKNOWLEDGMENTS
This work was supported in part by the Department of Energy under grant DE-SC0011726 and by the National Science Foundation under grant PHY-1607190.
FIG. 9. The 1PI transition amplitude A for the two particles with a near-threshold S-wave resonance. The connected 2 → 2 transition amplitude is the product of the 1PI amplitude A and a propagator for each of the 4 external legs. whose real and imaginary parts are We have denoted the real part of E pole by E X to distinguish it from the measured energy E X of the X in Eq. (8), which can be identified with the center of the resonance in the J/ψ π + π − channel. If Re[γ X ] > 0, the resonance is a bound state. It is a narrow bound state if |Re[γ X ]| is much larger than both Im[γ X ] and √ µΓ * 0 . If Re[γ X ] < 0, the resonance is a virtual state.
The behavior of the transition amplitude at complex energies E near the pole is This pole approximation for A(E) is a good approximation for the T-matrix element over some real range of the energy E only if the resonance is a narrow bound state. In this case, the residue of the pole in Eq. (A3) determines the vertex for the coupling of the bound state to a pair of particles: i √ 2πγ X /µ. This vertex can be used to calculate production rates of the bound state diagrammatically. In the case of the X, the resonance is in the S-wave channel for the superposition of charm mesons in Eq. (1). The vertices for the coupling of the X to D * 0D0 or to D 0D * 0 are therefore given by Eq. (10).

Short-distance production
If there is a reaction that can produce the two particles at short distances, the inclusive production rate of the two particles and their decay products can be determined by the optical theorem. The amplitude for the production of the two particles can be expressed as the product of A(E), which depends on the CM energy of the two particles, and a short-distance factor that is insensitive to E. We first ignore inelastic effects, and assume the T-matrix is exactly unitary. In the integral of the square of the amplitude over the momenta p 1 and p 2 of the two particles, it is convenient to change variables to the total momentum P = p 1 + p 2 and E: In the last step, we used the optical theorem for a 2-particle system that is exactly unitary. If the range of E is extended to negative values, the last expression in Eq. (A4) can take into account the production of bound states. We now consider inelastic effects. If they are taken into account in a way that ensures the positivity of Im[A(E)], the last expression in Eq. (A4) can be interpreted as a factor in the inclusive production rate of the resonance feature, which includes the entire contribution enhanced by the resonance. In addition to bound states, it includes final states from decays of the two particles as well as from their inelastic scattering.
If there is a near-threshold S-wave resonance, the resonance feature includes a peak above the scattering threshold called a threshold enhancement from the production of the pair of particles. The resonance feature may also include additional structure, such as a narrow peak below the threshold from a bound state or a peak near the threshold from a virtual state. The imaginary part of the simple model amplitude in Eq. (A1) at a real energy E can be expressed as The unitarity condition Im[A(E)] ≥ 0 requires the imaginary part of γ X to be positive. The first term in Eq. (A5) proportional to Im[γ X ] is the contribution from decays of the resonance into short-distance-decay channels, whose ultimate final states include particles with large momentum. In the case of X, they include J/ψ π + π − , J/ψ π + π − π 0 , J/ψ γ, ψ(2S) γ, and χ c1 (1P ) π 0 . The second term is the contribution from the constituent-decay channels. In the case of X, their ultimate final states are D 0D0 π 0 and D 0D0 γ. This contribution includes a threshold enhancement from production of a pair of constituents followed by the subsequent decay of D * 0 orD * 0 . The measured energy E X of the X in Eq. (8) can be identified with the center of energy of the Im[γ X ] term in Eq. (A5).
In the case of a narrow bound state, the short-distance production rate has a narrow peak below the threshold. The position of the peak, the center of energy E X of the resonance, and the real part E X of the pole energy in Eq. (A2a) are all approximately −Re[γ X ] 2 /2µ. By inserting the pole approximation in Eq. (A3) for the amplitude A(E) into Eq. (A5), the inclusive line shape near the peak can be approximated by where E X is given in Eq. (A2a). Inside an integral over the energy, this line shape can be further approximated by a delta function. Its coefficient can be expanded in powers of Im[γ X ]/Re[γ X ] and √ µΓ * 0 /Re[γ X ]. Up to relative corrections suppressed by three powers of 1/Re[γ X ], the line shape can be expressed as The square root reduces to |γ X | to leading order in Im[γ X ]/Re[γ X ]. If the resonance is not a narrow bound state, the real part E X of the pole energy in Eq. (A2a) is not a good approximation for the center of energy E X of the resonance. In the virtual-state limit where Re[γ X ] is negative and much larger in absolute value than Im[γ X ] and √ µΓ * 0 , both E X and E X are order Re[γ X ] 2 /µ, but E X is negative while E X is positive.
In Ref. [33], a simple model for the line shapes was used to illustrate different possibilities for the character of the resonance. The model for the transition amplitude is A(E) in Eq. (A1), which depends on Γ * 0 = 56 keV and the two adjustable parameters Re[γ X ] and Im[γ X ]. The inclusive production rate for the resonance feature is proportional to Im[A(E)] in Eq. (A5). The contribution from a specific short-distance decay mode i is γ i |A(E)| 2 , where γ i is an adjustable normalization factor proportional to the branching fraction into that decay mode. The total production rate of the resonance in that decay mode is obtained by integrating over the energy E. The integral over the energy range E min < E < E max depends logarithmically on the endpoints. To define a model with a finite production rate, we follow Ref. [33] in declaring the resonance to be the energy region between specified endpoints. Our model for the line shape in the short-distance decay mode i is therefore We also follow Ref. [33] in choosing the endpoints to be the D 0D0 π 0 threshold E min = −7.0 MeV and the D * + D − threshold E max = +8.2 MeV. We identify the measured energy E X of the X resonance in Eq. (8) with the center of resonance in the decay mode i = J/ψ π + π − . The center of energy E X for the line shape in Eq. (A8) can be defined by the condition If Im[γ X ] is fixed, Re[γ X ] can be adjusted to get a specified value of E X . The production rate in the decay channel i is proportional to γ i Emax E min |A(E)| 2 . This can be made independent of E X by adjusting the prefactor γ i as a function of E X . We choose a fixed value for the imaginary part of γ X : Im[γ X ] = √ µΓ * 0 = 7.4 MeV. We consider three cases with different values of the center of energy E X : • narrow bound state: E X = −0.30 MeV, which requires Re[γ X ] = 18.8 MeV. This is a little smaller than the value of γ X for a narrow bound state predicted from Eq. (11): 2µ|E X | = 24.1 MeV. We use this case to define a dimensionless normalization factor γ i = 1.
The line shapes for these three cases are illustrated in Fig. 10.

Bound-state wavefunction
In a quantum field theory, the Schrödinger wavefunction for a 2-particle bound state can be determined from the 2 → 2 transition amplitude for its constituents [38,39]. We take the particles with masses M 0 and M * 0 to have incoming momenta q 0 and q 1 , outgoing momenta q 0 and q 1 with q 0 + q 1 = q 0 + q 1 , and total energy E relative to the scattering threshold. The 1PI transition amplitude for a pair of particles with nonzero total momentum can be obtained from the 1PI transition amplitude in the CM frame in Eq. (A1) by replacing E by the Galilean-invariant combination of the total energy and the total momentum [31]: For simplicity of notation, we have set M * 0 +M 0 = M X . The connected 2 → 2 transition amplitude is the product of the 1PI amplitude A and a nonrelativistic propagator for each of the 4 external legs. The propagators for the incoming particles with energies E 0 and E 1 are i/(E 0 − q 2 0 /(2M 0 ) + i ) and i/(E 1 − q 2 1 /(2M * 0 ) + iΓ * 0 /2). If the external lines for one particle in the initial state and one particle in the final state are both amputated and if those lines are put on their energy shells, the residue of the pole can be factored into the product of a wavefunction that depends on the incoming momenta and a wavefunction that depends on the outgoing momenta [38,39]: The wavefunction is ψ(q 0 , q 1 ) = √ 8πγ X γ 2 X + µ 2 (q 0 /M 0 − q 1 /M * 0 ) 2 .

(A12)
This wavefunction is a function of the relative velocity of the constituents. Setting q 0 = q−k and q 1 = k, we reproduce the wavefunction in Eq. (28) for a bound state with momentum q and relative momentum k, except that Γ * 0 has been set to 0.

Resonance Feature from the Triangle Singularity
In our calculation of the loop amplitude F (W ) from the triangle singularity in Section IV, we assumed that the X is a narrow bound state with the sharp rest energy E X . This allowed the coupling of the X to a pair of charm mesons to be described by the momentumindependent vertex in Eq. (10). In the calculation of the cross section, the integral over the phase space includes an integral over a delta function at the sharp energy of the bound state. If the X is not a narrow bound state, the calculation of the cross section near the peak from the triangle singularity is more complicated. Whether or not the X is a narrow bound state, the distribution of the rest energy for the resonance feature is given by the factor of Im[A(E)] in Eq. (A4b). If X is a narrow bound state, the cross section for e + e − annihilation into X + γ near the D * 0D * 0 threshold is given in Eq. (18). It has a factor of |F (W )| 2 , where F (W ) is the loop amplitude in Eq. (22) for a bound state with the sharp rest energy E X . If X is not a narrow bound state, we should allow for the additional dependence on the energy E in the loop amplitude F (W, E). We should also replace the delta function in the rest energy of the bound state by the energy distribution proportional to Im[A(E)] in Eq. (A4b). The factor |F (W )| 2 in Eq. (18) can be replaced by where Im[A(E)] is given in Eq. (A5) and F (W, E) is given by the analytic expression in Eq. (22) with the coefficients a, b, and c replaced by The factor of 1/|γ X | in Eq. (A13) cancels the factor of | √ γ X | 2 from |F (W, E)| 2 , so the only dependence on γ X comes from Im[A(E)]. The variable q that appears in the coefficients b and c is the function of W − E that satisfies The variable q appears as a multiplicative factor in the analytic expression for F (W ) in Eq. (22). It also appears in the factor multiplying |F (W )| 2 in the cross section in Eq. (18). For those terms q, we can use the solution to Eq. (A15) with E = 0. The dependence of q on E is only essential in the argument of the logarithm in Eq. (22). In Fig. 11, we compare the cross section for Xγ near the peak from the triangle singularity for the three resonance cases itemized in Section A 2. The corresponding line shapes for the three resonance cases are illustrated in Fig. 10. The relative normalizations for the cross sections for the narrow bound state, the zero-energy resonance, and the virtual state are all determined by the condition that the three line shapes in Fig. 10 have the same area. The cross section for the narrow bound state has a shape very similar to that for the bound state with |E X | = 0.30 MeV in Fig. 6. The height of the peak in the cross section is a little larger for the zero-energy resonance than for the narrow bound state, and it is smaller for the virtual state. The peak for the zero-energy resonance and the virtual state is at a slightly lower energy and the full width at half maximum is smaller. However the three cross sections in Fig. 11 have roughly the same shape. Thus the peak from the triangle singularity cannot easily discriminate between a narrow bound state and other possibilities for the character of the resonance.