Doubly charmed baryon production in heavy ion collisions

We give an estimate of $\Xi_{cc}^{++}$ production rate and transverse momentum spectra in relativistic heavy ion collisions. We use Boltzmann transport equations to describe the dynamical evolution of charm quarks and diquarks inside quark-gluon plasma. In-medium formation and dissociation rates of charm diquarks are calculated from potential non-relativistic QCD for the diquark sector. We solve the transport equations by Monte Carlo simulations. For $2.76$ TeV Pb-Pb collisions with $0-10\%$ centrality, the number of $\Xi_{cc}^{++}$ produced in the transverse momentum range $0-5$ GeV and rapidity from $-1$ to $1$ is roughly $0.02$ per collision. We repeat the calculation with a melting temperature $250$ MeV above which no diquarks can be formed. The number of $\Xi_{cc}^{++}$ produced in the same kinematic region is about $0.0125$ per collision. We discuss how to study diquarks at finite temperature on a lattice and construct the anti-triplet free energy in a gauge invariant but path dependent way. We also comment on extensions of the calculation to other doubly heavy baryons and doubly heavy tetraquarks and the feasibility of experimental measurements.

Recently the LHCb Collaboration reported the observation of a doubly charmed baryon carrying two units of positive charge, Ξ ++ cc , with a mass m(Ξ cc ) ≈ 3621 MeV [1]. Though it is still unclear why the observed mass differs from the previous SELEX result [2], the existence of such doubly charmed baryons is now on a more solid ground. The particle is stable under strong interactions and only decays weakly. The structure of Ξ ++ cc can be thought of as an up quark bound around a deeply bound state (diquark) of two charm quarks [3]. Just as a pair of heavy quark and heavy anti-quark attract each other and can form a bound state in the color singlet channel, a pair of two heavy quarks also attract and can form a bound state, a heavy diquark, in the anti-triplet representation.
The peculiar properties of Ξ ++ cc have stimulated new theoretical and experimental research. Here we consider the production of Ξ ++ cc in high energy heavy ion collisions, where a hot nuclear environment, the quark-gluon plasma (QGP), is produced. Previous work was based on quark coalescence at hadronization and assumed that heavy quarks are thermally distributed [4,5]. Here we pursue out a more dynamical approach considering the formation of bound heavy diquarks within the quarkgluon plasma and the incomplete equilibration of the heavy quark spectrum.
In hadron-hadron collisions, it is difficult to produce a pair of heavy quarks in the color anti-triplet at leading order in a fragmentation process. On the other hand, the coalescence process involving two independently produced charm quarks is sensitive to the relative momentum between the heavy quark pair. In proton-proton collisions, the relative momentum is uncontrolled and likely large, suppressing the coalescence. Heavy ion collisions have two advantages for Ξ ++ cc production: First, the rapidity density of charm quarks produced in a single colli- * Electronic address: xiaojun.yao@duke.edu † Electronic address: mueller@phy.duke.edu sion is higher. Second, the deconfined QGP medium lasts roughly 10 fm/c, during which time the charm quarks can diffuse in the QGP via interactions with light quarks and gluons. This is confirmed by recent measurements from the STAR Collaboration, which shows that charm quarks participate in the collective flow of the QGP [6]. As a result, the relative momentum of a charm quark pair can be on the order of the QGP temperature. The coalescence probability into a charm diquark bound state is thus enhanced if the temperature of the QGP is not too high. After its formation the charm diquark also diffuses in the QGP because it carries color charge. At the same time, the charm diquark may dissociate by absorbing a real or virtual gluon. So the whole process is a dynamical in-medium evolution involving charm diquark formation, diffusion and dissociation. This is similar to the in-medium evolution of heavy quarkonia, such as the J/ψ, except that the heavy diquarks carry color while the quarkonia are color neutral. At the transition from the deconfined QGP phase to the hadronic phase, the charm diquarks hadronize into doubly charmed baryons by absorbing an up or down quark from the medium.
We will describe the in-medium dynamical evolution of charm quarks and diquarks by a set of coupled Boltzmann equations analogous to the transport equations for in-medium heavy quarks and quarkonia [7]. By connecting the transport equations with the initial production of charm quarks from the hard collision and the hydrodynamical background, we obtain an estimate of the yield and p T -spectrum of Ξ ++ cc in Pb-Pb collisions at 2.76 TeV. Finally, we study the static screening effect of the QGP on the production process.
The set of coupled Boltzmann transport equations for the charm quark and diquark distribution functions f (x, p, t) is given by where all the collision terms C, C ± depend on x, p, t.
Here we will focus on the ground charm diquark state cc(1S) because excited states are loosely bound and cannot survive at high temperature. In the following, by charm diquark we mean the cc(1S) state. The collision terms C c and C cc describe their scattering with thermal constituents of QGP. This process has been described as two-body scattering in the framework of the linearized Boltzmann equation [8][9][10]. Here we use the elastic scattering rate calculated and implemented in Ref. [11] to describe the in-medium diffusion. The diquark gain term C + cc is from the combination of a charm quark pair by gluon emission and the loss term C − cc is from dissociation by gluon absorption. The formation and dissociation of diquarks also change the charm quark distribution function, which are represented by C ± c . We calculate the diquark formation and dissociation rates in QGP to the lowest order in potential nonrelativistic QCD (pNRQCD) for the diquark sector [12,13]. The pNRQCD for the quarkonium sector has been used to study quarkonia dissociation rates inside QGP [14]. The effective field theory can be derived from QCD under the hierarchy of scales M M v M v 2 , T, m D where M = 1.3 GeV is the charm quark mass, v ∼ 0.4 is the relative velocity of cc inside the diquark, T is the QGP temperature, and m D is the Debye screening mass. If T or m D scales as M v, the Debye static screening of the color attraction is so strong that no diquark bound states can be formed inside QGP. So the above hierarchy of scales is relevant to the diquark formation. The pNRQCD is a systematic expansion in v or 1/M (NR expansion) and r, the relative distance between the charm quark pair inside the diquark (multipole expansion). Its Lagrangian is given by: where higher order interaction terms in 1/M and r are omitted. The Lagrangian of light quarks and gluons is just QCD with momenta k M v. The degrees of freedom are the anti-triplet T(R, r, t) and sextet Σ(R, r, t) where R denotes the center-of-mass (c.m.) position and r the relative coordinate. They are defined as where T l and Σ ν are the anti-triplet and sextet fields while t l and σ ν are the generators of the corresponding representations. They are given by The equations of motion of the anti-triplet and sextet are Schrödinger equations with the Hamiltonians expanded in powers of 1/M where D R is the covariant derivative associated with the c.m. position. By the virial theorem, −∇ 2 T,Σ . So the order of the relative kinetic term is accounted as 1/M 0 , not suppressed. The c.m. kinetic term is suppressed because momenta k ∼ M v have been integrated out in the construction and then D R M v. Higher-order terms of the potentials are also suppressed by 1/M which include relativistic corrections, spin-orbital and spin-spin interactions. We only work to order 1/M 0 since the charm quark mass is large. At this order, the Hamiltonians only contain the relative kinetic term and V (0) T,Σ . Inside the deconfined QGP, the potential is flattened and can be approximated by Coulomb interactions Since we keep track of the evolution of both the bound diquarks and unbound charm quarks in the Boltzmann equations, the potentials have no imaginary parts. The interaction between the anti-triplet diquarks and the medium can be decomposed into two parts: a part that only changes the c.m. motion and leaves the bound state intact and the other part that only modifies the relative motion and can destroy the bound state. The decomposition is explicit in the pNRQCD Lagrangian by the multipole expansion. At the order we are working, the c.m. motion part is fully described by the gauged kinetic term of the anti-triplet field, in the same way as the interaction between heavy quarks and the medium. The changes of the c.m. motions of diquarks are treated as diffusion in the Boltzmann equation, in the same way as the heavy quark diffusion (see C c and C cc in expression (1)). The change of the relative motion is described by terms of at least linear order in r. For example, the anti-triplet can interact with the sextet via a color dipole interaction where the chromoelectric field is given by q, ǫ * λ , a and t a F is the generator of the fundamental representation.
At leading order in r, the transition between unbound charm quark pairs and bound diquarks can only occur between an unbound sextet and a bound anti-triplet. The Feynman diagram of the transition via gluon absorption or emission is shown in Fig. 1. For simplicity, we only consider the interaction with on-shell gluons in the QGP. Transitions caused by virtual gluons (inelastic scattering with medium constitutes) are at next order in α s and neglected here. The scattering amplitude in Coulomb gauge is given by where k 1,2 are the c.m. momenta, p rel is the relative momentum between the unbound quark pair and q = |q| is the gluon energy. In the matrix element, |ψ 1S is the hydrogen-like 1S wave function for the bound diquark in the anti-triplet, and |Ψ p rel is the Coulomb wave function for the unbound sextet. The 1S binding energy is given by E 1S = −α 2 s M/9. According to the power counting explained above, the c.m. kinetic energies will be neglected. Throughout this paper we set α s = g 2 /(4π) = 0.4.
To calculate rates, we need to average and sum over certain quantum numbers. For convenience, we define where p rel and k 2 are the relative and c.m. momenta of the unbound charm quark pair with momenta p 1 and p 2 . The pre-factor 1 2 avoids double counting in the phase space of two charm quarks. The g-factors are given by where J = 1 is the diquark spin, s = 1 2 is the heavy quark spin, N c = 3 is the number of colors, d 6 = 6 is the sextet multiplicity and d3 = 3 is the anti-triplet multiplicity. For the formation process, one needs to average over the initial sextet multiplicity and only a fraction d 6 /N 2 c of unbound charm quark pairs are in the sextet, which can form a diquark by radiating out a gluon at the order of r and (1/M ) 0 . The formed 1S diquark is a color anti-triplet and thus has to be in the spin triplet because of the an-tisymmetric nature of fermions. So another spin factor 2J+1 (2s+1) 2 = 3 4 is inserted. For the dissociation process, one needs to average over the initial anti-triplet multiplicity. The phase space measure is relativistic for gluons and non-relativistic for charm quarks and diquarks, which is consistent with our field definitions. Formation from unbound anti-triplet pairs only happens at higher orders in r and 1/M . The gain and loss collision terms in the Boltzmann transport equations can be written as where the "δ−derivative" symbol is defined as where the δ in the second line denotes the standard functional variation and h(p 1 , p 2 , · · · , p n ) and a(p i ) are arbitrary independent functions. In C ± c two such "δ−derivatives" are involved because the initial or final states contain two charm quarks.
The rate of charm quarks combining Γ f and the dissociation rate of a diquark Γ d can be defined as The scattering amplitude and the rate are calculated in the rest frame of the diquark for dissociation and that of the unbound quark pair for formation, where the pN-RQCD is valid. The Bose distribution of medium gluons n (q) B is boosted into the rest frames, respectively. The two frames are not equivalent but since the gluon energy is small compared to M (T M ), the difference is suppressed by T /M . We test the implementation of the formation and dissociation rates in a static QGP box with a constant temperature. After evolving for a sufficiently long period, the system of charm quarks and diquarks reaches thermal equilibrium. The equilibrium test is similar to that for heavy quarks and quarkonia [7].
To solve the transport equations, an initial condition is needed. Due to the large mass, the charm quark can be thought of being produced from the initial hard scattering in heavy ion collisions, before the QGP is formed. The initial transverse momentum and rapidity distribution from the hard scattering is calculated from FONLL [15] with the nuclear parton distribution function (PDF) EPS09 [16]. The nuclear PDF contains a modification of the proton PDF due to nuclear many-body effects. The FONLL calculation is done with the renormalization and factorization scale m T = M 2 + p 2 T . The number of charm quarks produced in one collision event is determined by σT AA , the product of the cross section σ per binary collision calculated in FONLL, and the nuclear thickness function T AA derived from binary collision models. Here we will focus on collisions with 0−10% centrality, which corresponds to impact parameters from 0 to 5 fm roughly and T AA ≈ 23 mb [17].
The initial position of the charm quark produced is sampled using the Trento model [18], a binary collision model. The model assumes the heavy ion collision is a superposition of a number of nucleon-nucleon collisions and calculates the spatial probability distribution where two nucleons from the approaching nuclei scatter. The charm quark production is a short-distance process, implying that its initial position is roughly the same as the location where the two parent nucleons scatter.
Each binary collision also deposits a certain amount of energy and entropy into the system. The Trento model also gives the initial energy and entropy densities. These are then fed into a 2 + 1 dimensional viscous hydrodynamical simulator VISHNew [19,20], which numerically solves the hydrodynamical equations with the energy-momentum tensor for given initial conditions. Here e and p are the local energy density and pressure, and u µ is the local fourvelocity of the QGP. Π is the bulk stress with the bulk viscosity ζ, and π µν is the shear stress tensor with the shear viscosity η. Here the angle bracket means traceless symmetrization.
With the initial condition and hydrodynamical background given, we solve the transport equations by test particles Monte Carlo simulations. The hydrodynamical simulation is assumed to start at the co-moving time τ = 0.6 fm/c. Before this, we assume the charm quarks are just free-streaming without interactions. After τ = 0.6 fm/c, we consider three types of processes at each time step ∆t = 0.04 fm/c in the laboratory frame: diffusion, formation and dissociation.
First, for each charm quark and diquark, we determine their thermal scattering rate with medium constituents. The product of the rate and time step ∆t gives the scattering probability. Then we use random numbers to determine whether a certain process occurs. If so, we sample the momenta of the incoming medium constituent from a thermal distribution and obtain the momenta of outgoing particles by energy-momentum conservation. Finally, we update both particles' momenta and positions after one time step.
Second, for each diquark, we calculate its dissociation rate and probability within a time step as above. If the diquark is determined to dissociate, we replace it by two unbound charm quarks whose momenta are determined from energy-momentum conservation and whose positions are given by that of the diquark just before the dissociation.
Finally, for each charm quark with position y i and momentump i , whose neighboring charm quarks have positions y j and momentap j , we need to determine the diquark formation rate by using expressions (14), (18), (21). A problem appears, because the two quark distributions should be evaluated at the same position, but the product of two delta functions is ill-defined. We introduce a position dependence of the combination probability by means of a Gaussian function with a width chosen as the diquark Bohr radius a B = α s M/3. This ensures that the combination rate for a widely separated charm quark pair vanishes. The product of the local distributions in (14) is thus replaced with where the sum runs over all unbound charm quark pairs. For each charm quark i, the diquark formation rate in expression (21) involves a sum over j. If a diquark is formed, we replace the unbound charm quark pair by a diquark whose momentum is determined by momentum conservation and whose position is given by the centerof-mass position of the quark pair as indicated in (27). When particles reach the hadronization hypersurface determined by the local transition temperature T c ≈ 154 MeV, each diquark combines with a thermal up or down quark to form a doubly charmed baryon. Here we use a simple hadronization model: a massless up or down quark is sampled from a Fermi-Dirac distribution with the temperature T c , and its momentum is added to the diquark momentum to determine the baryon momentum. The baryon energy is fixed by the momentum and vacuum mass m(Ξ cc ). We assume all diquarks end up as the ground Ξ cc states because excited states decay to the ground state much faster than the weak decay of the ground state [21,22]. In this way, roughly half the diquarks end up as Ξ ++ cc . A more realistic hadronization model would include the effect of the baryon wave function on the coalescence probability.
We have simulated 40,000 nuclear collision events. In each event, the initial charm quark momentum is sampled over the range p T ∈ [0, 30] GeV and y ∈ [− 8,8]. At the end of each calculation, we accept Ξ ++ cc in the kinematic range p T ∈ [0, 5] GeV and y ∈ [−1, 1]. The p T spectra integrated over this rapidity range are shown in Fig. 2. The yield within this kinematic range is N (Ξ ++ cc ) ≈ 0.02 per collision.
So far, we have assumed that the diquark can be formed at any temperature. This cannot be true due to the Debye screening of the attractive color force inside the QGP. To understand the influence of Debye screening on Ξ ++ cc production, we repeat the calculation but assume a melting temperature T m = 250 MeV above which the charm diquark cannot be formed inside the QGP. The yield in the same kinematic range is then reduced to N (Ξ ++ cc ) ≈ 0.0125 per collision. The melting temperature of heavy diquarks can be studied from their free energies, in a similar way as quarkonia melting temperatures [23]. The free energy of a heavy quark pair could be studied on a lattice by calculating the correlations of two Polyakov loops at dif- ferent lattice locations, where each Polyakov loop corresponds to a static thermal heavy quark [24]. The free energy projected onto the color anti-triplet state can be used to study the binding energies and spectral functions of diquarks, from which one can obtain the melting temperature. The projections onto the anti-triplet and sextet states were first studied in Ref. [25]. In the appendix, we explain how to project onto the anti-triplet in a gauge invariant but path dependent way. We also show that under a weak coupling expansion, the free energy of a pair of heavy quarks in the anti-triplet is the sum of the free energies of two individual heavy quarks and their attractive potential energy. A previous gauge dependent lattice study can be found in Ref. [26].
The calculation presented here can be improved in several ways. First, one can include higher-order corrections to the in-medium processes. The in-medium potentials of the diquark can also be made temperature-dependent by performing matching calculations between lattice results of Wilson loops and pNRQCD. Furthermore, one can use more realistic hadronization models. Finally, effects of the initial charm quark momentum distribution modifications from the pre-equilibrium effects could be studied.
The calculation can be extended to the production of other doubly heavy baryons, such as Ξ bb and Ξ bc , and doubly heavy tetraquarks, among which the bbūd ground state with J P = 1 + is predicted to be stable [27][28][29][30][31]. The stability of heavy tetraquarks has been investigated previously in Ref. [32]. For Ξ bb , the only difference is that fewer bottom quarks are produced than charm quarks. This implies that the probability of having two bottom quarks come close and form a bottom diquark is much smaller. Thus, one expects a correspondingly smaller yield of Ξ bb . For Ξ bc , there exist extra dipole terms in the pNRQCD Lagrangian for transitions among anti-triplets (or sextets) [12], which means that an unbound pair of bottom and charm quarks in the anti-triplet channel can form a bound bc diquark via a dipole transition. For tetraquarks, the in-medium evolution of heavy quarks and diquarks proceeds in the same way, but the anti-triplet diquark hadronizes by coalescing with two light antiquarks. This process is analogous to the formation of an antibaryon containing a single heavy antiquark, while the formation of a doubly heavy baryon is analogous to the creation of a heavy meson. Heavy baryon (Λ c ) emission is known to be enhanced relative to heavy meson (D 0 ) emission in relativistic heavy ion collisions [33] as a consequence of quark recombination from the thermal quark-gluon plasma [34], compared with protonproton collisions. A similar enhancement of the production of doubly heavy tetraquarks, relative to the production of doubly heavy (anti-)baryons, can be expected. The measured ratio Λ c /D 0 ≈ 1 in Au+Au collisions at RHIC suggests that the yield of doubly heavy baryons and tetraquarks should also be approximately equal.
Finally, we discuss the feasibility of experimental measurements. The crucial factor is the yield-to-background ratio. Based on our calculations, the number of Ξ ++ cc produced at the LHC energies may be large enough. But at the same time, higher collision energies mean higher levels of background. Though a measurement is currently difficult, it is promising that the noisy background difficulty will be overcome in the future with detector upgrades such as the ALICE Inner Tracking System up-grade. With the high-resolution detectors, one can apply stricter topological cuts to reduce the level of background and increase the yield-to-background ratio. Just as the STAR Collaboration first measured the Λ c production in heavy ion collisions with the newly installed Heavy Flavor Tracker [33], measurements of doubly heavy baryons and even bound tetraquarks in heavy ion collisions may become possible in the future. Experimental measurements rely on the reconstruction from decay products of Ξ ++ cc . The decay properties of doubly heavy baryons have been intensely studied [3,[35][36][37][38][39][40][41][42].
In conclusion, we have used Boltzmann transport equations to describe the in-medium formation, dissociation, and diffusion of charm diquarks. Based on it, we estimate the production rate and p T spectra of the doubly charmed baryon Ξ ++ cc in central Pb-Pb collision at 2.76 TeV. It will be of great interest if experimental efforts are taken to try to measure Ξ ++ cc in heavy ion collisions. A measurement of the production rate would allow us to extract the melting temperature of the charm diquark in QGP from the above calculation. Comparison can be made with the melting temperature calculated from lattice results of the free energy of the anti-triplet. These experimental and lattice studies would provide valuable information to our understanding of QCD at finite temperature and properties of QGP. The free energy of a heavy quark pair in the anti-triplet can be defined as |s s|e −βH W k j ((0, β), (r, β))ψ j (r, β)ψ i (0, β)ψ † i (0, 0)ψ † j (r, 0)W † jk ((0, 0), (r, 0))|s .
In the static heavy quark limit [24], where T is the time ordering operator. The starting and ending points of the Wilson line along the Euclidean time direction are the same due to the periodicity of gauge fields at finite temperature and is denoted as the Polyakov line L(r). Then using the anti-commutation relation of heavy quark operators it can be shown where Ô T ≡ |s s|e −βHÔ |s and TrL is the Polyakov loop. Both the correlation terms in the above expression are gauge invariant because of the cyclic property of the trace and the periodicity of gauge fields. Schematic diagrams for the two correlation terms are shown in Fig. 3.
In a similar way, the sextet free energy can be defined as The free energy of an anti-triplet heavy quark pair is the sum of the free energies of two individual heavy quarks and their color attractive potential energy.
In a similar way, F QQ(6) (r) = 2F Q + 1 3 (A.14) The free energy of a sextet is the sum of the free energies of two individual heavy quarks and their color repulsive potential energy. Though up to order g 2 the anti-triplet and sextet free energies are independent of the Wilson line paths in the definition, they are generally dependent on the paths beyond the leading order [43].