Nonlocal Fractional Quantum Field Theory and Converging Perturbation Series

The main purpose of this paper is to derive a new perturbation theory (PT) that has converging series. Such series arise in the nonlocal scalar quantum field theory (QFT) with fractional power potential. We construct PT for the generating functional (GF) of complete Green functions (including disconnected parts of functions) $\mathcal{Z}\left[j\right]$ as well as for GF of connected Green functions $\mathcal{G}\left[j\right]=\ln \mathcal{Z}\left[j\right]$ in powers of coupling constant $g$. It has infrared (IR)-finite terms. We prove that the obtained series, which has the form of a grand canonical partition function (GCPF), is dominated by a convergent series, in other words, has majorant, which allows to expand beyond the weak coupling $g$ limit. Vacuum energy density in second order in $g$ is calculated and researched for different types of Gaussian part $S_{0}[\phi]$ of the action $S[\phi]$. Further in the paper, using the polynomial expansion, the general calculable series for $\mathcal{G}\left[j\right]$ is derived. We provide, compare and research simplifications in cases of second-degree polynomial and hard-sphere gas (HSG) approximations. The developed formalism allows us to research the physical properties of the considering system across the entire range of coupling constant $g$, in particular, the vacuum energy density.


Introduction
Why do we need a nonlocal quantum field theory (QFT for short) in physics? Nonlocal QFT is a perturbatively, in each order of the perturbation theory (PT for short) with respect to the small coupling constant g, nonrenormalizable theory. The set of such theories is much wider than the set of perturbatively renormalizable ones. In particular, a nonlocal QFT can be formulated in a spacetime of arbitrary dimension d. Hence, this research is in some sense "orthogonal" to two-dimensional d = 2 conformal and integrable QFTs. Thus, nonlocal QFT was created precisely in order to consider the spacetime dimension d as an arbitrary parameter .
Further, nonlocal QFT makes it possible to describe the high-energy physics of scalar particles (scalar nonlocal QFT) and in quantum electrodynamics (QED for short), as well as low-energy physics in quantum chromodynamics (QCD for short). For example, to describe the low-energy physics of light hadrons in QCD and quark confinement, the so-called Virton-Quark Model and the Quark Confinement Model of Hadrons are generally accepted [3,[6][7][8]. Nowadays, research in the field of nonlocal quantum theory of gravity and nonlocal QFT in curved spacetime is also gaining popularity [14,16,[23][24][25].

arXiv:2303.16011v3 [hep-th] 21 Aug 2023
This duality is well known, for example, in statistical physics and plasma physics, and its essence is as follows: S-matrix of nonlocal nonpolynomial Euclidean QFT, originally formulated in terms of the functional (path) integral, can be rewritten in terms of the grand canonical partition function (GCPF for short) of some "gas with interaction". This duality is used in both directions. On the one hand, this makes it possible to apply the methods of statistical physics in QFT, in particular, to use the cluster expansion when calculating the S-matrix. On the other hand, this makes it possible to apply QFT methods in statistical physics, in particular, to use the saddle-point method to calculate the functional (path) integral equal to the GCPF. Here let us note that this duality is used in the research of sin-Gordon and sinh-Gordon models, more precisely, their nonlocal versions.
Finally, within the framework of nonlocal QFT, it is possible to research the strong coupling limit in interaction constant g [3,[19][20][21][22]. The method used to research the strong coupling limit in nonlocal QFT originated in the theory of polarons (in polaronics). This method is based on finding the minorant (lower estimate) and majorant (upper estimate) for the quantity of interest to us, as functions of the coupling constant g. The Jensen and Hölder inequalities are often used to find these estimates. If both estimates tend to each other in the limit g → +∞, then each of them coincides with the quantity of interest to us (in this limit). In fact, this is the squeeze theorem. For example, the vacuum energy E (g) was obtained in the strong coupling limit in [3,19,20].
Why do we need a nonlocal quantum field theory in mathematics? The rigorous definition of the functional integral in a nonlocal QFT, the proof of the definition correctness (the existence and uniqueness of such an object), as well as the calculation of the latter, is a complicated but interesting mathematical problem [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. Let us note that even the calculation of such an object needs some definition. In what terms is the functional integral considered calculated? As quoted above, there is extensive mathematical literature devoted to the theory of integration in separable Hilbert and Banach spaces (HS and BS for short) with respect to Gaussian measures. At the same time, the "list of integrals" given in such literature is hardly sufficient for practical applications: results are usually given in cases where the integrand is a polynomial or again a Gaussian exponent, which corresponds to a Gaussian (free, without interaction) QFT. Our paper aims to narrow this gap a bit by adding another result to the "list of integrals".
Another complicated but interesting mathematical problem is the rigorous definition of the Dyson-Schwinger, Schwinger-Tomonaga and functional (nonperturbative, exact) renormalization group flow (FRG flow for short) equations in nonlocal QFT [42,44,45,[48][49][50][51][52][53]. These equations are formulated in terms of variational derivatives, therefore, they are differential. Again, as quoted above, there is extensive mathematical literature devoted to differential equations in infinite-dimensional spaces (separable HS and BS). However, as in the case of integrals, the "list of solutions" given in such literature is hardly sufficient for practical applications. Obtaining solutions to such equations in terms of functional Taylor series is associated with solving infinite hierarchies (systems) of integro-differential equations for various families of Green functions. Such hierarchies are usually solved approximately. Let us note that FRG flow equations formulated as nonlocal QFT appear in many physical sciences: QFT itself, condensed matter physics, critical behaviour theory, stochastic theory of turbulence, etc.
Finally, in the framework of nonlocal QFT research in mathematics, let us note one possible generalization to the case of locally convex topological vector space (LCTVS for short) [54]. Such a seemingly mathematical generalization can lead to far-reaching consequences in physics. The topology of LCTVS is given by a family (set of arbitrary cardinality) of seminorms, generalizing separable HS and BS. In fact, instead of one scalar product or norm, a family of such arises. Such a family, hypothetically, can be associated with an infinite set of different ultraviolet (UV for short) scales, which makes it possible to deselect one of them. It turns out, a theory in which there is no distinguished UV scale. And, since the selected UV scale is the main problem of nonlocal QFT, such a hypothetical theory may already be a fundamental theory of nature. At the time of writing this paper, such a study seems to be the subject of the future.
There are still a few questions left to answer in the Introduction section. Why a scalar field? On the one hand, from the theory of groups and Lie algebras' point of view, the quantum theory of a scalar field has the smallest number of symmetries. For this reason, in such a theory there are no Ward-Takahashi identities that can simplify the S-matrix calculation. In this sense, the scalar theory is the most complicated. On the other hand, to construct such a theory, it is the measure theory in infinite-dimensional spaces that is needed. In this sense, the scalar theory is the simplest, no additional symmetry relations need to be taken into account.
Why nonpolynomial self-interaction? In nonpolynomial theories that satisfy some general principles, the S-matrix in terms of the GCPF is a convergent series in the interaction constant g. Let us note that for polynomial theories, for example, for the theories ϕ 4 and ϕ 6 , such a series, being a PT series in the interaction constant g, is asymptotic. An asymptotic series arises as a consequence of the incorrectness of the mathematical transformations performed in the process of its derivation. In particular, the corollary of the monotone convergence theorem (MCT for short) for alternating series is violated. And although the latter is only a sufficient condition, but not necessary and sufficient, nevertheless, the emergence of an asymptotic series means that the permutation of summation and functional integration doesn't really take place. In nonpolynomial theories of a certain class, this problem doesn't arise.
Why fractional power self-interaction |ϕ| α , where α ∈ (1, 2)? The research of such a theory seems natural for a number of reasons. The first reason is that such α is greater than the power in the source term and less than the power in the Gaussian theory. Therefore, in such a theory, it is possible to apply the corollary of the MCT for alternating series. The second, less obvious, reason is that such a theory is implicitly related to all the polynomial theories at once. This will be demonstrated in our paper. And, finally, it can be expected that the method of calculating the S-matrix in such a theory can be applied in the future to a wider class of theories. In this sense, in addition to the research of the theory itself, the result of the paper is a new method.
Finally, why S-matrix and vacuum energy E ? The problem of QFT is considered solved if an explicit expression for the S-matrix is obtained. Here again, it is important to note that the very notion of an "explicit" expression needs to be defined. It is generally assumed that the corresponding functional integral must be calculated. The S-matrix contains complete information about scattering: the variational derivatives of the S-matrix determine the Green functions due to the interaction of particles. The Green functions make it possible to reconstruct the corresponding scattering cross-sections, which are experimentally observable quantities. For this reason, the research of the S-matrix is the most important. Further, the simplest quantity contained in the S-matrix is the vacuum energy E of the theory. It is this value that is most easily calculated. Let us note that the calculation of the functional integral leads to results that depend on the UV scale of the theory. In some cases, this dependency may be undesirably "strong". In this case, one has to invent some kind of subtraction scheme or change the original formulation of the theory. In this paper, we aim to calculate the functional integral in nonlocal QFT for given values of the theory parameters and to research how the results obtained depend on the corresponding parameters. This paper has the following structure. In section 2 we present a summary of the problem and obtained results. This section can be read independently of the rest of the paper. All the necessary definitions and notations are given within the section. A section 3 is devoted to physical background and motivation. A section 4 is devoted to: 1.
explanation of the Gaussian measure concept; 2.
derivation of PT series and proof of its convergence; 3.
exponentiation of the obtained series.
1. exact computation of first and second order in coupling constant g; 2.
expression for the second-degree polynomial approximation.
In the Discussion section, we give a final discussion of all the results obtained in the paper.
In the Conclusions section, we highlight further possible areas of research.

Summary of Problem and Obtained Results
In this paper, Euclidean quantum theory of scalar field ϕ in R d with action S[ϕ]: is considered. Here ϕ : R d → R -real-valued function (scalar field), g : R d → Rreal-valued function (IR coupling constant), j : R d → R -real-valued function (source, appearing in the expression below), L = G −1 -inverse propagator, defined on the space of scalar fields Φ. We don't specify this space in the present section, but in the following sections we will assume that Φ is a separable HS. Convention: hereafter we use the same letter for an operator and its integral kernel, for example, L and L(x, y), and we also use the same physical term for an operator and its integral kernel, for example, the term "propagator" for G and G(x, y). It is always clear from the context what exactly we are talking about. Also in the paper we will assume that L(x, y) = L(y, x), hence G(x, y) = G(y, x). Finally, the propagator G is a trace class operator and in all the final calculations, where possible, we will assume the translational invariance of the Gaussian theory, which implies G(x, y) = G(x − y). Complete (including disconnected parts) Green functions GF Z in terms of functional integral over primary (opposite to composite) field ϕ in such a theory: is derived in terms of convergent perturbation series in powers of coupling constant g, which is a function of x in the general case: In the expression (2) we use the notation "big differential" under the integral sign only formally: It is well known that there is no Lebesgue measure on infinite-dimensional separable HS and BS. However, together with the part of the integrand, the Gaussian exponent, this notation (up to a constant) means an integral over the Gaussian measure γ G with zero mean and covariance operator (propagator) G in the corresponding space. This is how the expression (2) should be understood in the framework of a rigorous mathematical theory. Details are provided in the following sections.
In the expression (3) the following notations are introduced: is the interaction scattering matrix (hereinafter -S ε -matrix) and Z 0 [j] is the Gaussian GF; 2.
G n is n × n matrix with elements (G n ) ab = G(x a , x b ) is the restriction of a propagator on finite-dimensional space R 2n ; 3.
φ(x) = Gj(x) = R d dy G(x, y)j(y) is the classical field (S ε -matrix argument) corresponding to the source j; 4. ε > 0 is the parameter for defining integrands on null sets (in final results ε → +0) and 1 n is n × n identity matrix with elements δ ab (Kronecker delta). Further in the paper we have proved the existence of the connected Green functions GF G I,ε [φ] = ln (Z [j]/Z 0 [j]) (index I means "interaction part"). At the same time, it is convenient to consider this functional depending on φ, but not on j. Let us note that another definition of the connected Green functions GF is also introduced in the literature [44,48]. At the same time, such a functional is considered depending on j. Convention: hereafter we choose a normalization of functional integration such that the value Z 0 [j = 0] = 1. This choice corresponds to the normalization of the Gaussian measure γ G in Φ.

Further, PT series for
, and the nth term of the series has the following form: In the expression (4) the following notations are introduced: 1. G C,n is the set of all connected undirected graphs with no loops and multiple edges on n vertices; 2.
ν ab (Γ) is an adjacency matrix of a graph Γ (the differential operator ∂ s ab is raised to the power of this quantity); 3.
s ab are auxiliary variables, using for extracting connected contributions to G I,ε,n [φ].
Further, we have calculated the first two orders in g(x) = gχ Q (x), where χ Q is the indicator function of d-dimensional cube centered at the origin with Vol Q = V, for the vacuum energy density w vac = E /V, namely (italic Γ is the Gamma function, 2 F 1 is the Gauss hypergeometric function and E = −G I,ε [0] is the vacuum energy): and applied these formulas for the cases of virton propagator [3,8] (exponential form factor, which is a smooth UV-cutoff function in momentum space with UV parameter µ) and Euclidean Klein-Gordon propagator with mass m.
For the Klein-Gordon propagator with mass m and sharp UV-cutoff function in momentum space with UV parameter µ when the modes are frozen at |k| > µ and the propagator in momentum representation G(|k| > µ) = 0, we have obtained for the vacuum energy density w vac the following expression: Besides, we have provided the following polynomial formula for G I,ε,n [φ], derived by Legendre polynomial P n expansion (with an arbitrary coupling constant g): In the expression (8) the following notations are introduced: ) is the n-particle quantum entangler -the main difficulty in calculating any GF; is the auxiliary vector quantity, introduced for the compactness of formulas; 3.
Convention: summation in the last line of the expression (8) is carried out over all l ab ≥ 0 satisfying the conditions n ∑ b=1 l ab = 2i + ν a (where l ab = l ba ) and l aa is even. This expression is non-zero for even n ∑ a=1 ν a and zero for odd.
This yields the simple approximation formula for w vac if one chooses the degree of approximation polynomial is equal to 2 and again g(x) = gχ Q (x): 7 of 63 Beyond that, we derived another approximation formula for w vac , based on HSG approximation of the matrix G n . This approximation is quite often used in statistical physics. Supposing for G n the ansatz: with n × n matrix J n of ones, Heaviside step function θ and parameters δ and γ, we arrive at the following expression for w vac : where v = π d/2 Γ(d/2+1) δ d is the "volume of one hard-sphere particle" and parameters δ and γ are determined by the equations: Let us make one remark about the equations (12). The HSG approximation is used in the statistical physics for potentials that decay rapidly at infinity in x. This is reflected in the equations (12) when integrating over R d . This usually also means that the product of a power function and the power of the propagator is again integrable. All the propagators considered in the paper satisfy the conditions of HSG approximation applicability. The other reasonable equations for parameters δ and γ are also possible, and these are chosen to make formulas of HSG approximation simpler for substituting the particular cases of propagators. All the obtained formulas were checked for the limiting case α = 2 and applied to the virton propagator (exponential form factor, smooth UV-cutoff function) as well as to Euclidean Klein-Gordon propagator with mass m and sharp UV-cutoff function. Finally, we present approximate plots of vacuum energy density w vac for all values of g for both described approximations and provide corresponding approximate formulas. These formulas are valid for strong coupling g as well as for weak, which is quite a rare case for quantum field models.
Using Weierstrass M-test, we proved the convergence of such a series for j = 0 and we have obtained the majorizing series explicitly. As for j ̸ = 0, we found the majorant depending on the regulator ε. There is still an open question about finding a more accurate majorant whose dependence on epsilon is negligible (for the case of non-zero source j). Therefore, we will remove ε in the quantities calculated for j = 0, such as vacuum energy and its density. And in the quantities depending on j we won't remove ε, which will be reflected in their indices, such as G I,ε , etc.
Performed considerations and research will help in the understanding of interacting quantum fields' general properties. Besides, they can potentially be applied for the strong coupling limit of nonlocal ϕ 4 and ϕ 6 theories, which will be the subject of further research and publication of the authors.

Definitions and Notations
The starting point of our research is the complete (including disconnected parts) Green functions GF Z in terms of functional integral over primary field ϕ. Traveling back, the explicit expression for Z is: We will also call this functional the GCPF. Let us note that here and after we follow the normalizing convention Z 0 [j = 0] = 1, which makes this definition of GF coherent with those via Gaussian Measure in the next section. The function (distribution) L(x, y) is the integral kernel of the operator L, which defines the quadratic part of a nonlocal QFT (Gaussian theory) action. Nonlocal QFT can be considered as a regularization of the Euclidean Klein-Gordon theory, but it is also of independent interest. Let us note that the integral kernel G(x, y) of the inverse operator G = L −1 is the UV finite propagator of a nonlocal QFT in any space dimension d. G(x, y) is also called the Green function of Gaussian theory. In all the final calculations, where possible, we will assume the translational invariance of the Gaussian theory, which implies G(x, y) = G(x − y), but at the same time, we will consider coordinate dependent coupling constant g. We will often choose g(x) = gχ Q (x), where χ Q is the indicator function of d-dimensional cube centered at the origin with Vol Q = V. We will also suppose that G(x) ≥ 0 and G(x) ≤ G(0) for all x. The functional integration is understood in the sense of Gaussian measure with covariance operator (propagator) G. In this paper, we assume thath = c = 1.
Generating functional Z is a regular functional, so it could be expanded in a functional Taylor series at j = 0 [44,48,52]: where the repeated index i means integration over some spacetime variable x i . This means the following: Coefficients D (n) 1,...,n = D (n) (x 1 , . . . , x n ) are called n-particle or n-point correlation functions of a functional Z since they represent n-point correlators of the theory: Correlators are the object of interest since they are the quantities which determine the cross-section of scattering of particles and therefore can be directly measured in particle physics as well as in statistical physics. In Gaussian theory one has a usual Gaussian integral, and we will denote the corresponding GF as Z 0 . Also recall the convention we follow: Z 0 [j = 0] = 1. In this case, for GF Z 0 we have a simple expression: It is also known that Z can be represented as an exponent of the other regular functional G, which is called the connected Green functions GF: and from statistical physics point of view Z is a GCPF, as well as G, is a grand thermodynamic potential (up to a temperature factor).
The functional G can also be expanded in a functional Taylor series, and its coefficients are called connected n-particle or n-point Green functions: The terminology goes from Feynman diagrams and will be clarified in the following while considering the Meyer cluster expansion. Of course, all D (n) (x 1 , . . . , x n ) can be expressed via G (n) (x 1 , . . . , x n ) by expanding the exponent. Throughout this paper, we denote connected n-particle Green functions with the mathcal font G (n) 1,...,n rather than the Green function of Gaussian theory G(x − y).
In statistical physics a finite limit G/V at V → ∞ exists. Therefore, we can write G as follows: where f is a volume density of G, which is simply the pressure for the case of homogeneous systems. It can be shown that where E is a vacuum energy of the considering QFT, and it is useful to consider its volume density The quantity w vac is a very useful quantity which in particular can tell whether the system has phase transition or not. This relation is literal. From Lehmann-Symanzik-Zimmermann reduction formula, one can obtain that the expression for the S-matrix n-particle correlation functions in momentum representation reads (the repeated variable means integration): where G kk ′ is the two-particle Green function of Gaussian theory in momentum representation. Of course, G −1 = L, but it is accepted to denote the corresponding operator as G −1 . For deriving the S-matrix, that is exactly the functional which functional Taylor coefficients are S (n) 1 ′ ,...,n ′ , one should simply change variables from sources to so-called "classical fields" (φ instead of ϕ), namely: so the S-matrix satisfies the following relation (up to a certain point, it is not necessary to introduce the regulator ε, mentioned in the section 2): as well as "normalized" GF Z I : and normalized GF G I of connected Green functions (this functional can be considered both depending on the source and depending on the classical field, which was also mentioned in the section 2): where G 0 [j] is the Gaussian theory connected Green functions GF. Such a functional is equal to the quadratic: Next, we are going to describe two approaches to the perturbative calculation of n-particle functions D (n) and G (n) . Let us consider the expansions of GFs Z and G in series: which could be power series in some parameters of the system such as coupling constant g. Then they induce the corresponding expansions of n-particle functions in powers of the same parameter: as well as for the S-matrix. In the next subsection we are going to introduce the conventional method for obtaining of such expansions.

Standard PT and Its Inapplicability
Usually, when considering QFT with GF Z: with analytic interaction potential V(ϕ), one can obtain PT series via transformation: is understood in the sense of power series in operator δ δj . Consequently, one has to demand the analyticity of the potential V(ϕ) in ϕ for this method to work. Let us note that it is this operator that leads to asymptotic series in QFT. Then, after introducing classical field φ, the final formula reads: The should be properly defined since it contains multiple variational derivatives in the same point, but in the final answer all the regularizations of this expression could be removed. Then, expanding both exponents, one can obtain a series in coupling constant g. As already mentioned above, these series are usually asymptotic. As one can see, this way relies significantly on the analyticity of the potential V(ϕ). It is invalid for the case of fractional power self-interaction potential.

Motivation of Study
At first glance, it can appear to be nonphysical to consider models with interaction potential of the form V(ϕ) = g(x)|ϕ| α (explicit dependence on x through g is not indicated in the potential). Though, there are several reasons to study it.
Firstly, it is a kind of exactly solvable model in QFT, where one can produce converging series in g. So it is an object of common interest because it is useful to learn the behavior of different QFT systems to: understand, what kind of properties they can own and, consequently, what one can look up in other models; 2.
develop new methods and approaches, which can be subsequently extended and applied to other systems to gain new results; 3.
since we have the converging series for such a theory, we can hope to construct their analytic continuations to points α > 2, in particular in the vicinity of α = 4, which is an object of common interest.
Secondly, there are reasonable expectations, that in theory with the potential V(ϕ) = g(x)|ϕ| α it is possible to complete explicit transition in obtained formulas from weak to strong coupling regimes. Moreover, there is the conjecture on the following duality between scalar field theories: This duality can be very useful for the research of nonlocal ϕ 4 theory in a strong coupling regime as well as a nonperturbative description of phase transitions. This duality can be obtained with the use of the Parseval-Plancherel identity to functional integral with the following substitution of saddle-point asymptotic in the interaction part, which turns out to be valid in renormalizable cases. Though, for careful checking of this duality one should obtain PT series of asymptotic expansions of ϕ 4 theory from the strong coupling regime of |ϕ| 4 3 , which demands the general studying of |ϕ| α theories for α ∈ (1, 2). Describing, checking and using of this duality is a subject of current studies of the authors, and will be the topic of the nearest publications.
On the one hand, it should be noted that the first of three substitutions in the expression (33) poses a danger: if the original Gaussian measure γ G was given by a trace class covariance operator G, the dual Gaussian measure γ L is determined by an operator L with an infinite trace value. Therefore, the integral over γ L is ill-defined. This is consistent with the conventional point of view that ϕ 4 is quantum trivial in d ≥ 4.
On the other hand, from the Gaussian measure in separable HS point of view, there are no evident reasons for GF to be in any sense "trivial" for any d since x is the index of a vector and the field ϕ is a vector in the corresponding HS Φ. The integration theory in Φ is insensitive to the dimension of the index and is constructed uniformly for any Φ. For the exponential of a polynomial, the integral over a "good" Gaussian measure is defined in rigorous mathematical sense [46,47]. For this reason, quantum triviality may turn out to be a consequence of the incorrectness of mathematical transformations in the calculation.
Be that as it may, the authors would like to find rigorous mathematical proof of the quantum triviality or to get to the conclusion that this triviality is a PT artefact only. And the present paper is the first step to ascertain the truth in this complicated question.

Gaussian Measure
Let us rewrite the expression for the complete Green functions GF Z as an integral over a Gaussian measure γ G with functions g and j in L 2 (R d ) (vector space of Lebesgue square-integrable functions): where we understand functional integral as integral over HS Φ = L 2 (R d ) -linear space of equivalence classes of functions from L 2 (R d ).
By definition, Gaussian measure γ G is a measure in HS Φ defined by some covariance operator (propagator) G : Φ → Φ, which is demanded to be trace-class. Let us introduce {e n } ∞ n=1 -a family of eigenvectors of G with corresponding eigenvalues {λ n } ∞ n=1 , so that Ge n = λ n e n , ∀n ∈ N. We also denote as Φ N := ⟨e i ⟩ N i=1 the linear span of first N eigenvectors, and as P N : Φ → Φ a projection operator on this linear span. Then, by definition, integral of any function f depending only on P N ϕ over Gaussian measure γ G reads: Then, under certain conditions, the integral of the generic function f : Φ → C can be computed using Dominated Convergence Theorem (DCT for short) and convergence of the All the described results on Gaussian measure can be found with derivation in [46,47].

Physical and Mathematical Approaches to Measure Theory in Infinite-Dimensional Spaces of Functions
It is useful to compare the mathematical view on functional integration with the physical one. In physics, we write the expression (30). The thing is that there is no translation invariant countably-additive Lebesgue measure on infinite-dimensional separable BS or HS. So if one wants to develop a measure theory in separable HS, he has to reject either countable additivity or translational invariance. So in mathematics, there are several approaches to integration in functional spaces, which are:
Ito integral (Wiener measure, which corresponds to significant restriction or sigmaalgebra to the one generated by cylindrical sets only) [63].
In this paper, as has been previously explained we use Gaussian measures language. In the way the physicists understand the functional measure, one can write informally: It means that informally Gaussian measure determined by some propagator G is a "physical" measure multiplied by an exponential of minus quadratic part of the action with "kinetic operator" G = L −1 .

About Continuity of Functionals
Let us make one important remark that will further ensure the correctness of all subsequent transformations. Consider the following functional (a function defined on a function): so it is positive, measurable and f [ϕ] ≤ e j(x)ϕ(x)dx ≤ e ∥j∥ 2 ∥ϕ∥ 2 which is clearly integrable with a finite value of the integral. Moreover, for proving that the sequence as N → ∞ almost everywhere it is necessary and sufficient to show that R d dx g(x)|ϕ(x)| α is a continuous functional on Φ, since exponent and R d dx j(x)ϕ(x) are continuous. Let the coupling constant (non-negative function) g be bounded and Lebesgue integrable. Since α ∈ (1, 2), the following is true (we denote the measure for which the coupling constant is the measure density by the same letter g): Thus the functional R d dx g(x)|ϕ(x)| α is continuous for all ϕ ∈ Φ. Here we have used a well-known consequence of Hölder inequality for the set R d with finite measure g(R d ). We would like to underline that this functional is not continuous when α > 2 because in this case ∥ϕ∥ g,α is not finite if ∥ϕ∥ 2 is finite (which is equivalent to ϕ ∈ Φ). As a result, one can apply DCT and find that: That's quite a good outlook that if one considers a functional integral from a Gaussian measure point of view, one can't compute a functional integral with α > 2 (e.g. α = 4, which is currently the most interesting in the physical case) by simply substituting P N ϕ instead of ϕ and taking a limit N → ∞. At least DCT is not applicable in this case. However, it should be remembered that DCT gives only sufficient conditions. Perhaps some other theorem may be useful in the case α > 2.

Reduction of Functional Integral to the Series from Finite-Dimensional Ones (GCPF Expansion)
In this section, we are going to derive PT type series, which will converge. Since the power of potential is not a natural number, the traditional method of carrying out an interaction action from integral in the form of a differential operator is not suitable. So, we follow the other way and in some sense get all Feynman diagrams summed into one integral, which appears to be exactly the coefficient of g(x 1 ) . . . g(x n ) in the obtained series.
The described approach of constructing PT repeats and develops methods described in [1- 3,21]. Let us proceed to the computation itself. Since the potential U(ϕ) = |ϕ| α is even function, the following is true: where two points are important: 1. we understand Fourier transform F (internal transform is performed over the variable ϕ, but not x) in the sense of distributions from S ′ (R) -the space of linear continuous functionals on the Schwartz space S(R); 2. this distributional Fourier transform is consistent with the following formula for the integrable function f (ϕ): Unfortunately, the potential U(ϕ) = |ϕ| α have Fourier transform only in generalized sense in S ′ (R). And as usual, careful computations with generalized Fourier transform are quite subtle, so we would like to avoid them. For that purpose, we approximate U(ϕ) with some smooth and integrable function U Λ (ϕ) from the Schwartz space S(R), so that U Λ (ϕ) → |ϕ| α pointwise, when Λ → ∞. Then both U Λ (ϕ) and F [U Λ (ϕ)] are usual functions in S(R) rather than distributions, and we are entitled to calculate their Fourier transform via the usual integral formula (42). We will make all the transitions for smooth and integrable function U Λ (ϕ). Also, we will assume that 0 ≤ U Λ (ϕ) ≤ |ϕ| α for all ϕ and that all the functions U Λ (ϕ) are even in ϕ, which will be useful in calculation of PT series majorant. Namely, we can write: since 0 ≤ U Λ [ϕ] ≤ |ϕ| α and the left-hand side integral converges absolutely, so DCT justifies this transition. After that, we are able to write: and this formula is much more convenient for the calculations. In the following calculations, we will extract this limit from functional integral also due to DCT and will consider regularized functional Z Λ [j]. And at the end of our computation, we will calculate the limit Λ → ∞. So we start by considering the "regularized" functional: and due to DCT described in previous section Further, we write the expression (45) in terms of the Fourier transform of potential and use that F [U Λ (ϕ)] is in S(R): Now let us expand the external exponent in a Taylor series in powers of the coupling constant g: Exchanging summation and integration, we arrive at the following expression: The correctness of this permutation will be proven independently in the next section. Finally, we have to integrate over ϕ. After introducing the "modified" sourcej (the second term in the right-hand side contains the Dirac delta function δ): we obtain usual Gaussian integral with linear exponent: Using the expression (49), we arrive at the following equality: Now we are going to use Parseval-Plancherel identity, but disappointingly the matrix (G n ) ab = G(x a , x b ) (the restriction of a propagator on finite-dimensional space R 2n ) is only semi-positive, which will be proved in section 4.2.1. However, det(G n ) ab = 0 only on null sets (sets of zero measure), which made it possible not to take care of it up to this moment. But in the following, we are going to approximate integrand with polynomials so we would like it to be bounded. So we introduce the Z Λ,ε instead of Z Λ , adding −ε n ∑ a=1 t 2 a term into the series nth term exponent for all n. And then we will calculate the limit ε → +0 in a final result, using the obtained majorant and its (in)dependence of ε. Hence, we have: At this point, it is evident that every term in the above series converges to the term with the same number in the series for Z Λ when ε → +0. At present, we can use Parseval-Plancherel identity: In the expression (52) the following notations are introduced: is "effective" discrete source, since it was obtained from j firstly by acting the "continuous" operator G, and then by its "discrete" inverse.
Let us note that both R n and G n are symmetric, hence diagonalizable in an orthonormal basis by orthogonal transformation. Moreover, all G n are semi-positive, so G n + ε1 n is positive and hence invertible. This proves that all written integrals converge absolutely, so we can interchange the order of integration in t a and x a variables.
In the section 4.3, we will show that the obtained series converges uniformly in Λ for generic j and also in ε for j = 0. The idea of the proof is that we can choose U Λ (ϕ) such, that 0 ≤ U Λ [ϕ] ≤ |ϕ| α . After that, we can apply Weierstrass M-test to the series for Z Λ,ε if we prove that the series above with U Λ (ϕ), replaced by |ϕ| α , converges. We will do it later in the paper.
Up to this moment, we finish on the expression for the "normalized" GF Z I,ε , removing the regulator Λ: , in accordance with the notations, introduced in the section 3. In particular, we can write the value Z I [0] for the case of j = 0, also removing the regulator ε: Here the matrix G n is degenerate on null sets, but it doesn't affect the integral because of the existence of majorant not depending on ε. This expansion of GF in perturbation series looks similar to the expansion of the GCPF of non-ideal gas with potential G(x a , x b ), but with additional inner integrals over ϕ a . Keeping this in mind, we will refer to Z I,ε,n as the n-particle canonical partition function. These inner integrals are in fact the key difference between statistical physics and QFT, warranting the complication of the last one.

Properties of G n -Matrices
It is useful to visualize the typical structure of matrix G n . We will do this for the translation invariant case (the general case is done in a similar way), since this is the case that will be considered in all the final results. The desired expression is: There are two limiting cases for this matrix: 1.
When all x a are equal, then:

2.
When all x a are infinitely faraway, then: Let us prove the semi-positivity of G n . Since the operator G : Φ → Φ is positive, the equality ⟨ j |G| j ⟩ ≥ 0 holds ∀j ∈ Φ. And for distributions this inequality also holds from the reasons of continuity. Namely, choosing a smooth sequence of functions approximating some distribution of the form: we get exactly: which completes the proof of semi-positivity of G n for all values of x a . For this reason, G n + ε1 n is a positive definite matrix for all n as well as its inverse, for any ε > 0.
The other way to prove this fact (in the translation invariant case) is to use Bochner's theorem, which claims the semi-positivity of G n from the fact that propagator G(x − y) is from Schwartz space and its Fourier transform F [G](k) ≥ 0 is non-negative. Though, this proof needs one more additional assumption that G(x − y) ∈ S(R d ).
By construction, the matrix G n is symmetric, therefore, diagonalizable. Let us denote its eigenvalues as 0 ≤ λ In particular, we found the bound for maximum eigenvalue of G n (and, simultaneously, the minimal eigenvalue of G −1 n , when G n is invertible): This bound will be used in the future computation of majorizing series.

Coordinate-Free Part
We always can exchange summation and integration when the series terms are positive due to MCT. And we have to check the absolute convergence of the obtained series to use DCT. So, it's sufficient to take the absolute value of every term, exchange sum and integral and prove that this series converges. At this moment, we have to find the upper bound of the series terms: where the notation of J n (α) is introducing for the convenience. Now we are going to estimate |ϕ 1 | α . . . |ϕ n | α with restriction ∥ϕ∥ 2 = r 2 . Our aim is to bound the integrand on every sphere in R n centered at the origin, and then use spherical coordinates. We start with the Lagrangian function (λ is the Lagrange multiplier): ln ϕ a − λ ∥ϕ∥ 2 − r 2 , and find its unconditional extremum: Taking the sum over a to get λ we receive: So, finally, the maximum is reached at: So, one can estimate (and these bounds are strict): In addition to the previous bound, we can change variables by the orthogonal transfor- where the factor det(G n + ε1 n ) cancels out with the Jacobian after changing of variables in the integral. Here we also use the fact that G n is almost everywhere invertible. Hereafter we use the notation (A n ) β ab for the matrix element of a power of matrix, namely: for n × n matrix A n and rational number β. Now we have to transform the part with source: n ∑ a,b=1 We denote for the shortness (which will not get confused with the source due to index 0): where the last transition follows from the Euclidean norm definition. We underline that this quantity depends on x.
In this paper we will mainly consider the calculation of vacuum energy, which corresponds to the case j = 0. Hence, we won't think about careful estimations for the case of nonzero source, and write the most rough estimation. Namely: where we used the fact that G n ≥ 0, and hence maximal eigenvalue of R n is no more than 1/ε. Here ∥ϕ(x a )∥ is a Euclidean norm of the vector (ϕ(x a )) n a=1 , which is bounded with n sup x∈R d |ϕ(x)|. Then we have: Next, we recall the inequality ∥G n + ε1 n ∥ ≤ n(G(0) + ε). The integral in the right hand side of the expression (60) can be easily calculated via series expansion in j 0 : As a first step for bounding the obtained series, we use log-convexity of Gamma function Γ for m ≥ 1: and rewrite: Moreover, for m → ∞ the equivalence relation holds: and it could be numerically checked, that for m ≥ 0: so the following upper bound holds: The last series can be calculated exactly in terms of the error function erf: so we finish at: Let us write the result: Finding the asymptotic of the expression (61) right hand side, we get to: where C n > 0 is a dimensionless constant which grows no faster than exponentially. As a result, we see that our series converges for α < 2. For α = 2 it also converges, but it could be proven in a much more simple way. Namely, it is easy to calculate GF Z for α = 2 exactly, which will be done further in the paper, and the corresponding expansion will converge. For α = 2 we will consider only the case g(x) = gχ Q (x). As the result, we have some at least asymptotic expansion in powers of g. So, since even asymptotic expansions in the predetermined system of functions (power functions of the coupling constant {g n } ∞ n=0 in our case) are unique, then these two expansions must coincide. This proves that for α = 2 our perturbation series expansion also converges. So we have to deal with the coordinate integrals to finish the proof.

Coordinate Part
The next aim is to bound the result of coordinate integration. Namely, remembering the notation g(R d ) = R d dx g(x), we have: Recalling the upper bound for j 0 , we obtain finally: so the obtained series for the GF Z I,ε converges for j = 0. As for j ̸ = 0, the convergence of series requires additional research and calculation of a more accurate majorant, since the obtained one increases infinitely, when ε → +0, therefore, it can be argued that the series converges only for ε ̸ = 0. Though, as it was mentioned above, we mainly restrict ourselves to the case j = 0.

Exponentiation of Series Using Meyer Cluster Expansion
We have just obtained the converging PT series for the GF Z I,ε . In terms of asymptotic series, there is often proved "the exponentiation of connected diagrams", which means roughly that the GF Z I,ε is an exponent of other regular (in the source j) functional. We are going to establish the same result in our case. It will be useful since after such exponentiation (for the coupling constant g(x) = gχ Q (x)) we will obtain only the first power of a volume V = Vol Q rather than all natural powers (in general, the first power of g(R d )). This fact will simplify significantly the extraction of physically measurable quantities, which will be done in the following sections of the paper.

Formulation and Applicability
We start with the following definition. Let (X, Σ, µ) be a measure space, that is a triple: X is a set, for example, R N with some natural number N, Σ is a σ-algebra on X and µ is a complex-valued measure (which shouldn't be confused with the UV-cutoff parameter µ) on a measurable space (X, Σ). Complex-valued measures are also called charges. We denote as |µ| the variation of µ (which shouldn't be confused with the total variation -the value of |µ| on X). In the case, where µ(dy) = g(y)ν(dy) and ν is a non-negative measure, |µ|(dy) = |g(y)|ν(dy). Given a complex-valued measurable symmetric function ζ on X × X, we introduce an abstract GCPF as follows: The term n = 0 in the sum is defined by 1. In the case of classical gas with interaction, n is the number of particles. We denote by G n the set of all (undirected, no loops) graphs with n vertices, and G C,n ⊂ G n the set of all connected graphs with n vertices. Next, let us introduce the following combinatorial function on finite sequences (x 1 , . . . , x n ): The product here is taken over edges of Γ. A sequence (x 1 , . . . , x n ) is a cluster if the graph with n vertices and two vertices i and j whenever ζ x i , x j ̸ = 0 are connected. The cluster expansion allows to express the logarithm of an abstract partition function as a sum over clusters. Namely, the following theorem holds [59].
Theorem 1 (Cluster expansion). Assume that |1 + ζ(u, y)| ≤ 1 for all u, y ∈ X, and that there is exist a non-negative function a on X such that for all u ∈ X, and X e a(x) |µ|(dx) < ∞. Then the following is true: Combined sum and integrals converge absolutely. Furthermore, we have for all y 1 ∈ X: Unfortunately, we cannot apply this theorem to the terms Z I,ε,n of the PT series for GF Z I,ε , since there appears det G n , spoiling the form of the terms in the GCPF expansion. So, we apply it to the terms Z I,Λ,ε,n of the PT series for in the form before using Parseval-Plancherel identity. Hence, we start from (G ε = G(0) + ε for short): Let us note that, for clarity from statistical physics point of view, we consider the case of translation invariant propagator G(x a − x b ), but all the formulas presented below can be easily transferred to the general case of propagator G(x a , x b ).
In our case the point y = (x, t), the function 1 + ζ(y a , y b ) = e − 1 2 G(x a −x b )t a t b , and the measure µ, appearing in the general theory, has the form: To use the theorem one have to find the function a from its formulation. But in the following we will do all the manipulations directly and see that in our case all the transitions are valid without explicit presenting such a function.

Rewriting Series in Terms of Exponent
Henceforth, we will keep in mind the measure definition (68). Hereby, we start from the expansion of GF Z I,Λ,ε : Using the definition of the function ζ, we rewrite the expression (69) as follows: Further, the sum over all possible graphs Γ in the right hand side of the expression (70) can be represented in the following way: In the expression (71) we have decomposed the graph Γ into the connected parts (graphs) (Γ 1 , . . . , Γ k ). Each Γ i ∈ G C,n is a connected graph with a set of vertices V i with the cardinality n i . All the sets V i form a partition of the set {1, . . . , n} : There are k! such sequences for each Γ, since the order of the sets Γ i doesn't matter. The sum over Γ can thus be realized by first summing over k, then over the partitions V 1 , . . . , V k , and then over connected graphs on the sets V i . After substituting this decomposition into the GF Z I,Λ,ε we can sum over the cardinalities n 1 , . . . , n k ≥ 1 with condition n 1 + · · · + n k = n. And the number of partitions of n elements into k subsets with n 1 , . . . , n k elements is given by the multinomial coefficient n! n 1 !...n k ! , which therefore should be inserted. As a result, the following chain of equalities for GF Z I,Λ,ε is valid: In accordance with the notations, introduced in the section 3, let us denote by G I,Λ,ε = ln Z I,Λ,ε the "normalized" GF of connected Green functions with regulations Λ and ε. Applying cluster expansion we made above, we have the following formula for it: To go further, recall the definition of the adjacency matrix ν ab (Γ) for a given graph Γ: We note that ν aa = 0, which is equivalent to our already existing requisition for the graph to have no loops. Using this definition, as well as the definition of the function ζ, we can write and transform the following product: We simply have presented the differences as integrals of derivatives (s ab are the auxiliary variables), using the Fundamental Theorem of Calculus. We have also expanded the sum in the exponent to all graph edges and added ν ab (Γ) multipliers that do not affect the terms for presenting edges and annihilate the contributions of the terms corresponding to the absent ones. Further, we understand the powers of differential operators ∂ s ab as follows: This means exactly that the edges (a, b), that are absent in the graph Γ, contribute only by 1, since for them ν ab (Γ) = 0. And as for the presenting edges, they give the desired difference after the action of the integral of the derivative. Such rewriting of the products provides an unification of the expressions and will be extremely useful for obtaining compact forms of the results that follow in the paper. Finally, in accordance with section 2, we recall the convenient notation: This quantity depends on ε and s ab , but we won't specify it in the notations for shortness.
Traveling back to the GF G I,Λ,ε in terms of the PT series G I,Λ,ε = ∞ ∑ n=1 G I,Λ,ε,n , we find for the n-particle contribution, which is the contribution of all connected graphs with n vertexes, the following expression: We underline that in the PT series for the connected Green functions GF G I,Λ,ε the index n starts from one rather than zero, as in the PT series for the complete Green functions GF Z I,Λ,ε in (69). It is not a mistake, and it arrives from the very statement (66) of cluster expansion theorem.
After substituting the measure µ(dx a dt a ) definition (68), we arrive at the following expression for the n-particle contribution: and, after applying Parseval-Plancherel identity to integrals over t a , the expression for the n-particle contribution becomes: This formula is not a very convenient one since s ab variables are included in a very complicated way because of det(G n,Γ ) and G −1 n,Γ . Though, this expression is a very important one. Indeed, it explains that it is sufficient to consider the terms Z I,Λ,ε,n in the initial (unexponentiated) GF Z I,Λ,ε , and then, in order to obtain G I,Λ,ε,n , it is enough to apply the operator O (everything that the operator depends on is omitted in the notation): and we will use this notation in the following. Summarising, we can get the results for the connected Green functions GF G I,Λ,ε from the ones for Z I,Λ,ε with the following steps: 1.
changing (G n ) ab → (G n,Γ ) ab in the quadratic part of the action; 2.
posterior applying the operator O. Equivalently, on the language of formulas: We additionally declared in the arguments of GFs the matrices G n and G n,Γ we use in the obtained perturbative expansions to avoid confusions. The notation Z I,Λ,ε,n [φ, (G n,Γ ) ab ] means that this functional must be expressed in terms of ϕ, which is fixed, and the replacement of (G n,Γ ) ab is made only in explicit dependence. Let us also notice the form of G I,Λ,ε,n without introducing new variables s ab : The expression (82) will be useful for the following HSG approximation.

Final Form of the Connected Green Functions GF
For the convenience, in accordance with section 2, we recall the notation Q n,Γ for the n-particle quantum entangler: and also we introduce the notation for the linear functional in vector ϕ a : Basically, χ a are the sources in the coordinate representation, which are acted by the operator G, after which we "pull them back" by the discrete restriction of G "modulated" as in the expression (77) inverse. Now we are going to substitute the introduced notations in the expression (79) as well as remove regulator Λ from it. We are eligible for the last action since the majorant for Z I,Λ,ε doesn't depend on Λ for generic j as well as on ε for j = 0. Hence, we have: for an arbitrary source j. Here, exactly as in the (53), we introduced the notation: And for j = 0 we can remove all the regulators and write: These formulas will be used intensively in the section 5.

A Short Way to Obtain Coefficients of Exponentiation
At the end of the section, consider again the coupling constant g(x) = gχ Q (x). If the expression Z (V, G, g) = e V· f (G,g) is proved, then the following equality is true: So, the thing is that we can extract f (G, g) as the coefficient in V if we get some expansion of Z (V, G, g) in powers of V. This follows from the uniqueness of (asymptotic) expansion in a prescribed system of functions (power functions of the volume {V k } ∞ k=0 in our case).

Calculation of PT Series Terms
Let us start by noting that the integrals of the form: are the particular cases of so-called Gelfand hypergeometric functions [56]. Unfortunately, they are too generic and their properties are excessively complicated. However, if we are going to get some practical results it is pointless to express anything via them. Besides, these integrals are enough sophisticated to be expressed simply in terms of even generalized Gauss and Lauricella hypergeometric functions.

Off-Diagonal Terms Expansion
The first and most evident way is to expand the quadratic exponent in its off-diagonal terms, and then compute the obtained integrals directly. However, this approach also fails because of the excessive complexity of the obtained terms. Let us demonstrate it. We start with the expansion for zero source: Here we have already removed the regulators Λ and ε in accordance with (54). Then, substituting the series: and calculating integrals over ϕ a , we arrive at the following expression: There are two main obstacles for applying this formula: 1. it is necessary to calculate the inverse of n × n matrix, depending on x a ; 2.
it is necessary to integrate the elements of the inverse matrix over x a .
These are the main reasons we will develop other computational techniques for the practical calculation of some concrete quantities.

Approximation of Generic Term via Polynomials
In this section we will approximate prefactor |ϕ 1 . . . ϕ n | α as a function f (ζ) := |ζ| α for ζ = ϕ 1 . . . ϕ n with polynomials. There is a couple of constructive ways, and we compare several of them. We will consider three types of polynomial approximations:
However, we will present the explicit formulas only for the Legendre polynomials, since they will be the most suitable for solving the problems, formulated in our paper. To show this, we will consider other families of polynomials. However, they will not be required to display our final results, so they will not be explicitly given. The complete guide about all the mentioned families of polynomials can be found in [55,57].

Background on Constructive Approximation with Legendre Polynomials
Legendre polynomials are defined as an orthogonal system with respect to the integral scalar product generated with the Lebesgue measure on the closed interval [−1, 1]. That is, P n is a polynomial of degree n, such that (we denoted the independent variable as t, since the domains of t and ζ do not coincide): The standardization P n (1) = 1 fixes the normalization of the Legendre polynomials with respect to the L 2 -norm on the closed interval [−1, 1]. There is an explicit formula for Legendre polynomial of degree n: where ( n k ) is the binomial coefficient, and they obey: The feature of the formula (93) is that due to the second binomial coefficient, all the terms for k with the different residue modulo 2 are zero. So, for even n polynomial P n consists of even monomials, and for odd n -of odd monomials. So, about a half of the terms in the sum in (93) are zero, and one should bear this in mind during practical calculations. The first few even Legendre polynomials according to this definition are: Further, for the absolutely continuous function f with derivative f ′ of bounded variation on [−1, 1], the expansion in Legendre polynomials converges uniformly and absolutely: due to Jackson's theorem [55]. Here, the coefficients c n are determined by the formulas: There also exists a bound for the residue of the series truncated on Nth term, but it is not very useful for our case. It is so because it is an estimate for the absolute value of error, so the consequent estimation of GF will be excessively rough, since we have to take the sign into account. As it will be explained in the next subsection, we will approximate a function f (t) = |t| α , for α ∈ (1, 2) and t ∈ [−1, 1], choosing some "normalized" variable t. For the considering case the conditions of Jackson's theorem are satisfied for α ∈ (1, 2), so we have the uniform and absolute convergence of Legendre series to f (t). We write for the finiteterms approximation, using that f is even (for the convenience of writing some further formulas, we denote the coefficients of the Legendre polynomial as u i ), the following expression: Let us note that because of the absence of the odd Legendre polynomials in the approximation, we will use the numeration for the approximation coefficients as in (98) rather than in (97). And ultimately, to avoid any misconceptions, we will write for the coefficients with the "new" numeration c q ( f ) rather than c q . Let us also note that the total degree of polynomial approximation h N is deg h N = 2N. The coefficients c q ( f ) for f (t) = |t| α have the form: and this formula is fair also for q = 0, as a computation shows. So for h N we get: In the following subsections we will use only Legendre polynomials. Therefore we can expand f in a series, putting N → ∞. However, for the comparison of the Legendre approximation with the Bernstein one (which is only approximation rather than expansion), we will keep N finite until the end of the different polynomial approximations analysis.

Construction of Approximation for the Terms
According to theorems on polynomial approximation, e.g. the Stone-Weierstrass theorem, it is possible to approximate the function by polynomials on compacts only. Moreover, if we try to approximate (even pointwise) a function |ζ| α for all ζ ∈ R via any kind of polynomials, e.g. Hermite or Laguerre polynomials (multiplied by the corresponding weight functions so that such expansions converge), we obtain a diverging series after term-by-term integration. And the thing is that there is a problem in the radial (in terms of spherical coordinates in R n ) direction. As a possible solution, before approximating by polynomials, one should evaluate the integral (88) in the radial direction, introducing spherical coordinates in R n (hyperspherical coordinates, if n > 3, but for brevity we will use this term for n > 1). In the remaining integral over hypersphere we can approximate |t| α , where t -some convenient function on the hypersphere with values in [−1, 1], by polynomials in t. There will be no such a problem since on a hypersphere t ∈ [−1, 1], so all the approximations and expansions will converge. After such an approximation one can go back to R n where it is easier to calculate the integrals using the source trick. It is a brief summary of the material in the current subsection.
Let us denote the integral, which we are going to work with, as (n > 1): The integrand has the symmetry under ϕ → −ϕ, since Q n,Γ and |ϕ a | α are even functions, hence (in terms of the hyperbolic cosine): In this integral we can introduce hyperspherical coordinates and proceed as: where we choose ϕ a = rξ a and ⃗ ξ ∈ S n−1 is the unit normal vector. Integration over r gives: In the expression (104) we follow the commonly accepted notation of 1 F 1 for the confluent hypergeometric function. This function is defined in terms of the following series (and this series converges for any finite value of z and thus defines an entire function of z): where (a) n is a rising factorial. Then we expand 1 F 1 and obtain the following expression: (1/2) 2k Γ n(α + 1) 2 One can easily notice that this expansion can be obtained by expanding cosh in the formula (103). As in the section 4.3 on majorant calculation, we can bound the following relation, using (56) and (58): for all x a . And now we can specify the function t: Having defined an appropriate function t, consider finite-degree even polynomial approximation (98) of the function f . Let us note the following: the expression (98) describes a general polynomial approximation such that the polynomials h N are defined for all integer N > 0 and deg h N = 2N. and we will often use the same notations, as for Legendre polynomials, but for generality and simplification of formulas, we won't specify the concrete form of (real) coefficients u i and c q . So, here {P 2q (t)} ∞ q=0 is some general set of polynomials suitable for pointwise approximation of the function f , in other words, h N (t) → f (t) for all t ∈ [−1, 1]. We will also suppose for generality, that the coefficients c q ( f ) can also depend on N, which is actual for approximation with Bernstein polynomials.
In the following we will denote as I n,N the integrals obtained from I n by substitution polynomial approximation h N (t) instead of f (t) in the integrand. Let us note that I n,N → I n due to the DCT, when N → ∞. This fact follows from the uniform convergence of Legendre series to the considered function (as discussed in subsection 5.2), which means the uniform boundedness of partial sums. This convergence will also be checked numerically in the subsection 5.2.4.
Thus we arrive at the following expression: Recalling the identity for n, s > 0: we can return to R n due to homogeneity: It is worth noting that without the ratio Γ n(α+1) 2 /Γ n 2 + ni + k we would receive cosh after summation instead of confluent hypergeometric function 1 F 1 . At the same time, we don't write out the last one, since it is convenient to continue the transformations of the series itself.
We can rewrite the integrand, using the multinomial formula (all the β a ≥ 0): for some integer constants m a ≥ 0. There is a common method for their calculation consisting in introducing auxiliary variables η a , called "sources", in similar fashion to path integrals, and further differentiation over them. Namely, we use the identity, following from DCT: where the partial derivatives ∂ i := ∂ ∂η i . If all η a = 0, we get the desired equality. Keeping the previous expression in mind, we obtain: (1/2) 2k Starting from this formula and below, we found its useful to introduce the notation G I,ε,n [j] N . We specify the value N to the right of j, so as not to clutter up the number of indices to the left of j. By definition, G I,ε,n [j] N is obtained from G I,ε,n [j] by substituting the approximation I n,N instead of I n . And owing to the convergence I n,N → I n , when N → ∞, it is also true that G I,ε,n [j] N → G I,ε,n [j]. Recall that in the expression (111) s ab ν ab and ε enter in three ways: 1. Q n,Γ (x, φ) in the exponent; 2.
(G n,Γ ) ab in the exponent.
It is interesting that in fact the formula (111) provides the expansion of the connected Green functions GF for fractional-power interaction theory in terms of contributions of power-interaction theories with some "weights". However, it is not the same if one would like to expand |ϕ| α in the action (1) itself. In such an approach there would be no "damping factor" of the form Γ n(α+1) 2 /Γ n 2 + ni + k , and the series would be divergent due to its absence. And the reason for this factor to appear is that before the approximation we reduced the integral from non-compact R n to a compact S n−1 . As a result, we won't lose the convergence on any step and therefore get the convergent series rather than asymptotic. And this achievement was the primary aim of the present paper.

Combinatorial Expansion
We want to derive analytical formulas for the quantities: for any integer m a ≥ 0, for a = 1, . . . , n. The authors are sure that such formulas were derived a long time ago, but we could not find them in the literature in a suitable form. So let us derive them from the basics for our own usage. In fact, the formulas we will obtain are nothing else but the general formulas for the symmetry coefficients (up to some factor) of Feynman diagrams.
As already mentioned about the expression (109), in this expression it is not hard to recognize the moments of n-dimensional Gaussian distribution with the covariance matrix G n,Γ , which we will denote by brackets of "averaging": The index of the angle brackets indicates the covariance matrix we use for the computation of this quantity. We will refer to the quantity (113) as correlator or multidimensional Gaussian moment with integer powers. We have already introduced another notion of correlator in the section 3 (referring to the path integral), but it will be clear from the context what notion we mean in every particular case. We will also refer to the variables ϕ a in correlators as fields. For now, we will consider covariance matrix G n rather than G n,Γ for the simplification of notations (in this subsection we will also omit the index n for the same reason). Though, we won't suppose any of its special properties except it is symmetric and has equal terms on the diagonal which we will denote as G aa = G(0) (a generalization without this property is straightforward). Then the result for G n,Γ can be obtained with the simple substitution G ab → (G n,Γ ) ab .
There is the Isserlis-Wick theorem which claims that the multidimensional Gaussian moment with integer powers can be expressed in terms of the covariance matrix elements as the sum over all possible pairings: We will denote the general element of P m as { {a 1 , a 2 }, . . . {a m−1 , a m } }. And these indices are exactly ones from the formula (114). So one can determine the set of pairings with the array {a k } m k=1 , and we will use this array to determine the indices of i a . And now we are going to calculate the number of all equal terms in the sum over all pairings in the expression (114) assuming, that i 1 , . . . , i m 1 = 1; i m 1 +1 , . . . , i m 2 = 2; . . . and i m n−1 +1 , . . . , i m n = n, which corresponds to our object of interest (113).
We start from the case of n = 2, so the index a takes only values 1 and 2. Then there are only two different types of terms G ij : 1.
when i ̸ = j, then G ij = G 12 . As a result, every term in (114) has the form G(0) q (G 12 ) p for some integers p, q ≥ 0. So we have to count the number of pairings {a k } m k=1 , corresponding to every possible p and q. Suppose we consider ϕ m 1 1 ϕ m 2 2 G for m 1 > 0 and m 2 > 0. From the conservation of the total points number one can deduce that m = m 1 + m 2 = 2q + 2p. So, for given m 1 and m 2 the number q is uniquely determined with p.
We receive the combinatorial factor from the following considerations: there are m 1 fields with index 1 and m 2 fields with index 2. It is possible to visualize them with the Figure 1. We will denote all the enterings of the fields with the points, and for every pairing we will draw the segment ending in the points constituting a pair. We will also refer to these segments as edges, being inspired by the graph theory. In the following we will mean under the notion "configuration of pairings" the set of all pairings which give the equal contributions in the sum (114). In fact, this is a set of pairings (subset of P m ), consisting of { {a 1 , a 2 }, . . . {a m−1 , a m } } for some a k which possess the same number of pairings between different clusters of points and inside every cluster. For n = 2 it is equivalent that they have the same p, but the number of "degrees of freedom" n(n − 1)/2 rises for greater n. One can draw an analogy with statistical physics: we identify microstate with the particular set { {a 1 , a 2 }, . . . {a m−1 , a m } } and macrostate is what we call the configuration of pairs.
Then we have m 1 + m 2 points in total, and they all have to be connected by such segments in pairs. And these points can be naturally distributed between two clusters, corresponding to the indices of the fields. For the convenience, we also introduce the notations l a := m a − p for the number of points which will be paired with the points from the same cluster. From the Isserlis-Wick theorem (114) it follows, in particular, that for odd m 1 + m 2 the considering correlator is zero, so we will consider only even case. So, the term G(0) q (G 12 ) p corresponds to p segments with the ends in different clusters, and q segments, which have ends in the same cluster. To draw all the segments which such condition one have to: choose p points from every cluster which will give rise to the segments with ends in the different clusters; 2.
choose the way we draw the lines between the chosen p points in every cluster; 3.
choose the way the remaining m 1 − p and m 2 − p points in first and second clusters correspondingly will be connected with the segments inside their clusters.
The first number is ( m 1 p )( m 2 p ) by combinatorial definition, according to the combinatorial rule of product, where ( n k ) = n! k!(n−k)! is a number of combinations of k elements from a subset of cardinality n. The second one is p!, since for the first point (for any fixed enumeration) in the first cluster we have p variants of choosing a pair, for the second one -(p − 1), and so on. And the third number equals to l a ! 2 la /2 (l a /2)! for each of two clusters, so in total these two contributions should be multiplied. The explanation for the last formula is that we have to choose firstly one pair, which is possible in l a (l a − 1)/2 ways, then the another one, for which we have (l a − 2)(l a − 3)/2 ways, and so on. Finally, we have to divide the product of such ways' numbers by the permutation number of points' pairs, which is (l a /2)!. We also note that the numbers l a have to be even for any possible configuration, and the number of pairs inside the cluster is l a /2.
In total for the number M p of pairings with contribution G(0) q (G 12 ) p , using the same combinatorial rule of product, we have the following expression: We have specially written this formula in such notations to make its generalisation for n > 2 more clear. At this moment, for the two-dimensional Gaussian moment with integer powers m 1 and m 2 the expression reads: where the summation is carried out over all p, such that l a = m a − p are even, since only these configurations of pairings are possible. For the considering case, when m 1 + m 2 is even, the two cases are realisable: m 1 and m 2 are both even or both odd. In the first case the summation condition means that the summation is carried out over only even p, and in the second -over only odd p. The sense of these facts becomes clear from Figure 1. Let us generalize the resulting formulas to the n-dimensional case. Acting in the same spirit, one can verify the validity of the following expression: (117) In accordance with the notations, introduced in the section 2, the summation in the expression (117) is carried out over all l ab ≥ 0, satisfying the conditions n ∑ b=1 l ab = m a (where l ab = l ba ) and l aa is even. This expression is non-zero for even n ∑ a=1 m a and zero for odd. The formula (117) truly generalises (116) and we will use it in our further calculations. Let us derive the expression (117). It can be done in the same manner as for n = 2, using the same graphical interpretation. We start from introducing slightly different notations. In general case we have n clusters, consisting of m 1 , . . . , m n points, correspondingly. At present, the contribution to the sum in (114) of a given pairings configuration is prescribed with the numbers l ab (indices a, b ∈ {1, . . . , n}): 1.
the numbers l ab (a ̸ = b) of edges connecting ath and bth clusters; 2.
the numbers l aa of edges connecting the vertices inside the ath cluster.
We will suppose that in the notation l ab always a ≤ b, but for the clarity of formulas (summation conditions in (117) Further, it is necessary to count how many pairings configurations are described with the set of numbers l ab for a < b and l aa . The logic is similar to the case n = 2. We have to multiply: 1.
the number of ways to choose from every ath cluster l ab points going to bth cluster for all b. This is a typical problem of multiple choices from m a objects firstly l a1 , then l a2 and so on until l an objects (including l aa objects). As a result, this number equals to the multinomial coefficient ( m a l a1 ... l an ); 2.
for every ath cluster, the number of ways to form pairs inside of it from the remaining l aa points, which is l aa !
The expression (120) gives us exactly the result (117) after substituting the explicit form of M {l ab } . Let us note again, that all l aa have to be divisible by 2, otherwise such a configuration is impossible. This condition is reflected in the summation condition in (117).

Numerical Analysis of the Polynomial Approximations
In this subsection we are going to use different families of polynomials for the approximation and compare the results. Let us compare pointwise approximations of |t| α by Bernstein, Chebyshev and Legendre polynomials with each other. The results are represented on the Figure 3 for α = 4/3. For other α ∈ (1, 2) there is no significant difference. From the Figure 3 one can draw a conclusion that Legendre and Chebyshev polynomials are better for the approximation of integrands than Bernstein ones, since for the same degree 2N in the "majority of points" they give smaller error. Moreover, errors of Legendre and Chebyshev polynomial approximations oscillate and change sign, rather than Bernstein polynomial approximation, so one should expect that there will be some "cancellation of errors" in integrals, which will additionally raise the precision of the approximation. Now we are going to compare Bernstein, Chebyshev and Legendre polynomial approximations for integrals. We would like to simplify the general formula (111), since we don't want to waste too many resources for choosing the best variant of the listed three. Namely, we will: 1.
assume that the source j = 0 to escape additional parameter, as well as ε = 0, since it is possible to remove this regulator for zero source;

2.
replace all the graphs Γ in the sum inside the operator O (80) with the complete graph K n , which is the graph without loops, where every two vertices are connected with an edge. This means that ν ab (K n ) = (J n ) ab − δ ab (the matrix J n is the matrix of ones) for all a, b ∈ {1, . . . , n}, and the sum over graphs transforms into the contribution of the K n , multiplied by the total number of connected graphs on n vertices, which we will denote by |G C,n | (the explicit value of this number in terms of n is not needed); 3.
assume that G(x) ≡ G(0) = const for all x. From this follows that G n becomes proportional to a matrix of ones, i.e. G n = G(0)J n . Besides, (G n,K n ) aa = G(0) and (G n,K n ) ab = s ab G(0) for a ̸ = b.
This gives us the following rough approximation of connected Green functions GF, keeping in mind the definition (113): Now we can calculate coordinate integrals, since all the dependence of x a is left only in coupling constants g(x a ) as well as the result of the integro-differential operator over the variables s ab action. To that end, we should substitute the combinatorial formula (117). And we arrive at the following expression: and the last factor gives nothing else but the additional restriction of the summation condition l ab > 0 for a < b. However, we ignore this condition, despite the fact that this will lead to a worse approximation, since this will also lead to a decrease in the resources expended. Folding combinatorial sum back into correlator for the dimensionless covariance matrix (which is exactly the matrix of ones J n ), we get to the approximation: (123) From the formula (123) it is clear that for the comparison of different approximations it is useful to consider the following quantities H(n, N): where we have removed the common (depending on n only) factors for all polynomial approximations from (123), but left the factor 1/n! to avoid factorial growth. Recall that the degree of polynomial approximation in all the cases is equal to 2N. So in fact all the parameters of connected Green functions GF polynomial approximations are roughly encoded in function H, explicit formula for which is: Further, for different kinds of polynomial approximations there will be different coefficients u i and c q ( f ). In the following, we will mark the family of polynomials we use for the approximation in the index of H. Thus, we have: for Legendre polynomial approximation, and: for Bernstein polynomial approximation, as well as: for Chebyshev polynomial approximation. The plots of H L,B,Ch (n, N) for n = 2, . . . , 5 and N = 1, 2, 4 and α = 4/3 are presented in the Figure 4. The picture for other values of α ∈ (1, 2) is qualitatively the same. From the plots one can see that Legendre and Chebyshev approximations for integrals also converge much better than Bernstein approximation. At the same time, formulas for approximation by Legendre polynomials are simpler than for approximation by Chebyshev polynomials. This, in particular, explains our choice of the Legendre polynomials. Let us also note that using the second-degree approximation gives an error which is about 20% and the error of fourth-degree is about 10%. However, it is important to note that this error occurs in the model system considered in this subsection. As the results of the following subsections show, in a more realistic case, the errors can be larger.

Final Results of Polynomial Approximations for Connected Green Functions GF
According to all points described above, we will focus on Legendre polynomial approximation. And exactly this kind of approximation we will mean in the following under the term "polynomial approximation". We will write these expressions in compact forms using the already introduced notations as well as few new ones. Namely, we will additionally introduce the notations for the arguments of Gamma functions and a new operator O g : Then the expression (111) will take a form: where we have rewritten correlators, using dimensionless covariance matrix G n,Γ /G(0) for the convenience. And for zero source j similarly: In this form the structure of expressions is much more clear, so in the following we will use these formulas for representing connected Green functions GF.
For the reference we also write down the extended form of the expressions (127) and (128). Substituting all the introduced notation, we receive big, big formulas: In accordance with the notations, introduced in the section 2, the summation in the last line of the expression (129) is carried out over all l ab ≥ 0, satisfying the conditions n ∑ b=1 l ab = 2i + β a (where l ab = l ba ) and l aa is even. The expression (129) is the main polynomial approximation formula for the following computations. And for j = 0 we can also remove the regulator ε and write: As it follows from the discussion in the section 3, vacuum energy E = −G[0]. So in the following computations of vacuum energy density this formula will be our starting point. Finally, let us note that due to the convergence of Legendre series we have obtained in fact the expression for G I,ε,n [j] in terms of repeated series. In other words, in formulas (127-130) one can proceed to the limit N → ∞: Moreover, at least for j = 0 the sums over n and q can be permuted because of the absolute convergence of Legendre series. Hence, at least for G I [0] we have obtained the expression in terms of double series which is very convenient for calculations. Therefore: and the same for G I,ε [j]. These formulas are the main result of the presented paper. In the next subsections we will rewrite the result of applying the operator O g to the correlator to make this formula more similar to usual formulas from Feynman diagram approach.

Two Types of Graphs in the Formulas for Connected Green Functions GF
To go further, let us provide an improvement of graphical interpretation of the terms from Isserlis-Wick theorem (114), developed in the subsection 5.2.3. Our new graphical interpretation will be more convenient for our purposes as long as we deal with all the combinatorics.
We are going to consider the clusters (in terms of subsection 5.2.3) as vertices, and to draw the edges (corresponding to pairings) between the clusters themselves (in particular, starting and ending clusters can coincide). Therefore, we are going to describe the particular pairings' configuration in terms of the unoriented graph with possible loops and multiple edges, possibly disconnected. It is illustrated on the Figure 5, from which one can recognise the ordinary Feynman diagrammatics. Informally, we took the graphical interpretation from the subsection 5.2.3 and "constricted" every cluster to a point. And the pairings inside every cluster became loops on the corresponding vertices.
Considering a correlator ϕ m 1 1 . . . ϕ m n n (where we won't specify the covariance matrix in the notation for the brevity) we will enumerate vertices with the indices of corresponding fields ϕ a or with the coordinates x a (which come from the covariance matrix with the same indices). In such notations the degree of a vertex (the number of edges that a given vertex belongs to) with the number a equals to m a and l ab is the number of edges between the vertices a and b. Let us note that the adjacency matrix of the described graph is (l ab ) n a,b=1 . In the following we will call these graphs as Isserlis-Wick graphs. In the literature they are also known as Feynman graphs.
Therefore, in this paper we have two types of graphs: 1. Isserlis-Wick graphs, that are the subject of discussion in this subsection. We will denote them as Γ IW ; 2.
Meyer graphs, which came from the Meyer cluster expansion in the subsection 4.4, and still remaining in the operator O. We will denote them as Γ M ; And they are essentially different. While the Isserlis-Wick graphs can have loops, multiple edges or be disconnected, the Meyer graphs, in opposite, have no loops, multiple edges and have to be connected. However, there is a link between these two types of graphs, explaining also why physicists have only one notion of Feynman diagrams, serving simultaneously for all the purposes. Namely, for the Meyer graphs that give a non-zero contribution to the expression (130) the following inclusion holds: This inclusion is non-trivial for the sets of edges, since the sets of vertices are the same by definition. Let us prove this inclusion. Consider the formula (130). Next, consider concrete Γ M . Let it contain an edge between vertices a and b. In this case ν ab (Γ M ) = 1. Then differentiation with respect to s ab is non-trivial if and only if the corresponding value l ab > 0. But this means that there is such an edge in Γ IW defined by the given set l ab . Repeating this argument for all the edges in Γ M , we obtain the required inclusion.
In the same manner, if ν ab (Γ M ) = 0 for distinct a and b, then for every non-zero contribution has to be l ab = 0 for the same a and b. This means that if there is at least one edge in Γ IW , then there is an edge in Γ M .
The following useful statement is also true: the summation over Meyer and Isserlis-Wick graphs reduces to the summation over connected Isserlis-Wick graphs, at least for the zero source. Let us prove this statement. With the direct considering of cases, it can be shown that: However, for this product to be non-zero each of its terms have to be non-zero. And for every set {l ab } as follows from (134) there is only one possible adjacency matrix ν ab (Γ), satisfying this condition, and it has the following form: for all nonequal a and b. This expression can also be interpreted as the condition for Meyer graph to have the edge if and only if there exists at least one edge between these vertices in Isserlis-Wick graph. Therefore, for every {l ab } (every Isserlis-Wick graph) in the sum over connected (Meyer) graphs ∑ Γ∈G C,n there is no more than one term stays alive, and it is completely determined with the condition (136). Though, in the sum over Meyer graphs there could be no such a graph, since Meyer graphs are always connected but the condition (136) can be satisfied for disconnected graphs also. So the connectivity of Isserlis-Wick graph is the necessary and sufficient condition for the term to give the non-zero contribution.
Summarising, the operator O leaves only connected Isserlis-Wick graphs, and we finish at the expression: where the index "C" means that the summation in the correlator is carried out only over the connected Isserlis-Wick graphs. Equivalently, in the combinatorial sum (117) only those {l ab } should be taken into account, that correspond to the connected graphs. In calculations using computer algebra systems, one can check the connectivity of a Isserlis-Wick graph using built-in functions, recalling that (l ab ) n a,b=1 give the adjacency matrix. In conclusion, the expression (137) is an explicit representation of (connected Green functions) GF for fractional power interaction as a weighted sum of integer-power interactions contributions. However, all our effort were for deriving the explicit formulas for weights and obtaining the expression in terms of the converging series rather than asymptotic, which turned out to be quite a complicated task.

Simple Approximate Formula for General Term
In this subsection we will not use the general results on Meyer and Isserlis-Wick graphs, derived in subsection 5.2.6 since in the particular case that we are going to consider in this subsection, it is more convenient to carry out all the calculations in a different, simpler and more intuitive way. However, for the higher degrees of polynomials their use is unavoidable because of the calculations complexity.
As one can see from the model system, considered in the subsection 5.2.4, approximation with the Legendre polynomial of the second degree gives an error which is about 20%, and this is quite an accurate result. Remarkably, in this case one can finish up with a not very complicated formula for the vacuum energy. To obtain such a formula we will start from the expression (128) and assume g(x) = gχ Q (x): . If one writes down this formula in more details, there will appear the following coefficients: which we denote for the further use. So, we want to calculate G I,n [0] N analytically for N = 1, which corresponds to the second-degree polynomial approximation. We consider the terms with n = 1, 2 and n > 2 separately. For the first two orders the calculation is straightforward: for n = 1 and: for n = 2. Further it will turn out that this expression will occasionally satisfy the formula for n > 2, so we will put this term into the sum over n in the following. One important remark: in the calculation of such coordinate integrals one has to change variables and the integration domain will deform. Though, we are interested mainly in finding results in the thermodynamic limit (finding the leading contribution in V, when V → ∞). And in this limit the integration domain deformation is not important. We will describe the thermodynamic limit in subsection 6.1, and it will be more convenient to discuss such integrals there. Let us proceed to the case n > 2. Because of the operator O, the constant terms in all s ab will not contribute due to the differentiation. In more detail: it is so, since the sum in (80) is carried out over the connected graphs, so they have at least one edge for n > 1. As a result, in O at least exists one derivative, which will annihilate the constant term. This means, due to (117), we can consider only i ̸ = 0, so there will remain only a single term in a sum N ∑ q=0 q ∑ i=0 , namely i = q = N = 1: At present, for the evaluation of the Isserlis-Wick correlators for the second degrees of fields, it is useful to directly sort all the pairings rather than use the obtained combinatorial formulas. Namely, one can sort all the configurations of pairings for ϕ 2 1 . . . ϕ 2 n G n,Γ G(0) by the number r of vertices, pairing with themselves. Other n − r vertices have to constitute some number of closed chains, i.e. graphs whose vertex degrees are equal to 2, in terms of recently proposed graphical interpretation. This is the case since the degree of every vertex in any Isserlis-Wick graph is two, and all such graphs are disjoint unions of loops and chains. Then we can write: The summation indices i a in the right-hand side of the expression (142) take values in {1, . . . , n} and enumerate all different closed chains (one closed chain as a graph can be prescribed with different sequences of indices, what will be discussed in more details further). In the expression (142) we have written only the terms with the unique chain, and denoted the others as ". . .". We won't write them down explicitly, since in the following they will give only zero contribution. Further, substituting the last expression into the expression (141), we arrive at the following result: Because of the operator O, each term in the sum n ∑ r=0 survives if and only if ν ab (Γ) = 1 for a of b belonging to all the vertices which both lie in the chain. From this follows that the graph Γ has to be contained in this chain. But since Γ is connected, this condition can't be satisfied for the disconnected Isserlis-Wick graphs. This means, that in the sum in (143) survive only the terms corresponding to the connected graphs referring to pairings. In particular, all the terms denoted as ". . ." disappear, since they have at least two non-trivial chains and hence they are disconnected. And only the term with r = 0 in the remaining sum gives the non-zero contribution in ∑ Γ∈G C,n . One can note that this term corresponds to a cyclic chain, schematically written as Moreover, due to the invariance of every term under permutations of x a in coordinate integration, we can suppose that i a = a. And the number of different chains from n vertices equals to 2 n n!/2n, since n! is a total number of permutations of x a , and we have to divide it by the number of permutations, assigning the same chain, which is 2n. The thing is, we can make cyclic permutations of finite set of all x a , which don't change a chain (and there are n of such transformations) as well as reflections of finite set of all x a (which should not be confused with the transformation x a → −x a ), for which we have exactly 2 ways. In total it gives us the factor 2n. Moreover, we have to multiply the obtained number n!/(2n) by 2 n since we have two possible choices of fields in each vertex. Though this derivation is valid only for n > 2, this formula also gives the right result for n = 2. The illustration of different enumerations of vertices in a given chain is also presented on the Figure 6.
Calculating the integrals over all s ab except for s 1,2 , . . . , s n,1 , we obtain: where we have used that after permutation of x a there will be n! equal terms, so n! in the numerator and denominator cancel each other out. The remaining integrals over s ab are easily calculated: Let us change the variables to simplify this expression: and the last argument is the sum of all the previous: The Jacobian of such a transformation is equal to 1, since it is a linear map with uppertriangular matrix. Here we also assume the limit V → ∞, so we won't change the integration domain (since it is "big enough yet"), referring to the independent consideration of thermodynamic limit in the subsection 6.1. So we receive the result for n > 2: The coefficients A n,1,α,1 can also be simplified: so finally the expression for vacuum energy density (21) reads: where we assume the thermodynamic limit V → ∞.
The expression (144) is the simple approximate formula we worked for, and in the sections 6 and 7 we will apply it for the calculations of vacuum energy density for particular cases of propagators G. Unfortunately, the analysis of similar formula even for fourthdegree approximation polynomial, i.e. N = 2, becomes much more complicated due to significant variety of Isserlis-Wick graphs with the desired properties. Let us make a remark that in this subsection we have got in a different way all the results from 5.2.6 for the particular case, namely we have checked that all the non-zero contributions correspond to the connected Isserlis-Wick graphs.
Let us also note that in the limiting case α = 2: which coincides with the exact result for α = 2. We will obtain this result in the subsection 6.3. This is a simple test of the expression (144), since the function f (t) = t 2 belongs to the linear span of zero-degree and second-degree Legendre polynomials.

Hard-Sphere Gas Approximation
The obtained general formulas from the subsection 5.2.5 are still difficult for analytical calculations since of the coordinate integrals. So we have to make some assumptions to obtain simpler and hence more useful formulas. Namely, the formulas become significantly simpler if we use for G n the HSG approximation, inspired by statistical physics. In this subsection, we will also consider the coupling constant in the form g(x) = gχ Q (x), emulating a system with fixed coupling constant in finite volume V. In accordance with the expression (10) in the section 2, we keep all the diagonal elements equal to G(0), since they are constant, and "cut off" all the off-diagonal elements for |x a − x b | > δ. For the graphical illustration on the Figure 7, we will denote the right hand side of the expression (10) as G n,HSG . Finally, according to the section 2, we will denote the "volume of one hard-sphere particle" as v. Let us note, that all the quantities in HSG approximation will be expressed in terms of v and γ, so the obtained formulas won't have an explicit dependence on the definition of these parameters.
We start from the formula (82), which precedes the Parseval-Plancherel identity, in order to avoid the appearance of G −1 n,Γ and 1/ det G n,Γ . Substituting the notion of measure (68) and suggesting that j = 0, we arrive at the following expression: Applying for the matrix G n HSG approximation (10) and introducing auxiliary variables s ab , we obtain the following result: where we have used the fact, that: so the dependence of the coordinates factorises. Let us note that all the dependence of coordinates remains only in Heaviside functions. Now we can calculate the integrals over x a approximately with a commonly accepted approximation analogical to one in the statistical physics of hard-sphere gas. We can place the first particle in the full volume which gives us V and the following particle should be no further from the first than δ: n ∏ a=1 Q dx a n ∏ a<b θ(δ − |x ab |) ≈ Vv n−1 .
Traveling back to expression (146), we get: The expression (147) is much simpler than expression (146), because the integrals over x a in (147) have already been "calculated". Well, it is still a problem to calculate directly the remaining integrals over t a , so we use the developed technique of polynomial approximations. Exactly like in the subsections 4.4.2 and 5.2, we get to the following expression ("taking out of brackets" all the transformations that lead us to it): where we also used the convergence of the constructed polynomial approximations and wrote the result as a series as well as removed all the regulators. Finally, substituting the explicit expression for the correlators (117) and using the definition of vacuum energy density w vac (21), we arrive at the result: One can also rewrite the last combinatorial sum, introducing l aa = 2q a :  (150) is useful for decreasing of the computational time for numerical calculations using the formula (149). Summarizing, we went ahead with some relatively simple approximate expression for vacuum energy density w vac for a generic nonlocal propagator G. It can be simply realized using any soft suitable for symbolic computations, such as Wolfram Mathematica. Its main advantage is that unlike in general formulas for GF one doesn't need to compute numerically or analytically coordinate integrals. This fact makes us suppose that HSG approximation is good enough for investigating the qualitative and partly quantitative properties of nonlocal theories.

PT Calculation of System's Physical Characteristics
Using the results obtained one can calculate the following characteristics of the quantum scalar field in nonlocal theory: 1.
quantum scalar field vacuum energy density; 2.
quantum scalar field Green functions in terms of functional integral over primary field ϕ, for example, two-particle Green function; 3.
quantum scalar field composite operators Green functions, which are also often called "form factors".
In this paper, we will focus on the vacuum energy density computation. We will use the first few exactly calculated terms of PT series in coupling constant g as well as approximate expressions for general terms and explore what physical results can be obtained from them. Derivation for the non-zero source general case remains the same but includes more technical details and large calculations. Everywhere in sections 6 and 7 we will suppose g(x) = gχ Q (x) in all the obtained formulas. We start by considering special cases n = 1, 2 of the general formula (81) (everywhere in this section we refer to graphs only as the Meyer graphs 5.2.6): 1.
for n = 1 there are no s ab variables, since a < b, and the only graph is the single-point, therefore: 2. for n = 2 there is one s ab variable, namely s 12 , and the only one connected graph without loops Γ, for which ν 12 = 1. So the matrix G 2,Γ has a form: and the formula (81) in this case reads as follows: Let us note, that we consider only the zero source case, for which the index I of GFs can be omitted due to the definitions (25) and (26). Now we are going to calculate the first orders of PT from section 4 exactly. And it turns out that for n = 1, 2 it is possible to do without using the constructed polynomial approximations or HSG approximation. For n = 1 the result can be expressed in terms of elementary, and for n = 2 in terms of Gauss hypergeometric functions.

First Orders of Complete Green Functions and Connected Green Functions GFs
In this subsection we will provide exact expressions for Z n for n = 1, 2. Then we get from them the corresponding G n . We start with the expression (54): The matrix G n can be degenerate on null sets, though is doesn't affect the convergence of integrals because of the obtained majorant in 4.3.
The first order calculation is straightforward. Applying the scaling transformation of ϕ and calculating the integrals over ϕ and then over x, we arrive at the following result: The second order calculation demands more mental and technical effort. First, we write the matrix R 2 := (G 2 ) −1 , inverse to the matrix G 2 : values of this function on the eigenvalues of A, and we won't present its derivation here, referring to textbooks in functional analysis. Expanding the logarithm of an operator, using equivalent for analytical functions and bounded operators definition of function of an operator, we obtain the following series: Now we rewrite the powers of operators in terms of integrals, using the notion of operator's G integral kernel G(x − y) and also substitute the coupling constant g(x) = gχ Q (x), which cuts off the integration domains: Similarly to the subsection 5.2.7, we change variables and reduce the number of integrals by one, which leads to the appearance of the factor V: Finally we can write down the expression for the vacuum energy density (21): From this expression one can see that it coincides with Legendre polynomial approximation for α = 2, obtained in the subsection 5.2.7.
6.3.2. Comparison of PT Formulas for the First Orders for α = 2 with the Exact Result Let us recall the formulas (154) and (160) for G 1 and G 2 and substitute α = 2 into them. Taking into account, that: we receive the following compact expressions: which coincides with directly obtained result for α = 2. So, our formulas pass this simple test, which witnesses indirectly their rightness.

Research of First PT Terms in Vacuum Energy Density
In this subsection we are going to calculate the first PT terms for vacuum energy density. Formerly we have obtained the general formulas (154) and (164), and now we are going to substitute different propagators into these formulas: virton propagator [3], which is purely non-local, and Euclidean Klein-Gordon propagator (with sharp cut-off), which admits a local limit. Then we will consider a local theory as a certain limit of a nonlocal one, being its regularization. However, let us note that the nonlocal theory is of important independent interest.

Research of Second-Degree Legendre Polynomial Approximation Formula
Let us start from the general formula (144) for second-degree Legendre polynomial approximation (in the thermodynamic limit, in which the vacuum energy density is of greatest interest): In the following, we will consider only the terms with n ≥ 2, since the case n = 1 is already simple enough. Now we are going to substitute different propagators into this formula.
1. Nonlocal Case. For the virton propagator the integrals in the expression (182) are Gaussian, therefore, can be calculated explicitly: Thus we arrive at the following expression: for the entire function: and the parameter: The summation over n indeed starts from one, but it is convenient to separate the linear term into two parts.

2.
Local Case. Again we restrict ourselves to the cases d = 2, 3. We have already calculated G(0) for Klein-Gordon propagator. For integrals, that we have for n ≥ 2, the following chain of equalities is true: .
We removed the UV regulator again. This formula is valid for d = 2, 3. It is convenient to introduce the following dimensionless parameters: In contrast to the nonlocal case, in the local one we define the parameter z in terms of m. In the following subsection we will explore the way to scale parameters of a nonlocal theory to obtain a non-trivial local limit. Looking ahead, we announce that keeping exactly these parameters constant is a necessary condition of the non-trivial local limit existence. So, in these terms the expression for the vacuum energy density reads: The numerical plots of (183) and (187) are presented in the Figure 9. Their discussion will be given in the subsection 7.5. Let us note, that these formulas are not very friendly for numerical calculations, though, their asymptotics may be of separate interest. In the present paper, however, we restrict ourselves to numerical results.

Local Theory as a Limit of Nonlocal One
In this subsection we are going to research briefly whether it is possible to scale nonlocal theory in some way to get a local one. We start from the writing the generic term of the polynomial approximation for G I [0]. Everywhere in this subsection the term "converging integral" means the converges of integral when the region of integration extends from Q to R d . And the same for the term"diverging integral". We also will use the notations from the formulas (117) and (128).
According to these formulas, we have the following n-particle term for the G I [0] (to simplify the analysis, only one term from the combinatorial sum is written out and up to a some numerical coefficient depending on the term): We have written this expression from the dimension considerations and translational invariance of the considering propagator. Let us assume that all the coordinate integrals except one converge, and in this case m is the only dimensional parameter in the integrals. So for the vacuum energy density we have the following expression: for some function Ψ (formal power series). Therefore, in order to get a non-trivial local limit one should scale m, g and G(0) so that the arguments of Ψ stay constant. However, it is only a necessary condition but not sufficient, since for the degrees of the polynomial approximation higher than two there appear diverging integrals and one have to solve the renormalization problem. It is well-known that in d = 2 all the power theories are renormalizable and in d = 3 only the theories with power which ≤ 5 are renormalizable. This means that for d = 3 it is impossible to gather all the "infinities" coming from the monomials with degrees higher than 4 into some new parameters. So, we can't write a better approximation for d = 3 than the fourth-degree approximation.
The conclusion is that being applied to local theories in d > 2 our method has a limited precision of degree 4 in d = 3 and degree 2 in d ≥ 4. Though, even the second-degree approximation can provide some tool for (rough) quantitative research.

Hard-Sphere Gas Approximation
The formulas for vacuum energy density (149) and (150), obtained using HSG, are convenient and simple for numerical calculations. As a practise in statistical physics and numerical simulations show, HSG approximation gives worthy results, consistent with the experiment. Unfortunately, we have to deal with the power series in system parameters, so it requires accurate work with precision of calculations to distinguish the errors from real results. We make the truncations of both series in the expression (149), namely n 0 for index n and N for index q, and use a numerical experiment to explore the truncations influence on the results. We assume α = 4/3, and one can check that there will be no qualitative differences for others α ∈ (1, 2).
Recall that in the HSG approach we have approximated the matrix G n by the Heaviside step function and parameters γ and δ in accordance with the expression (10). Let us introduce new variables: where ε vac is the vacuum energy, contained in one "hard-ball" (specific energy). Further, we are going to plot the dependency ε vac (γ, z), since for HSG approximation γ and z are the most natural parameters. Figure 11. Plot of specific energy ε vac (γ = 0.5, z) (left) and surface of ε vac (γ, z) (right) for n 0 = 3, N = 3 and α = 4/3. they could be compared with virton and Klein-Gordon cases separately, but after the proper change of parameters to γ and v according to (12).

Analysis of HSG Approximation
In the case of HSG approximation we have found the inflexion point. It was observed from the third-order in the coupling constant g and the sixth-order polynomial approximation, which we encoded in the plot as (3, 3) (three particles and sixth degree). It slightly changed with increasing the approximation order that was checked up to (4, 3) (four particles and sixth degree). Thus we expect that this result stays relevant for higher orders of approximation.
Both the virton and Klein-Gordon Gaussian theory Green functions satisfy the HSG approximation conditions. Though, the error of this approximation demands separate research, we can extract some predictions from this fact. For example, we can guess that the next orders of polynomial approximation addition (without HSG) should result in arising of an inflexion point in virton and Klein-Gordon cases.

Analysis of PT First Orders
In the general propagator case there is no point of inflexion in orders g and g 2 (165) and the function w vac is convex upwards (concave function). The existence of the inflexion point will be determined with the next PT orders.
We have also calculated the order g 3 for the fourth-degree approximation in the virton theory, but we won't present it here because of this paper size. However, this contribution is positive for g > 0, so it produces an inflexion point referring to the virton propagator. Recalling the behaviour of polynomial approximations in the HSG approach, we expect that a fourth and higher PT orders as well as higher degrees of approximations won't affect the existence of such an inflexion point. In the Klein-Gordon case the coordinate integrals are much more complicated, so we postpone its research for future papers.

Analysis of Virton and Klein-Gordon Second-Degree Approximation Formulas
At the same time, we do not observe an inflexion point in the case of second-order polynomial approximation for the virton (183) and Klein-Gorgon (187) vacuum energy densities. Numerically its absence follows from plots, and analytically from the asymptotics of the corresponding series. Namely, we have: for some positive constants A, B, a, b. We don't present the derivation of these asymptotics in this paper since they are used only for reference. From the written asymptotics it follows that both functions are convex upwards, which (combined with numerical plots) proves the absence of the inflexion points in the approximations of the second degree. Moreover, their convexity differs from the one predicted with the PT first three orders (at least, for the virton case) and the HSG consideration. Recall that we have deduced from these approaches that the right behaviour at infinity is downward convex (convex function). Thus we come to the conclusion that these formulas are not applicable for strong coupling limit.
Moreover, the vacuum energy density in Klein-Gordon theory behaves non-physically in the large coupling constants region because it becomes negative for sufficiently large z. It can be seen both from plots and asymptotics. Moreover, from the asymptotic, it follows that such negative values are reached for any ξ if z is sufficiently large.
Summarising the discussion of the second-degree polynomial approximation formulas, they fail to reach a strong coupling limit. Though, they can describe the intermediate values of the coupling constant (rather than very small only), and therefore they can be the computational alternative for PT first orders. The plots of vacuum energy density for second-degree polynomial approximation should be cut on some value of g.

Discussion
Let us start the discussion section from a brief overview of the work carried out in the paper. This paper starts from the theory with fractional power self-interaction R d dx g(x)|ϕ(x)| α with α ∈ (1, 2) for the nonlocal quadratic part. We rewrite the interaction potential through the successive application of straight and inverse Fourier transform (41). After that, in GF we expand the exponent with potential (50) and face the problem: we need to calculate complicated Gaussian expectation value of fractional potential Fourier transform. The solution to this problem arises from the Parseval-Plancherel identity (52), which brings us back to the original potential. Proof of the obtained series convergence for zero sources is presented in subsection 4.3.
Our next goal is to introduce the connected Green functions GF (18) using Meyer theorem (66). It generalizes the standard technique commonly used in statistical physics and results in the general expression for the connected Green functions GF in terms of coupling constant powers series (79).
The main part of our work focuses on the research of the "fractional Gaussian moments" (101) expansion in terms of integer moments. We perform careful approximation with the Legendre polynomials after restriction of the integration domain on a compact (hypersphere) S n−1 in every n-particle contribution, resulting in (128), rewrite the obtained sum over graphs in a form convenient for direct calculations, and end with the expression for connected Green functions GF in terms of weighted sum of theories with integer powers contributions, obtaining this way an immediate generalisation of Feynman technique, but with converging series.
In the final part we focus on the explore of the obtained PT applications for the cases of virton and Klein-Gorgon propagators as well as for the HSG case. Thus, we provide comprehensive thorough research of various generating functionals. There is a brief list of the obtained results: 1.
Construction of non-asymptotic PT series for generating functional Z in powers of coupling constant g (53), proof of its convergence in section 4.3 and derivation of similar formulas for GF G = ln Z (82). We assume that these PT series have a strong coupling limit; 2.
Derivation of the general formulas (53) and (85) for Z n and G n in fractional potential scalar field theory (1) with arbitrary nonlocal quadratic part of the action. These expansions are the generalization of common PT where we explicitly have a sum over Wick graphs (Feynman diagrams) with additional weights; 3.
Construction of calculable and converging approximations for Z n and G n (111) with any precision, obtained with the polynomial approximation of function f (t) = |t| α in [−1, 1]. It results in the fact that we can extend our attention to non-integer potential theories and rewrite it through a weighted sum of integer potential theories. As a consequence for the probability theory, we provide the series representation for multidimensional Gaussian distribution moments of fractional order, which hasn't been described in the literature yet, as far as we know; 4.
Application of HSG approximation for Z n [0] and G n [0] (148), obtaining vacuum energy density (149) and clarifying its dependence on parameters of the theory (189); 5.
Finding the inflection point for the vacuum energy density in HSG approximation (Figure 10). We expect that this result remains if one adds the next orders of approximation and also in the exact result.
There were made multiple checks and comparisons of the obtained results.

Conclusions
Traveling back to the very beginning, to the motivation of the research, presented in this paper: the research of nonlocal field theories with interaction potential g|ϕ| α may shed light on the strong coupling behaviour of nonlocal ϕ 4 theory providing such a definition