What are the fundamental properties of nuclear matter? How does it behave under extreme conditions? – These are some questions that ab initio calculations of hadrons and nuclear matter using numerical lattice techniques are trying to answer. Focussing on one, we ask: Is there room for a bound dibaryon state in nature?
The dibaryon (a 6quark state) was predicted over 30 years ago in a model calculation [1]. So far, however, its detection has eluded experimental physicists, and theoretical studies remain inconclusive. It is therefore interesting to take the step from the well determined modelindependent numerical calculations of lattice Quantum Chromodynamics (QCD) for single baryons (3quark states, such as the proton) to dibaryons – thus paving the way to ever more complicated nuclei. We study the spectrum of (di)baryons in terms of correlation functions that determine the properties of a particle state, such as its mass and whether it is a loosely bound system of individual particles or a tightly bound state (Fig. 1). The interpolating operators involved in correlation functions are related to the particle's wavefunction, which is mostly determined by the individual (valence) quarks it consists of. This is the realm of QCD and our goal is the numerical ab initio lattice QCD calculation of dibaryon properties.
Quantum Chromodynamics (QCD)
At the core of all nuclear matter lies QCD. It is the fieldtheoretical description of the fundamental strong force, which determines the interactions between quarks and gluons, i.e. the constituents of hadrons.
QCD has two distinctive properties: asymptotic freedom at high energies, and confinement at low energies. Asymptotic freedom corresponds to a small coupling constant in the theory, and in this regime perturbation theory gives a highly precise description of the physical phenomena. Nuclear matter, however, is a lowenergy phenomenon and a different approach must be used to explore hadronization, i.e. the formation of bound states and matter from individual quarks. Lattice QCD is a firstprinciples approach that has exactly this capability.
Lattice QCD
In lattice QCD simulations, the properties of baryons and nuclei can be calculated by studying correlation functions of local operators. These expectation values are linked to functional integrals of QCD that are evaluated stochastically using Monte Carlo methods. In addition, every 'measurement' must be performed on large ensembles of field configurations generated by a Markov process and distributed according to the QCD action. The extremely large number of variables in these simulations makes it impossible to use conventional computing resources. At the same time, lattice QCD simulations are well suited for massively parallel computers, as they involve local interactions and require the inversion of large sparse matrices. Lattice QCD is thus well suited for high performance computing (HPC).
To set up a lattice QCD calculation, spacetime has to be discretized on a hypercubic lattice. This induces systematic effects that need to be considered. These are the finiteness of the introduced lattice spacing between the spacetime points at which measurements are made, and the finite extent of the lattice. Effects due to the breaking of symmetry (e.g. rotational) down to a subgroup introduced by the hypercubic lattice must also be considered. Moreover there is the added complication that the simulations become increasingly expensive as the physical limit is approached, specifically when tuning the light quarks to their physical values. To handle these effects, multiple series of measurements are usually made at several unphysical quark masses, several different lattice spacings and lattice volumes. Through closely monitoring the lattice QCD results, a safe and reliable extrapolation to the physical limit may be performed. In the case of the single proton, this has been performed extensively so that a reliable extrapolation to the physical point with controlled systematics can be made (Fig. 2). However, in the much more complicated case of nuclear matter, i.e. nuclei of mass number A ≥ 2, this type of calculation is still in its infancy. As a consequence, the following initial calculation is on a lattice volume of 64 × 323 (time × space3) points, for a pion mass of 450 MeV and a spatial extent of 2 fm. A larger and more physical ensemble of (mπ,=320 MeV, 96 × 483, and 3 fm) will shortly be analyzed on JUQUEEN.
Dibaryons
Hadrons are composite particles, consisting of quarks bound by the strong force. There are 6 quarks in total. However we simulate only the 3 lightest quarks: the up (u), down (d) and strange (s) quarks (in ascending order of mass). There are two families of hadrons that have been observed in nature. The first family is that of quarkantiquark states called mesons, with the pion being its most prominent member. The second is that of 3quark states referred to as baryons, of which the stable proton (quark content uud), the long lived neutron (udd, with a lifetime of ~15 minutes) and Λbaryon (uds) are members. In addition, the dibaryon is a proposed 6quark state (Fig 1). The Hdibaryon is a state with strangeness quantum number –2 and quark content udsuds. This state was predicted in the 1970s, and an indication that such a state exists would be that it is energetically favored over two distinct baryons, i.e. if the Hdibaryon mass is less than twice the single Λbaryon mass.
Numerics
Our preliminary findings indicate that, in order to resolve the mass difference between two individual Λbaryons and a true dibaryon state, a statistical precision of < 1% for the masses of the individual and dibaryon states is required. This means that an extremely large number of individual lattice QCD measurements are required, of order 30,000 (corresponding to approximately 3.8 Mcorehrs on JUQUEEN for our chosen parameters).
To set up a calculation of the Hdibaryon, a QCD interpolating operator that describes it must be used. This is achieved by forming a combination of quark fields so that the resulting interpolating operator has the correct quantum numbers to be related to the Hdibaryon wavefunction. These quark fields then have to be connected ('contracted') to form a measurable correlation function. Numerically, any single lattice QCD measurement of a hadron mass thereby consists of three main parts: the calculation of quark propagators (matrix inversions); techniques used to obtain a good overlap with the state of interest, which requires some tuning of the operator (Fig. 3); and the contraction of the quark propagators to form a hadronic correlation function. For a baryon, the first two are the most expensive. However, due to the larger number of combinations of quark propagator contractions that enter the dibaryon correlation function, stemming from the additional 3 quarks, all three steps are equally significant for a dibaryon measurement, thus making it more computationally demanding.
The number of combinations of quark propagators required for a (di)baryon calculation is equal to the factorials of their respective quark content N = Nu!Nd!Ns!, which is 2 for the proton (uud), 8 for the Hdibaryon (udsuds), and 36 for a deuteron (uududd). For other nuclei of mass number A ≥ 2 this rapidly increases. Consequently, the Hdibaryon system is much more complex and demanding than the single proton. To improve the efficiency of the calculation, we employ two techniques: a 'blocking' algorithm and a 'truncated' solver.
Blocking Algorithm
To tackle the large number of contractions, it is convenient to use the socalled 'blocking' algorithm [3] to handle and simplify the potentially large number of combinations of quark propagators (Fig. 4). This involves a precontraction of the quark propagators at the source of the correlation function, so that the precontracted 'block' can be reused to efficiently combine the quark propagators at the sink. This removes the need for a new source contraction for each combination.
Truncated Solver
Quark propagators are calculated through the inversion of a large sparse matrix using an iterative solving procedure. The truncated solver method [4] utilizes a less precise stopping condition (tolerance of 103 compared to the usual 1010) for the inversion. There is of course an associated bias with respect to a high precision solve. However this can be corrected for, which allows many more measurements to be made in a given time compared to the usual highprecision method and results in at least a 20 % gain in statistical precision.
The truncated solver and the aforementioned blocking algorithm are two very new and modern technical developments that are key to making this calculation possible. With these improvements, the calculation may be achieved using existing HPC resources in a sensible timeframe and has become feasible for lattice ensembles large enough to model the dibaryon without significant distortion due to the finiteness of the lattice spacing and with lattice parameters approaching those of the physical situation.
Results
We show results gathered on a total of ~16,000 measurements [5] obtained in 1.9 Mcorehrs of computing time on JUQUEEN. The mass of the Hdibaryon was determined using the exponential decay of the correlation function. For large time separations the ground state dominates the correlation function and excited states are exponentially suppressed. The ground state mass is most easily computed from taking the logarithmic derivative of the correlation function to determine an 'effective' mass. The ground state mass can then be determined by performing a fit in the plateau region of the effective mass (seen in Fig. 5) where it is independent of the time separation. However, this method is only suited for a determination of the ground state. A more sophisticated method is required for a study of the excited states. For instance, we employ the method of solving a generalized eigenvalue problem (GEVP). The choice of interpolating operator used to measure the Hdibaryon is not unique, and one is free to choose any wavefunction that gives the correct quantum numbers. Depending on how well a given operator describes the actual physical state of the Hdibaryon, the correlation function will have a better or worse overlap with the true Hdibaryon wavefunction. Solving the eigenvalue problem requires a set of such differently overlapping operators and provides the mass eigenvalues of the operator matrix, which allows the corresponding state's mass to be retrieved. Our analysis uses a set of 4 trial wavefunctions to extract the ground and first excited states of the Hdibaryon mass spectrum. The preliminary results in Fig. 5 show that with the present statistical accuracy of ~ 4 %, the Hdibaryon mass is consistent with that of two individual particles. However for a truly conclusive result the statistical error on the dibaryon data must be significantly reduced.
Conclusions
The study of nuclear matter from ab initio lattice QCD calculations is still in its infancy and requires toptier HPC resources for its study. Our first analysis shows a tendency of the Hdibaryon to be unbound. However, the statistical error needs to be substantially reduced. Additionally, the calculation has to be repeated on a series of lattice ensembles with ever more physical, and hence expensive, parameters. Even at this early stage the results do however give insights into the nature of light nuclei, specifically the Hdibaryon [6], that are beyond the reach of other nonab initio approaches. With the continued investment of HPC resources on platforms such as JUQUEEN we will be able to answer the question: is there room for a bound dibaryon state in nature?
The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time for a GCS Large Scale Project on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC). This work was supported by the DFG via SFB 1044 and grant HA 4470/31.
References

[1] Jaffe, R.L.
Phys. Rev. Lett. 38 (1977) 195.

[2] von Hippel, G.M., Jäger, B., Rae, T.D., Wittig, H.
JHEP 1309 (2013) 014.

[3] Doi, T., Endres, M.G.
Comput. Phys. Commun. 184 (2013) 117.

[4] Blum, T., Izubuchi, T., Shintani, E.
Phys. Rev. D88 (2013) 094503.

[5] Francis, A., Miao, C., Rae, T.D., Wittig, H.
PoS (LATTICE 2013) 440, arXiv:1311.3933, 2013.

[6] WalkerLoud, A.
PoS (LATTICE 2013) 013, arXiv:1401.8259, 2013.
contact: Hartmut Wittig, wittig[at]kph.unimainz.de