The Yakubovsky equations provide a non-perturbative framework for investigating few-body bound and scattering systems in different sectors of physics, including atomic, nuclear, and particle physics. These equations have been utilized extensively in both momentum [1–24] and configuration [25–45] spaces. The Yakubovsky equations fundamentally rely on two-body (2B) transition matrices, denoted as t(ϵ), which are derived from the solution of the Lippmann-Schwinger (LS) equation, considering either positive (scattering states) or negative (bound states) 2B subsystem energies ϵ. Notably, solving the LS equation for positive energies can prove to be a numerically demanding task due to the presence of singularities. Contrarily, for negative energies, the LS equation must be computed for the 2B subsystem energies, which are determined by the Jacobi momenta of the third and fourth particles.

In this study, we utilize a version of the Yakubovsky equations for four-body (4B) bound states that directly incorporates 2B interactions, eliminating the need for 2B t − matrices. This t − matrix-free approach to the Yakubovsky equations has been previously solved in configuration space by Lazauskas et al. [29, 39, 40, 44]. In our work, we present the t − matrix-free coupled Yakubovsky equations in momentum space. Here, the Yakubovsky components linked to the 3 + 1 and 2 + 2 chains are derived as a function of the Jacobi momentum vectors, directly including 2B interactions. To determine the 4B binding energies, we solve the coupled Yakubovsky integral equations using separable potentials with Yamaguchi and Gaussian form factors, as well as the non-separable Malfliet-Tjon potential, with all potentials projected into the s − wave channel. Our numerical findings highlight that the t − matrix-free version of the Yakubovsky integral equations, when utilizing 2B interactions, aligns perfectly with results derived from the conventional Yakubovsky integral equations that employ 2B t − matrices. In a related research, this t − matrix-free approach is successfully employed for relativistic three-body (3B) bound states [46]. This led to a version of the relativistic Faddeev equation that directly employs 2B boosted interactions [47, 48], eliminating the need for 2B boosted t − matrices [19, 20, 24].

This paper has been structured into several sections. Section 2. provides a brief overview of the Yakubovsky equations for 4B bound states, bypassing the use of 2B t − matrices. Herein, we formulate the coupled Yakubovsky equations in momentum space, considering Jacobi momentum vectors and presenting a simplified form for s − wave interactions. Section 3. Outlines our numerical results for 4B binding energies as calculated from the t − matrix-free Yakubovsky equations, alongside a comparison with results obtained from conventional t − matrix-dependent Yakubovsky equations. Lastly, Section 4. presents our conclusion and discusses prospects for succeeding studies.

The conventional form of Yakubovsky equations that describe the bound state of four identical particles interacting through pairwise interactions reads as follows [4]. ψ 1 = G 0 t P 1 + P 34 ψ 1 + ψ 2 , ψ 2 = G 0 t P ̃ 1 + P 34 ψ 1 + ψ 2 , (1) where ψ ₁ and ψ ₂ denote Yakubovsky components of the 4B wave function, corresponding to 3 + 1 and 2 + 2 type chains, respectively. G 0 = ( E − H 0 ) − 1 represents 4B free propagator, while P = P ₁₂ P ₂₃ + P ₁₃ P ₂₃, P ̃ = P 13 P 24 and P ₃₄ serve as the permutation operators. The 2B t − matrix is determined by the LS equation t = V + V G 0 t . (2) The coupled Yakubovsky Eq. 1 can be restructured to yield another form of Yakubovsky equations as follows ψ 1 = G 0 V P 1 + P 34 ψ 1 + ψ 2 + G 0 V ψ 1 , ψ 2 = G 0 V P ̃ 1 + P 34 ψ 1 + ψ 2 + G 0 V ψ 2 , (3) where the 2B interactions V are being directly utilized as input to the Yakubovsky equations, consequently eliminating the need for the 2B t − matrices. The representation of the t − matrix-free form of the coupled Yakubovsky Eq. 3 in momentum space leads to the following coupled 3D integral equations ψ 1 u 1 , u 2 , u 3 = G 0 u 1 , u 2 , u 3 ∫ d 3 u 2 ′ V s u 1 , 1 2 u 2 + u 2 ′ × { ψ 1 u 2 + 1 2 u 2 ′ , u 2 ′ , u 3 + ψ 1 u 2 + 1 2 u 2 ′ , 1 3 u 2 ′ + 8 9 u 3 , u 2 ′ − 1 3 u 3 + ψ 2 u 2 + 1 2 u 2 ′ , − u 2 ′ − 2 3 u 3 , 1 2 u 2 ′ − 2 3 u 3 } + 1 2 ∫ d 3 u 1 ′ V s u 1 , u 1 ′ ψ 1 u 1 ′ , u 2 , u 3 , ⟨ v 1 , v 2 , v 3 | ψ 2 ⟩ = 1 2 G 0 v 1 , v 2 , v 3 ∫ d 3 v 3 ′ V s v 1 , v 3 ′ × 2 ψ 1 v 3 , 2 3 v 2 + 2 3 v 3 ′ , 1 2 v 2 − v 3 ′ + ψ 2 v 3 , − v 2 , v 3 ′ + ∫ d 3 v 1 ′ V s v 1 , v 1 ′ ψ 1 v 1 ′ , v 2 , v 3 ) , (4)

with the symmetrized 2B interaction defined as V _s (a, b) = V(a, b) + V(a, − b) and 4B free propagators characterized through the following expressions G 0 u 1 , u 2 , u 3 = E − u 1 2 m − 3 u 2 2 4 m − 2 u 3 2 3 m − 1 , G 0 v 1 , v 2 , v 3 = E − v 1 2 m − v 2 2 2 m − v 3 2 m − 1 . (5) The Jacobi momenta u _i and v _i (i = 1, 2, 3) correspond to 3 + 1 and 2 + 2 chains, respectively [11]. The coupled Yakubovsky Eq. 4 can be simplified for the s − wave interactions as ψ 1 u 1 , u 2 , u 3 = 4 π G 0 u 1 , u 2 , u 3 × ∫ 0 ∞ d u 2 ′ u 2 ′ 2 ∫ − 1 1 d x ′ V u 1 , Π 1 u 2 , u 2 ′ , x ′ × ψ 1 Π 1 u 2 ′ , u 2 , x ′ , u 2 ′ , u 3 + 1 2 ∫ − 1 1 d x ψ 1 ( Π 1 u 2 ′ , u 2 , x ′ , Π 2 u 2 ′ , u 3 , x , Π 3 u 2 ′ , u 3 , x ) + ψ 2 Π 1 u 2 ′ , u 2 , x ′ , Π 4 u 2 ′ , u 3 , x , Π 5 u 2 ′ , u 3 , x } ] + ∫ 0 ∞ d u 1 ′ u 1 ′ 2 V u 1 , u 1 ′ ψ 1 u 1 ′ , u 2 , u 3 , ψ 2 v 1 , v 2 , v 3 = 4 π G 0 v 1 , v 2 , v 3 × ∫ 0 ∞ d v 3 ′ v 3 ′ 2 V v 1 , v 3 ′ × ∫ − 1 1 d x 3 ′ ψ 1 v 3 , Π 6 v 2 , v 3 ′ , x 3 ′ , Π 7 v 2 , v 3 ′ , x ) + ψ 2 v 3 , v 2 , v 3 ′ ] + ∫ 0 ∞ d v 1 ′ v 1 ′ 2 V v 1 , v 1 ′ ψ 2 v 1 ′ , v 2 , v 3 , (6)

where the shifted momentum arguments are defined as [49]. Π ̃ 1 = 1 4 u 2 2 + u 2 ′ 2 + u 2 u 2 ′ x 1 2 , Π 1 = u 2 2 + 1 4 u 2 ′ 2 + u 2 u 2 ′ x 1 2 , Π 2 = 1 3 u 2 2 + 64 9 u 3 2 + 16 3 u 2 u 3 x ′ 1 2 , Π 3 = u 2 2 + 1 9 u 3 2 − 2 3 u 2 u 3 x ′ 1 2 , Π 4 = 1 4 u 2 2 + 4 9 u 3 2 + 2 3 u 2 u 3 x ′ 1 2 , Π 5 = u 2 2 + 4 9 u 3 2 − 4 3 u 2 u 3 x ′ 1 2 , Π 6 = 2 3 v 2 2 + v 3 2 + 2 v 2 v 3 x 1 2 , Π 7 = v 3 2 + 1 4 v 2 2 − v 2 v 3 x 1 2 . (7)

For comparison purposes, the representation of the conventional t − matrix-dependent form of the coupled Yakubovsky Eq. 1 in momentum space can be found in Appendix A. Upon comparing the t − matrix-free and conventional t − matrix-dependent forms of the Yakubovsky equations - specifically, Eqs. 4 and (A1) - it is evident that the t − matrix-free form incorporates an extra term. This term involves integration over the 2B interaction and the Yakubovsky components without interpolations on momenta or angles. Despite this, its numerical solution proves to be more straightforward and cost-effective than that of the conventional t − matrix-dependent that necessitates solving the LS equation to compute the 2B t − matrices for all required 2B subsystem energies, which depend on the magnitude of the second and third Jacobi momenta.

The numerical solution of the coupled Yakubovsky integral Eq. 6 for the calculation of 4B binding energy demands solving an eigenvalue equation, where the physical binding energy corresponds to an eigenvalue equal to one. The Lanczos iterative method is implemented for solving such eigenvalue equation [24, 49–52]. We employ Gauss-Legendre quadratures to discretize the continuous momentum and angular variables with a hyperbolic mapping for Jacobi momenta and a linear mapping for angle variables [53]. This allows us to properly capture the behavior of the Yakubovsky components of the 4B wave function at both small and large momenta. In each iteration step of solving the coupled Yakubovsky integral equations, to accurately perform multi-dimensional interpolations on the shifted momentum arguments given in Eq. 7, we employ the Cubic-Hermite spline method due to its combination of high accuracy and computational speed [54].

Our numerical analysis presents a comparison between 4B binding energy obtained from the t − matrix-free approach and the conventional t − matrix-dependent formulation of the coupled Yakubovsky equations, namely, Eqs. 6 and A4. For our numerical analysis, we employ two models of one-term separable potential with the following general form V p , p ′ = λ g p g p ′ , (8) where λ represents the potential strength. The potential form factor g(p) for the Yamaguchi-type potential is defined as g(p) = 1/(p ² + β ²) [55], while for the Gaussian potential, it takes the form g(p) = exp(−p ²/Λ²) [56]. Furthermore, to provide a comprehensive validation of our formalism and code, we also incorporate an s-wave non-separable Malfliet-Tjon (MT) potential, comprises two attractive and repulsive terms [57]. V p , p ′ = ∑ i = 1 2 λ i 2 π p p ′ ln μ i 2 + p 2 + p ′ 2 + 2 p p ′ μ i 2 + p 2 + p ′ 2 − 2 p p ′ . (9) In Table 1, we provide the parameters for the potentials employed in our calculations, which include MT model V (MT-V), Yamaguchi potential model IV (Y-IV), and a Gaussian potential. The strength of the Gaussian potential was adjusted to reproduce the deuteron binding energy of −2.225 MeV, with a form factor parameter Λ = 0.7 fm⁻¹.

Table 2 presents the convergence of the 4B binding energy as a function of the number of mesh points for the Jacobi momenta magnitudes u _i and v _i . The table provides a side-by-side comparison of results obtained using both the t-matrix-free and conventional t-matrix-dependent forms of the coupled Yakubovsky equations. Specifically, our results show that the t − matrix-free Yakubovsky equations yield a 4B binding energy of −30.08 MeV for the MT-V potential and −36.27 MeV for the Yamaguchi-IV potential. These values closely align with the −30.07 MeV [58] and −36.3 MeV [59] obtained by other groups using the t − matrix-dependent Yakubovsky calculations. The comparison validates the potential of the t − matrix-free adoption of Yakubovsky formulation to solve 4B bound state problems efficiently, matching the precision of the conventional method but potentially offering more straightforward computational requirements. Moreover, the convergence behavior remains consistent across different numbers of mesh points for the magnitude of Jacobi momenta, emphasizing the computational robustness of the t − matrix-free formulation.

The Yakubovsky approach is a powerful method to study few-body bound and scattering systems. However, the solution of these equations can be computationally demanding due to inherent singularities in 2B t − matrices when dealing with scattering problems and the need to calculate them for 2B subsystem energies dictated by second and third Jacobi momenta when 4B bound state problems are considered. This study utilizes a version of the coupled Yakubovsky equations for 4B bound states that directly incorporates 2B interactions in momentum space, avoiding the use of the 2B t − matrices. The efficacy of this approach is validated through the calculation of 4B binding energies in momentum space using, both the separable potentials with Yamaguchi and Gaussian form factors, and the non-separable Malfliet-Tjon potential. Our findings align well with results from the conventional form of the coupled Yakubovsky integral equations incorporating 2B t − matrices. The extension of calculations to include more general interactions, beyond just the s − wave, is currently in progress.