We give an estimate of Ξ þþ 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 nonrelativistic 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 Ξ þþ 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 Ξ þþ 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 antitriplet 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. DOI: 10.1103/PhysRevD.97.074003 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 antiquark 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 antitriplet 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 quark-gluon 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 antitriplet 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...