Optimized Two-Baryon Operators in Lattice QCD

A set of optimized interpolating operators which are dominantly coupled to each eigenstate of two baryons on the lattice is constructed by the HAL QCD method. To test its validity, we consider heavy dibaryons $\Omega_{3Q}\Omega_{3Q}$ ($Q=s,c$) calculated by (2+1)-flavor lattice QCD simulations with nearly physical pion mass. The optimized two-baryon operators are shown to provide effective energies of the ground and excited states separately stable as a function of the Euclidean time. Also they agree to the eigenenergies in a finite lattice box obtained from the leading-order HAL QCD potential $V({\boldsymbol{r}})$ within statistical errors. The overlapping factors between the optimized sink operators and the state created by the wall-type source operator indicate that $V( \boldsymbol{r})$ can be reliably extracted, no matter whether the spacetime correlation of two baryons is dominated by the ground state or the excited state. It is suggested that the optimized set of operators is useful for variational studies of hadron-hadron interactions.


I. INTRODUCTION
Accurate determination of hadron-hadron interactions from QCD is one of the most challenging problems in nuclear and particle physics [1,2]. To achieve the goal, two theoretical methods have been proposed in lattice QCD (LQCD); the Lüscher's finite volume method [3] and the HAL QCD method [4][5][6]. The former focuses on the temporal behavior of the hadronic correlations, from which the scattering phase shifts are extracted through the Lüscher's finite volume formula. On the other hand, the latter considers the spacetime behavior of the hadronic correlations, from which physical observables are extracted through the equal-time Nambu-Bethe-Salpeter (NBS) amplitude, or the NBS wave function in short.
Although these two methods are like looking at the two sides of the same coin and are theoretically related with each other, it has been realized that the "naive" plateau fitting in the Lüscher's finite volume method without variational analysis sometimes leads to misleading results for two-baryon systems. It originates from the fact that separating the ground state from excited states of two baryons below inelastic threshold is exponentially difficult on the lattice as the splitting in energy becomes zero in the large volume limit [7]: Detailed account of this issue is given in Refs. [8,9] (see also, Refs. [10][11][12]).
In contrast, the time-dependent version of the HAL QCD method [6] does not need to identify each level, since the same NBS integral kernel (i.e. the energyindependent non-local potential) governs all the elastic scattering states simultaneously. In other words, any superposition of the NBS wave functions of different energies below inelastic threshold leads to the same nonlocal potential. By using this nice property, one can construct a set of optimized interpolating operators and make a firm connection between the Lüscher's finite volume method and the HAL QCD method for two baryons, as demonstrated in the (2+1)-flavor LQCD simulations with heavy pion mass, m π 510 MeV, and the lattice volumes, (3.6, 4.3, 5.8 fm) 3 [9]. Such optimized operators can also be used to check the validity of the derivative expansion of the non-local potential in the HAL QCD method.
The purpose of this paper is to further explore the idea of optimized operators proposed in Ref. [9] and study whether it is applicable to (2+1)-flavor LQCD near the physical pion mass, m π 146 MeV, with a large lattice volume, (8.1 fm) 3 . We consider two heavy dibaryons Ω 3Q Ω 3Q in the 1 S 0 channel with Q = s and c, where Ω 3Q implies the spin-3/2 baryon composed of three valence quarks with flavor Q. Both systems are in the unitary regime with large scattering lengths as demonstrated recently by the (2 + 1)-flavor LQCD simulations [13,14]. The statistical errors in these systems are relatively small in comparison to the baryons with light valence quarks, so that one can make quantitative analysis on the effect of the optimized operators as well as on the uncertainty due to derivative expansion of the non-local potential. This paper is organized as follows. After a brief review of the HAL QCD method in Sec. II, we introduce arXiv:2201.02782v2 [hep-lat] 6 Apr 2022 a general framework to define the optimized two-baryon operators through the HAL QCD potential in Sec. III. Details of our lattice setup are given in Sec. IV. Numerical results and discussions on eigenfunctions, effective energies and the overlapping factors for Ω 3Q Ω 3Q are presented in Sec. V. Sec. VI is devoted to summary and concluding remarks. The analyses for the higher excited states are given in Appendix. A.

II. THE HAL QCD METHOD
The equal time NBS amplitude for two Omega baryons with energy E is defined in the Euclidean spacetime as whereΩ is a local interpolating operator for Ω ≡ Ω 3Q (Q = s and c) whose explicit form can be found in Refs. [13,14]. Z Ω is the wave function renormalization factor and |2Ω, E is the 2Ω eigenstate with the center of mass energy E = 2 m 2 Ω + k 2 . Using Eq.(1) and the Haag-Nishijima-Zimmermann reduction formula for composite particles [15], the interaction between baryons can be identified as the energyindependent non-local potential U (r, r ) in the equaltime NBS equation [4,5]: In the time-dependent HAL QCD method [6,16], the spacetime correlation function R(r, t) is introduced as a linear superposition of the NBS wave functions below the inelastic threshold; = n a n Ψ En (r)e −(∆En)t + O(e −(∆E * )t ). ( Here E n is the eigenenergy of the n-th elastic scattering state in a finite box, and ∆E n = E n − 2m Ω with ∆E n ≤ ∆E * : The inelastic threshold is denoted by ∆E * ∼ Λ QCD ∼ 300 MeV. The overlapping factors a n = 2Ω, E n |J (0)|0 depend on the choice of the source operator J (0) at t = 0. If t (∆E * ) −1 ∼ 0.7 fm, effects from the inelastic states are exponentially suppressed as O(e −(∆E * )t ), so that one can rewrite Eq.(2) into the timedependent HAL QCD equation, with H 0 = −∇ 2 /m Ω . The advantage of this equation is that we do not need to separate each eigenenergy E n to extract the potential U . Since the potential is spatially localized in QCD by construction, its finite volume correction is exponentially suppressed for large lattice volume. It has been also demonstrated in Refs. [9,17] that observables (such as the phase shifts) do not depend on the choice of source operators J (0) (either wall-type source or smeared source) by taking the ΞΞ system as an example with the lattice volumes (3.6, 4.3, 5.8 fm) 3 and m π 510 MeV, as long as the non-locality of the potential is well approximated by the derivative expansion.
We note that it is practically useful to make a derivative expansion of the non-local potential as Then the leading-order (LO) potential is obtained as

III. OPTIMIZED SINK OPERATORS
Once the HAL QCD potential is obtained from the LQCD simulation of R(r, t) in Eq.(4), one can calculate eigenfunctions and eigenenergies in a finite lattice box by solving Eq.(2). This enables us to construct a set of optimized interpolating operators which couple strongly to each eigenstate, as originally proposed in Ref. [9]. We now apply this idea to find optimized sink operators in the present systems.
Let us first rewrite Eq.(2) in a three-dimensional lattice box by introducing the Hamiltonian H with the discretized Laplacian H 0 and the non-local potential U (r, r ); where ε n is related to ∆E n as From the n-th NBS wave function Ψ En (r) for Hermitian matrix U (r, r ), 1 one may construct an optimized sink operator as a projection to each n-th state, This is expected to couple primarily to |2Ω, E n for large t. Two-baryon temporal correlation after the projection 1 For non-Hermitian U (r, r ), Ψ † En (r) in Eqs. (9,10) should be replaced by m N −1 nm Ψ † Em (r) with Nnm = r Ψ † En (r)Ψ Em (r) being the norm kernel [5,16].
which can be used to extract effective two-baryon energy for each n as This quantity should have plateau structure as a function of t and approach to ∆E n in Eq.(8) for large t. Note that the projection to each eigenstate is possible only if we have spatial information of the correlation function R(r, t). For similar attempts based on the spatial information of the correlation function in different physical contexts, see Refs. [18][19][20].
In the case where only the LO potential V (r) in Eq. (6) is available, eigenfunction ψ n (r) from the LO Hamiltonian H LO ≡ H 0 +V is an approximation of the exact n-th NBS wave function Ψ En (r). Hence resultant quantities in Eqs. (8)-(11) are approximate ones. Then, whether ∆E eff n (t) has a plateau and approaches to Eq.(8) at large t provide a confidence test on the truncation of the derivative expansion of U .
For later purpose, let us introduce the unprojected temporal correlation, R(t) = r R(r, t), which can be decomposed as where b n = a n r Ψ En (r). Apparently, R(t) is a superposition of different states and approaches to the ground state only when t which becomes unrealistically large for large volume and/or heavy hadrons. (Here L is the number of lattice sites in one spatial direction and a is the lattice spacing.) Note that R(t) may create a "fake plateau" for intermediate values of t due to excited state contaminations [8,9,21].

IV. LATTICE SETUP
Numerical data used in this study are obtained from the (2+1)-flavor gauge configurations with Iwasaki gauge action at β = 1.82 and nonperturbatively O(a)-improved Wilson quark action with stout smearing at nearly physical quark masses [22]. The relativistic heavy quark action [23] for the charm quark [24] is used to remove cutoff errors associated with the charm quark mass up to next-to-leading order. We have a 0.0846 fm (a −1 2.333 GeV) and L = 96, which leads to La 8.1 fm. The pion, kaon, Ω 3s , and Ω 3c masses read (m π , m K , m Ω3s , m Ω3c ) (146, 525, 1712, 4796 MeV). The correlation functions are calculated by the unified contraction algorithm [25]. For the source operator J (0), we use the wall-type with the Coulomb gauge fixing.
In order to increase statistics, forward and backward propagations are averaged, the hypercubic symmetry on the lattice (4 rotations) are utilized, and multiple measurements are performed by shifting the source position along the temporal direction. The total measurements for 2Ω 3s (2Ω 3c ) amounts to 307,200 (896), where 2Ω 3Q is a shorthand notation for Ω 3Q Ω 3Q . Note that current statistics for 2Ω 3s are twice as those in Ref. [13]. The statistical errors are evaluated by the jackknife method throughout this paper. For more numerical details, see Refs. [13,14]. 1. (Color online). The one-dimensional projection of the LO potential in the A1-rep, V (r) ≡ V (r)| r=|r| , for 2Ω3s [13] at Euclidean time t/a = 17 (red circles), and for 2Ω3c [14] at Euclidean time t/a = 26 (green squares).

FIG.
In Fig. 1, we show the one-dimensional projection of the LO potential on the three-dimensional lattice box, V (r) ≡ V (r)| r=|r| , obtained by the time-dependent HAL QCD method in the A 1 representation of the cubic group SO(3, Z) (A 1 -rep in short) for 2Ω 3s [13] at t/a = 17, and for 2Ω 3c [14] at t/a = 26. These Euclidean times are chosen such that they are large enough to suppress contaminations from excited states in the single-baryon correlator and simultaneously small enough to avoid exponentially increasing statistical errors. The potential for 2Ω 3c is shorter ranged with smaller repulsive core than that for 2Ω 3s . In both systems, the LO potentials V (r) are localized around the origin and are approximately spherical functions as shown in Fig. 2 for V (x, y, z = 0). Also, V (r) in Fig. 1 can be well-fitted by the three-range Gaussians, V fit (r) = i=1,2,3 α i exp(−β i r 2 ) [13,14]. Due to the large cancellation between the medium-range attraction and the short-range repulsion, only one loosely bound state appears in the infinite volume (L = ∞) for  [13] at Euclidean time t/a = 17 (left), and for 2Ω3c [14] at Euclidean time t/a = 26 (right).

A. Eigenfunctions of HLO on the lattice
Using the LO potential V (r), the discretized Hamiltonian H LO ≡ H 0 + V can be diagonalized on the threedimensional lattice box with the size (La) 3 (8.1fm) 3 and L = 96. Under the periodic boundary condition, we thus obtain the eigenenergies ∆E n and associated eigenfunctions ψ n (r) in the A 1 -rep. Note that ψ n here represents an approximation of the exact n-th NBS wave function Ψ En associated with the LO potential V as discussed in Sec. III. The eigenfunctions are normalized r |ψ n (r)| 2 = 1 with a convention ψ n (0) > 0. Shown in Fig. 3 are the first four states (n = 0, 1, 2, 3), ψ n (x, y, z = 0), together with ∆E n . The eigenfunctions are distorted by the boundary condition and have only discrete rotational symmetry. Therefore those in the A 1rep not only contain component with angular momentum l = 0 but also components with l = 4, 6, · · · . Such a mixing becomes prominent as r and/or n increase. Shown in Fig. 4 are the one-dimensional projection of the above eigenfunctions, ψ n (r) ≡ ψ n (r)| r=|r| on the three-dimensional lattice box. Note that the number of nodes of ψ n (r) is equal to n as expected from the quantization condition given by the periodic boundary condition. The characteristic size of the ground state is smaller for the 2Ω 3c than that for 2Ω 3s . Shown together by the black solid lines are the bound-state eigenfunctions ψ inf. (r) in the infinite (L = ∞) and continuous (a = 0) space obtained by solving the Schrödinger equation with H ∞ LO = −∇ 2 /m Ω + V fit (r). We find that ψ inf. (r) and ψ 0 (r) are indistinguishable for r < 3 fm in both cases.

B. Effective energies on the lattice
Let us now utilize the eigenfunctions ψ n (r) to define optimized two-baryon sink operators S n (t) in Eq.(9) and evaluate temporal correlators R n (t) in Eq.(10) to derive the effective energy ∆E eff n (t). Shown in Fig. 5 by the open squares and open diamonds as a function of t/a are ∆E eff n (t) for n = 0 and 1, respectively. To be consistent with the Euclidean time employed to extract V (r), we choose t/a = 17 ± 3 (26 ± 3) for 2Ω 3s (2Ω 3c ). The colored bands in Fig. 5 are ∆E n obtained from H LO as discussed in Sec.V A. Open pentagons are the "naive" ground-state energy, ∆E(t) = 1 a ln R(t) R(t+1) , obtained from the unprojected temporal correlations, R(t) = r R(r, t).
By comparing ∆E eff n (t), ∆E n and ∆E(t) in Fig. 5, we find the following points: (i) Each effective energy (∆E eff 0 (t) and ∆E eff 1 (t)) has its own plateau for both 2Ω 3s and 2Ω 3c . In addition, good agreements between ∆E eff 0,1 (t) and ∆E 0,1 (t) are found.
(ii) The above observations indicate that different eigenstates are properly separated by the projection, and the uncertainty from the derivative expansion of U is within the statistical errors. Otherwise, eigenfunctions ψ 0,1 (r) obtained from V would be different from the exact NBS wave functions Ψ E0,1 (r).
(iii) Although ∆E(t) looks stable, its magnitude does not agree with the ground state energy ∆E 0 . This is the typical "fake plateaux" behavior due to the contamination from the excited states in R(t). This phenomenon for two-baryon systems was extensively studies in Refs. [8,9,21].
The same analyses for n = 2 and 3 are given in Appendix A where the agreement between ∆E eff n (t) and ∆E n is also observed within the statistical errors.

C. Overlapping factors
The coefficients a n in Eq.(10) and b n in Eq. (12) represent the overlapping strengths of the projected and un- They are normalized as r |ψn(r)| 2 = 1 with a convention ψn(0) > 0. The red, green, blue, and purple wireframes correspond to ψn(x, y, z = 0) with n = 0, 1, 2, and 3, respectively. The corresponding eigenenergy ∆En is shown at the top of each panel.
projected sink operators to the state created by walltype source, respectively. By using the eigenfunctions obtained in Sec. V A, these coefficients can be calculated as a n = r ψ † n (r)R(r, t)e (∆En)t and b n = a n r ψ n (r) [9]. Note that the source and sink operators are not Hermitian conjugate with each other, so that the overlapping coefficients are not necessarily positive. To see the relative magnitude of the coefficients, we plot |a n /a 0 | and |b n /b 0 | versus ∆E n in Fig. 6. These ratios for 2Ω 3s (2Ω 3c ) are evaluated at the center of the plateau, t/a = 17 (t/a = 26), in Fig. 5.
From the red squares for 2Ω 3s in Fig. 6, one finds that the wall-type source couples primarily to the n = 0 state, secondly to the n = 1 state with |a 1 /a 0 | 0.4, and negligibly to the higher states |a 2,3 /a 0 | < 0.1. On the other hand, the red squares for 2Ω 3c indicate that the walltype source couples primarily to the n = 1 state with |a 1 /a 0 | 3. This is due to the fact that the size of the n = 0 state for 2Ω 3c is rather small, so that its coupling to the extended wall-type source is weak. As we mentioned in Sec. II, the advantage of the HAL QCD method is that the decomposition of R(r, t) into each eigenstate is not required to derive the potential, since the potential is independent of n as long as the system is in the elastic region. Present results show that the interaction potentials can be extracted reliably in the time-dependent HAL QCD method with the derivative expansion for both 2Ω 3s and 2Ω 3c , regardless whether the ground state or the (first) excited state dominates the two-baryon correlation function.
Let us now turn to the discussion on the unprojected temporal correlation R(t) which produces the fake plateau, i.e. the black pentagons in Fig. 5. From |b 1 /b 0 | in Fig. 6, one may evaluate the value of t/a above which the true ground state saturation is reached by 14   ∆E(t). Let us demand the systematic error in the effective energy due to the first excited state contamination is bounded, in the case of b 1 /b 0 > 0 and ∆E 0 < 0. Taking ε = 0.1 and using ∆E 0,1 in Fig. 3, and b 1 /b 0 0.2 (6) for 2Ω 3s (2Ω 3c ), we need to have t/a > 580 for 2Ω 3s and t/a > 1500 for 2Ω 3c , respectively. However, R(t) for such large t would suffer from exponentially large statistical errors [7].

VI. SUMMARY AND CONCLUDING REMARKS
In this paper, we explored the idea of the optimized interpolating operators originally proposed in [9] on the basis of the time dependent HAL QCD method. To reduce the statistical errors, we considered heavy dibaryons Ω 3Q Ω 3Q (Q = s, c) in the 1 S 0 channel and extracted the leading-order HAL QCD potential V (r) localized in space. It was then used to obtain the eigenfunctions and eigenenergies on a finite lattice box with the periodic boundary condition. The eigenfunctions ψ n (r) were then used to construct the projected two-baryon sink operators S n (t) such that they couple predominantly to each n-th state. The temporal correlations R n (t) with such optimized sink operators are used to extract corresponding effective energies ∆E eff n (t) which were found to have plateau structure for each n. Moreover, they agree quantitatively with ∆E n calculated by solving the Schrödinger equation on the lattice with V (r). Such a feature provides an indirect evidence that the leading-order poten-tial V is a good approximation of the non-local potential U for 2Ω 3Q systems within the statistical errors. It also implies that the baryon-baryon interaction potential can be reliably extracted from the time-dependent HAL QCD method, no matter whether the two-baryon correlation function is dominated by the ground state or the (first) excited state, while the analysis of the unprojected temporal correlation R(t) leads to a wrong effective energy associated with a fake plateau.
Although we applied our projection only to the sink operators in this paper, one may apply the idea to the source operators as well to further improve the stability and accuracy of ∆E eff n (t). Such an approach can also provide an optimized operator basis (S n ) for the conventional variational method [26] to study hadron-hadron interactions through the matrix correlation, C nn (t) = 0|S n (t)S † n (0)|0 . Further Investigation along this line will be reported elsehwhere.