Dynamical implications of gluonic excitations in meson-meson systems

# Dynamical implications of gluonic excitations in meson-meson systems

## Abstract

We study meson-meson interactions using an extended basis that allows calculating coupling of an ordinary meson-meson system to a hybrid-hybrid one. We use a potential model matrix in this extended basis which at quark level is known to provide a good fit to numerical simulations of a system in pure gluonic theory for static quarks in a selection of geometries. We use a combination of resonating group method formalism and Born approximation to include the quark motion using wave functions of a potential within a cluster. This potential is taken to be quadratic for ground states and has an additional smeared (Gaussian) for the matrix elements between hybrid mesons. For the parameters of this potential, we use values chosen to 1) minimize the error resulting from our use of a quadratic potential and 2) best fit the lattice data for differences of and configurations of the gluonic field between a quark and an antiquark. At the quark (static) level, including the gluonic excitations was noted to partially replace the need for introducing many-body terms in a multi-quark potential. We study how successful such a replacement is at the (dynamical) hadronic level of relevance to actual hard experiments. Thus we study effects of both gluonic excitations and many-body terms on mesonic transition amplitudes and the energy shifts resulting from the second order perturbation theory (i.e. from the respective hadron loops). The study suggests introducing both energy and orbital excitations in wave functions of scalar mesons that are modelled as meson-meson molecules or are supposed to have a meson-meson component in their wave functions.

## 1. Introduction

Given the availability of both the numbers generated by lattice simulations of quantum chromodynamics and continuum models of the hadronic systems, an effective use of the numbers could be to improve the models through constraints of getting a least chisquare difference with the numbers for the corresponding discrete quarks and antiquarks configurations. Such lattice-improved models can then be reasonably used for all spatial configurations to eventually give dynamical predictions for experimentally measurable quantities like meson masses, meson-meson bindings and cross-sections and shifts (polarization potentials) to meson masses arising through meson-meson loops, etc. For one pair of quark and antiquark, a well established such use of lattice results is substituting in a Schrdinger equation a Coulombic-plus-linear quark-antiquark potential supported by lattice QCD calculations (see ref. [1, 2, 3] and others) for the ground state of the gluonic field between a quark and antiquark. Now that lattice results for excited state of the gluonic field are also available for years, even some dynamical uses of excited state gluonic field potentials have been worked out [4, 5, 6].

Uses worked out by others are either limited to numerical calculations without an explicitly written excited state gluonic field potential or the potential used originates from flux tube [7] or string models [8, 9]. Each of these approaches has its usefulness. What we add to this series of works is ourselves writing an analytical quadratic plus exponentially falling expression for the excited state gluonic field potential between a quark and antiquark and fit its parameters to the lattice data for the excited state gluonic field values available for discrete quark antiquark separations in [10]. This is reported in our previous work [6] as well. This work of us actually suggests and evaluates few other expressions for the excited state gluonic field potential as well. But the dynamical applications in it are present for a system whose valence quark contents are limited to one quark and one antiquark.

Through the present paper, we extend work on the dynamical implications of gluonic excitations to multiquarks that can be composed to more than one hadronic subclusters. For this extension, we need quark-level potentials that can model the more complicated gluonic field of this multiquark system. For this we combine the modelling of the spatial distribution of this gluonic field reported in ref. [11] (along with its fits to the continuum limits of the corresponding lattice simulations) with a realistic three colour structure. We have to do this combination because using all three colours the direct lattice simulations of all the Wilson loops relevant to a two quark two antiquark system are perhaps limited to ref. [12] and those mentioned in refs.[9-11] of ref. [13]. These works use a basis for a system that is truncated to the ground state of the glounic field. In comparison ref. [11] extends the basis to include the gluonic excitations and its description of the spatial distribution of the glounic field is more complete. But the lattice simulations in ref. [11] were carried out in a two-colour approximation to save computer time.

Our purpose here is to take advantage of the relatively complete basis and spatial distribution models of ref. [11] but using all three colours in the quark-level potential we use for our calculations of the dynamical hadron-level implications for a meson-meson system. In changing a two colour based model to a full three colour one, the number of colour basis states that interact (for any the spatial configuration) remain the same (see eqs. A.1 and A.3 of ref. [14], along with eqs. 5.1 to 5.4 and fig. 5.1 of ref. [15]) and we had to essentially only replace some colour overlap factors with proper values as elements of matrices of the same order. But the inter-quark elementary potential of ref. [11] turned out to be problematic for including the quark motion for the hadron level implications and thus we had to replace their numerically fitted ground state quark antiquark potential by a constant plus quadratic confining potential term and the additional potential for the gluonic excitation by one of the form . As written in ref. [6], only can be used in solving the integrals of our present work analytically although few other forms for the gluonic potential are also suggested in ref. [6] with less as compared to . Only then we are able to perform a full meson level dynamical calculations for transition amplitudes from one set of quark-antiquark clusterings (mesons) to the other. Using these amplitudes we are also able to study certain properties of the polarization potentials for a meson-meson system.

We had to use simpler interquark potentials that can be symbolically integrated at a later stage, after necessary multiplications by wave functions of positions, to complete the adiabatic-approximation-based treatment of a system mentioned below. Being not limited by such demands of later integrations, the form and then parameters values of the continuum model proposed in ref. [11] were chosen to simply minimize

 χ2A=1N(G)N(G)∑i=1(Ei−Mi)2/△Ei, (1)

where is the number of data points for geometry G. The geometries in ref. [11] numerically worked on were (quarks at the corners of) squares, rectangles, tetrahedra and some other less symmetric geometries and (linear). For each data point , the lattice energy was extracted by solving the following eigenvalue equation (see eq. 2 of ref. [16])

 WTiklaTil=λ(T)iWT−1iklaTil (2)

for that approaches to as (Euclidean time) , and then subtracting the energy of two separated clusters from the to get the lattice-generated binding energy for the data point. The values of the , and thus of , depend only on the numerical values of the elements of the matrix of the Wilson loops. The values of and depend on the number of Wilson loops evaluated; for the system these were taken to be 1 and 2. Two of the corresponding Wilson operators (whose vacuum expectation values are the Wilson loops ) are shown for example in Fig. 1.5 of  [14]. Knowing the Wilson loops, the procedure of getting can be found for example in eqs. 4, 11, 12 and 15 of ref. [16]. The arguments for continuum limits being achieved before extracting are given in ref. [17].

The in eq.(1) are obtained by subtracting the energy of two separated clusters from the eigenvalues of a matrix obtained through a model of the system. For this, the model has to give a basis and an operator whose representation with respect to the basis gives a potential matrix . are obtained by setting the determinant of equal to zero, with being the energy of two separated clusters and the (overlap) matrix of an identity operator in the basis. Searching for the model, the simplest way to extend a two-particle potential model to a few-body is to use the potential for each pair of particles in the few-body system and simply add up such two-body potentials. This approach has been successful in atomic and many-nucleon systems; the corresponding two-body interaction being described by Coulombic and Yukawa potential, for example. For a hadron (or a system of hadrons) composed of many quarks, antiquarks and the gluonic field, the lowest order perturbative Feynman amplitudes are of this sum-of-pair-wise form. Though Feynman diagrams themselves become impractical for typical hadronic energies because of larger couplings, models have been tried which simply replace the two-body Coulombic potential (essentially a Fourier transform of the Gluonic propagator) by more general Coulombic-plus-linear form; see ref. [18]. This approach is not free of problems; for example it leads to inverse power van der Waals’ potentials [19] between separated colour-singlet hadrons which are in contradiction with experimental data. But this model has many phenomenological successes and it is worthwhile inquiring if

1) it provides a basis and operator to generate a potential matrix, and

2) how good is the chisquare if the eigenvalues of the resulting matrix are used as in eq.(1).

The answer provided by ref. [11] and earlier related works is that the model does generate a matrix of the required kind. But the resulting chisquare, defined by eq.(1), is too bad; see Fig. 4 of ref.[20]. To refine the model we can improve the basis beyond the defined as  [11, 21]

 |1⟩=(q1q3)(q2q4),|2⟩=(q1q4)(q2q3), and|3⟩=(q1q2)(q3q4), (3)

