Study of virtual annihilation and retardation effects in relativistic two-, four-, and six-body wave equations in scalar QFT

The variational method within the Hamiltonian formalism of QFT using Darewych reformulated model has been used for scalar particles and antiparticles interacting through a scalar mediating field. We have investigated the relativistic effects such as virtual annihilation interactions and retardation effects for relativistic two-, four-, and six-body wave equations in scalar QFT. Approximate ground-state solutions have been studied for different strengths of coupling, for both massive and massless mediating fields where the virtual annihilation terms or retardation effects in the wave equations have been included or eliminated.


Introduction
The study of few-or many-body systems in relativistic quantum field theory (QFT) is of fundamental interest and has played a major role in the progress of modern physics. It contributes to the development of our idea about the nature of particle interactions. For example, we can mention that the thorough study of the spectrum of the hydrogen atom was an early test of the theory of quantum electrodynamics (QED). The Bethe-Salpeter (BS) formalism [1,2] is the traditional tool for analyzing relativistic bound state, but this method is hard to implement and is all but intractable for systems that contain more than two particles. The study of two-body bound states using BS equation in the ladder approximation has been known as the Wick-Cutkosky model (WC) [3,4].
The relativistic quantum-mechanical treatment of bound states of a particle in a given field is principally based on the Dirac and Klein-Gordon equations. The Dirac equation is suitable for particles, like electrons, with spin 1/2 and the Klein-Gordon equation is used to solve problems of particles without spin like p þ or p À : The Klein-Gordon equation is also called the scalar equation of relativistic quantum mechanics.
An important problem in quantum field theory has been the description of relativistic n-body scalar particles wave equations or n-body scalar particles and antiparticles wave equations (spin-0 bosons), including the interactions, from the underlying Lagrangian. In previous works by Emami-Razavi [5,6], it has been demonstrated that, by employing a reformulated model in QFT proposed by Darewych [7,8], formulas concerning relativistic n-body wave equations for scalar particles and/or antiparticles can be obtained (see also [9]). It has also been proven that using Darewych's formalism [7,8] relativistic n-fermion wave equations (particles and antiparticles; spin-1/2) in quantum electrodynamics can be obtained [10,11]. One of the appealing features of having n-body wave equations formula is that we can explicitly write the relativistic wave equations for an arbitrary number of n including the interactions. Therefore, one can easily write, for example four-or six-body wave equations, and thereupon one can perform different calculations (employing various numerical methods) to study relativistic effects such as retardation effects or virtual annihilation interaction terms.
In this paper we study the effects of inclusion or exclusion of retardation terms or virtual annihilation interactions for relativistic two-, four-and six-body scalar particles and antiparticles which interact through a scalar mediating field which can be massless (l ¼ 0) or massive (l=m ¼ 0:15). The case of two-body system has been investigated before [12,13]. However, we remind it again here to better follow the four-and six-body cases. In this work the coupling constant that shows the strength of the interaction is presented by the symbol a which is related to another constant g by the relation a ¼ g 2 =16pm 2 : a can vary from a relatively small coupling (such as a $ 0:1) to strong coupling numbers such as a $ 1: Our calculations show that even though the retardation effects or virtual annihilation interactions can be small effects relative to Yukawa interactions at low coupling for the two-, four-, and six-body systems, for large values of a (0:4\a 1) the magnitude of the contribution of virtual annihilation terms and retardation effects augments for the ground-state energy solutions of the systems under study.
The presentation of the manuscript is as follows. The first section is the introduction. The second section is about variational methods and the reformulated scalar model. The relativistic wave equations, virtual annihilation interactions and retardation effects are reminded in the following section. In the subsequent section, the results of approximate ground-state energy solutions for various cases have been discussed. Concluding remarks have been presented in the last section. Appendix contains some formulas and examples concerning relativistic wave equations.

