A post-Gaussian approach to dipole symmetries and interacting fractons

We use a post-Gaussian variational approach to non-perturbatively study a general class of interacting bosonic quantum field theories with generalized dipole symmetries and fractonic behaviour. We find that while a Gaussian approach allows to carry out a consistent renormalization group (RG) flow analysis of these theories, this only grasps the interaction terms associated to the longitudinal motion of dipoles, which is consistent with previous analysis using large $N$ techniques. Remarkably, our post-Gaussian proposal, by providing a variational improved effective potential, is able to capture the transverse part of the interaction between dipoles in such a way that a non trivial RG flow for this term is obtained and analyzed. Our results suggest that dipole symmetries that manifest due to the strong coupling of dipoles, may robustly emerge at low energies from short distance models without that symmetry.


Introduction
Fractons are exotic quasiparticles arising in lattice models which have raised a remarkable interest in both high energy physics and condensed matter [1][2][3][4][5][6][7][8].Continuous models exhibiting fractonic behaviour are realized through a specific class of field theories possessing a generalized notion of global symmetry known as subsystem symmetries and/or multipole symmetries.These are sensitive to the details of the underlying lattice and not compatible with Lorentz invariance [9].
Subsystem symmetry defines symmetry operators acting independently on subspaces of the total system.The most simple examples of these models involve a single free real scalar field.As a fact, any interaction term allowed by this symmetry is irrelevant under renormalization, making these theories essentially free at low energies [7,8].
A closed related symmetry is multipole symmetry, and concretely dipole symmetry, the type of symmetry in which we will focus in this work.In models with dipole symmetry, a U(1) charge and the associated dipole moment are conserved [10].Contrarily to models with subsystem symmetries, models with dipole symmetries can be made rotationally invariant.Interacting models consistent with these symmetries are built using complex scalar fields.The resulting bosonic theories are characteristically non-Gaussian and strongly correlated as, in contrast to standard field theories, the only symmetry allowed terms with spatial derivatives in the Lagrangian have four or more powers of the fundamental fields [10].This makes these theories very hard to analyze through mean field or perturbative techniques.
References [11,12] have studied large N versions of models with conserved dipole moment in static and out equilibrium settings respectively.The large N approach allows to study the physics of these systems at a finite interaction strength.While they obtained the Green's function (charge-charge propagator) of the models, the large N factorization of higher correlation functions does not allow to compute the dipole propagator from which it is possible to fully access the non-Gaussian structure of the dipole-dipole interaction.
Regarding renormalization group (RG) analysis of these models, it was pointed out in [13] that microscopic models with dipole symmetries may look quite unrealistic, highly anisotropic, and fine-tuned.There, authors found that adding as a perturbation, dipole-symmetry-disallowed kinetic terms to the type of large N models commented above, the one-loop β-function shows that the UV model, which is not dipole-symmetric but a more realistic and less fine-tuned model, flows to the dipole-symmetric point at the deep infrared limit, regarding the dipole-symmetry as an emergent phenomenon at low energies.In [14] a model of weakly interacting scalars with subsystem symmetries and non-relativistic fermions in 2+1 dimensions was analized in perturbation theory.It was shown that the 1-loop β-function of the fracton-fermion coupling does not flow due to an emergent symmetry of the effective Lagrangian.In [15], authors, using renormalization group analysis, investigated the emergence of exotic non-Fermi liquids in two dimensions, in a model where fermions interact with a bosonic fracton system.In [16], a RG procedure was applied to a lattice model with fractonic behaviour and authors found that all screening effects generating non-transverse fractonic couplings vanish, thus indicating that dipole interactions might only exist transverse to the dipole direction.Indeed, they also found that the coefficient of the transverse fractonic term did not get renormalized.
Given these results, it would be desirable to invoke a nonperturbative analysis of models with dipole symmetries which, in addition to their manifest theoretical interest, are expected to be realized both experimentally and in nature [17][18][19].
In this work, we carry out this analysis through a variational approach based on the functional Schrödinger picture in QFT [20,21].Being the Schrödinger representation not manifestly Lorentz invariant poses no particular problem here because we deal with theories whose symmetries are incompatible with Lorentz invariance.We use a method that allows to extend predictions for these systems beyond mean field theory and the Gaussian variational approch in a systematic way.Our method is based on nonlinear canonical transformations (NLCT) [22][23][24] from which it is possible to build extensive variational wavefunctionals which are certainly non-Gaussian.
Our results suggest that, in spite of all of the particularities of the models that we have addressed, our approach, provides useful insights into the physics of field theories with fractonic behaviour.Concretely, we find that a Gaussian approximation, while allowing to carry out a consistent RG flow analysis of the theory at 1-loop based on the Gaussian effective potential (GEP), only grasps the interaction terms associated to the longitudinal motion of dipoles.This agrees with the large N approach taken in [11,12] and is not surprising, as the Gaussian variational approach coincides with the large N solution for an interacting model.
Remarkably, our proposal of non-Gaussian ansatz wavefunctional based on NLCT, by providing (-a variational approximation of-) the connected part of the four point function between charges (a.k.a dipole two point funcion) by means of an improved post Gaussian effective potential, is able to capture the transverse part of the the interaction between dipoles in such a way that a non trivial RG flow for this term is obtained and analyzed.This is something that cannot be derived from a pure Gaussian ansatz or large N techniques.Namely, our non-perturbative calculation of the RG flow through its β-functions, suggests that dipole symmetries that manifest due to the strong coupling of dipoles, arise at low energies from less symmetric (in the dipole sense) UV models.
The paper is organized as follows: in Section 2 we review some background on fractonic field theories with a dipole symmetry.In Section 3 we carry out an analysis of these type of models through a variational Gaussian wavefunctional representing the ground state of the theory.Once a self-consistent solution of the variational parameters is provided, we carry out an analysis of the RG flow of the theory based on the Gaussian effective potential.Section 4 presents a post-Gaussian variational wavefunctional builded through non-linear canonical transformations (NLCT).Through the NLCT technique, we compute the connected part of the dipole-dipole interaction that contributes to improve the energy expectation value and the effective potential, from which we study the RG flow of the model through the β-functions for the dipole-dipole interaction terms.We provide a final account of the work in Section 5.

Interacting Fracton Model
In this Section we review some background on field theories with dipole symmetry.A complex scalar field φ has been used to describe these theories, invariant both under global phase rotations φ → e iΛ φ corresponding to conservation of charge, and also invariant under linear phase rotations, φ → e i d•x φ for constant d, corresponding to conservation of dipole moment [10].This is a special case of theories that poses a generalized notion of global symmetry known as multipole symmetries (in this case, dipole) [9].

Dipole symmetry
We consider a translational and rotationally invariant continuum quantum field theory in (D +1) dimensions of a charged scalar field φ(t, x) which, in addition, it is also invariant under U (1) phase rotations, and dipole transformations.That is, under a U (1) phase rotation Λ and dipole transformation d, the scalar field transforms as (2.1) In order to build an action invariant under this transformation it is useful to find covariant operators, that is, operators that are invariant under this transformation up to a phase factor [10].The result is that, in contrast to standard field theories, these theories do not possess any covariant operators featuring spatial derivatives, that is, acting on only a single φ operator.Instead, as it has been shown in [10], the lowest order covariant spatial derivative operator contains two factors of φ, taking the form: which transforms covariantly as D ij (φ, φ) → e 2(iΛ−d•x) D ij (φ, φ); see [10].Using these covariant operators, we can then write down a lowest order action for this theory as: where λ is a dipole-dipole coupling constant and U ( φ φ) is an interaction potential term for the isolated charges.We note that the usual spatial kinetic term ∂ i φ ∂ i φ is forbidden by dipole symmetries, and the simplest terms with spatial derivatives involve at least four powers of the charged field written above.
For the sake of convenience in the rest of this paper, we note that dipole symmetry is simpler to analyze in momentum space.Namely, a dipole transformation acts on φ(t, k) as and on the conjugate field by the opposite shift, that is, as translations in momentum space.Thus, dipole-invariant interactions must be invariant under translations in momentum space and therefore they must only depend on momentum differences.As an illustration, one may write the interaction D ij (φ † , φ † )D ij (φ, φ) in momentum space as [25] where we use the notation Dk ≡ d D k/(2π) D and the vertex is defined as with k mn = k m − k n .

Model of interacting Fractons with dipole-symmetry
In this work, following [10] and [11] we consider a more "sophisticated" model than (2.5) given by the action is the traceless part of D ij (φ, φ) and γ is a complex parameter.A mass term is included which have no interplay with the spatially dependent phase rotation.The constants λ T (T is for tensor) and λ S (S is for scalar) are arbitrary couplings describing the longitudinal motion of dipoles(S) and the transverse one (T ) respectively.This model posses a characteristic non-Gaussian form and thus it is hard to deal with it.
In this work we use variational approaches through the functional Schrödinger picture in QFT [20,21].The functional Schrödinger picture works with the Hamiltonian of the theory in (2.7), which, written in momentum space reads as [11] where the conjugate momentum of the field φ(p) is π(p) ≡ −iδ/δφ(−p), in such a way that [φ(p), π(q)] = iδ D (p + q) and the interaction vertex between dipoles is given by The Hamiltonian (2.9) is bounded from below for λ T ≥ 0 and λ T + D λ S ≥ 0 [26].In this theory, the only allowed processes are those in which the total dipole moment is conserved.Hence, isolated charges have restricted mobility and thus show fractonic behaviour.Remarkably, the non-Gaussian dipole-dipole interactions do not forbid a fractonic charge to move completely as far as the (restricted) charge mobility arises from processes in which the total dipole moment is conserved.Namely, the locality of dipole-dipole interactions in (2.9) forces that such processes must be mediated by a propagating dipole between the two charges [5,10].
Given this qualitative expected behaviour, in this work we will mainly be focused on providing some quatitative responses on i) the characterization of the dipole dispersion relation and ii) using standard renormalization group (RG) analysis, investigate how and to what extent it can be considered that the dipole-symmetry is an emergent phenomenon at low energies.In the next Sections, we investigate on this topics by invoking variational nonperturbative techniques.

Gaussian Variational Approach
In this Section we carry out an analysis of the theory given by the Hamiltonian in (2.9) through a variational approach using the functional Schrödinger picture in QFT [20,21].Being the Schrödinger representation not manifestly Lorentz invariant, poses no particular problem here as we tackle with theories whose multipole and subsystem symmetries are incompatible with Lorentz invariance.Therefore we use a Gaussian variational wavefunctional representing the ground state of the theory.By the variational method we provide a self-consistent solution of the variational parameters and carry out a RG flow analysis of the theory based on the Gaussian Effective Potential [27,28].

Gaussian wavefunctional
The variational Gaussian wavefunctional that we take as an ansatz for the ground state of the theory (2.9) is given by where the variational kernel G(k) satisfies and N = [det(2πG)] −1/2 .As it will be commented below, the real space representation of G(k), G(x − y) = φ(x) φ(y) , may be interpreted as a the one point function of a dipole creation operator, that is, an operator creating a charge at y and anti-charge at x and thus it can be used as a dipole order parameter [11].

Ground state energy density
The energy density of the ground state, using Wick's theorem, is given by where Vol ≡ δ D (0) and in the last line we have used that being the forward scattering limit of V (4) (k 1 , k 2 , k 3 , k 4 ) [11].This is in accordance with the Gaussian ansatz being a consistent truncation of the Dyson series.In fact, this imposes that the model is consistent only when the forward limit is positive, i.e, λ S > 0. With this, one may write after defining the self-energy Σ(k) as The variational parameter G(k) is thus obtained by which yields the truncated Dyson equation (gap equation) and thus Several comments are due here.
1.The variational principle and the Gaussian form of the variational wave functional (3.1) implies a gap equation that determines the Green function G(k) in the static case.
2. The variational Green function G(k) acts as order parameter for dipole breaking: if the position space Green's function φ(x) φ(0) δ D (x), then the dipole symmetry acts on it.In other words: G can be understood as the one-point function of an operator that creates a dipole, with a charge at 0 and an anticharge at x. Thus, being G(k) nonzero, it is understood as a dipole condensate in analogy with the usual charge condensate associated with ordinary symmetry breaking [11].
3. We note that the dipole symmetry acts on the momentum space field by φ(k) → φ(k − d), i.e. by a shift of momentum, and on the conjugate field by the opposite shift.Dipole transformations then act on (3.10) 4. For infinite volume, the theory in Eq (2.7) posses a continuous non-compact dipole symmetry R D .This implies the existence of a family of solutions Ψ d G φ, φ with the same energy, labeled by allowed dipole transformations.In a lattice regularization, d will be valued in the space of momenta, that is the Brillouin zone.Thus, the number of allowed dipole transformations, and so Ψ d G φ, φ solutions related by dipole symmetry, would be then equal to the number of lattice sites N sites ∼ Vol • V BZ , with V BZ the (-dimensionless-) volume of the Brillouin zone [11].
5. As commented above, the dipole symmetry is lattice sensitive, and thus the theory (2.7) exhibits features of UV/IR mixing characteristic of fracton models.Nevertheless, as stated in [25], the UV-sensitivity showed by these models, is mild in comparison with models exhibiting subsystem symmetries.As discussed in [29], for the later, the UV/IR mixing refers to the low-energy mixing among small and high momenta.In other words, depending on the chosen direction, the low-energy modes can have very high momenta (see Appendix A).As a consequence it has been suggested that renormalization group is not applicable to models exhibiting fractonic behaviour [30].However, it has been argued that this can be addressed by conforming the integration of the high-energy modes to the symmetries of the fracton models [31].
6.In a theory without interactions, G(k) ∝ 1/m, that is, there are no propagating charges which amounts to a momentum dependent self-energy Σ(k) = 0. Instead, interactions imply Σ(k) = 0 and thus allow restricted charge mobility.The precise nature of this mobility can be established by finding self consistent solutions for Σ(k).Unlike the case of subsystem symmetries, Σ(k) in (3.6) is rotationally invariant.Indeed, following [11], we note that the vertex V (4) (k, q) ∼ ||k − q| 2 + 2γ| 2 is a polynomial in |k| of degree 4, and thus, according to (3.6), so Σ(k) is too.Therefore, we make a rotationally invariant ansatz with parameters (a 0 , a 1 , a 2 ) to be self-consistently determined.In Appendix B, we find these solutions numerically.
7. From Eqs (3.4) and (3.5), one realizes that the Gaussian approximation to the ground state only "sees" the S-term of the dipole-dipole interaction, as only captures the forward scattering limit of (2.10).This result agrees with the large N approach taken in [11,12] and is not surprising, as it is well known that the Gaussian ansatz coincides with the large N solution for an interacting model.An obvious question is if there is a systematic way to improve the ansatz going beyond the Gaussian regime and if this procedure is capable to grasp the transverse part of the interaction vertex.This will be addressed in the next Section but before, we will analyze further the Gaussian solution.

A more general Gaussian wavefunctional
For subsequent developments, it is worth to consider a more general Gaussian ansatz.To this end, we consider a wavefunctional that accounts for the existence of a nonzero charge condensate, and which can be obtained from Ψ G φ, φ in (3.1) by with Operationally, expectation values w.r.t Ψ SG amount to expectations values w.r.t Ψ G under the substitutions Under this ansatz, the energy density reads with Σ σ (k) given by We note that in (3.16) where in the last equality, we used that V (4) (0, 0) = λ S γ 2 and we imposed For σ = 0, this single condition imposes that Σ(k) to be gapless at k = 0. Contrarily, if Σ(k) is gapped, the system has to be such that σ = 0 (U (1) is unbroken).In other words, this single restriction allows to consider different phases of the model as it is explained in Appendix B, [11].With this, the energy density now reads as and the variational procedure yields the truncated Dyson equation (gap equation),