and the operator beyond

 H=−4∑i=1[mi+ˆP2i2mi]+∑i

with being is the potential energy of a pair with the gluonic field between them in the ground state. Or both the basis and the operator can be improved. What ref. [11] does is to improve directly the matrix (representation) after writing down the underlying basis. The authors do this in a number of ways. One model, termed model II, uses the same basis but multiplies the off-diagonal elements of the overlap and potential energy matrices (that is, the representations of the identity and potential operators

 ∑i

respectively) by a few-body gluonic field overlap factor with as the tension of the string connecting a quark with an antiquark, the area of a surface bounded by external four lines connecting two quarks and two antiquark and approximately; theoretical arguments suggest should be the area of the corresponding minimal surface, though in ref. [11] half of a sum of four triangles was used for numerical convenience. At the quark level, this model II was noted to much reduce the chisquare of eq.(1). This model has been worked out in  [22, 23, 24] till meson-level transition amplitudes. The dynamical calculations require a kinetic energy term as well. As this is taken to be apart from some technical considerations of hermicity, proportional to the overlap matrix and hence its off-diagonal elements are also multiplied by the overlap factor. Thus provides one parametrization that connects QCD simulations with hard experiments.

But model II is not the best model of ref. [11]; the paper continues to then improve the basis by including the gluonic excitations as well. That is, it extends the basis by including the states

 |1⋆⟩=(q1q3)g(q2q4)g,|2⋆⟩=(q1q4)g(q2q3)g, and|3⋆⟩=(q1q2)g(q3q4)g. (5)

Here denotes a state where the gluon field is excited to the lowest state. (The excited states of gluonic field can, for example, be seen in the QCD numerical simulations; see ref.[25, 10] and others). When the overlap, potential and kinetic energy matrices are written in this extended basis, their order increases to rather than previous . If in addition, we introduce many body terms in this extended model, new kind of gluonic field overlap factors (, ) appear in the off diagonal terms resulting in what ref. [11] terms model III giving the least chisqure in ref. [11]; (see eq.(15) below); our truncation to matrices is explained before this equation. As mentioned above, the purpose of our present paper is to work out this improved model III of ref. [11] till the meson-meson scattering amplitudes and energy shifts. As this improved model III includes the gluonic excitations, it consider transitions from three ground state quark states to the ones having gluonic excitations. And by adding to it the quark motion (wave functions) to reach the hadron level, we are now able to study transitions from ground state meson-meson systems to hybrid-hybrid systems.

A worth-mentioning aspect we have studied is the hadron-level implications of the differences of the gluonic-excitation-including model III and the sum-of-pair-wise approach. The quark-level work in ref. [11] can be interpreted to mean that with suitable extensions and modifications something like a sum-of-pairs approach can be a good approximation for a tetraquark system as well. It states ”At the shortest distances, upto about 0.2fm, perturbation theory is reasonable with the binding being given mainly by the , , and states inter-acting simply through the two-quark potentials with little effect from four-quark potentials”. However, ref. [11] states, ”for large inter-quark distance (greater than 0.5 fermi), quark-pair creation can no longer be neglected. However, in the intermediate energy range, from about 0.2 to 0.5 fm, the four-quark potentials act in such a way as to reduce the effect of the , , and states so that the binding is dominated by the , , and states, which now interact among themselves again simply through the two-quark potentials with little effect from four-quark potentials.” This suggests that models involving only two-quark potentials could be justified provided excited gluon states (such as , , ) are included on the same footing as the standard states , ,. We have checked if such features survive at the experimentally meaningful hadron level, by comparing the dynamical implications of

(1) a model extended to the gluonic excitations but otherwise sharing many features with the sum-of-pair-wise approach, with

(2) a model that includes explicit many-body terms but does not include gluonic excitations.

Thus we report if after including the gluonic excitations a sum-of-two-body potential model can replace to some extent many-body potential terms in a tetraquark system even at a hadronic level. Specifically, we have calculated in both kind of models meson-meson transition amplitudes from () a ground state meson-meson clustering to () a different ground-state clustering and to a clustering of gluonic-excited mesons. are elements of the meson-meson scattering theory -matrix, can be termed as phase shifts, transition potentials or meson-meson coupling, and their absolute squares give meson-meson differential cross sections [18]. Moreover, using these transition amplitudes in the second order perturbation theory, we study shifts (in both kind of models) in a ground state meson-meson energy due to coupling to a different ground-state clustering and to an excited state meson-meson system (i.e. to hybrid loops). These energy shifts are what are also termed polarization potentials  [26].

To reach the hadronic level, we have included the quark motion through quark wave functions. To find the quark position dependence in multiquark systems, a number of methods are used such as variational method [27, 28], Born-order diagrams [18], and resonating group method [29]. Variational approach is used by Weinstein and Isgur to optimize a meson-meson wave function in a quadratic  [27] and later coulomb plus linear [28] potential and a hyperfine term combined with a sum-of-two-body approach. Then they projected the meson-meson state onto free meson wave functions to estimate a relative two meson wave function which gives the equivalent meson-meson potential and obtained the meson-meson phase shifts. The similar results can be obtained by using Born-order quark exchange diagrams  [18] in a non-relativistic potential model to describe low energy scattering of mesons.

In the present paper, we have used a formalism of the resonating method as used in [22]. In the resonating group method, the dependence on the internal co-ordinates of the system is specified before solving the problem to integrate out the degrees of freedom corresponding to the internal coordinates of clusters of the system. At a later stage, because of the complexity of the calculations we use a Born approximation to specify the dependence on the vectors connecting the centers of masses of our mesonic clusters as well. Moreover, we have not included in our basis an explicit diquark-diantiquark state. In the weak coupling limit such a state is a linear combination of the meson-meson states and thus cannot be included in a basis. Away from the weak coupling this can be included. But inclusion of a third clustering states in the basis did not affect the lattice simulations reported in ref. [11]. Later works  [13, 12] using were done with only two clusterings. Thus we have expanded our two-quark two antiquark wave function in a basis that is limited to four meson-meson states: ground and excited states for each of the two possible clusterings.

In this exploratory work, we have taken all the constituent quark masses to be same as that of a charm quark and we have taken all the spin overlaps to be 1 without calculating them. Actually we have not studied the effects of the quark degree of freedom, not its orbital angular momentum as well; we want to study only dynamical effects of including the gluonic excitations in our basis for the same quark configurations. In section 2, the potential model in the extended basis in the pure gluonic theory is introduced for the static quarks. Basically, in this section we tell where does the model of ref. [11] fits in our full scheme that incorporates the quark motion through a resonating group method formalism meaning pre-specifying quark-antiquark wave functions within clusters. The coupled integral equations for the remaining inter-cluster wave function are written in section 3. In section 4, these integral equations are solved to calculate the transition amplitudes and energy shifts we are studying. The numerical results for meson-meson system transition amplitudes and energy shifts in both the above mentioned kinds of models with concluding remarks are given in section 5. The partial wave analysis results of transition from ground state to excited state gluonic field are also reported in section 5.

## 2. q2¯¯¯q2 Potential Model (In the Extended Basis)

Using the adiabatic approximation, the total state vector of a system containing two quarks, two anti-quarks and the gluonic field between them can be written as sum of product of quark position dependence function and the gluonic field state . (The gluonic state is defined as a state approaching to colour state in the limit of quark anti-quark separation approaching to zero. Here .) The function can be written as

 ϕK(r1,r2,r¯3,r¯4)=ϕK(Rc,RK,yK,zK),

with . is the overall center of mass co-ordinate of the whole system.

With the notation of Fig. 1, the relative co-ordinates and are defined as

 R1=12(r1+r¯3−r2−r¯4),y1=r1−r¯3,z1=r2−r¯4 R2=12(r1+r¯4−r2−r¯3),y2=r1−r¯4,z2=r2−r¯3 R3=12(r1+r2−% r¯3−r¯4),y3=r1−r2,z3=r¯3−r¯4,

being the vector joining the centers of mass of the mesonic clusters and ; similarly about and . Now using the resonating group method, the quark position dependence function can be written as a product of function of known dependence on ,, and of unknown dependence on . i.e.

 ϕK(Rc,RK,yK,zK)=ψc(Rc)χK(RK)ψK(yK,zK). (6)

Thus, the two quarks two antiquarks state vector can be written as

 ∣Ψ(q1q2q¯3q¯4)⟩=∑k∣k⟩gψc(Rc)χk(RK)ψk(yK,zK), (7)

where

 ψk(yK,zK)=ξk(yK)ξk(zK),

and being the normalized solutions of the Schrdinger equation for quadratic confining potential (written in eq.(10)) for a pair of quark-anti-quark within a cluster. We take, for the zero relative orbital momentum () of a quark w.r.t. the antiquark of the cluster,

 ξK(yK)=1(2πd2)34exp(−y2K4d2),ξK(zK)=1(2πd2)34%exp(−z2K4d2). (8)

Here is the size of meson (detail is written after eq.(10)) and being the constitute quark mass. In our case, is the mass of c-quark and equal to GeV as used in  [30].

After writing the form of the wave vector, we describe our Hamiltonian, starting with the limit when each gluonic field overlap factor . In this limit, the Hamiltonian whose representation matrices in the basis {} would become those in ref. [11] is

 H=−4∑i=1[mi+ˆP2i2mi]+∑i

where is the potential energy of meson for the ground state gluonic field and is the difference between ground state and excited state gluonic field potential. We take the kinetic energy in the non-relativistic limits. This limit is also used in a recent work by Vijande ref. [31] that deals with multiquark system (two quarks and two antiquarks) to study the spectrum using a string model for the potential. In ref. [11], potential energy matrix elements are written so that the potential energy for each pair is equal to for the matrix elements of the Hamiltonian between gluonic ground states, and it is equal to for the matrix elements between the gluonic-excited states. We have modeled these two forms by taking for the ground state matrix elements and for the elements between gluonic-excited states. For the elements between ground and excited state gluonic field, the value that results from their parameter being fitted to 4 (in their Table 1) is surprisingly and not any value between 0 and 1. A possibility is that this is a result of them taking the area , we mentioned above in our introduction, in the form of average of the sum of triangle areas instead of the theoretically motivated minimal surface area. Thus we have somewhat explored between 0 and 1 and in addition to which we have mainly studied.

In the above equation, (operating on particle) has eight components with . Each component is equal to , where are the Gell-Mann matrices. is used in the superscript to avoid the possible confusion with index . We used the potential with the colour structure of one gluon exchange in the form given in ref. [11]. With the use of ground state potential in the realistic coulombic plus linear form, it becomes impossible for us to solve the integral equations appearing below in eqs.(21-24). Therefore we used the parametrization of the static pairwise two quark potential as

 vij=Cr2ij+¯¯¯¯C,withi,j=1,2,¯¯¯3,¯¯¯4. (10)

In this simple harmonic potential, the parameters is related to size () of wavefunction () through the relation with , and for consistency of the diagonal term of the integral eq.(21) GeV [22] with GeV being the mass of a charmonium meson. The parameter is chosen in such a way to reduce the error resulting from a use of this quadratic potential instead of the realistic one. The error may be both in the wave functions of the distance between a quark and antiquark within clusters and those of the vectors joining the centers of masses of the two clusters. As for the first dependence, we found that the maximum of the overlap integral of each of the wave functions and of the quadratic potential and that of a more realistic coulombic plus linear potential () is at . (For the parameters of the realistic potential we used values  [30] and  [30] for mesons composed of charm quarks.) This overlap is shown in Fig. 2. A similar work was done in ref. [32] for lighter quarks. They found that the overlap wave function of SHO (quadratic potential) and that of coulombic plus linear can be made as large as with the suitable adjustment of parameters.

For the additional term in the potential for the gluonic excitation, the usual flux tube ( [7]) or string based analytical expressions become impractical for us, as mentioned in the introduction. Thus for that we tried an of the form of

 △v⋆ij=Ae−Br2ij. (11)

This gaussian gluonic potential () is a smeared form of as written in appendix of ref. [28]. From the Fig. 3 of [10], we get the potential energy difference between ground and excited states for different r values (). We choose A and B for which becomes minimum. is defined as

 χ2=n∑i=1(εi−Aexp[−Br2ij])2,

with being the number of data points. This gives

 A=1.8139GeV,B=0.0657GeV2.

For finding the wave function corresponding to our total potential , we used the variational method with an wave function

 ξ⋆K(yK)=ny2Kexp(−py2K). (12)

The normalization of this w.r.t gives

 n=(4234p74)(1512π34).

This leaves us with one variational parameter chosen to minimize the expectation value of the two body Hamiltonian in the excited state gluonic field wave function. This gave . For this value of , the overlap of wave function of the quadratic potential plus and that of coulombic plus linear plus within a hybrid cluster became . Both wave functions are shown in in Fig. 3. Having much reduced the errors in the in-cluster factors of the total wave function, the question remains how much the inter-cluster factors of the (terms of the) total state vector are affected by our use of convenient but not realistic potentials. For the inter-cluster wave functions, eventually we use below in eq. (25) plane wave forms which get their justifications from the validity of Born approximation for our problem regardless of potential expressions we use. This plane wave form has only one usual parameter (the wave number) and eq.(C10) below relates its value for the ground as well as excited state inter-cluster waves to the very good values of and that almost give realistic ground state and excited state wave functions within clusters. But the relations between the inter-cluster wave numbers and the and assume a quadratic confinement and this may affect our numerical results but hopefully not at least the qualitative features we are pointing out. Perhaps it is worth mentioning here that properties of systems were calculated using quadratic confinement in ref. [27], and then with the realistic potential in ref. [28] and both the works favoured the existence of meson-meson molecules.

Now we combine our Hamiltonian and all the wave functions we have mentioned in the Schrdinger equation for the meson-meson system, which means that the overlap of with an arbitrary variation of state vector vanishes where is the state vector of the whole system. In , we considered only the variation in (see eq.(7)), as in the resonating group method. Thus we wrote

 ⟨δΨ∣H−Ec∣Ψ⟩=∑k,l∫d3Rcd3RKd3yKd3zKψc(Rc)δχk(RK)ξk(yK)ξk(zK)g⟨k∣H−Ec∣l⟩gψc(Rc)χl(RL)ξl(yL)ξl(zL)=0 (13)

for and The arbitrary variations ’s for different values of are linearly independent and hence their co-efficient in eq.(13) should be zero. With the trivial integration performed to give a finite result, this leads to

 ∑l∫d3yKd3zKξk(yK)ξk(zK)g⟨k∣H−Ec∣l⟩gχl(RL)ξl(yL)ξl(zL)=0, (14)

where

 g⟨k∣H−Ec∣l⟩g=g⟨k∣KE+V+4m−Ec∣l⟩g.

Elements of and matrices are defined below in eq.(17) and eq.(18). In ref. [11] it is stated that, in the lattice QCD simulations, it was found that the energy of the lowest state was always the same in both a and description, provided or had the lowest energy. In addition the energy of the second state was, in most cases, more or less the same. Two level approximation is also used in later works [13, 12] of the tetraquark system. Considering this, we include only two topologies (1 and 2), meaning four states. According to the model of ref. [11] the overlap matrix N in this truncated 4-basis is given by

 N={g⟨k∣l⟩g}=⎡⎢ ⎢ ⎢ ⎢⎣1f/30−fa/3f/31−fa/300−fa/31−fc/3−fa/30−fc/31⎤⎥ ⎥ ⎥ ⎥⎦. (15)

And the potential matrix is

 V={g⟨k∣V∣l⟩g}=⎡⎢ ⎢ ⎢ ⎢⎣V11V12V11⋆V12⋆V21V22V21⋆V22⋆V1⋆1V1⋆2V1⋆1⋆V1⋆2⋆V2⋆1V2⋆2V2⋆1⋆V2⋆2⋆⎤⎥ ⎥ ⎥ ⎥⎦. (16)

Here,

 V11=−43(v1¯3+v2¯4)V12=V21=49f(v12+v¯3¯4−v1¯3−v2¯4−v1¯4−v2¯3)V22=−43(v1¯4+v2¯3)V21⋆=V2⋆1=−fa18(√2(v⋆1¯3+v⋆2¯4)−16√2(v⋆1¯4+v⋆2¯3)−√2(−v⋆12−v⋆¯3¯4))V1⋆2=V12⋆=−fa18(√2(v⋆1¯4+v⋆2¯3)−16√2(v⋆1¯3+v⋆2¯4)−√2(−v⋆12−v⋆¯3¯4))V1⋆1⋆=16(v⋆1¯3+v⋆2¯4)V1⋆2⋆=V2⋆1⋆=−118fc(−(v⋆1¯3+v⋆2¯4+v⋆1¯4+v⋆2¯3)+10(v⋆12+v⋆¯3¯4))V2⋆2⋆=16(v⋆1¯4+v⋆2¯3)V1⋆1=V11⋆=V2⋆2=V22⋆=0, (17)

with , being defined above (after eq.(9)). The coefficients of and , resulting from the operator, are given in Table 1 in the Appendix A. The kinetic energy matrix of the two quarks two anti-quarks is taken to be

 KE={g⟨k∣KE∣l⟩g}=\emphN(f)12k,l(−12m4∑i=1∇2i)\emphN(f)12k,l. (18)

The kinetic energy in the same form is also used in ref. [22].

The gluonic field overlap factor f, as written in introduction, is suggested by ref. [11] as

 f=exp[−bskfS], (19)