Variational methods and the reformulated model
From a practical point of view, the variational method is independent of the magnitude of the coupling constant strength, contrary to perturbation theory which usually depends on expansion of some small parameters. The variational methods in Hamiltonian QFT is similar to Schrödinger type of representation of few-body systems, and can easily be applied to systems of more than two bodies [14].
The variational methods have been moderately used in QFT since the construction of realistic trial states can be a difficult task. The variational approach to the treatment of relativistic bound states has received increasing attention since the early papers on the variational methods in QFT (see, for examples, [15,16]). We should note that the variational approach is basically as good as the trial states we use [17][18][19].
As pointed out before, the variational method is only as good as the trial states that we use. In the functional formulation we are restricted mainly to wave functionals of Gaussian type, due to the difficulty of handling analytically non-Gaussian functional integrals. Another possible approach that we can employ is to expand the trial state on a Fock-space basis. For example, the Tamm's work [20] is along these lines, though it was not formulated variationally. In this paper, we have applied the approach given by Darewych [7] to particle-antiparticle bound states in the scalar Yukawa model. In this model, we have assumed that particles or antiparticles have no spin. This means that they are scalar particles or antiparticles. Moreover, the mediating fields which act between particles-antiparticles are also scalar mediating fields.
In Ref. [6] simple trial states have been employed to obtain the n-body wave equations in the scalar QFT by the variation method. Thereupon, one can easily write the relativistic wave equations for the two-, four-, and six-body systems. The equations have the Schrödinger non-relativistic limit, with Yukawa potentials in the event of a massive mediating scalar field or Coulombic interparticle potentials on condition of having a massless mediating scalar field. Ground-state energy solutions of the equations have been approximately calculated for various strengths of the coupling. A comparison of various cases have been considered where virtual annihilation terms and/or retardation effects have been included or excluded for our different systems under study.
For the sake of completeness and having a self-contained paper it is better to recall some of the main points from previous works [12,13]. The Lagrangian density for scalar particles which interact via a mediating scalar field is where g indicates the coupling constant. We are working in the unit h ¼ c ¼ 1: m is the mass of particle or antiparticle. / and / Ã present the field associated with the particle and antiparticle, respectively. vðxÞ designates the mediating field which can be massless (i.e., l ¼ 0), (cf. [3]), or massive ðl 6 ¼ 0Þ.
The Lagrangian density obtained in the reformulated model is (see [7] or [17])

A r h i v e o f S I D
where dx 0 ¼ dt d 3 r, qðxÞ ¼ Àg/ Ã ðxÞ/ðxÞ; v 0 ðxÞ denotes the free field and Dðx À x 0 Þ indicates a covariant Green function in a way that we have As it has been mentioned, for example, in [7] or [12,13] the theories based on (2.1) and (2.2) are equivalent in the sense that they lead to the same invariant matrix elements in different orders of perturbation theory. The Lagrangian density (2.2) can be written as: where we have suppressed the free part of the Lagrangian in the above equation. As a matter of fact, one should note that the trial states, jw n i, that we employ (see Eq. (3.1)) are insensitive to H v (Eq. 2.7) and H I 1 (Eq. 2.8 ) in the sense that hw n j :Ĥ v : jw n i ¼ 0 and hw n j :Ĥ I 1 : jw n i ¼ 0 for the current choice of trial states.
where dk ¼ d Nþ1 k and k 2 ¼ k m k m . To indicate our notation, we quote the usual decomposition of the fields in N þ 1 dimensions: The momentum-space operators, A y ; A; B y ; B obey the usual commutation relations, the non-vanishing ones are The Hamiltonian is normal ordered (specified by :Ĥ:) since we are not interested in vacuum-energy matters in this work.
In the Hamiltonian formalism we wish to find solutions of the equation where jW trial i is a suitable trial state with adjustable parameters or functions. We set the time t ¼ 0 for the calculations.