Renormalization and β-function
Following [28], we define the Gaussian Effective Potential (GEP) The GEP usually contains divergences only because it is written in terms of bare parameters m and λ S .Once reexpressed in terms of renormalized parameters m R and λ R , V eff [σ] becomes manifestly finite.This reparametrization of the theory, or "renormalization", can not change the physical content of the theory.That is, choosing different renormalized parameters m R and λ S would lead to a different looking but equivalent V eff [σ].As noted in [32], the GEP is exactly renormalization-group (RG) invariant, which means that all ways of defining the renormalized parameters are equivalent, and our task reduces to merely finding the most convenient one.Noteworthily, this is not the case for the effective potential computed at one-loop, for which the RG invariance is spoiled by a "renormalization-scheme-dependence problem", just as in perturbation theory.
Renormalized mass.From the Gaussian effective potential (3.23), one may define the renormalized mass m R as which gives where we have introduced the convenient notation The renormalized parameter m R , as defined above, is very convenient because it turns out to be the mass of a one-particle excitation in the σ = 0 vacuum [28].It is a finite quantity once the momentum integral is dimensionally regularized as explained in Apendix B.
Renormalized coupling.The renormalized coupling constant is defined by and is given by, where λ B = 2γ 2 λ S = 2 V (4) (0, 0) and With the renormalized parameters (3.25) and (3.28), it would be possible to write a finite expression for V eff [σ].This would be quite convenient in order to study the phase diagram of the theory and its critical points.We leave for a future investigation the use of the GEP to study the phase diagram of the models under consideration while focusing the scope of this work on knowing how the coupling constant changes with the energy scale.
To this end, now following Coleman and Weinberg [27] and Barnes and Ghandour [28], we define the renormalized coupling as where the mass scale M is nonzero but arbitrary.The result yields As M is merely an arbitrary substraction point [27], it is non-physical and thus, no physical quantity can depend on it.In particular, for V eff , one may write the RG-flow equation This implies that we may define a β-function as which reads as with Λ(M ) ≡ M 2 I 3 (M ), and It is convenient to write β = β(λ M ).To this end we replace λ B with λ M everywhere except for first order terms.This substitution is valid to 1-loop, since it induces changes only for higher loops [20].That is, we may invert Eq. (3.31) to solve for λ B in terms of λ M to obtain Then substituting into (3.34)one obtains the β-function at 1-loop To debunk the exact behaviour of λ M on the energy scale M , requires to exactly determine Λ(M ).This can only be done numerically, using the methods described in Appendix B. In Figure  the case shown in the figure, there is no fixed point λ F for which β(λ F ) = 0.This means that (3.37) is always positive and thus implies that the S-part of the interaction vertex grows towards the UV from an arbitrary bare inital value.

Post-Gaussian Variational Approach
In this Section we present a post-Gaussian variational wavefunctional builded through nonlinear canonical transformations (NLCT) [22][23][24].The goal is the following: as shown above, the Gaussian ansatz only captures the longitudinal/scalar part of the interaction but obviates the transverse/tensor part.Here, using a concrete class of post-Gaussian wavefunctionals generated through the NLCT technique, we compute the connected part of the dipole-dipole interaction (which amounts to a connected part of a four point function between charges).From this, it is possible to directly compute many interesting physical properties of the model (2.7) which are hidden for both the Gaussian ansatz and the large N approach.Before going into this, it is worth to highlight some relevant issues related with applying a variational method in QFT [33][34][35].To study the the most relevant features of the ground state of the theory, the choosen variational states must posses nice calculability properties.That is to say, in order to first optimize and then compute physical properties of the theory, one needs to evaluate expectation values of operators such as n-point functions.Here the problem is obvious as evaluations of expectation values of these quantities with respect to general non-Gaussian states are highly limited.In practice, this fact has restricted the use of variational wavefunctionals to the Gaussian case.This poses a severe limitation to explore truly nonperturbative effects.For instance, Gaussian wavefunctionals, which amount to a set of decupled modes in momentum space, cannot describe the prototypical interaction between the high and low momentum modes of an interacting QFT, not to say the interplay between these modes in theories with the strong UV/IR mixing effects of those that we are considering in this work.For this reason, a feature that one should require for a variational wavefunctional is to include parameters describing the interplay between the high energy and low energy modes and how the former affects the low energy physics of our theory.
In the following, such a class of variational ansatz is used to study the theory in (2.7).

Non-Gaussian states through NLCT
Following [22][23][24], a class of extensive non-Gaussian variational wavefunctionals can be nonperturbatively built as1 where The functional structure of (4.2) implies that the expectation value of an operator O, that in general depends on φ, φ, π, π, with respect to Ψ N G reduces to the computation of a Gaussian expectation value of the transformed operator O → U † O U. The relevant poit here is that concrete and suitable forms of B, while leading to non-Gaussian wavefunctionals (4.1), automatically truncate after the first non-trivial term, the infinite nested commutator series which inevitably arises as one applies the Hadamard's lemma.As a result of this truncation, the calculation of any expectation value Ψ N G | O | Ψ N G reduces to the computation of a finite number of Gaussian expectation values (calculability).In addition, the exponential nature of U asures an extensive volume dependence of observables such as the energy of the system.In other words, the non-Gaussianities captured by a trial wavefunctional of the form (4.1) persits in the thermodynamic limit.Finally, being U unitary, the normalization of the state is preserved to the one imposed to the gaussian wavefunctional.
In this work, we build a post-Gaussian class of wavefunctionals through an operator B of the form2 with p,q 1 ,q 2 ,q 3 ≡ Dp Dq 1 Dq 2 Dq 3 .Here, α is a variational parameter that keeps track on the deviation of the wavefunctional and any observable from the Gaussian case.Following [22][23][24], it is understood that an efficient truncation of the commutator series in (4.3) is such a one that terminates after the first non-trivial term when (4.3) is applied to the canonical fields of the theory.To this end, it is straightforward to see that the variational function f (p, q 1 , . . ., q m ), that is symmetric w.r.t. the exchange of q i 's3 , must be constrained to observe f (p, q 1 , . . ., q m ) = 0 , p = q i , and f (p, q 1 , . . ., q m ) f (q i , k 1 , . . ., k m ) = 0 , ( for m = 1, 2, 3.These constraints force that the action of U on the canonical field operators of the model (2.9), φ(k), φ(k) and π(k), π(k) is given by and with f (p, q, k, r) π(p) φ(q) φ(r) δ D (k − p − q − r) . (4.9) The quantities in (4.7) and (4.9) are nonlinear field functionals that shift the degrees of freedom of the canonical fields of the theory by a nonlinear polynomial function of other degrees of freedom.
Being U unitary, one can check that Eqs.(4.6)-(4.9)ensure that i.e, the canonical commutation relations (CCR) still hold under the nonlinear transformed fields.For this reason, the transformations above are known as nonlinear canonical transformations (NLCT).
A remarkable property of the wavefunctionals (4.1) generated through NLCTs is that non-Gaussian corrections to Gaussian correlation functions can be obtained in terms of a finite number of Gaussian expectation values.Let us illustrate this with an example that is relevant when com-puting the expectation value of the Hamiltonian (2.9).Under the NLCT implemented through B in Eq. (4.4) we have where • • • refers to an evaluation w.r.t the non-Gaussian wavefunctional (4.1) and • • • refers to a Gaussian expectation value w.r.t (3.1).For reading convenience, we have schematically compressed the details of the arguments of the f -variational functions.It must be noted that the Gaussian expectation values appearing above contain an equal amount of fields φ as conjugate fields φ.As it will be shown, this is a convenient feature in order to variationally cover the structure of the vertex

Post-Gaussian improved ground state energy density
Now we are interested in evaluating E N G ∼ H with H in (2.9).In order to this, apart from the four-point function stated above, we must compute where we used (4.6) and (4.9).When taking into account the truncation constraints (4.5) one obtains the same result as for the Gaussian ansatz, that is Remarkably, while under dipole transformations, the nonlinear field shift in Eq (4.7) does not posses a clear transformation, the correlators above clearly state that the dipole one function still behaves as With this, now we compute the dipole-dipole interaction term in the Hamiltonian (2.9).Before that, we note that when evaluating the interaction term with the non-Gaussian ansatz, the O(α 0 ) term amounts to the Gaussian expectation value4 with O(α, α 2 ) refering to non-Gaussian corrections.To be more explicit, denoting the Gaussian energy expectation value by E G , one may write where ∆ E (1) and ∆ E (2) refer to O(α) and O(α 2 ) non-Gaussian corrections respectively.Interestingly, these correction terms are completely related to the connected part of the four-point function φ(k 1 ) φ(k 2 )φ(k 3 )φ(k 4 ) .That is, substracting the disconnected terms in Eq. (4.11) one obtains In words, the connected part of the four point function completely determines the the non-Gaussian corrections to the ground state energy or an improved version of the effective potential in (3.23).
∆ E (1) correction.This term corresponds to the evaluation of χ (1) 2 with, χ By symmetry, χ 1 + χ 2 ≡ χ (1) so we focus on the first one, which explicitly reads with f (k 3 , p, q, r) ≡ f (k 3 , p, q, r) δ D (k 3 − p − q − r).Using the truncation constraints (4.5) and Wick's theorem, the result reads as with Γ 1 (p, q, r) = c(p, q, r) V (4) (−r, −q, p + q + r, −p) (4.20) where we have used the compact notation c(p, q, r) ≡ f (p + q + r, p, q, r) .(4.21) Therefore, the O(α)-contribution to the post-Gaussian improved energy expectation value ∆ E (1)  can be written as It is worth to remark that in (4.22), the "variationally dressed vertex" Γ 1 (p, q, r), captures a range of values of (2.10) larger than the forward scattering limit of the Gaussian term.In other words, ∆ E (1) grasps features of the tensor part of the dipole-dipole interaction that are inaccessible to the Gaussian approach.∆ E (2) correction.This term corresponds to the evaluation of After a lengthy albeit straightforward calculation, the result reads as ∆ (4.25)As before, the variationally dressed vertex Γ (4) 2 (p, q, r, s) covers a wider range of values of the vertex interaction (2.10) than the forward scattering Gaussian limit, thus capturing the λ T part of the dipole-dipole interaction.
The post-Gaussian improved energy density E N G in Eqs (4.15), (4.22) and (4.24) as said above, is builded in terms of a field transformations that shift some degrees of freedom of the canonical fields through a nonlinear polynomial function of other degrees of freedom.This is implemented by the structure of the variational parameters c(p, q, r).Concretely, in Appendix C we adopt a functional structure for c(p, q, r) such that the IR modes of φ are nonlinearly shifted by its UV modes.In this sense, it is expected that an optimized set of c(p, q, r) parameters can conform, in a variational sense, the integration of the UV modes to the symmetries of the model (2.7).
This integration helps to build the "effective variational vertices" Γ 1 (p, q, r) and Γ 2 (p, q, r, s) involving three and four propagators G(k) respectively.In this sense, noting that (in the forward scattering limit) V (4) can be understood as a vertex involving two propagators G(k), Γ 1 in (4.20) is the combination c − V (4) and Γ It is important to remark at this point that general equations for the optimal values of the variational parameters G(k), Σ(k) and c(p, q, r) can be obtained, given a fixed α5 , by deriving E N G w.r.t.them and then equating to zero.This yields a set of non-linear coupled equations that must be self-consistently solved.Finding analytical expressions for this solutions is difficult.In this regard, our aim here is to provide expressions that explicitly show the relation between the variational parameters and the coupling constants of the model under consideration.With The variational parameter c(p, q, r) (white arrow) combines with V (4) to build the "3-legged" variational vertex Γ (4) 1 and the correction to energy density ∆ E 1 Eq (4.22).Down: Diagramatic representation of (4.25) (left).Two variational parameters c(p, q, r) (white arrows) combine with V (4) to build the "4-legged" variational vertex Γ this aim, from here in advance we choose α in such a way that with the coupling constants λ's not necessarily small.With this choice, the optimization equations greatly simplify and i) G(k) and Σ(k) can be reduced to the Gaussian solution Eq. (3.22) and ii) the optimization of the variational parameters c(p, q, r) is decoupled from them and can be carried over the non-Gaussian correction ∆E (1) + ∆E (2) .The simplification occurs as a consequence of the structure of the optimization equations (discussed in Appendix C) and is never due to any additional assumptions.It is clear that if one chooses α not fulfilling (4.26) for solving the equations, G(k), Σ(k) have not the Gaussian structure anymore.In the limit posed by (4.26), it is obvious that the main contribution for a post-Gaussian improvement of the ground state energy of the theory in (2.7) is given by ∆ E (1) .