with  [7]. In ref. [13] the gluonic field overlap factor is used in the Gaussian form as

 f=exp[−kfbs∑i

employed in for interpreting the results in terms of the potential for the corresponding single heavy-light meson. In ref. [13], the simulations that are fitted by using are for the configurations when the gluonic field is in the ground state i.e. overlap matrix is a matrix. In ref. [11], simulations are reported with 2-colour approximation. But in ref. [13], lattice simulations performed for are reported. Our overlap, potential and kinetic energy matrices are also written in , so we use the multiplying sum of area form of (written in eq.(20)) with (as used in  [13]) for numerical convenience and have not used the minimal area form. When we are observing the dynamical effects for the ground state, our overlap, potential and kinetic energy matrices are matrices and we use with . But when we incorporate the excited state gluonic field our overlap, potential and kinetic energy matrices become matrices, in the upper left block of these matrices, the form of remains the same but the value of is changed to according to the conclusion of ref. [11]. In the other blocks , are also used and defined in refs. [11] [21] as

 fa=(fa1+bsfa2S)exp[−bskaS],
 fc=exp[−bskcS].

If we take as a function of area as defined in  [11], it becomes unmanageable to solve the integral equations (21-24) and hence we have taken to be a constant but have tried a variety of its values to explore how much our conclusions depend on its value. As for , the fit in ref. [21]  [11] of the model to the lattice data favours which implies that i.e. the excited configurations interact amongst themselves in the way expected from perturbation theory. Thus we have used .

## 3. Coupled Integral Equations

Using ”N”,”V” and ”KE” elements in eq.(14), we got four integral equations for four different values of k or l. Then we do and integrations. All the integrations required are in the Gaussian form or modified Gaussian form (with a polynomial in the integrand multiplying the Gaussian exponential) and we integrate analytically. For , in eq.(14), is independent of and and, thus, can be taken out of integrations. After the integration, the result is dependent co-efficient of . For , and are replaced by their linear combinations with one of them as identical to and other one independent of it as . The jacobian of transformation is equal to 8. Then we integrate the equation w.r.t . Integration leaves the following four equations:

 Δ1(R1)χ1(R1)+∫d3R2F(R1,R2)χ2(R2)+∫d3R2E1(%R1,R2)χ⋆2(R2)=0, (21)
 Δ2(R2)χ2(R2)+∫d3R1F(R1,R2)χ1(R1)+∫d3R1E1(%R2,R1)χ⋆1(R1)=0 (22)
 F4(R1)χ⋆1(R1)+∫d3R2E2(R1,R2)χ2(R2)+∫d3R2E3(% R1,R2)χ⋆2(R2)=0 (23)
 F4(R2)χ⋆2(R2)+∫d3R1E2(R2,R1)χ1(R1)+∫d3R1E3(R2,% R1)χ⋆1(R1)=0. (24)

The symbols are defined in the appendix B. We have eventually replaced by , and . With trivial integration on , we have eq.(14) that is independent of . Now, after the integration on , the above four integral equations (21-24) depend only on and . So every quantity which we want to calculate depends on and . In eq.(21-22), the first two terms in each equation, containing and , are for the ground state. It is noted that in these terms (observing the definitions of the symbols in appendix B), there is no dot product of vectors and . So the results from these terms should not depend on the angle between and , called . The third term in each of eq.(21-22) is due to the gluonic ground and excited states. In these terms dot product of two vectors ( and ) appear, so the results from these terms depend on .

## 4. Solving the Integral Equations

Now taking the three dimensional Fourier transform of eq.(21,23) with respect to and eq.(22,24) with respect to , we get formal solutions , , , and as shown in appendix C. Because of the coupling to the gluonic excitations, it become difficult to solve the integral equations for non trivial solutions for , and analytically as done in  [22, 23]. In  [22], the meson wave functions, including the gluonic field overlap factor, is , separable. So there the integral equations can be solved analytically by replacing and . But in our present work, the meson-meson wave functions are not separable in ,. So we use the Born approximation (as used in  [18] for meson-meson scattering) to solve the integral equations. Our results given below also justify our use of the Born Approximation. Using this approximation, we use the solutions () of eqs.(21-24) in absence of interactions (meaning )

 χi(Ri)=√2πexp(ıPi.Ri),andχ⋆i(Ri)=√2πexp(ıPis.Ri) (25)

for . Here the coefficient of is chosen so that it makes as Fourier transform of . Similarly the coefficient of is chosen. Using this approximation, the integration on and can be performed to get (written in eq.(C11)).

can be calculated (as in ref. [22]) by considering the coefficient of containing the from eq.(C11). As in this equation, there is no coefficient having , so it gives . can be calculated by considering the coefficient of containing the from eq.(C11) in the following eq.

 T12=Mπ2P1[coef.of1△1(P1)containingχ2(R2)], (26)

with being the mass of meson. Similarly can be calculated by considering the coefficient of containing the from eq.(C11)

 T12⋆=Mπ2P1s[coef.% of1△1(P1)containingχ⋆2(R2)]. (27)

The relation between off-diagonal transition and scattering matrix element is written as

 S=I−2ιT, (28)

where , , and represent scattering, transition, and identity matrices respectively. The eq.(28) can also be written as

 Sij=δij−2ιTij. (29)

for . Using the transition matrix elements, the contribution to the energy shift of meson-meson system () through states (both with and without gluonic excitations) can be calculated by using the stationary state perturbation theory, i.e.

 Ei=E0i+Tii+∑i≠j∫∞0|Tij|2E0i−E0jdPj, (30)

with the initial state and intermediate state . We have considered initial states where the gluonic field should is in ground state, so , but intermediate . Here , , is the energy of a ground state (1 or 2) of meson-meson system, and may be the energy of the other meson-meson ground state or that of an gluonic-excited meson-meson state.

## 5. Results and Conclusions

1-The transition amplitude , from one meson-meson ground state to other, is calculated by using eq.(26) with  [13](without the incorporation of gluonic excited states). Its dependence on the center of mass kinetic energy is shown in Fig. 4 below. As it is noted that magnitude of transition element is less than 1, so this result shows the validity of Born approximation. In result 3 we compare the transition amplitude of this many body ground state gluonic field model at with the transition amplitude obtained from a model that is extended to gluonic excitations along with changing .

2- The transition matrix element , for transition from ground state to excited state gluonic field with , depend on the parameters , and . also depends on (the angle that makes with and makes with ). We take parameter as a constant as discussed earlier in section 2. For , we take different values of to see the effects of on . Fig. 5 shows the dependence of on at . By taking different values of , was calculated. We found that the results are slightly different for different . This slight angle dependence is not directly reported here, but can be found by a linear combination of the corresponding spherical harmonics with coefficients for each value of energy read from Figs. (8-10) that report the partial wave amplitudes that result from this angle dependence. For , the dependence of on is shown in Fig. 6. And for , the dependence of on is shown in Fig. 7. These graphs show that the magnitude of transition amplitude is increasing with the increase of