Helium nuclei, deuteron and dineutron in 2+1 flavor lattice QCD
We calculate the binding energies for multi-nucleon bound states with the nuclear mass number less than or equal to 4 in 2+1 flavor QCD at the lattice spacing of fm employing a relatively heavy quark mass corresponding to GeV. To distinguish a bound state from attractive scattering states, we investigate the volume dependence of the energy shift between the ground state and the state of free nucleons by changing the spatial extent of the lattice from 2.9 fm to 5.8 fm. We conclude that He, He, deuteron and dineutron are bound at GeV. We compare their binding energies with those in our quenched studies and also with several previous investigations.
pacs:11.15.Ha, 12.38.Aw, 12.38.-t 12.38.Gc
Lattice QCD has a potential ability to quantitative understand the nature of nuclei, whose characteristic feature is a hierarchical structure in the strong interaction. The nuclear binding energy is experimentally known to be about 10 MeV per nucleon, which is much smaller than the typical energy scale of hadrons. A measurement of the binding energies is therefore the first step for direct investigation of nuclei in lattice QCD. A key ingredient in the study is a systematic change of the spatial volume of the lattice to distinguish a bound state from an attractive scattering state.
We carried out a first attempt to measure the binding energies of the He and He nuclei in quenched QCD with a rather heavy quark mass corresponding to GeV, thereby avoiding heavy computational cost Yamazaki et al. (2010). We followed this work by a renewed investigation of the bound state for the two-nucleon channel in quenched QCD at the same quark mass, which found that not only the deuteron in the S channel but also the dineutron in the S channel is bound Yamazaki et al. (2011). Independently, NPLQCD collaboration reported the possibility that a bound state is formed in both channels at GeV in 2+1 flavor QCD Beane et al. (2012a). They later confirmed the bound states for the helium nuclei and the two-nucleon channels at GeV in 3-flavor QCD taking a different choice for the quark and gluon actions Beane et al. (2012b).
In this article we report on our investigation of the dynamical quark effects on the binding energies of the helium nuclei, the deuteron and the dineutron. We perform 2+1 flavor lattice QCD simulation with the degenerate up and down quark mass corresponding to GeV. Four lattice sizes are employed to take the infinite spatial volume limit: , , and , whose spatial extent ranges from 2.9 fm to 5.8 fm with the lattice spacing of fm Aoki et al. (2010a).
For the helium nuclei our main interest lies in the magnitude of the binding energies, since all studies carried out so far, both in quenched and in unquenched QCD and for several quark mass values, agree on the bound state nature for helium nuclei. Much more intriguing is the two-nucleon system, for which there are two ways to study. One is a direct investigation Fukugita et al. (1994, 1995); Beane et al. (2006, 2010); Yamazaki et al. (2011); Beane et al. (2012a, b) in which one calculates the two-nucleon Green’s functions directly in lattice QCD, and the other is an indirect calculation by means of the two-nucleon effective potential extracted from the two-nucleon wave function in lattice QCD Ishii et al. (2007); Aoki et al. (2010b).
So far only the former method has reported the binding energies of the two-nucleon systems. In quenched QCD the bound state nature has been confirmed for both channels at GeV in our recent work Yamazaki et al. (2011). On the other hand, unquenched studies show a complicated situation. A somewhat early study in 2+1 flavor QCD with a mixed action Beane et al. (2006) reported a positive energy shift (repulsive interaction) in both channels at GeV. More recently, however, deep bound states were observed at GeV in 3-flavor QCD Beane et al. (2012b). We hope to shed light on this situation with our own investigation in 2+1 flavor QCD.
This paper is organized as follows. In Sec. II we explain simulation details including the simulation parameters and the interpolating operators for the multi-nucleon channels. Section III presents the results of the binding energies for the helium nuclei, the deuteron and the dineutron. We compare our results with those in the previous studies. Conclusions and discussions are summarized in Sec. IV.
Ii Simulation details
ii.1 Simulation parameters
We generate 2+1 flavor gauge configurations with the Iwasaki gauge action Iwasaki (2011) and the non-perturbative -improved Wilson quark action at with Aoki et al. (2006). The lattice spacing is fm, corresponding to GeV, determined with GeV Aoki et al. (2010a). We take four lattice sizes, , , and , to investigate the spatial volume dependence of the ground state energy shift between the multi-nucleon system and the free nucleons. The physical spatial extents are 2.9, 3.6, 4.3 and 5.8 fm, respectively. Since it becomes harder to obtain a good signal-to-noise ratio at lighter quark masses for multi-nucleon systems Lepage (1989); Fukugita et al. (1995), we employ the hopping parameters which correspond to GeV and GeV and the physical value for the strange quark mass. These values are chosen based on the previous results for and obtained by the PACS-CS Collaboration Aoki et al. (2009, 2010a).
We employ the domain-decomposed Hybrid-Monte-Carlo (DDHMC) algorithm Lüscher (2003, 2005) for the degenerate light quarks and the UV-filtered PHMC (UVPHMC) algorithm Ishikawa et al. (2006) for the strange quark employing the Omelyan-Mryglod-Folk integrator Omelyan et al. (2003); Takaishi and de Forcrand (2006). The algorithmic details are given in Ref. Aoki et al. (2009). We summarize simulation parameters in Table 1 including the block sizes in DDHMC and the polynomial order in UVPHMC. We take for the trajectory length of the molecular dynamics in all the runs. The step sizes are chosen such that we obtain reasonable acceptance rates presented in Table 1. We generate the gauge configurations in a single run except for the case for which we carry out two runs. The total trajectory length is about 2000 for all the volumes, except 4000 for the case of smallest volume.
ii.2 Calculation method
We extract the ground state energies of the multi-nucleon systems and the nucleon state from the correlation functions
with being appropriate operators for He, He, two-nucleon S and S channels, and the nucleon state (see the next subsection for actual expressions).
We carry out successive measurements in the interval of 10 trajectories. The errors are estimated by jackknife analysis choosing 200 trajectories for the bin size for all volumes, except for the largest volume for which we use 190. The numbers of configurations are listed in Table 2. We attempt to extract as much information as possible from each configuration by repeating the measurement of the correlation functions for a number of sources at different spatial points and time slices. For the and lattices, we calculate the correlation functions not only in the temporal direction but also in the three spatial directions exploiting the space-time rotational symmetry. We found that this procedure effectively increases statistics by a factor of four. This factor is included in the number of measurements on each configuration given in Table 2
We are interested in the energy shift between the ground state of the multi-nucleon system and the free nucleons on an box,
with being the lowest energy level for the multi-nucleon channel, the number of nucleon and the nucleon mass. This quantity is directly extracted from the ratio of the multi-nucleon correlation function divided by the -th power of the nucleon correlation function
where the same source operator is chosen for the numerator and the denominator. We also define the effective energy shift as
which is useful to check the plateau region in later sections.
ii.3 Interpolating operators
We use an interpolating operator for the proton given by
where and and are the Dirac index and the color indices, respectively. The neutron operator is obtained by replacing by in the proton operator. To save the computational cost we use the non-relativistic quark operator, in which the Dirac index is restricted to the upper two components.
The He nucleus has zero total angular momentum, positive parity and zero isospin . We employ the simplest He interpolating operator with zero orbital angular momentum , and hence with being the total spin. Such an operator was already given long time ago in Ref. Beam (1967),
with being up/down spin of each nucleon, and are obtained by replacing in by for the isospin. Each nucleon in the sink operator is projected to zero spatial momentum.
We also calculate the correlation function of the He nucleus whose quantum numbers are , and . We employ the interpolating operator in Ref. Bolsterli and Jezak (1964),
with the zero momentum projection on each nucleon in the sink operator.
The two-nucleon operators for the S and S channels are given by
In the spin triplet channel the operators for the other two spin components are constructed in a similar way. We take average over the three spin components.
The quark propagators are solved with the periodic boundary condition in all the spatial and temporal directions using the exponentially smeared source of form
after the Coulomb gauge fixing. We choose the smearing parameters depending on the volume (see Table 2) in order to obtain reasonable plateaux of the effective energy for the ground states in the multi-nucleon channels as well as for the nucleon. For the source operators explained above we insert the smeared quark fields of Eq. (12) for each nucleon operator located at the same spatial point . Each nucleon in the sink operator, on the other hand, is composed of the point quark fields, and projected to zero spatial momentum.
ii.4 Difficulties for multi-nucleon channel
There are several computational difficulties in the calculation of the correlation functions for the He and He channels. One is a factorially large number of Wick contractions for the quark-antiquark fields. A naive counting gives for a nucleus composed of protons and neutrons, which quickly becomes prohibitively large beyond three-nucleon systems, e.g., 2,880 for He and 518,400 for He. To overcome the difficulty, we use the reduction techniques proposed in our exploratory work Yamazaki et al. (2010). After the reduction, only 1107 (93) contractions are required for the correlation function in the He (He) channel. Other reduction techniques for the large number of the Wick contractions have been proposed for the multi-meson Detmold and Savage (2010) and multi-baryon Doi and Endres (2012); Detmold and Orginos (2012) channels.
Another difficulty in studying a multi-nucleon bound state is the identification of the bound state nature in a finite volume, because an attractive scattering state yields a similar energy shift due to the finite volume effect Lüscher (1986); Beane et al. (2004); Sasaki and Yamazaki (2006). To solve the problem we need to investigate the volume dependence of the measured energy shift Yamazaki et al. (2010, 2011): For a scattering state, the energy shift decreases in proportion to at the leading order in the expansion Lüscher (1986); Beane et al. (2007), while for a bound state the energy shift remains at a finite value in the infinite spatial volume limit. In order to distinguish a non-zero constant from a behavior in the energy shift, we employ four spatial extents from 2.9 to 5.8 fm.
We first show the effective nucleon mass on the (5.8 fm) box in Fig. 1 as a typical result. The plateau of the effective mass is clearly observed. A fit result of the correlation function with an exponential form is also drawn in the figure with the one standard deviation error band. We list the nucleon mass together with the pion mass in Table 2.
iii.2 He nucleus
The effective energy shift defined in Eq. (4) is plotted in Fig. 2. The signal is clear up to , beyond which the statistical error increases rapidly. The energy shift is extracted from of Eq. (3) by an exponential fit over the range of –14. The fit result is denoted by the solid lines with the statistical error band in Fig. 2. The systematic error in the fit is estimated from the variation of the fit results with the minimum or maximum time slice changed by . Results with similar quality are obtained on other volumes. We summarize the values of with the statistical and systematic errors in Table 3.
Figure 3 shows the volume dependence of as a function of . The inner bar denotes the statistical error and the outer bar represents the statistical and systematic errors combined in quadrature. The negative energy shifts are obtained in all the four volumes. We extrapolate the results to the infinite volume limit with a simple linear function of ,
The systematic error is estimated from the variation of the results obtained by alternative fits which contains a constant fit of the data and a fit of the data obtained with a different fit range in . The non-zero negative value obtained for the infinite volume limit shown in Fig. 3 and Table 3 leads us to conclude that the ground state is bound in this channel for the quark masses employed. The binding energy MeV, where the first error is statistical and the second one is systematic, is consistent with the experimental result of 28.3 MeV and also with the previous quenched result at GeV Yamazaki et al. (2011). Note that the error is still quite large.
A recent work in 3-flavor QCD at GeV reported a value 110(20)(15) MeV for the binding energy of He nucleus Beane et al. (2012b). This is about three times deeper than our value. Whether this difference can be attributed to the quark mass dependence in unquenched calculations needs to be clarified in future.
iii.3 He nucleus
Figure 4 shows the effective energy shift of Eq. (4). The quality of the signal is better than the He channel in Fig. 2. An exponential fit of in Eq. (3) with the range of –14 yields a negative value, which is denoted by the solid lines with the statistical error band in Fig. 4. The systematic error in the fit is estimated in the same way as in the He case.
As listed in Table 3 we find non-zero negative values for the energy shift for all the volumes. The volume dependence is illustrated in Fig. 5 as a function of with the inner and outer error bars as explained in the previous subsection. We carry out a linear extrapolation of Eq. (13). The systematic error is estimated in the same way as in the He channel. The energy shift extrapolated to the infinite spatial volume limit is non-zero and negative, see Fig. 5 and Table 3, which means that the ground state is a bound state in this channel. The value of MeV is roughly three times larger than the experimental result, 7.72 MeV, though consistent with our previous quenched result at GeV Yamazaki et al. (2011).
In 3-flavor QCD MeV was reported Beane et al. (2012b) at a heavier quark mass corresponding to GeV. Here again future work is needed to see if a quark mass dependence explains the difference from the experiment.
iii.4 Two-nucleon channels
iii.4.1 Present work
In Fig. 6 we show the time dependence for of Eq. (4) in the S channel. The signals are lost beyond . We observe negative values beyond the error bars in the plateau region of –14. We extract the value of from an exponential fit for of Eq. (3) in the range of –14. The systematic error of the fit is estimated as explained in the previous subsections.
Figure 7 shows the result for in the S channel on the (5.8 fm) box. The value of is again negative beyond the error bars in the plateau region, though the absolute value is smaller than in the S case. The energy shift is obtained in the same way as for the S channel.
The volume dependences of in the S and S channels are plotted as a function of in Figs. 8 and 9, respectively. The numerical values of on all the spatial volumes are summarized in Table 4, where the statistical and systematic errors are given in the first and second parentheses, respectively. There is little volume dependence for , indicating a non-zero negative value in the infinite volume and a bound state, rather than the dependence expected for a scattering state, for the ground state for both channels.
The binding energies in the infinite spatial volume limit in Table 4 are obtained by fitting the data with a function including a finite volume effect on the two-particle bound state Beane et al. (2004); Sasaki and Yamazaki (2006),
where and are free parameters, is three-dimensional integer vector, and denotes the summation without . The binding energy is determined from
where we assume
The systematic error is estimated from the variation of the fit results choosing different fit ranges in the determination of and also using constant and linear fits as an alternative fit forms. We obtain the binding energies =11.5(1.1)(0.6) MeV and 7.4(1.3)(0.6) MeV for the S and S channels, respectively. The result for the S channel is roughly five times larger than the experimental value, 2.22 MeV. Our finding of a bound state in the S channel contradicts the experimental observation. These features are consistent with our quenched results with a heavy quark mass corresponding to GeV Yamazaki et al. (2011).
iii.4.2 Comparison with previous studies
A number of studies have been performed for the two-nucleon channel after the first work of Ref. Fukugita et al. (1995). It is therefore instructive to summarize the results and make a comparison with each other. Table 5 tabulates, in chronological order, the results for for the S and S channels together with the pion mass and the spatial extent in physical units. The numbers are plotted in Figs. 10 and 11 for the S and S channels, respectively, as a function of .
The early studies in Refs. Fukugita et al. (1995); Beane et al. (2006); Aoki et al. (2010b) employed a single volume, and we do not observe a common feature or trend among them. The positive values for in Ref. Beane et al. (2006) means repulsive interaction for both channels, which is not seen in other studies. The results for in Ref. Aoki et al. (2010b) is an order of magnitude smaller compared to other groups, probably due to significant contamination from excited states.
The four recent studies Yamazaki et al. (2011); Beane et al. (2012a, b) have made a systematic investigation of the spatial volume dependence. Our quenched and 2+1 flavor results show qualitatively the same feature that the binding energy for the S channel is much larger than the experimental value and the bound state is observed in the S channel. The 2+1 flavor results from Ref. Beane et al. (2012a, b) at GeV give non-zero negative values for in both channels on the (3.9 fm) box, which are consistent with our results as shown in Table 5. Unfortunately, the extrapolation to the infinite spatial volume limit introduces large errors so that becomes consistent with zero within the error bars. The most recent study Beane et al. (2012b) worked at a heavier quark mass of GeV in 3-flavor QCD, and found large values for the binding energies: 25(3)(2) MeV for the S channel and 19(3)(1) MeV for the S channel Beane et al. (2012b). While all recent studies are consistent with a bound ground state for both S and S channels when quark masses are heavy, quantitative details still need to be clarified.
Iv Conclusion and discussion
We have calculated the binding energies for the helium nuclei, the deuteron and the dineutron in 2+1 flavor QCD with GeV and GeV. The bound states are distinguished from the attractive scattering states by investigating the spatial volume dependence of the energy shift . In the infinite spatial volume limit we obtain
While the binding energy for the He nucleus is comparable with the experimental value, those for the He nucleus and the deuteron are much larger than the experimental ones. Furthermore we detect the bound state in the S channel as in the previous study with quenched QCD, which is not observed in nature. These findings and the enhanced binding energies at GeV in 3-flavor QCD Beane et al. (2012b) tell us that a next step of primary importance is to reduce the up-down quark mass toward the physical values. A possible scenario in the two-nucleon channels is as follows. The binding energy in both channels diminishes monotonically as the up-down quark mass decreases. At some point of the up-down quark mass the binding energy in the S channel vanishes and the bound state evaporates into the attractive scattering state, while the binding energy in the S channel remains finite up to the physical point. This is a dynamical question on the strong interaction, and only lattice QCD could answer it.
Numerical calculations for the present work have been carried out on the HA8000 cluster system at Information Technology Center of the University of Tokyo, on the PACS-CS computer under the “Interdisciplinary Computational Science Program” of Center for Computational Sciences in University of Tsukuba, on the T2K-Tsukuba cluster system and HA-PACS system at University of Tsukuba, and on K computer at RIKEN Advanced Institute for Computational Science. We thank our colleagues in the PACS-CS Collaboration for helpful discussions and providing us the code used in this work. This work is supported in part by Grants-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology (Nos. 18104005, 18540250, 22244018) and Grants-in-Aid of the Japanese Ministry for Scientific Research on Innovative Areas (Nos. 20105002, 21105501, 23105708).
- Yamazaki et al. (2010) T. Yamazaki, Y. Kuramashi, and A. Ukawa (PACS-CS Collaboration), Phys.Rev. D81, 111504 (2010).
- Yamazaki et al. (2011) T. Yamazaki, Y. Kuramashi, and A. Ukawa, Phys. Rev. D84, 054506 (2011).
- Beane et al. (2012a) S. Beane et al. (NPLQCD Collaboration), Phys. Rev. D85, 054511 (2012a).
- Beane et al. (2012b) S. Beane et al. (NPLQCD Collaboration) (2012b), eprint arXiv:1206.5219 [hep-lat].
- Aoki et al. (2010a) S. Aoki et al. (PACS-CS Collaboration), Phys. Rev. D81, 074503 (2010a).
- Fukugita et al. (1994) M. Fukugita, Y. Kuramashi, H. Mino, M. Okawa, and A. Ukawa, Phys. Rev. Lett. 73, 2176 (1994).
- Fukugita et al. (1995) M. Fukugita, Y. Kuramashi, M. Okawa, H. Mino, and A. Ukawa, Phys. Rev. D52, 3003 (1995).
- Beane et al. (2006) S. R. Beane, P. F. Bedaque, K. Orginos, and M. J. Savage, Phys. Rev. Lett. 97, 012001 (2006).
- Beane et al. (2010) S. R. Beane et al. (NPLQCD), Phys. Rev. D81, 054505 (2010).
- Ishii et al. (2007) N. Ishii, S. Aoki, and T. Hatsuda, Phys. Rev. Lett. 99, 022001 (2007).
- Aoki et al. (2010b) S. Aoki, T. Hatsuda, and N. Ishii, Prog.Theor.Phys. 123, 89 (2010b).
- Iwasaki (2011) Y. Iwasaki (2011), eprint arXiv:1111.7054 [hep-lat].
- Aoki et al. (2006) S. Aoki et al. (CP-PACS/JLQCD Collaborations), Phys. Rev. D73, 034501 (2006).
- Lepage (1989) G. P. Lepage (1989), eprint CLNS-89-971, C89-06-04.
- Aoki et al. (2009) S. Aoki et al. (PACS-CS Collaboration), Phys. Rev. D79, 034503 (2009).
- Lüscher (2003) M. Lüscher, JHEP 0305, 052 (2003).
- Lüscher (2005) M. Lüscher, Comput. Phys. Commun. 165, 199 (2005).
- Ishikawa et al. (2006) K.-I. Ishikawa et al. (PACS-CS Collaboration), PoS LAT2006, 027 (2006).
- Omelyan et al. (2003) I. P. Omelyan, I. M. Mryglod, and R. Folk, Comput. Phys. Commun. 151, 272 (2003).
- Takaishi and de Forcrand (2006) T. Takaishi and P. de Forcrand, Phys. Rev. E73, 036706 (2006).
- Beam (1967) J. E. Beam, Phys. Rev. 158, 907 (1967).
- Bolsterli and Jezak (1964) M. Bolsterli and E. Jezak, Phys. Rev. 135, B510 (1964).
- Detmold and Savage (2010) W. Detmold and M. J. Savage, Phys. Rev. D82, 014511 (2010).
- Doi and Endres (2012) T. Doi and M. G. Endres (2012), eprint arXiv:1205.0585 [hep-lat].
- Detmold and Orginos (2012) W. Detmold and K. Orginos (2012), eprint arXiv:1207.1452 [hep-lat].
- Lüscher (1986) M. Lüscher, Commun. Math. Phys. 105, 153 (1986).
- Beane et al. (2004) S. R. Beane, P. F. Bedaque, A. Parreno, and M. J. Savage, Phys. Lett. B585, 106 (2004).
- Sasaki and Yamazaki (2006) S. Sasaki and T. Yamazaki, Phys. Rev. D74, 114507 (2006).
- Beane et al. (2007) S. R. Beane, W. Detmold, and M. J. Savage, Phys. Rev. D76, 074507 (2007).
- Bratt et al. (2010) J. Bratt et al. (LHPC Collaboration), Phys. Rev. D82, 094502 (2010).
|# config.||bin size||# meas.||[GeV]||[GeV]|
|[MeV]||fit range||[MeV]||fit range|
|[MeV]||fit range||[MeV]||fit range|
|Ref.||quark action||# flavor||[GeV]||[fm]||[MeV]|
|Fukugita et al. (1995)||Wilson||0||0.72||2.7||29.8(6.9)||14.7(4.3)|
|Beane et al. (2006)||Mixed (DW on Asqtad)||2+1||0.35||2.5||16(19)||16(13)|
|Aoki et al. (2010b)||Wilson||0||0.38||4.4||0.97(37)||0.68(26)|
|Yamazaki et al. (2011)||Wilson-clover||0||0.80||3.1||10.2(2.2)(1.6)||6.1(2.3)(2.2)|
|Beane et al. (2010)||Aniso. Wislon-clover||2+1||0.39||2.4||1.6(2.6)(4.3)||3.9(1.7)(2.6)|
|Beane et al. (2012a)||Aniso. Wilson-clover||2+1||0.39||3.0||22.3(2.3)(5.4)||10.4(2.6)(3.1)|
|Beane et al. (2012b)||Stout Wilson-clover||2+1||0.81||3.4||25(3)(2)||16(3)(1)|