Renormalization and β-functions
In this subsection, we obtain an improved effective potential with the aim of studying the RG flow of the model (2.7) through the β-functions for the dipole-dipole interaction.
The renormalization of post-Gaussian variational calculations is rather involved.In [39] it was showed that contributions generated by the NLCT-wavefunctionals, which are related to 1particle irreducible diagrams not taken into account by the Gaussian approximation, require new counterterms that in principle could be determined systematically by taking the derivatives of the effective potential.Nevertheless, this approach is preconditioned by an analytic optimization of the energy expectation value, that, in general, is not feasible (see Appendix C).
Fortunately, it is possible to make advance by using simplified functional forms for the variational parameters c(p, q, r).This approach can be taken as a compromise route to gain further understanding of the problem [22][23][24]39]; while it leads to a sub-optimal approximation, a particular ansatz for the NLCT c(p, q, r), motivated by the optimization of a part of the energy expectation, amounts to a renormalized post-Gaussian effective potential (see Appendix C).
The following discussion thus takes for granted that the (sub-)optimal values of the variational parameters of the ansatz have been found.
NLCT-Improved Effective Potential.Here we consider the the non-Gaussian corrections to the Gaussian effective potential in (3.23) generated by the wavefunctional with Ψ SG defined in (3.12).These corrections are given by with and 2 (p, q, r, s) G(p)G(q)G(r)G(s) 1 (0, q, r) G(q)G(r) , According to the definitions in (4.20) and (4.25), we observe that while ∆V 0 eff [σ] and ∆V 2 eff [σ] capture both the T -part and the S-part of (2.10), the term ∆V 4  eff [σ] only captures the S-part.It is thus convenient for our interests to split V eff into a T -term and S-term.To this end, we decompose (2.10) as After this splitting, one can write the post-Gaussian improved effective potential as with and Here, V G eff [σ] refers to the Gaussian effective potential in (3.23) and the T and S labels refer to the vertex splitting in (4.31).
RG flow and β-functions.In order to obtain the running of the coupling constants with the energy scale, following [27] as before, we define the renormalized couplings as At this point, we note the following: by working in the scaling limit where α λ T, S 1, it is reasonable to consider that a good approximation for the renormalized S-coupling and its RG-flow β-function is given by the Gaussian term, that is to say, λ M is given by Eq. (3.31) and where Λ S (M ) is the one appearing in (3.37).It is obvious that the non-Gaussian terms provide α-corrections to the S-interaction that might be interesting to compute.However, we foresee that a bigger interest relies on seeing how the non-Gaussian terms in (4.33) provide a non-vanishing β-function for the T -interaction part between dipoles.Therefore, here we will focus on V T eff [σ] that, conveniently rearranged, reads as In addition, we note that when plugging (4.38) into (4.36), in the limit (4.26), the contributions coming from O(α 2 ) terms in (4.38) can be effectively discarded and keep only considering the O(α) terms in way such that, (4)T 1 T (q − r, 2p + q + r) .
With this we write, where Let us first focus on I 2 (M ).Specifically, let us note the result of the second derivative w.r.t σ 2 of the arguments inside the momentum integrals yields ≡ V (4) (k, 0).At this point we remark that terms with higher number of "propagators" G(k) and "vertices" V (S) k correspond to higher loop terms in the dipole-dipole interaction captured by the non-Gaussian ansatz.Therefore, it is sensible that we keep in our calculation, only those terms with the lowest number of them.In this sense, the terms multiplied by the arbitrary substraction point M involve a high number of G(k)'s, and we will assume that they can be safely discarded w.r.t the other terms.Under the same assumption, I (1) 0 (M ), which contain terms with two "vertices" V (S) k and seven "propagators" G(k), may be discarded as well.With this, one can approximate where λ S and λ T are the bare parameters in (2.7) and T (q, 2 p + q) + (p ↔ −p) T (q − r, q + r) V (S) q G(r)G(q) 3 + (q ↔ r) . (4.46) A renormalized T -coupling can be defined by taking the substraction point to be M = 0, as that, interestingly, depends on both bare parameters of the model.It is remarkable also that it is possible to define this renormalized coupling only for α = 0. Finally, we compute the RG flow of λ M , which after a lenghty albeit straightforward computation yields and where we use the compact notation T (q, 2 p + q) + (p ↔ −p) , h 1 2 (p, q, 0) = V (S) p c(p, q, 0) V T (q, 2 p + q) + (p ↔ −p) , h 2 1 (0, q, r) = V (S) q c(0, q, r) V T (q − r, q + r) , h 2 2 (0, q, r) = V (S) r c(0, q, r) V T (q − r, q + r) . (4.51) Let us summarize our results as where we note that Λ S (M ) in the first equation is a rescaled version of the one appearing in (3.37).Eqs (4.52) represent the main results of this work for several reasons: 1. First, they show that our post-Gaussian ansatz, by providing (-a variational approximation of-) the connected part of the four point function (a.k.a dipole two point funcion) is able to capture i) the transverse component of the dipole-dipole interaction and consequently, ii) a nonvanishing β-function for this part of the interaction.This is something that cannot be derived from a pure Gaussian ansatz or large N techniques.
2. As it was commented in Section 2, the Hamiltonian (2.9) is bounded from below for λ T ≥ 0 and λ T + D λ S ≥ 0 [26].Thus, our non-perturbative calculation of the β-functions shows an intringuing feature.Noting that, being Λ T (M ) a priori, a non negative function of M 2 , the result in (4.52) suggests that, with respect to the transverse part of the dipoledipole interaction, the model is asymptotically free, that is, the renormalized coupling λ M decreases as we move into the UV.This is a highly non-trivial result.This kind of dipoleasymptotic freedom suggests that symmetries that manifest only due to the strong coupling of dipoles, may arise at low energies from less symmetric, in the dipole invariance sense, UV models.Of course, it would be customary to check, by carrying out explicit optimizations of the ansatz, if Λ T (M ) is always positive or posses zero crossings that would imply more involved behaviours.
3. Regarding the last point, there is an interesting limit for the model in Eq. (2.7) that can be reached by taking λ S → 0 + while holding λ S λ T = λ fixed.In this limit, the renormalized coupligs are well defined and read as In addition, both β-functions vanish, with The latter implies that the tensor part of the dipole-dipole coupling does not flow for any arbitrary initial bare value of λ T .From this point of view, in this limit, the model amounts to a well-defined quantum field theory with fractonic behaviour, a result that seems to be consistent with those in [14] and [16].The RG flow diagram plotted in Fig. 4.2, shows the limit commented above in the region where the strenghth of the flow is close to zero (dark region).
4. It is worth to frame these results within the context of the robustness in QFT exposed in [7].There, it is commented that usually, in condensed matter systems, the UV models posses not many global internal symmetries G U V , being them typically, Z 2 , U (1), or the trivial one.It is assumed that by a fine tuning of the short distance parameters, one may find a low-energy theory with an enhanced emergent global symmetry G IR , for instance, a dipole symmetry.The interesting question is whether the low-energy model will preserve the emergent symmetry G IR once the short distance parameters are slightly deformed.Our results in Eqs.(4.52) suggest that a small deformation of the short-distance system can not wreck the G IR dipole symmetry at long distances and therefore, we can say that this symmetry is robust.

Conclusions
In this paper we have studied the ground state properties and the RG flows of a continuous model of interacting fractons using a nonperturbative variational approach based on the functional Schrödinger picture in QFT.
We have found that a Gaussian approximation, while allowing to carry out a consistent RG flow analysis of the theory at 1-loop, only grasps the interaction terms associated to the longitudinal motion of dipoles.
In order to investigate the true non-Gaussian features of these models, we have proposed a non-Gaussian ansatz wavefunctional based on NLCT.Without losing sight on the fact our results have been obtained for a concrete NLCT, we argue that those are illustrative and general enough (in the NLCT sense) as to provide a useful kind of understanding on the essence of the nonperturbative physics associated to the models under consideration.Concretely, by providing a variational improved effective potential, the post-Gaussian ansatz is able to capture the transverse part of the the interaction between dipoles in such a way that a non trivial RG flow for this term is obtained and analyzed.This is something that cannot be derived from a pure Gaussian ansatz or large N techniques.Our non-perturbative calculation of the RG flow suggests that dipole symmetries that manifest due to the strong coupling of dipoles, robustly arise at low energies from UV models without that symmetry.In other words, these symmetries account for an emergent phenomena rather than being well defined properties of microscopic models.Some interesting future directions would be to perform explicit optimizations of the improved effective potential in order to study the phase transition between the broken and unbroken phases of the theory (see [40] for a recent work on this line).Indeed, as the the theory is non-Gaussian and strongly coupled close to the critical point transition, even a mean field approximation of the transition is not possible.It would be interesting to clarify if the phase transition is a second order continuous one or first order.Other interesting future research lines consist in extending these kind of nonperturbative analysis to dipole-symmetric fermionic models, possibly obtained by generalizing [11], and to apply similar techniques to study interacting models possesing subsystem symmetries.
Finite Temperature We have shown that for the Gaussian ansatz, the energy density reads as The procedure to obtain the free energy density at a finite inverse temperature β is simple and can be extensively discussed in [23,41].For the Gaussian approach to the model in (2.7) this amounts to, which tends to the ground state values as β → ∞ as well as F 0 asymptotes to (B.1) .Operationally speaking, the expectation value of any quantity at finite temperature can be computed using the replacements where the equality allows a simple separation into IR and UV contributions: the first term is the zero-temperature result, and the second term, which involves the relativistic Bose statistical factor, decreases exponentially at large momentum.Noteworthily, the self-energy at finite temperature reads Σ(k) explicit solutions From (B.4), we write a regularized version of the the self-energy (B.4) as where the vertex V (4) (k, q) = λ S 4 ||k − q| 2 + 2γ| 2 is a polynomial in |k| of degree 4 and so Σ(k) is too.The regularization above effectively removes contributions of high energy modes as coth(βx) ∼ 1 as βx → ∞, that is, the prescription removes the contributions of high energy momentum modes.As the theory in Eq (2.7) is rotationally invariant, we make a rotationally invariant ansatz Σ Plugging this ansatz into (B.5) and picking the terms in the vertex that yield rotationally invariant terms, one obtains with the constraint σ a 0 = 0 and λ = λ S /4.This amounts to a set of coupled equations for the parameters (a 0 , a 1 , a 2 , σ) that can be solved numerically.In this work we have been mainly interested in solving the equations for σ = 0 and β → ∞.We show some of these solutions in

C Appendix C. Optimization of the post-Gaussian ansatz
Optimization of c(p, q, r) The post-Gaussian energy expectation value is given in (4.15) and here we conveniently rewritte it as = p,q,r Γ 1 (p, q, r) + Γ (4) 1 (−p, q, r) G(p) G(q) G(r) , 2 (p 1 , p 2 , q 1 , q 2 ) G(p 1 )G(p 2 )G(q 1 )G(q 2 ) . (C.1) As commented in the main text, we choose α in such a way that with the coupling constants λ T and λ S not necessarily small.With this choice, the optimization equations greatly simplify as G(k) and Σ(k) can be approximated at O(α) to the Gaussian solution in Eq. (3.22).This leaves the optimization of the variational parameters c(p, q, r) as a decoupled process that can be carried out by doing δ E N G δ c(p, q, r) = −2α δ χ(1) δ c(p, q, r) + 4 α 2 δ χ(2) δ c(p, q, r) = 0 .
As a result one obtains the integral equation 2 r V (4) (−r, −p, q, p − q + r) [c(−q, −p, r) + c(−q, p, r)] G(r) = V (4) (p, q) , (C.4) which can also be written as 4 q,r V (4) (−r, −p, q, p − q + r) [c(−q, −p, r) + c(−q, p, r)] G(q) G(r) = Σ(p) , (C.5) where, noting that α amounts to a normalization value for the parameters c(p, q, r) that we use as a tracking parameter, solutions to the equation above refer to the fully meaningful quantity c(p, q, r) = α c(p, q, r).Analytical solutions to this equation may be challenging to obtain, even approximately.In this situation, it thus difficult to gain any insight on the structure of these parameters.Fortunately, it is possible to make advance by using a simplified functional form for the variational parameters c(p, q, r).
Functional form of f (p, q 1 , q 2 , q 3 ) A suitable way of accomplishing (4.5) is decomposing f (p, q 1 , q 2 , q 3 ) = g(|p|, |q 1 |, |q 2 |, |q 3 |) L(p) H(q 1 ) H(q 2 ) H(q 3 ) , (C.6) where g(|p|, |q 1 |, |q 2 |, |q 3 |) is scalar function to be determined by energy minimization and we have imposed that L(p)•H(p) = 0, that is, the domains of momenta where L and H are different from zero must be disjoint, up to sets of measure zero.Refs [22,23] provide a useful functional form for L and H as with ∆ 0 < ∆ 1 being two variational and coupling dependent momentum cutoffs and Γ(x) ≡ θ(1 − |x|) with θ(x) the Heaviside step function.On very general grounds, f (p, q 1 , q 2 , q 3 ) can be understood as separating the Fourier components of the field φ into non-overlapping domains of "high"-H and "low"-L-momenta.The variational parameters ∆ determine the size of these nonoverlapping regions in momentum space.The effect of the functional form of f (p, q 1 , q 2 , q 3 ) on, for instance, Eq. (4.6) is to shift the "low"-L-momentum modes of φ by a nonlinear function of the "high"-H modes, providing therefore, once optimized, a variational integration of the effects of the high energy modes into the low energy physics of the problem.With this, and noting that the main contribution for a post-Gaussian improvement of the ground state energy is given by χ(1) in (C.1), a (sub-)optimal way of finding the variational parameters f (p, q 1 , q 2 , q 3 ) consist in [22,23,38]
(3.1) we show the behaviour of Λ(M ) for D = 3.To this end, Σ M (k) in (3.17) is numerically determined for different values of M .The Figure (3.1)shows a positive function of the energy scale M .For
4.25)  amounts to c − V(4) − c .This suggests the interpretation that each c(p, q, r) creates a new "insertion" point into the variational vertex to which a new propagator can be attached.This is pictorially shown in Fig 4.1 where the different terms in Eq. (4.15) are presented diagramatically.

Fig. 4 . 2 :
Fig. 4.2: The RG flow implemented by Eqs (4.52).The figure plots the vector field ( β, β)| M =0 and it has been generated assuming that Λ T (0) ∼ Λ S (0).The arrows indicate the direction of the flow and the background color refers to the strengh of the flow, from weak or zero (dark) to a bigger values (orange).In the dark central region the model does not flows and is a well defined QFT.The rest of the diagram shows that there is no perturbation of the UV parameters that could ruin the fractonic IR behaviour (robustness).

(D. 5 )
Being U = e B unitary, one can check that Eqs.(D.4)-(D.5)ensure that the canonical commutation relations (CCR) still hold under the nonlinear transformed fields.