Relativistic wave equations, virtual annihilation terms and retardation effects
For the sake of completeness and to have a self-contained manuscript it is better to recall some main points of the nbody scalar particle and antiparticle relativistic wave equations in Ref. [6]. In the reformulated scalar QFT variational model simple trial states have been used (see Eq. 3.1) where n ¼ n 1 þ n 2 (n 1 : number of particles and n 2 : number of antiparticles, n 1 ¼ n 2 or n 2 ¼ n 1 À 1Þ. One should note that since we are studying here the two-, four-, and six-body systems the total number of particles and antiparticles in the present paper is even. We have for (n/2) particles and (n/2) antiparticles the following simple trial state [6]: where G is normalizable to unity and it should be a wellbehaved function with adjustable parameters. The matrix element for the rest-plus-kinetic energy of the n-body system is The matrix elements related to the interactions of our system have nðn À 1Þ=2 terms (all of them present attractive nature, i.e., gravity-like). We have also n 2 =4 terms (when n is an even number) analogous to virtual annihilation terms (interaction among the particle and antiparticle pairs). The matrix element for the interaction of the system under study has the following form [6]: As an example we can quote the following. In case we consider three particles and three antiparticles (n ¼ 6Þ, we will have ten interparticle interactions (labeled G 1 ; G 2 . . .; G 15 ; all of them attractive, i.e., gravity-like) and we will have six terms that correspond to virtual annihilation terms. These interactions occur among particle and antiparticle pairs (labeled A 1 ; A 2 ; . . .; A 9 ). The virtual annihilation interactions are repulsive contact potentials if l\2m.
As it can be seen, for example in [21], the description of interactions of two or more than two bodies with equal or comparable masses is much more demanding than the case of particles with non-identical masses. We cannot generally differentiate between the test particles or antiparticles and the source of the field for the case of equal or comparable masses. Furthermore, we should consider all particles on an equal footing. Moreover, the interactions among particles or antiparticles do not propagate with infinite velocity. Therefore, there are further complications for the treatment of interactions among particles or antiparticles. These phenomena lead to the terms corresponding to ''retardation effects''.
The matrix element for the interaction related to our nbody wave equations (retardation effects are included) is [6] (for j\k): (3.4) indicates that the terms with indices j and k are left out. We can write down: One should note that if n ¼ 2, we will have Q i¼1...n ðj;kÞ d N ðp 0 i À p i Þ ¼ 1. We should also consider the statistical factors [22] in our systems under study due to the fact we have identical particles or antiparticles. This means that we should not double count the states of identical particles [22]. For the case where n is an even number, the statistical factor is

A r h i v e o f S I D
1=ðn=2Þ! Â 1=ðn=2Þ!: For examples, if we have four-and six-body systems consisting of two particles-two antiparticles and three particles-three antiparticles, respectively, we will have the following statistical factors for four-body and six-body systems, respectively: 1=2! Â 1=2! and 1=3! Â 1=3!: In Eq. (3.4), i, j and k indicate indices. Even indices are considered for antiparticles and odd indices are considered for particles. The relativistic n-body wave equation for scalar particles-antiparticles can be obtained from the matrix elements as below [6] (for j\kÞ: : ð3:5Þ The first terms on the right-hand side of the above equation are related to attractive interactions among the n bodies. The second terms correspond to virtual annihilation interactions among odd-even or even-odd indices (particlesantiparticles pairs). The Eq. (3.5) has the Schrödinger-like structure, with positive-energy solutions only. This can easily be seen by putting the right-hand side of this equation to zero (i.e., no interactions). In Appendix we explicitly quote some examples of the two-, four-, and sixbody wave equations. We should mention that if the terms corresponding to retardation effects are eliminated we should substitute in the kernel of the equations the terms like 1=½l 2 À ðp 1 À p 2 Þ 2 by 1=½l 2 þ ðp 1 À p 2 Þ 2 ; and the terms like 1=½l 2 À ðp 1 þ p 2 Þ 2 by 1=½l 2 À 4m 2 : The non-relativistic limit of the relativistic n-body wave equations can be obtained (which corresponds to p 2 =m 2 \\1) by the Fourier transformation (from momentum space to coordinate space), ð3:6Þ Hence, the resulting n-body Schrödinger equation is . . .; r n Þ þ ðVðr 1 ; r 2 ; . . .; r n Þ À ÞWðr 1 ; r 2 ; . . .; r n Þ ¼ 0; ð3:7Þ where ¼ E À nm.
The non-relativistic potential Vðr 1 ; r 2 ; . . .; r n Þ can be explained as follows. The first terms of Eq. (3.8) correspond to sum of attractive Coulomb potentials for the case of massless mediating field ðl ¼ 0Þ or Yukawa potentials for the case of massive mediating field ðl 6 ¼ 0Þ among the n bodies, and repulsive (if l\2m) contact potentials, due to virtual annihilation interactions among the particle and antiparticle pairs (among odd-even or even-odd indices of j and k). The n-body system potential in coordinate space in the non-relativistic limit can be written as below: Approximate ground-state energy solutions for two-, four-, and six-body systems with and/ or without virtual annihilations and retardation effects The simplest case for the systems under study is the twobody case. As mentioned before, the case of two-body system has been investigated before [12,13]. However, it is good to recall some key points of the two-body particle and antiparticle system since it will help to better understand the four-and six-body cases.
As pointed out before, the scalar Yukawa [or Wick-Cutkosky (WC)] model [3,4] has been investigated by several researchers in different formalisms and approximations. Figure 1 of Ref. [12,13] (which is plotted here again as Fig. 1) contains various results of WC or BS formalism as well as other calculations. Moreover, the results of Di Leo and Darewych [23] correspond with the solutions obtained in [12,13] (the two-body system in the present formalism) with x p 1 ¼ x p 0 1 (i.e., no retardation). In other words, the retardation term is not included in [23]. Furthermore, Refs. [3,4] and [23] do not contain the virtual annihilation interaction terms which are repulsive interactions. The effect of this repulsive interaction is to increase the ground-state energy solutions of the two-body system by a small amount or to decrease the binding energy of the two-body system [B 2 ¼ E 2 À 2 in units of m (i.e., m ¼ 1)] by a small amount. One should note that the magnitude of contribution of virtual annihilation interactions augments with increase of the coupling constant a: This means that at very high coupling (a $ 1) the contribution of this repulsive interaction becomes more important. Figure 1 of Refs. [12,13] (also plotted here as Fig. 1) shows the retardation effects for the two-body case. The inclusion of this effect causes to decrease the ground-state energy solutions of the two-body system by a small but not negligible amount. One should note that the magnitude of contribution of retardation effects augments with increasing of the coupling constant a: A comparison of the present two-body equation [12,13] with BS [1,2] or WC results [3,4] proves that for the present case we have stronger binding energy solutions particularly at strong coupling. One should keep in mind that the obtained results compared to BS or WC model are from different sets of equations. One of the differences between variational method Eq. (3.5) in this manuscript and BS equation is clear even in the no interaction case (a ¼ 0), that is, the former equation has no negative energy solutions. Thereupon, perhaps it is not surprising that there are differences between the variational method equations here and BS and WC models, particularly at strong coupling a: Figure 1 also contains the results obtained in [24] showing that the two-body eigenstates can be written down for Hamiltonian of the present model without free fields, provided that an empty vacuum state j0i, annihilated by both negative and positive-energy components of the field operator /ðxÞ; is employed. The use of such an empty vacuum state did give both negative and positive groundstate energy solutions. Moreover, The two-body system in Ref. [24] could be solved analytically for the massless mediating field case and also the results suggest a critical value of the coupling strength around a $ 1.
All the results related to different formalisms and approaches for the two-body system mentioned above have been presented in Fig. 1 (Fig. 1 of Refs. [12,13]).
The four-body system (consisting of two particles and two antiparticles) and six-body system (consisting of three particles and three antiparticles) cannot be solved analytically. Therefore some approximate solutions should be sought. In the Appendix, four-, and six-body wave equations have been written explicitly using the Eq. (3.5) for n ¼ 4 and n ¼ 6; respectively. One can obtain variational expression for the energy E n of the n-body system (for our case: n ¼ 4 or n ¼ 6) of the systems under study by substituting Gðp 1 ; p 2 ; . . .; p n Þ with analytical functions that contain adjustable constants or parameters to compute the energy of four-or six-body system (see Eqs. 3.2 and 3.4). E n ¼ 1 hw n jw n i hw n j :Ĥ : jw n i: ð4:1Þ We consider ground-state energy solutions and using the variational principle (Eq. 2.16) and taking simple forms for the wave function G we can obtain some approximate solutions. We have, namely where p 0 and m describe adjustable parameters. The calculations have been performed in the rest frame (with zero total momentum). For the four-or six-body wave equations, the multidimensional integrals that we have cannot be solved analytically, hence we should seek some numerical solutions using the Monte Carlo method [25]. The parameter p 0 indicates the scale of the approximate wave functions. We Curves from lowest to highest: Feshbach-Villar formalism [24]; with retardation and without virtual annihilation [12,13]; with retardation and with virtual annihilation [12,13]; without virtual annihilation interactions and retardation terms [

A r h i v e o f S I D
locate the minimum of the variational trial energy as a function of p 0 . Moreover, varying the parameter m does not change the numerical solutions sensibly (see, for example, [12,13]); therefore, we leave it at the value m ¼ 2 which is the non-relativistic, hydrogenic ground-state value. For the four-body system the numerical results have been presented for different cases, ''without'' and ''with'' retardation effects and virtual annihilation interactions. The virtual annihilation interactions are mostly not included in other works. One should note that the solutions corresponding to four-body relativistic wave equations when the retardation effects are included and virtual annihilation interactions are eliminated have been done previously in [14] and they are added here for the comparison purposes with other cases in Tables 1 and 2. The numerical solutions for the variationally obtained trial wave functions are given in Table 1 for the massive-exchange case with l ¼ 0:15 m and for the massless-exchange case l ¼ 0 in Table 2. These solutions have been also plotted in Fig. 2a, b, respectively. Energies, E 4 and p 0 have been provided in units of m (i.e., m ¼ 1) for given coupling a.
Binding energy B 4 is: B 4 ¼ E 4 À 4 in units of m (i.e., m ¼ 1); the binding energy augments when we increase the coupling a: This is because of the fact that for present scalar mediating field, the dominant one-quantum exchange interactions are attractive. The parameter p 0 corresponding to our wave functions also augments when we increase the coupling strength a: This means that the wave function becomes more peaked in coordinate space at higher values of the coupling. Similar results have been presented in Table 2 and Fig. 2b for the massless-exchange case (l ¼ 0). In this case, Yukawa potentials transform to Coulombic potentials in the non-relativistic limits. Although the qualitative behaviors of EðaÞ for l ¼ 0 are very close to that for l ¼ 0:15 m, the binding energy, for a specific a is stronger for l ¼ 0 than for l ¼ 0:15 m, and it sets in at a 0 ! 0.
The numerical results in Tables 1 and 2 generally show that the retardation effects and virtual annihilation interactions are somehow minor effects but not negligible in contrast to Yukawa interactions for the scalar four-body particle-antiparticle systems. Moreover, the contribution of these retardation effects and virtual annihilation terms increases at high coupling and they become more noticeable for the range of the values of coupling strength 0:5 a 1: Furthermore, for a given value of coupling a we can observe (from Tables or plotted figures) that for the case where virtual annihilation terms are included and retardation effects are excluded we have the highest values of the energy solutions (the top curves in the Fig. 2a,  b, respectively). In addition, for the case where virtual annihilation terms are excluded and retardation effects are included we have the lowest values of the energy solutions (the lowest curves in the Fig. 2a, b, respectively). We can observe these differences of the ground-state energy values with respect to inclusion and/or exclusion of virtual annihilation terms and retardation effects particularly at higher values of the coupling strength 0:5\a 1: We should mention that the four-body system similar to the case of two-body system in [24] (that was solved analytically for the massless mediating field case) may suggest a critical value of the coupling constant around a $ 1: This means that the values of energy solutions obtained for a ! 1 may not be reliable or no energy solutions exist beyond a critical value of the coupling a $ 1: Hence, we did not consider the energy results for a larger than 1. Note also that the error estimates in the numerical Similar to the case of four-body system (two particles and two antiparticles) we have the case of six-body problem (system consisting of three particles and three antiparticles). Tables 3 and 4 show the numerical results corresponding to the effects of inclusion or exclusion of virtual annihilation interactions and/or retardation effects for the six-body system for the massive-exchange case with l ¼ 0:15 m and for the massless-exchange case l ¼ 0, respectively. The wave function parameters p 0 are also available in Tables 3 and 4 for each case at a given coupling strength a: Note that as for the previous cases of two-, and fourbody systems the virtual annihilation terms have small contributions for the ground-state energy results. The virtual annihilation term is a ''repulsive'' interaction that increases the energy of the system by a small amount, particularly for coupling values around 0:1 a 0:4: However, the contribution of this interaction becomes more important (not negligible) at higher coupling strength values such as 0:5\a 1: Moreover, the contribution of retardation effects become more significant at higher values of coupling a: In addition, for a given value of coupling a if we consider that virtual annihilation interactions are included in the calculations, we will have two different cases. The first case is when retardation term is added and the second one is when retardation term is neglected. Similar to the case of twoand four-body systems the inclusion of retardation terms do decrease the energy solutions of the six-body system. This means that the binding energy of the systems in the example mentioned above will increase for both cases: the massive-exchange case with l ¼ 0:15 m and for the massless exchange case l ¼ 0, respectively. Note that binding energy B 6 is: B 6 ¼ E 6 À 6 in units of m (i.e., m ¼ 1).
The numerical solutions have been presented in Tables 3  and 4. These values have been also plotted in Fig 3a, b for the massive-exchange case with l ¼ 0:15 m and for the massless-exchange case l ¼ 0, respectively.
The case of six-body system, for particles only and no virtual annihilation interaction [5], was done before and it has been included in the Tables 3 and 4 as well as the calculation of the case of combined three particles and three antiparticles with virtual annihilations and with retardation effects [9]. The corresponding graphs are also put in Fig. 3a, b for comparison purposes. However, the other two cases when the retardation terms are included or excluded have not done previously and correspond to the present work.
We should note that for the six-body system under study, as for the two-and four-body cases mentioned before, we may have a critical value of the coupling around a $ 1: Note also that the error estimates in the numerical results obtained by Monte Carlo method is growing with increase of the coupling strength values. Furthermore, the error estimates augment notably with the increasing number of particles and antiparticles. This means that, for example for a ¼ 1; the error bars in the graphs related to Table 2 Four-body E 4 ground-state energies and wave function parameter, p 0 , for the massless-exchange case, l ¼ 0

A r h i v e o f S I D
the six-body system are much larger than the error bars corresponding to the four-body system investigated earlier in this section. Binding energy of the six-body system B 6 is: B 6 ¼ E 6 À 6 in units of m (i.e., m ¼ 1). As one can observe from the results of six-body system the binding energy of the system augments with increase of coupling strength a similar to the previous cases of two-and four-body systems. Although the qualitative behaviors of EðaÞ for l ¼ 0 are very close to that for l ¼ 0:15 m; the binding energy for a particular value of a is stronger for l ¼ 0 than for l ¼ 0:15 m; and it sets in at a 0 ! 0:

Concluding remarks
Approximate ground-state solutions of the two-, four-, and six-body systems have been presented for different strength of coupling, for both massive and massless scalar mediating field to compare the effects of virtual annihilation interaction terms and retardation effects in the relativistic equations. The numerical solutions of the two-, four-, and six-body systems have been discussed considering that ''the virtual annihilation interactions'' and ''the retardation effects'' are neglected or included.

A r h i v e o f S I D
The results for our various systems under study show that at low coupling the retardation effects and virtual annihilation interactions are minor effects compared to Yukawa interactions in our particle-antiparticle systems. In addition, at high coupling the results show that we have a noticeable increase in the effects of virtual annihilation interactions or retardation terms in the two-, four-, and sixbody systems.
As mentioned before, the effect of retardation terms is as follows. The binding energies of our systems under study have been increased by a small amount for low values of the coupling and by a noticeable amount for high values of the coupling. However, as we can see, for example in the work of Gross [26], the inclusion of retardation effects reduces the binding energy of the two-body system. The argument that we can use to express the difference between these results with respect to the inclusion of the retardation effects is as follows. In a general point of view, we should consider the types of interactions that have been used. As a matter of fact, retardation terms that did arise in the equations are relativistic effects. The relativistic equations in the current model have only positive-energy solutions. In other words, the present relativistic equations are Salpeterlike and not Klein-Gordon like. The relativistic results lie above non-relativistic solutions, as we can see, for example, in [5] or [12,13]. This fact, that the relativistic solutions are above the non-relativistic results, is the characteristic of formalisms, like the present model, which has no negative energy answers. One should also note that for the two-body case, the relativistic solutions which are obtained by employing an ''empty'' vacuum (Klein-Gordon-like equations) lie below the non-relativistic solutions. As one can see from eq. (50) of Darewych [24] for the massless-exchange case (the same as Todorov equation [27] using a quasipotential approach), at small values of the coupling a, the energy solutions are: E 2 ¼ 2 À 1 4 a 2 À 5 64 a 4 : To sum up, the retardation effect augments the binding energy of the systems mostly due to the characteristic of the present model or formalism.
From the practical point of view, it is of interest to remind what follows. Wheeler [28] speculated that two positronium atoms can combine to form the four-body system, namely positronium molecule (Ps 2 : e À e þ e À e þ ), in 1946. He also postulated the possible existence of larger atoms such as the six-body systems Ps 3 (e À e þ e À e þ e À e þ ). The first calculation of the binding energy of Ps 2 was performed [29] later on. In 2007 the positronium molecule was observed by Cassidy and Mills [30]. Similarly, the existence of a scalar ''fundamental'' particle, named Higgs boson, was theorized in 1960 by Higgs [31] and some other scientists. Observation of this new particle was confirmed in 2012 [32]. Hence, it would be also of interest to investigate the existence of bound states and related relativistic effects of the few-body Higgs bosons such as fourbody Higgs systems or even larger structures. For example, the existence of two-Higgs-boson bound states has been studied in references [33,34].
Therefore, from experimental and theoretical points of view it will be of fundamental interest to study the existence of bound states of scalar particles (or other exotic systems [35]) and the corresponding relativistic effects. Moreover, since the bound-state solutions, for example, of the scalar Higgs bosons are mainly found in the very strong coupling regime [34] the use of variational method seems to be more appropriate for the computations of relativistic effects of such ''scalar'' particles and/or antiparticles (the perturbation theory does depend on the coupling constant and the calculations may not be reliable) particularly at strong coupling.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix
For the two-body system (one particle and one antiparticle), choosing the trial state j w 2 i ¼ Z d N p 1 d N p 2 Gðp 1 ; p 2 ÞA y ðp 1 ÞB y ðp 2 Þj0i; ð6:1Þ we obtain the following wave equation: