# Cooperating Cracks in Two-Dimensional Crystals

Shizhe Feng1,2, Xiaodong Zheng3, Pengjie Shi2,
Thuc Hue Ly4, Jiong Zhao3,∗, Zhiping Xu2, Corresponding authors:
<EMAIL_ADDRESS>(Z.X<EMAIL_ADDRESS>(J. Z.) 1Department of
Mechanics, Sichuan University, Chengdu 610065, China
2Applied Mechanics Laboratory, Department of Engineering Mechanics, Tsinghua
University, Beijing 100084, China
3Department of Applied Physics, The Hong Kong Polytechnic University, Kowloon
999077, Hong Kong, China
4Department of Chemistry, City University of Hong Kong, Kowloon 999077, Hong
Kong, China

###### Abstract

The pattern development of multiple cracks in extremely anisotropic solids
such as bilayer or multilayer two-dimensional (2D) crystals contains rich
physics, which, however, remains largely unexplored. We studied crack
interaction across neighboring 2D layers by transmission electron microscopy
and molecular dynamics simulations. Parallel and anti-parallel (‘ _En-Passant_
’) cracks attract and repel each other in bilayer 2D crystals, respectively,
in stark contrast to the behaviors of co-planar cracks. We show that the
misfit between in-plane displacement fields around the crack tips results in
non-uniform interlayer shear, which modifies the crack driving forces by
creating an antisymmetric component of the stress intensity factor. The cross-
layer interaction between cracks directly leads to material toughening, the
strength of which increases with the shear stiffness and decreases with the
crack spacings. Backed by the experimental findings and simulation results, a
theory that marries the theory of linear elastic fracture mechanics and the
shear-lag model is presented, which guides the unconventional approach to
engineer fracture patterns and enhance material resistance to cracking.

Fracture is a catastrophic process that leaves cracks behind as materials
fail, which could also be harnessed to fabricate structures with specific edge
shapes in a positive perspective [1, 2]. The problem of cracking has been
rationalized in the theory of fracture mechanics since the seminar
contribution of Griffith in 1921 [3]. However, the main focus of studies has
been placed on single cracks nucleating and growing in a solid, although in
practice, cooperative dynamics of interacting cracks are often of vital
importance by defining the fracture patterns and energy dissipation processes
[4, 5]. Model problems of two interacting cracks have been studied in the
literature [6]. In-plane parallel cracks propagating in the same direction
kink away from each other depending on the mixing of fracture modes [7]. _En-
Passant_ or anti-parallel cracks repel and then attract while approaching each
other [8, 9, 10], which is attributed to shielding and amplification of the
crack-tip stress fields [11, 12].

Two-dimensional (2D) crystals and their layered assemblies have drawn notable
attention for promising applications in materials and devices [13]. One of the
unique features of these materials or structures is the high contrast between
the basal-plane and interlayer mechanical properties [14]. The weak van der
Waals (vdW) and/or electrostatic forces transfer the mechanical load across
the 2D layers such as graphene, hexagonal boron nitride (h-BN), and metal
dichalcogenides (TMDs), compete with the elasticity of 2D crystals and
modulate the lattice deformation [15]. It is thus of both fundamental and
practical interest to explore the behaviors of multiple cracks propagating in
different layers, which are modulated by the interlayer interactions. $In\
situ$ scanning electron microscopy (SEM) studies show bifurcated fracture
patterns in multilayer graphene, where mechanical energy dissipation in post-
fracture events of interlayer slip improves the fracture toughness [16].
Complex crack-tip behaviors were reported to be governed by highly variable
interlayer friction. Cracks advancing across a bilayer MoS2 can either
propagate, be blocked, or branch due to the combined effects of lattice
defects and interaction between the layers. Cracks propagating in one layer
can even break the other layer with small flaws located in their paths [17].
These experimental phenomena imply the significance of interlayer coupling in
defining fracture patterns in layered materials.

Here, we carried out an $in~{}situ$ scanning transmission electron microscopy
(STEM) study of the interaction between cracks advancing in neighboring layers
of 2D crystals. The underlying mechanisms governed by interlayer strain
transfer are explored by atomistic simulations and elucidated continuum
mechanics analysis, which offer an unconventional approach to enhance crack
resistance and material toughness.

Figure 1:  Attraction between cracks propagating in a bilayer of 2D crystals.
(a) Schematic illustration of the experimental setup. (b, c) Cracks in a
twisted ${\rm MoS_{2}}$ bilayer characterized by the atomic-scaled high-angle
annual dark field-scanning transmission electron microscope (HAADF-STEM). TL:
top layer. BL: bottom layer. Inset: fast Fourier transform (FFT) patterns
denoting the twist angle is $8^{\circ}$. (d) Geometric phase analysis (GPA)
shows asymmetric strain distribution around the crack tips in (c).

In the experiments, molybdenum sulfide (MoS2) bilayers were synthesized
through atmospheric pressure chemical vapor deposition (APCVD). Twisted
bilayers were obtained for the weak interaction between the 2D MoS2 layers
that cannot maintain epitaxial growth. Precracks were introduced to the
bilayers by exposure to ultrahigh intensity _e_ -beam ($>1$ pA/nm2), the
propagation of which is characterized by atomic-scaled high-angle annual dark
field-scanning transmission electron microscopy (HAADF-STEM). The
characterization was operated on the Thermo Fisher Spectra 300 TEM with a
field-emission gun at 300 kV [18]. The areas of observation for the crack
dynamics are completely outside of the initial area exposed to intensified _e_
-beams, therefore experiencing low dose exposure and free from noticeable
e-beam damage. The precracks in the two respective layers of a twisted MoS2
bilayer are often misaligned due to the anisotropic fracture toughness defined
by the lattice symmetry. However, these cracks tended to approach each other
during propagation (Fig. 1). Geometric phase analysis (GPA) [19] on the in
situ STEM results shows asymmetric strain distributions around the crack tips,
which are attributed to the interlayer shear among the stacking area between
cracks (Fig. 1(d)).

To rationalize the interaction between cracks across 2D layers, we carried out
molecular dynamics (MD) simulations to probe the dynamics of a pair of cracks
using bilayer graphene as a representative model, which can be extended to
other 2D crystals with a honeycomb lattice by modifying the force field
parameters (see Supplemental Material for details [20]). Graphene is chosen
for the study instead of MoS2 in the experiments because the corresponding
force fields of the former are more mature and verified.

Figure 2:  Fracture patterns of parallel cracks in bilayer graphene simulated
by molecular dynamics (MD). (a) Crack paths with edge precracks spaced by $H$,
obtained from simulations without and with interlayer van der Waals (vdW)
interaction. (b) The atomic shear strain obtained from bond-distortion
analysis [19]. The strain distributions at the crack tips in TL and BL show
reflected patterns about the $x$-axis. The dashed lines denote the axis of
symmetry of the atomic shear strain around the crack tip as the crack starts
to advance. (c) Crack paths predicted at different interlayer shear stiffness
$G_{\rm int}$ and $H=5$ nm. $G_{0}$ is the interlayer shear stiffness of the
graphene bilayer [21]. (d, e) The dependence of $\alpha=d_{\rm m}/H$ on $H$
(d) and $G_{\rm int}/G_{\rm 0}$ (e), where $d_{\rm m}$ is the distance
traveled by the cracks before they ‘merge’.

Our simulation results show that, firstly, interlayer shear partially releases
the strain energy and modifies the fracture toughness (Fig. S1 [20]). The
critical mode-I (tensile) stress intensity factor (SIF) or fracture toughness,
$K_{\rm Ic}$, is enhanced by increasing the crack spacing, $H$, and the
interlayer shear stiffness $G_{\rm int}$ (Figs. S1(a) and S1(b) [20]). By
comparing the crack configurations and stress concentration around the crack
tips with or without the interlayer interaction (Fig. S1(c) [20]), we find
that interlayer shear adds a restoring force, hinders crack opening, and
alleviates the in-plane stress concentration at the crack tips. This cross-
layer coupling directly leads to effective toughening.

On the other hand, the fracture patterns show that two parallel cracks in
adjacent layers attract each other and eventually ‘merge’ into a single path
in the $x$-$y$ plane (Fig. 2). The fracture paths display lattice kinks, which
are attributed to the lattice discreteness and fracture-resistance anisotropy
[22, 23, 24]. By ignoring the kinking features, the strength of attraction is
assessed through the crack deflection angle ($\theta_{\rm c}$) from the
distance traveled by the cracks before they merge ($d_{\rm m}$) (Fig. 2(a)).
The value of $\theta_{\rm c}$ at the coarse scale ($\sim d_{\rm m})$ is
related to the geometric factor $\alpha=d_{\rm m}/H$ as ${\rm
tan}(90^{\circ}-\theta_{\rm c})\sim\alpha$. The simulation results show that
$\alpha$ increases linearly with $H$ (Fig. 2(d)), suggesting that $\theta_{\rm
c}$ is inversely related to $H$. The attraction is thus weakened as the
overlap area between the cleaved edges of the two cracks increases. The
attraction between cracks is absent by turning off the interlayer interaction.

The strain distribution in the bilayer is analyzed by bond-distortion analysis
[19]. As the crack starts to propagate, the symmetry to the crack direction
(the $x$ direction) is broken and gradually recovered for configurations with
large values of $H$ (Fig. 2(b)). A mode-II SIF (shear, $K_{\rm II}$) emerges
at the crack tips even under mode-I loading as the result of the interlayer
interaction. The control of crack paths through $K_{\rm I}$ and $K_{\rm II}$
can be rationalized, for example, through the criterion of maximum energy
release rate (MERR) criterion [25]. The crack deflection at the coarse scale
depends on the competition between the driving force and the material
anisotropy in fracture resistance defined by the crystal symmetries as
evidenced by the kinking features. Consistently, experimental observation
identifies asymmetric stress at the crack tips of parallel cracks, modifying
their advancing directions (Fig. 1(d)). In addition, the effect of interaction
strength is studied by rescaling all the energy-related coefficients
($\epsilon,C,C_{6}$) of ILP potential [26]. The results show that the strength
of attraction between cracks increases with $G_{\rm int}$ (modified from the
original value of $G_{0}$), and the ratio $\alpha=d_{\rm m}/H$ decreases as a
result (Figs. 2(c) and 2(e)). These MD simulation results uncover the
underlying mechanisms with an atomic-level picture.

The stress focusing at the crack tip yields highly non-uniform displacement
fields. For bilayers containing one crack each, the difference between in-
plane displacement fields in the two layers results in interlayer shear. The
interlayer slip $\Delta{\bf u}$ of the graphene bilayer can be characterized
in the unit of the C-C bond length $b=0.142$ nm. Figs. 3(a)-3(c) shows
localized slips between two parallel cracks spaced by $H=5$ nm. The slip
amplitudes $\Delta u_{x}$ and $\Delta u_{y}$ are concentrated around the crack
tips, which are anti-symmetric and symmetric around the axis $x=0$,
respectively. Linear elastic fracture mechanics (LEFM) analysis (Supplemental
Material [20]) shows that the spatial distribution of $\Delta{\bf u}$ agrees
well with the MD simulation results (Figs. 3(d)-3(f)). The amplitude of slip
is slightly lower in the simulations due to the self-adjustment of lattice
registries [27, 28], which is not included in LEFM.

Figure 3:  Interlayer slip distribution. (a-f) Interlayer displacements,
$\Delta u_{x}$ and $\Delta u_{y}$, and the curl,
$\lvert\nabla\times\mathbf{u}\rvert$, predicted by MD simulations (a-c), and
by continuum theory analysis (d-f) for precracks with $H=5$ nm, respectively.
The lines in panel B show the distribution of $\Delta u_{x}$, $\Delta u_{y}$,
and $\lvert\nabla\times\mathbf{u}\rvert$ along the $y$-axis at the crack tip.

The discrete nature of 2D crystals results in the formation of Moiré patterns
in the bilayers, the distortion of which indicates the lattice strain. Both
experimental and simulation results show that crack opening results in
relative twisting of the bilayers. By ignoring the lattice registry at the
interface, we can further assume that the interlayer interaction follows a
linear constitutive shear model. This approximation remains valid before
interlayer sliding is activated and $|\Delta{\bf u}|<b$. The interlayer shear
stiffness $G_{\rm int}$ can then be treated as a constant, and the interlayer
shear stresses are $\tau_{xz}=G_{\rm int}\Delta u_{x}/h,\tau_{yz}=G_{\rm
int}\Delta u_{y}/h$ (Fig. 4(a)), where $h=0.335$ nm is the interlayer
distance. A shear-lag analysis (Supplemental Material [20]) is then used to
calculate the stress intensity factors (SIFs, $K_{\rm I}$ and $K_{\rm II}$).

The direction of crack deflection $\theta_{\rm c}$ at the coarse scale can
then be determined from MERR (Supplemental Material [20]). By using the
interlayer shear modulus of graphene ($G_{\rm 0}=3.48$ GPa) [21], we find that
the value of ${\rm tan}(90^{\circ}-\theta_{\rm c})$ decreases with $H$ for
$H<2$ nm, and then increases linearly with $H$ (Fig. 4(c)). The analytical
prediction agrees well with the MD simulation results summarized in Fig. 2(d).
Specifically, for cracks that align closely with $H<2$ nm, the overlap area
between the crack tips is small even if the amplitudes of relative in-plane
displacement and interlayer shear are large. As a result, the attraction is
weak, and no deflection is observed in the MD simulations. For $H>2$ nm, both
$K_{\rm I}^{\rm int}~{}(<0)$ and $K_{\rm II}^{\rm int}~{}(>0)$ decrease with
$H$, so does the ratio between overall SIFs, $\beta$, which results in crack
deflection with an angle diminishing as $H$ increases. In addition, the
normalized value of ${\rm tan}(90^{\circ}-\theta_{\rm c})$ decreases with
$G_{\rm int}$ in a logarithmic form, in consistency with the MD simulation
results (Fig. 4(d)).

Figure 4:  Theoretical prediction of crack interaction. (a) The shear-lag
model captures the load transfer between neighboring layers. (b) Calculations
of the stress intensity factors (SIFs, $K_{\rm I},K_{\rm II}$) under a
concentrated force ($P,Q$) at $z_{0}=x_{0}+iy_{0}$ measured from the crack
tip. The incomplete curve around the crack tip shows the distribution of
$G(\theta)$ under $K_{\rm I},K_{\rm II}(\neq 0)$. The direction of crack
propagation is determined by $\theta_{\rm c}$ = ${\rm max}G(\theta)$ according
to the criterion of maximum energy release rate (MERR). (c) The strength of
attraction, ${\rm tan}(90^{\circ}-\theta_{\rm c})$, plotted as a function of
$H$. (d) The normalized strength of attraction, ${\rm
tan}(90^{\circ}-\theta_{\rm c})/{\rm tan}\left[90^{\circ}-\theta_{\rm
c}(G_{0})\right]$, plotted as a function of $G_{\rm int}/G_{0}$ and compared
to the MD simulation results.

The interlayer load transfer by shear between 2D crystals depends on the twist
angle. MD simulations with a twist angle of $\theta=12^{\circ}$ show no
obvious cracks deflection from their intrinsic directions of propagation
defined by the loading condition and the crystal orientation-dependent
fracture toughness [24] (Figs. 5(a) and 5(b)). Compared to the bilayers in AB
stacking (Fig. 2(b)), shear strain at the crack tip in the twisted bilayer is
symmetric about the crack direction (Fig. 5(c)), suggesting that interlayer
sliding has negligible contribution to the in-plane lattice distortion. The
effect of twist angles on the fracture of bilayers can be further discussed by
tuning the shear stiffness $G_{\rm int}$ in the shear lag model. For graphene,
the shear stiffness of twisted bilayers is significantly lower than that in
AB-bilayers [29], and the mode-mixing factor $\beta=K_{\rm II}/K_{\rm I}$ is
negligible. The loading conditions and intrinsic anisotropic fracture
toughness of 2D crystals determine the crack direction. For MoS2 with a higher
interlayer shear modulus of $G_{\rm int}=16.4$ GPa [30] and lower in-plane
modulus of $E=270\pm 100$ GPa [31], the shear responses were reported to be
independent of the twist angle, possibly due to the low in-plane stiffness
that allows more significant in-plane lattice relaxation favoring local
registry match [32].

Our discussion also extends to crack configurations with different crack
lengths and directions of propagation in 2D bilayers. For parallel cracks with
different lengths in bilayer AB-stacking graphene, attraction is prevalent in
parallel cracks (with the same direction of propagation, Fig. S3 [20]). The
shorter crack always propagates ahead of the longer one, and the two cracks
chase each other and take the lead alternately since the critical SIF is
$K_{\rm Ic}\sim\sigma_{y}\sqrt{\pi a_{0}}$, where $\sigma_{y}$ is lower for
short cracks. For anti-parallel or ‘ _En-Passant_ ’ cracks in the bilayers,
the cracks avoid encountering each other during propagating, demonstrating
repulsion between the cracks (Fig. S4 [20]). The repulsive strength relies on
the distances between the two crack tips (Fig. 5(d)), where $H$ and $S$ are
the initial distances between the crack tips in the $x$\- and $y$-directions,
respectively. The deflection angle $\theta_{\rm c}$ is measured from MD
simulations, showing the repulsion becomes stronger as the cracks approach
each other firstly and then vanishes as the values of $H$ or $S$ become close
to $0$ (Fig. 5(e)). Moreover, the interlayer slips calculated from MD
simulations and continuum theory analysis show opposite directions of
$\Delta{\bf u}$ around the crack tips, in contrast to the parallel cracks
(Figs. 3 and S5 [20]). By integrating the SIFs over the crack-tip domain, we
find that $\beta<0$ at the crack tip is preserved. This result validates the
repulsion between anti-parallel cracks in the bilayer and its enhancement as
$H$ and $S$ decrease (Fig. 5(f)).

Figure 5:  Interaction between twisted graphene bilayer cracks and _En-
Passant_ cracks in a bilayer. (a, b) Fracture patterns of parallel cracks
obtained from MD simulations with/without interlayer interaction ($H=50$ nm).
The insets show the geometries, where the angles between the zigzag directions
and the $x$-axis in TL and BL are $24^{\circ}$ and $-24^{\circ}$,
respectively. (c) Shear strain at the crack tip in TL of the twisted bilayer
(corresponding to the simulations for panel A, with vdW interaction), which is
symmetric about the crack direction. (d) Fracture patterns of anti-parallel
cracks advancing toward each other. Here $H$ and $S$ are the initial distances
between the crack tips in the $x$\- and $y$-directions, respectively. (e) The
value of ${\rm tan}\theta_{\rm c}$ calculated from the fracture patterns as
functions of $H$ and $S$, showing repulsion in _En-Passant_ cracks (${\rm
tan}\theta_{\rm c}>0$). (f) Theoretical predictions of the mode-mixing factor
$\beta=K_{\rm II}/K_{\rm I}$ at the crack tip. The negative values of $\beta$
suggest anti-clockwise deflection, in consistency with the MD simulation
results.

For comparison, the interaction between co-planar cracks is also studied using
finite element analysis (FEA, Fig. S6 [20]), which can be inferred from the
mode-mixing factor $\beta=K_{\rm II}/K_{\rm I}$. Positive or negative $K_{\rm
II}/K_{\rm I}$ values result in clockwise or anti-clockwise crack deflection,
and thus attraction or repulsion, respectively. For parallel cracks with
identical lengths, the cracks tend to repel each other (Figs. S6(a) and S6(b)
[20]), while for those with length contrast, the interaction depends on their
crack-tip distances (Figs. S6(c) and S6(d) [20]). The longer crack advances by
experiencing repulsion and slight attraction by the neighboring shorter one as
the longer crack extends, while the shorter one tends not to advance due to
stress shielding. For _En-Passant_ edge cracks, the calculated SIFs show that,
when the cracks extend towards each other, the shear component $K_{\rm II}$
transits from negative to positive after the two crack tips encounter (Figs.
S6(e) and S6(f) [20]), explaining the attraction-after-repulsion behaviors.
Co-planar cracks interact by experiencing crack-tip stress shielding or
amplification that depends on the crack configurations. The cracks curve as a
result to minimize the shear stress and maximize the energy release rate [9].
However, for cracks in neighboring layers of a bilayer, the direction of
interlayer slip depends on the crack directions only. The in-plane shear
around the crack tips is always in the same direction, and as a result, only
attraction (repulsion) is observed for parallel (anti-parallel) cracks.

In summary, we demonstrate the control of crack paths through interaction
across neighbouring 2D layers by atomic-scale _in situ_ electron microscopy
and molecular dynamics simulations. The interaction that defines the crack
paths is attributed to interlayer shear resulting from non-uniform deformation
around the crack tips. The phenomenon is a direct consequence of the layered
anisotropy of 2D crystals and is well predicted by continuum mechanics
analysis by incorporating the linear elasticity of fracture mechanics and the
shear-lag model. The findings can guide the material design for engineered
fracture patterns and modified crack resistance. For example, by harnessing
the attraction between parallel cracks in neighboring layers, the ‘merging’ of
precracks can improve fracture resistance, provided that the interlayer
interaction can be strengthened.

###### Acknowledgements.

We wish to acknowledge the financial support by the National Natural Science
Foundation of China through grants 11825203, 52090032, and 52173230, the
National Key Basic Research Program of China grant No. 2022YFA1205400, the
China Postdoctoral Science Foundation through the grant 2023M731967, the Hong
Kong Research Grant Council General Research Fund through grants 15302522,
15301623. The computation was performed on the Explorer 1000 cluster system of
the Tsinghua National Laboratory for Information Science and Technology. S.F.
and X.Z. contributed equally to this work.

## References

  * Lawn  _et al._ [1977] B. R. Lawn, T. Wilshaw, and J. Rice,  _Fracture of Brittle Solids_ (Cambridge Solid State Science Series, 1977).
  * Feng  _et al._ [2023] S. Feng, W. Wang, S. Wang, X. Cui, Y. Yang, F. Xu, L. Liu, and Z. Xu, Controlling and visualizing fracture of 2D crystals by wrinkling, J. Mech. Phys. Solids 174, 105253 (2023).
  * Griffith [1921] A. A. Griffith, VI. The phenomena of rupture and flow in solids, Phil. Trans. Roy. Soc. A 221, 163 (1921).
  * Goehring  _et al._ [2015] L. Goehring, A. Nakahara, T. Dutta, S. Kitsunezaki, and S. Tarafdar,  _Desiccation Cracks and Their Patterns: Formation and Modelling in Science and Nature_ (John Wiley & Sons, 2015).
  * Quinn  _et al._ [2021] G. D. Quinn, J. J. Swab, and P. Patel, On radial, circumferential, and spiral cracks in fractured glass plates, Int. J. Appl. Glass Sci. 12, 189 (2021).
  * Broberg [1999] K. B. Broberg,  _Cracks and Fracture_ (Elsevier, 1999).
  * Ramesh and Simha [2008] M. Ramesh and K. Simha, Cylindrical crack vis a vis parallel crack pair: Shielding, energetics and asymptotics, Int. J. Fract. 152, 75 (2008).
  * Kranz [1979] R. L. Kranz, Crack-crack and crack-pore interactions in stressed granite, Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 16, 37 (1979).
  * Dalbe  _et al._ [2015] M.-J. Dalbe, J. Koivisto, L. Vanel, A. Miksic, O. Ramos, M. Alava, and S. Santucci, Repulsion and attraction between a pair of cracks in a plastic sheet, Phys. Rev. Lett. 114, 205501 (2015).
  * Schwaab  _et al._ [2018] M.-É. Schwaab, T. Biben, S. Santucci, A. Gravouil, and L. Vanel, Interacting cracks obey a multiscale attractive to repulsive transition, Phys. Rev. Lett. 120, 255501 (2018).
  * Xia and Hutchinson [2000] Z. C. Xia and J. W. Hutchinson, Crack patterns in thin films, J. Mech. Phys. Solids 48, 1107 (2000).
  * Ghelichi and Kamrin [2015] R. Ghelichi and K. Kamrin, Modeling growth paths of interacting crack pairs in elastic media, Soft Matter 11, 7995 (2015).
  * Novoselov  _et al._ [2016] K. Novoselov, A. Mishchenko, A. Carvalho, and A. Castro Neto, 2D materials and van der waals heterostructures, Science 353, aac9439 (2016).
  * Wang and Zheng [2007] L.-F. Wang and Q.-S. Zheng, Extreme anisotropy of graphite and single-walled carbon nanotube bundles, Appl. Surf. Sci. 90, 153113 (2007).
  * Feng and Xu [2021a] S. Feng and Z. Xu, Pattern development and control of strained solitons in graphene bilayers, Nano Lett. 21, 1772 (2021a).
  * Jang  _et al._ [2017] B. Jang, B. Kim, J.-H. Kim, H.-J. Lee, T. Sumigawa, and T. Kitamura, Asynchronous cracking with dissimilar paths in multilayer graphene, Nanoscale 9, 17325 (2017).
  * Jung  _et al._ [2018] G. S. Jung, S. Wang, Z. Qin, F. J. Martin-Martinez, J. H. Warner, and M. J. Buehler, Interlocking friction governs the mechanical fracture of bilayer MoS2, ACS Nano 12, 3600 (2018).
  * Han  _et al._ [2023] W. Han, X. Zheng, K. Yang, C. S. Tsang, F. Zheng, L. W. Wong, K. H. Lai, T. Yang, Q. Wei, M. Li, _et al._ , Phase-controllable large-area two-dimensional In2Se3 and ferroelectric heterophase junction, Nat. Nanotechnol. 18, 55 (2023).
  * Feng and Xu [2021b] S. Feng and Z. Xu, Strain characterization in two-dimensional crystals, Materials 14, 4460 (2021b).
  * [20] See the Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.000.000000 for details.
  * Liu  _et al._ [2012] Y. Liu, B. Xie, Z. Zhang, Q. Zheng, and Z. Xu, Mechanical properties of graphene papers, J. Mech. Phys. Solids 60, 591 (2012).
  * Liu  _et al._ [2010] Y. Liu, A. Dobrinsky, and B. I. Yakobson, Graphene edge from armchair to zigzag: The origins of nanotube chirality?, Phys. Rev. Lett. 105, 235502 (2010).
  * Kim  _et al._ [2012] K. Kim, V. I. Artyukhov, W. Regan, Y. Liu, M. F. Crommie, B. I. Yakobson, and A. Zettl, Ripping graphene: Preferred directions, Nano Lett. 12, 293 (2012).
  * Shi  _et al._ [2023] P. Shi, S. Feng, and Z. Xu, Non-equilibrium nature of fracture determines the crack paths, arXiv:2307.16126 (2023).
  * Feng  _et al._ [2022] S. Feng, K. Cao, Y. Gao, Y. Han, Z. Liu, Y. Lu, and Z. Xu, Experimentally measuring weak fracture toughness anisotropy in graphene, Commun. Mater. 3, 28 (2022).
  * Ouyang  _et al._ [2018] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nanoserpents: Graphene nanoribbon motion on two-dimensional hexagonal materials, Nano Lett. 18, 6009 (2018).
  * Zhang  _et al._ [2022] S. Zhang, Q. Xu, Y. Hou, A. Song, Y. Ma, L. Gao, M. Zhu, T. Ma, L. Liu, X.-Q. Feng, and Q. Li, Domino-like stacking order switching in twisted monolayer–multilayer graphene, Nat. Mater. 21, 621 (2022).
  * Yoo  _et al._ [2019] H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, K. Watanabe, G.-C. Yi, M. Kim, M. Luskin, E. B. Tadmor, E. Kaxiras, and P. Kim, Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene, Nat. Mater. 18, 448 (2019).
  * Hod  _et al._ [2018] O. Hod, E. Meyer, Q. Zheng, and M. Urbakh, Structural superlubricity and ultralow friction across the length scales, Nature 563, 485 (2018).
  * Zhao  _et al._ [2013] Y. Zhao, X. Luo, H. Li, J. Zhang, P. T. Araujo, C. K. Gan, J. Wu, H. Zhang, S. Y. Quek, M. S. Dresselhaus, _et al._ , Interlayer breathing and shear modes in few-trilayer MoS2 and WSe2, Nano Lett. 13, 1007 (2013).
  * Bertolazzi  _et al._ [2011] S. Bertolazzi, J. Brivio, and A. Kis, Stretching and breaking of ultrathin MoS2, ACS Nano 5, 9703 (2011).
  * Sun  _et al._ [2022] Y. Sun, Y. Wang, E. Wang, B. Wang, H. Zhao, Y. Zeng, Q. Zhang, Y. Wu, L. Gu, X. Li, and K. Liu, Determining the interlayer shearing in twisted bilayer MoS2 by nanoindentation, Nat. Commun. 13, 3898 (2022).


#  Vacuum zero point energy of self-interacting quantum fields in dS
background

Hassan Firouzjahia<EMAIL_ADDRESS>Haidar Sheikhahmadia
h.sh.ahmadi<EMAIL_ADDRESS>a School of Astronomy, Institute
for Research in Fundamental Sciences (IPM), P. O. Box 19395-5531, Tehran, Iran

###### Abstract

We consider self-interacting scalar fields with a conformal coupling in the dS
background and study the quantum corrections from bubble loop diagrams.
Incorporating the perturbative in-in formalism, we calculate the quantum
corrections in the vacuum zero point energy and pressure of self-interacting
fields with the potential $V\propto\Phi^{n}$ for even values of $n$. We
calculate the equation of state corresponding to these quantum corrections and
examine the scaling of the divergent terms in the vacuum zero point energy and
pressure associated to the dimensional regularization scheme. In particular,
we show that for quartic self-interacting scalar field the conformal
invariance is respected at two-loop order at the conformal point.

##  I introduction

Quantum field theory in dS background is a rich topic which has important
implications both observationally and theoretically DeWitt:1975ys ;
Birrell:1982ix ; Fulling:1989nb ; Mukhanov:2007zz ; Parker:2009uva .
Observationally, there are strong evidences that the early universe
experienced a period of inflation in which the background was nearly a dS
spacetime. In the simplest realization, inflation is driven by a scalar field,
the inflaton field, with a near flat potential Weinberg:2008zzc ;
Baumann:2022mni . While the inflaton field slowly rolls on top of its
potential, its quantum fluctuations are stretched on superhorizon scales,
which provide the seeds of large scale structure and the perturbations on CMB
Kodama:1984ziu ; Mukhanov:1990me . The basic predictions of the models of
inflation are that these primordial perturbations are nearly scale invariant,
Gaussian and adiabatic which are well supported by cosmological observations
Planck:2018vyg ; Planck:2018jri . In addition, numerous observations indicate
that the late universe is undergoing a phase of accelerating expansion. The
origin of dark energy as the source of the recent cosmological acceleration is
not known but a cosmological constant associated with the quantum zero point
energy of fields may be a good option Weinberg:1988cp ; Sahni:1999gb ;
Peebles:2002gy ; Copeland:2006wr ; Martin:2012bt .

On the theoretical sides, while there is no compelling theory of quantum
gravity at hand, but understanding quantum fields in curved backgrounds,
including the dS background, may shed some lights for the pursuit of a theory
of quantum gravity. Understanding important issues such as regularizations and
renormalization of cosmological correlations of quantum perturbations in dS
background can play important roles as well. More specifically, similar to
quantum field theories in flat spacetime, physical quantities such as the
energy momentum tensor, energy density and pressure suffer from infinities in
a curved spacetime. Therefore, it is an important question as how one can
regularize and renormalize the infinities to extract the finite physical
quantities. Furthermore, the fact that there is no unique vacuum in a curved
spacetime adds more complexities for the treatment of regularization and
renormalization in a curved spacetime Hawking:1975vcx ; Unruh:1976db ;
Unruh:1983ms ; Jacobson:2003vx ; Cozzella:2020gci ; Firouzjahi:2022rtn .

In this work we study the quantum fluctuations of a self-interacting scalar
field with non-minimal coupling to gravity in a dS background. The free
quantum fields with the conformal coupling in dS background Bunch:1978yq ;
Bunch:1978yw ; Avis:1977yn ; Ford:1977in ; Anderson:1985cw ; Anderson:1987yt ;
Buchbinder:1987fg ; Odintsov:1988wz ; Odintsov:1990qq ; Armendariz-
Picon:2023gyl ; Zhang:2019urk ; Ye:2022tgs ; Akhmedov:2022whm ;
Akhmedov:2020qxd ; Akhmedov:2019esv  and the self-interacting scalar field
with a quartic potential $V\propto\lambda\Phi^{4}$ Onemli:2002hr ;
Onemli:2004mb ; Brunier:2004sb ; Miao:2005am ; Prokopec:2008gw ;
Janssen:2008px ; Kahya:2009sz ; Miao:2010vs ; Glavan:2021adm ; Glavan:2020gal
; Woodard:2023vcw ; Beneke:2023wmt  are vastly studied in the literature. In
this paper we extend those works to more general self-interacting potentials
with the emphasize on vacuum zero point energy associated to the bubble
diagrams. More specifically, we study potentials of the form
$V\propto\lambda\Phi^{n}$ for even values of $n=4,6,...$ and calculate the
expectation values of the vacuum zero point energy and pressure associated to
the bubble diagrams. We employ dimensional regularization scheme
tHooft:1972tcz ; tHooft:1974toh ; Bollini:1972ui ; Deser:1974cz ;
Dowker:1975tf ; Barvinsky:1985an  to regularize the quantum infinities (for
earlier works concerning dimensional regularization scheme in dS spacetime see
for example Onemli:2002hr ; Onemli:2004mb ; Brunier:2004sb ; Miao:2005am ;
Prokopec:2008gw ; Janssen:2008px ; Kahya:2009sz ; Miao:2010vs ; Glavan:2021adm
; Glavan:2020gal ; Woodard:2023vcw ; Beneke:2023wmt ). This paper extends our
earlier work Firouzjahi:2023wbe  in which the vacuum zero point energy
associated to the bubble diagram for a free field with a conformal coupling in
a dS background were studied. In the presence of self-interaction, we have to
consider multiple bubble loop diagrams as depicted in Fig. 1. For an even
value of $n$, we have to consider a bubble diagram with $\frac{n}{2}$ loops to
obtain the order $\lambda$ corrections in the expectation values of physical
quantities.

Figure 1:  The Feynman bubble diagrams for $\lambda\Phi^{n}$ theory at the
leading order. For $n=3$ one is dealing with a nested (double time) in-in
integral while for even values of $n=4,6,8...$ one has a single time in-in
integral. For each even value of $n$, one has to consider a diagram with
$\frac{n}{2}$ loops to calculate the order $\lambda$ corrections in
expectation values such as $\langle\rho\rangle_{\Omega}$.

##  II The Setup

We consider a real scalar field $\Phi$ in a dS spacetime which is non-
minimally coupled to gravity with the conformal coupling $\xi$ and the self-
interacting potential $V(\Phi)=\lambda\Phi^{n}/n$ in which $\lambda$ is the
self-interacting constant. As will become evident soon, we consider even
values of $n$ with $n=4,6,...\,$. For the special case $n=4$, the coupling
constant $\lambda$ is dimensionless while for higher values of $n$ it has the
dimension $M^{4-n}$. In addition, as we are interested in the conformal limit,
and also to simplify the analysis, we assume the field is massless, $m=0$.
However, as we shall see below, one can easily restore the mass in our
formalism, though the equations will become more complicated. In our analysis
below, we treat the contribution of the self-interaction as a perturbation to
the free theory and calculate the vacuum zero point energy and pressure to
first order in coupling constant $\lambda$.

With the above discussions in mind, the action involving the scalar field is
given by

$S=\int
d^{D}x\sqrt{-g_{{}_{D}}}\left({-\frac{1}{2}\xi\Phi^{2}R}-\frac{1}{2}\nabla^{\mu}\Phi\nabla_{\mu}\Phi-{\frac{\lambda}{n}\Phi^{n}}\right)\,,$
(1)

where $D$ refers to the dimension of the spacetime and $g_{{}_{D}}$ stands for
the determinant of the metric. Since we employ dimensional regularization to
regularize the quantum infinities, we keep the spacetime dimension general and
set $D=4-\epsilon$ with $\epsilon\rightarrow 0$ as in conventional dimensional
regularization approach. We work in the test field limit where the background
geometry is the solution of the Einstein field equations with no backreactions
from the scalar fields. In order for this approximation to be consistent, the
vacuum zero point energy and pressure associated with the fluctuations of
$\Phi$ should be much smaller than the corresponding background quantities.

In the absence of the self-interaction the theory is classically conformal
invariant in four dimensional spacetime when $\xi=\frac{1}{6}$. However, as it
is well-known, this classical symmetry is anomalous under quantum corrections.
Furthermore, the addition of the potential can break the conformal invariance
at the classical level since the coupling constant $\lambda$ induces scale
into the theory when $n\neq 4$. While the theory is still classically
conformal invariant for $n=4$, but the quantum corrections break conformal
invariance in this case as well. We consider the arbitrary values of even $n$
while the special case of $n=4$ was studied extensively in the literature, see
for example the works of Woodard and collaborators Onemli:2002hr ;
Onemli:2004mb ; Brunier:2004sb ; Miao:2005am ; Prokopec:2008gw ;
Janssen:2008px ; Kahya:2009sz ; Miao:2010vs ; Glavan:2021adm ; Glavan:2020gal
; Woodard:2023vcw ; Beneke:2023wmt .

The background spacetime has the form of the FLRW metric,

$ds^{2}=a(\tau)^{2}\big{(}-d\tau^{2}+d{\bf x}^{2}\big{)}\,,$ (2)

where $\tau$ is the conformal time related to cosmic time via
$dt=a(\tau)d\tau$ in which $a(\tau)$ is the scale factor. In our limit of a
fixed dS background, $aH\tau=-1$, in which $H$ is the Hubble expansion rate
associated to the dS background. Since the dS spacetime is maximally
symmetric, the Ricci tensor and Ricci scalar are given as follows,

$\displaystyle R_{\mu\nu}=(D-1)H^{2}g_{\mu\nu},\quad\quad R=D(D-1)H^{2}\,.$
(3)

The scalar field equation is given by

${\Box\Phi-\xi R\Phi-{\lambda}{\Phi}^{n-1}=0}\,.$ (4)

Note that to simplify the above field equation, we use the convention that the
coupling constant has the form $\lambda/n$ instead of $\lambda/n!$ which is
usually used in QFT textbooks.

To study the quantum fluctuations, we introduce the canonically normalized
field $\sigma(\tau)$

$\displaystyle\mathrm{\sigma(\tau)}\equiv a^{\frac{D-2}{2}}\Phi(\tau)\,,$ (5)

in terms of which the action takes the following diagonal form,

$\displaystyle S=\frac{1}{2}\int d\tau d^{D-1}{\bf
x}\left[\sigma^{\prime}(\tau)^{2}-(\nabla\sigma)^{2}+\left(\frac{(D-4)(D-2)}{4}\big{(}\frac{a^{\prime}}{a}\big{)}^{2}+\frac{D-2}{2}\frac{a^{\prime\prime}}{a}-({\xi
R}+{\frac{\lambda}{a^{D}}\sigma^{2}})a^{2}\right)\sigma^{2}\right]$ (6)

where a prime indicates the derivative with respect to the conformal time.

To quantize the field we expand it in terms of the creation and annihilation
operators in $D-1$ dimension Fourier space as follows,

$\sigma\left(x^{\mu}\right)=\int\frac{d^{D-1}\mathbf{k}}{(2\pi)^{\frac{(D-1)}{2}}}\left(\sigma_{k}(\tau)e^{i\mathbf{k}\cdot\mathbf{x}}a_{\mathbf{k}}+\sigma^{*}_{k}(\tau)e^{-i\mathbf{k}\cdot\mathbf{x}}a_{\mathbf{k}}^{\dagger}\right)\,,$
(7)

in which $\sigma_{k}(\tau)$ is the quantum mode function and $a_{\mathbf{k}}$
and $a_{\mathbf{k}}^{\dagger}$ satisfy the following commutation relation in
$D-1$ spatial dimension,

$\left[a_{\mathbf{k}},a_{\mathbf{k^{\prime}}}^{\dagger}\right]=\delta^{D-1}(\mathbf{k}-\mathbf{k^{\prime}})\,.$
(8)

Correspondingly, the equation of motion of the free mode function (with
$\lambda=0$) from the action (6) is given by

$\sigma_{k}^{\prime\prime}(\tau)+\left[k^{2}+\frac{1}{\tau^{2}}\Big{(}{{D(D-1)}}\xi-\frac{D(D-2)}{4}\Big{)}\right]\sigma_{k}(\tau)=0\,.$
(9)

Note that the above equation is similar to the Mukhanov-Sasaki equation
associated to the inflaton perturbations in an inflationary background
Mukhanov:1990me . Notice that for $\xi=\frac{1}{6}$ in $D=4$ the second term
in the big bracket vanishes and the mode function reduces to its simple flat
form. This corresponds to the conformal limit which we consider in our
analysis below. However, note that in a general $D$-dimensional spacetime, the
conformal limit corresponds to the special value

$\displaystyle\xi=\xi_{D}\equiv\frac{D-2}{4(D-1)}\,.$ (10)

Imposing the Bunch-Davies (Minkowski) initial conditions for the modes deep
inside the horizon, the solution of the mode function from Eq. (9) is given in
terms of the Hankel function,

$\Phi_{k}(\tau)=a^{\frac{{2-D}}{2}}{\sigma_{k}}(\tau)={(-H\tau)^{\frac{{D-1}}{2}}}{\left({\frac{\pi}{{4H}}}\right)^{\frac{1}{2}}}{{e^{\frac{i\pi}{2}(\nu+\frac{1}{2})}}}H_{\nu}^{(1)}(-k\tau){\mkern
1.0mu}\,,$ (11)

where

$\nu\equiv\frac{1}{2}{\mkern 1.0mu}\sqrt{(D-1)^{2}{-4D(D-1)\xi}}\,.$ (12)

From the above expression we see that $\nu$ can be either real or pure
imaginary depending on the values of $\xi$. In our limit of interest where
$\xi$ is near its conformal value, $\xi$ is real. In particular, for
$\xi=\frac{1}{6}$ with $D=4$, we obtain $\nu=\frac{1}{2}$. As the in-in
integrals become non-trivial, and in order to prevent complications associated
to an imaginary $\nu$ in the mode functions, in our analysis below we assume
that $\nu$ is real. This imposes the bound $0<\xi\leq\frac{D-1}{4D}$ which for
$D=4$ corresponds to $\xi\leq\frac{3}{16}$.

The energy momentum tensor in the presence of self-interaction is given by,

$\displaystyle{{T_{\mu\nu}}}$ $\displaystyle=$
$\displaystyle(1-2\xi){\partial_{\mu}}\Phi{\partial_{\nu}}\Phi+(2\xi-\frac{1}{2}){g_{\mu\nu}}{g^{\alpha\beta}}{\partial_{\alpha}}\Phi{\partial_{\beta}}\Phi$
$\displaystyle+$
$\displaystyle\xi({R_{\mu\nu}}-\frac{1}{2}{g_{\mu\nu}}R){\Phi^{2}}+2\xi({g_{\mu\nu}}\Phi\Box\Phi-\Phi{\nabla_{\nu}}{\nabla_{\mu}}\Phi){-\frac{\lambda}{n}{g_{\mu\nu}}{\Phi^{n}}}\,.$
(13)

Employing the field equation (4) one can eliminate $\Box\Phi$ and using Eq.
(3), $T_{\mu\nu}$ is further simplified to

$\displaystyle\begin{array}[]{l}{T_{\mu\nu}}={\partial_{\mu}}\Phi{\partial_{\nu}}\Phi+\frac{{{g_{\mu\nu}}}}{2}(4\xi-1){\partial^{\alpha}}\Phi{\partial_{\alpha}}\Phi\\\
\,\,\,\,\,\,\,\,\,\,+\frac{\xi}{2}(D-1)\Big{(}2+(4\xi-1)D\Big{)}{H^{2}}{g_{\mu\nu}}{\Phi^{2}}-\xi{\nabla_{\mu}}{\nabla_{\nu}}{\Phi^{2}}+{{g_{\mu\nu}}{\lambda}(2\xi-\frac{1}{n}){\Phi^{n}}}\,.\end{array}$
(16)

Similarly, the trace of the energy momentum-tensor $T\equiv T^{\mu}_{\mu}$ is
obtained to be

$\displaystyle\begin{array}[]{l}T=2\Big{(}(D-1)\xi+\frac{2-D}{4}\Big{)}\Big{(}\partial^{\alpha}\Phi\partial_{\alpha}\Phi+D(D-1)\xi
H^{2}\Phi^{2}\Big{)}+{{\lambda}\Big{(}2\xi(D-1)-\frac{D}{n}\Big{)}{\Phi^{n}}}\,.\end{array}$
(18)

The energy density $\rho=T_{00}$ is:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{(1+4\xi)}{2}\dot{\Phi}^{2}+\frac{(1-4\xi)}{2}\nabla^{i}\Phi\nabla_{i}\Phi+\frac{H^{2}}{2}\Big{[}(1-4\xi){{\big{(}D(D-1)\xi)}}-2(D-1)\xi\Big{]}\Phi^{2}$
(19) $\displaystyle-$
$\displaystyle\xi\nabla_{0}\nabla_{0}\Phi^{2}-\lambda(2\xi-\frac{1}{n})\Phi^{n}\,.$

Finally, the pressure $P$ is given by

$\displaystyle P=\frac{1}{{D-1}}{\bot^{\mu\nu}}{T_{\mu\nu}}\,,$ (20)

in which $\perp^{\mu\nu}\equiv g^{\mu\nu}+u^{\mu}u^{\nu}$ represents the
projection operator and $u^{\mu}=(1,0,0,0)$ is the comoving four velocity.
Consequently, we obtain

$\displaystyle P=\frac{1}{{D-1}}(T+\rho)\,.$ (21)

In our analysis below, we will be mainly interested in vacuum expectation
values such as $\langle\rho\rangle$, $\langle P\rangle$ and $\langle
T\rangle$. This was studied for a free theory with a non-zero mass in
Firouzjahi:2023wbe  and here we extend these analysis in the presence of the
self-interaction $\lambda\Phi^{n}$. However, it is important to note that the
expectation value is with respect to the full vacuum which we denote by
$|\Omega\rangle$ so $\langle\rho\rangle\equiv\langle\Omega|\rho|\Omega\rangle$
and so on. Because of the interaction term $\lambda\Phi^{n}$ the vacuum
$|\Omega\rangle$ is different than the vacuum associated to the free theory
which is denoted by $|0\rangle$. To prevent confusion we define
$\langle\Omega|\rho|\Omega\rangle\equiv\langle\rho\rangle_{\Omega}$ while
$\langle 0|\rho|0\rangle\equiv\langle\rho\rangle_{0}$ and so on.

##  III Dimensional Regularizations and In-In Formalism

In this section we calculate vacuum expectation values such as
$\langle\rho\rangle_{\Omega}$ using dimensional regularization scheme in $D$
dimension. As in Firouzjahi:2023wbe , let us define

$\displaystyle\rho_{1}\equiv\frac{1}{2}\dot{\Phi}^{2}\,,\quad\quad\rho_{2}\equiv\frac{1}{2}g^{ij}\nabla_{i}\Phi\nabla_{j}\Phi\,,\quad\quad\rho_{3}\equiv\frac{1}{2}H^{2}\Phi^{2}\,.$
(22)

Then,

$\displaystyle\langle\rho\rangle_{\Omega}$ $\displaystyle=$
$\displaystyle(1+4\xi)\langle\rho_{1}\rangle_{\Omega}+(1-4\xi)\langle\rho_{2}\rangle_{\Omega}+\Big{[}(1-4\xi)D(D-1)\xi-2(D-1)\xi\Big{]}\langle\rho_{3}\rangle_{\Omega}$
(23) $\displaystyle-$
$\displaystyle\xi\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{\Omega}-\lambda(2\xi-\frac{1}{n})\langle\Phi^{n}\rangle_{\Omega}\,.$

Note that the first three terms in the first line in Eq. (23) are formally the
same as in Firouzjahi:2023wbe  except that in Firouzjahi:2023wbe  the
expectation values were with respect to the vacuum of the free theory. The two
terms in second line in Eq. (23) are new. First, we have a term with the
specific $\lambda$ coupling. Second, the contribution
$\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{\Omega}$ which is non-trivial. In
the analysis of Firouzjahi:2023wbe  it was noticed that
$\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{0}=0$. This is because
$\langle\Phi^{2}\rangle_{0}$ is a constant so it is easy to understand that
$\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{0}=\nabla_{0}\nabla_{0}\langle\Phi^{2}\rangle_{0}=0$.
However, in the presence of the interaction, we notice that
$\nabla_{0}|\Omega\rangle\neq 0$ so one can not simply take the derivative
outside the expectation value, i.e.
$\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{\Omega}\neq\nabla_{0}\nabla_{0}\langle\Phi^{2}\rangle_{\Omega}$.

Out of the five contributions into $\langle\rho\rangle_{\Omega}$ in Eq. (23),
the last term $\lambda\langle\Phi^{n}\rangle_{\Omega}$ is the easiest term to
calculate. This is because it has a factor of $\lambda$ and since we are
interested in first order corrections in $\lambda$, we can simply assume the
vacuum in this case is the free vacuum and

$\displaystyle\lambda\langle\Phi^{n}\rangle_{\Omega}\simeq\lambda\langle\Phi^{n}\rangle_{0}\,+{\cal
O}(\lambda^{2})\,.$ (24)

Note that the assumption that $n$ is even was necessary to obtain the above
result to leading order in $\lambda$. For odd values of $n$, the linear term
in $\lambda$ vanishes and one has to go to higher orders of $\lambda$ to
calculate $\lambda\langle\Phi^{n}\rangle_{\Omega}$. This brings additional
complexities involving in-in integral as we shall see in next section.

Since $\Phi$ is a Gaussian free field in the absence of interaction, one can
easily see that for even values of $n$

$\displaystyle\langle\Phi^{n}\rangle_{0}\simeq(n-1)!!\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{\frac{n}{2}}\,,$
(25)

in which $(n-1)!!=(n-1)(n-3)...1\,$. For example, for $n=4$, we have
$\langle\Phi^{4}\rangle_{0}=3\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{2}$. As
a result, the last term in Eq. (23) reads

$\displaystyle\lambda(2\xi-\frac{1}{n})\langle\Phi^{n}\rangle_{\Omega}\simeq(2\xi-\frac{1}{n})\lambda\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{\frac{n}{2}}+{\cal
O}(\lambda^{2})\,.$ (26)

The quantity $\langle\Phi^{2}\rangle_{0}$ can be viewed as the coincident
limit of the Feynman propagator. It plays important roles in our analysis
below in which the expectation values of the physical quantities in the
presence of interaction can be expressed in terms of
$\langle\Phi^{2}\rangle_{0}$.

###  III.1 Free Theory

Here we briefly review the results of Firouzjahi:2023wbe  for the free theory
which will be used in our following analysis as well. In the free theory with
$\lambda=0$, the vacuum $|0\rangle$ is dS invariant so from Eq. (23) we obtain

$\displaystyle\langle\rho\rangle_{0}=(1+4\xi)\langle\rho_{1}\rangle_{0}+(1-4\xi)\langle\rho_{2}\rangle_{0}+\Big{[}(1-4\xi)D(D-1)\xi-2(D-1)\xi\Big{]}\langle\rho_{3}\rangle_{0}\,.$
(27)

As we shall see below, all three components of $\langle\rho_{i}\rangle_{0}$
are expressed in terms of $\langle\Phi^{2}\rangle_{0}$ so let us calculate
this quantity. Performing the dimensional regularization analysis, we have

$\displaystyle\langle\Phi^{2}\rangle_{0}={\mu^{4-D}}\int\frac{d^{D-1}{\bf
k}}{(2\pi)^{D-1}}\left|\Phi_{k}(\tau)\right|^{2}\,,$ (28)

in which $\mu$ is a mass scale to keep track of the dimensionality of physical
quantities. We decompose the integral into the radial and angular parts as
follows

$\displaystyle{{\rm{d}}^{D-1}}{\bf{k}}={k^{D-2}}\;{\rm{d}}k\;{{{\rm{d}^{D-2}}}}\Omega{\mkern
1.0mu}\,,$ (29)

where ${\rm{d}^{D-2}}\Omega$ represents the $D-2$-dimensional angular part
with the volume

$\displaystyle\int\mathrm{d}^{D-2}\Omega=\frac{2\,\pi^{\frac{D-1}{2}}}{\Gamma\left(\frac{D-1}{2}\right)}\,.$
(30)

Combining all numerical factors and defining the dimensionless variable
$x\equiv-k\tau$ we obtain

$\displaystyle\langle\Phi^{2}\rangle_{0}={\frac{{\pi}^{\frac{3-D}{2}}{\mu}^{4-D}{H}^{D-2}}{{2}^{D}\Gamma\left(\frac{D-1}{2}\right)}}\int_{0}^{\infty}dx~{}x^{D-2}\left|H_{\nu}^{(1)}(x)\right|^{2}\,.$
(31)

Performing the integral, this yields 111We use the Maple computational
software to calculate the integrals.

$\displaystyle\langle\Phi^{2}\rangle_{0}=\frac{{{\mu^{4-D}}{\pi^{-\frac{D}{2}-1}}}}{2^{D}}\Gamma\Big{(}{\nu+\frac{D}{2}-\frac{1}{2}}\Big{)}\Gamma\Big{(}{-\nu+\frac{D}{2}-\frac{1}{2}}\Big{)}\Gamma\Big{(}{-\frac{D}{2}+1}\Big{)}\cos\big({\pi{\mkern
1.0mu}\nu}\big{missing}){H^{D-2}}\,.$ (32)

In particular, note that for the conformal limit with $\nu=\frac{1}{2}$, the
above expression vanishes. With $\langle\Phi^{2}\rangle_{0}$ at hand and using
Eq. (26), the last term for $\langle\rho\rangle_{\Omega}$ in Eq. (23) is
calculated accordingly.

Following similar steps to calculate $\langle\rho_{i}\rangle_{0}$ we obtain
the following relations Firouzjahi:2023wbe

$\displaystyle\langle\rho_{1}\rangle_{0}=(D-1)\xi\langle\rho_{3}\rangle_{0}\,,\quad\langle\rho_{2}\rangle_{0}=-(D-1)\langle\rho_{1}\rangle_{0}=-(D-1)^{2}\xi\langle\rho_{3}\rangle_{0}\,,$
(33)

where from Eq. (22),
$\langle\rho_{3}\rangle_{0}=\frac{H^{2}}{2}\langle\Phi^{2}\rangle_{0}$ with
$\langle\Phi^{2}\rangle_{0}$ given in Eq. (32). In the conformal limit where
$\langle\Phi^{2}\rangle_{0}=0$, we see that $\langle\rho_{i}\rangle_{0}=0$ and
correspondingly $\langle\rho\rangle_{0}=0$.

Similarly, calculating $\langle T\rangle_{0}$ one can show that $\langle
T\rangle_{0}=-D\langle\rho\rangle_{0}$, and correspondingly Firouzjahi:2023wbe

$\displaystyle\langle P\rangle_{0}=-\langle\rho\rangle_{0}\,.$ (34)

This is the expected result indicating the local Lorentz invariance for the
free theory in which one expects to locally have $\langle
T_{\mu\nu}\rangle=-\langle\rho\rangle g_{\mu\nu}$.

It is important to note that in the above results, $D$ is general so to
perform the regularization we consider $D=4-\epsilon$. One sets
$\epsilon\rightarrow 0$ at the end with the understanding that the singular
terms in physical quantities with inverse powers of $\epsilon$ are cancelled
by appropriate counter terms as in standard QFT analysis.

It is useful to look at the results in some limits of interests. In the case
of massless field with no conformal limit, $\xi=0$, one obtains

$\displaystyle\langle\rho\rangle_{0}^{\mathrm{reg}}=\frac{3H^{4}}{32\pi^{2}}=-\frac{1}{4}\langle
T\rangle_{0}^{\mathrm{reg}}\,,\quad\quad\quad(\xi=0)\,.$ (35)

On the other hand, as we noticed before, for the special case of conformal
point with $\xi=\frac{1}{6}$, $\langle\rho\rangle_{0}^{\mathrm{reg}}=\langle
T\rangle_{0}^{\mathrm{reg}}=0$. However, if one restores the mass so the
theory is not conformally invariant, one obtains Firouzjahi:2023wbe

$\displaystyle\langle\rho\rangle_{0}^{\mathrm{reg}}=-\frac{H^{4}}{96\pi^{2}}\beta^{2}+\frac{H^{4}}{64\pi^{2}}\Big{[}\ln\Big(\frac{H^{2}}{4\pi\mu^{2}}\Big{missing})-\frac{1}{2}\Big{]}\beta^{4}+{\cal
O}(\beta^{6})\,,\quad\quad\quad(\xi=\frac{1}{6})\,,$ (36)

in which $\beta\equiv m/H$.

###  III.2 In-in formalism

To calculate the first four terms in Eq. (23) to first order in $\lambda$, we
need to implement the in-in formalism which perturbatively relates the vacuum
expectation values of the interacting theory to the vacuum expectation of the
free theory as follows Weinberg:2005vy

$\langle{\cal
O}(\tau)\rangle_{\Omega}=\Big{\langle}0\Big{|}\bar{T}e^{i\int_{\tau_{0}}^{\tau}H_{I}\left(\tau^{\prime}\right)d\tau^{\prime}}{\cal
O}_{I}(\tau)Te^{-i\int_{\tau_{0}}^{\tau}H_{I}\left(\tau^{\prime}\right)d\tau^{\prime}}\Big{|}0\Big{\rangle}\,,$
(37)

where $T$ and $\bar{T}$ stand for time ordering and anti-time ordering
respectively. The subscript $I$ in the right hand side of the above equation
indicates that all quantities are calculated in the interaction picture, i.e.
with the mode function of the free theory given by Eq. (11). The initial time
is $\tau_{0}=-\infty$ while the final slicing $\tau$ is the time when the
measurement on the quantum operator ${\cal O}$ is made. As we work in an
unperturbed dS background, we have $-\infty<\tau^{\prime}\leq\tau<0$. Since
the time integrals in Eq. (37) are non-trivial, we shall restrict ourselves to
the case where the upper limit $\tau\rightarrow 0$, i.e. the measurement is
being made towards the future boundary of dS. In an inflationary setup with
deviations from an exact dS background, the final slicing $\tau\rightarrow 0$
corresponds to the time of end of inflation. Finally, $H_{I}$ represents the
interacting Hamiltonian, which in our case is

$\displaystyle H_{I}=\frac{\lambda}{n}a^{D}\int d^{{D-1}}{\bf
x}\,\Phi(x)^{n}\,.$ (38)

Note that the factor $a^{D}$ appears because of the volume element
$\sqrt{-g}$.

To leading order in $\lambda$, the correction in $\langle{\cal
O}(\tau)\rangle_{\Omega}$ is given by

$\displaystyle\langle{\cal O}(\tau)\rangle_{\Omega}=\langle{\cal
O}(\tau)\rangle_{0}+2\mathrm{Im}\int_{-\infty}^{\tau}d\tau^{\prime}\big{\langle}{\cal
O}(\tau)H_{I}(\tau^{\prime})\big{\rangle}_{0}\,.$ (39)

The first term above is the contribution in the absence of interaction which
were calculated in Firouzjahi:2023wbe  for ${\cal O}=\rho,P$ etc. Our goal
below is to calculate the above integrals for four different operators
$\rho_{1},\rho_{2},\rho_{3}$ and $\nabla_{0}\nabla_{0}\Phi^{2}$ as appearing
in Eq. (23). In order to isolate the contribution of the free theory, we
denote the last term in Eq. (39) by ${\Delta}\langle{\cal O}(\tau)\rangle$.

Let us start with $\langle{\rho_{3}}(\tau)\rangle_{\Omega}$, yielding

$\displaystyle{\Delta}\langle\rho_{3}(\tau,{\bf x_{0}})\rangle$
$\displaystyle=$
$\displaystyle\frac{H^{2}}{2}{\Delta}\langle\Phi^{2}(\tau,{\bf x_{0}})\rangle$
(40) $\displaystyle=$ $\displaystyle\frac{\lambda
H^{2}}{n}\imaginary[\int_{-\infty}^{\tau}d\tau^{\prime}a(\tau^{\prime})^{D}\int
d^{{D-1}}x\big{\langle}\Phi(\tau,{\bf x_{0}})^{2}\Phi(\tau^{\prime},{\bf
x})^{n}\big{\rangle}_{0}\Big{]}\,.$

Note that ${\bf x_{0}}$ is an arbitrary reference point in background where
the measurement is made. However, because of the spatial translation
invariance, we can set ${\bf x_{0}}=0$ so we do not specify ${\bf x_{0}}$ in
the rest of the analysis below. From the above expression we see that for odd
values of $n$, the expectation values vanish in the light of Wick theorem.
Therefore, for odd values of $n$ one needs to go to second order of
perturbations, leading to order ${\cal O}(\lambda^{2})$ corrections. This in
turn requires nested time integrals (i.e. a double time integral) which are
more complicated than the single time integral over $\tau^{\prime}$ which we
encounter for even value of $n$ as given in Eq. (40). For this reason, as
mentioned before, we restrict our analysis to even values of $n=4,6,...\,$.

There are two different types of contributions when performing the
contractions in Eq. (40). The first type is in the form $\langle\Phi(\tau,{\bf
x_{0}})^{2}\rangle_{0}\langle\Phi(\tau^{\prime},{\bf x})^{n}\rangle_{0}$. With
a bit of efforts one can show that this contribution has no imaginary
component so this contribution vanishes. The second type contracts each term
of $\Phi(\tau,{\bf x_{0}})$ with a term in $\Phi(\tau^{\prime},{\bf x})^{n}$.
There are total $n(n-1)!!$ possibilities for these contractions. After
performing the Wick contractions we obtain (for further details of Wick
contractions see Appendix A),

$\displaystyle\Delta\langle\rho_{3}(\tau)\rangle={(n-1)!!}\lambda
H^{2}\,\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{\frac{n-2}{2}}\,\mu^{4-D}I_{3}(\tau)\,,$
(41)

in which

$\displaystyle
I_{3}(\tau)\equiv\int_{\tau_{0}}^{\tau}d\tau^{\prime}a(\tau^{\prime})^{D}\int\frac{d^{D-1}\bf{q}}{(2\pi)^{(D-1)}}\mathrm{Im}\Big{[}\Phi_{q}(\tau)^{2}{\Phi_{q}^{*}}(\tau^{\prime})^{2}\Big{]}\,.$
(42)

Following similar steps for $\rho_{1}$ and $\rho_{2}$ we have

$\displaystyle\Delta\langle\rho_{1}(\tau)\rangle={(n-1)!!}\lambda
H^{2}\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{\frac{n-2}{2}}\,\mu^{4-D}I_{1}(\tau)\,,$
(43)

with

$\displaystyle
I_{1}(\tau)\equiv\int_{-\infty}^{\tau}d\tau^{\prime}{\frac{a(\tau^{\prime})^{D}}{a(\tau)^{2}}}\int\frac{d^{D-1}\bf{q}}{(2\pi)^{(D-1)}}\mathrm{Im}\Big{[}\Phi_{q}^{\prime}(\tau)^{2}{\Phi_{q}^{*}}(\tau^{\prime})^{2}\Big{]}\,,$
(44)

and

$\displaystyle\Delta\langle\rho_{2}(\tau)\rangle={(n-1)!!}\lambda
H^{2}\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{\frac{n-2}{2}}\mu^{4-D}I_{2}(\tau)\,,$
(45)

with

$\displaystyle
I_{2}(\tau)\equiv\int_{-\infty}^{\tau}d\tau^{\prime}{\frac{q^{2}a(\tau^{\prime})^{D}}{a(\tau)^{2}}}\int\frac{d^{D-1}\bf{q}}{(2\pi)^{(D-1)}}\mathrm{Im}\Big{[}\Phi_{q}(\tau)^{2}{\Phi_{q}^{*}}(\tau^{\prime})^{2}\Big{]}.$
(46)

In addition, we have to calculate
$\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{\Omega}$ as well, which is given
by

$\displaystyle\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{\Omega}=\frac{2}{a^{2}}\Big{[}\langle\Phi^{\prime
2}\rangle_{\Omega}+\langle\Phi\Phi^{\prime\prime}\rangle_{\Omega}+\frac{1}{\tau}\langle\Phi\Phi^{\prime}\rangle_{\Omega}\Big{]}\,.$
(47)

Calculating each term as above, we obtain

$\displaystyle\langle\nabla_{0}\nabla_{0}\Phi^{2}\rangle_{\Omega}={(n-1)!!}\lambda
H^{2}\big{(}\langle\Phi^{2}\rangle_{0}\big{)}^{\frac{n-2}{2}}\mu^{4-D}I_{4}(\tau)\,,$
(48)

in which

$\displaystyle
I_{4}(\tau)\equiv\int_{-\infty}^{\tau}d\tau^{\prime}{\frac{a(\tau^{\prime})^{D}}{a(\tau)^{2}}}\int\frac{d^{D-1}\bf{q}}{(2\pi)^{(D-1)}}\mathrm{Im}\Big{[}\Big{(}\Phi_{q}^{\prime}(\tau)^{2}+\Phi_{q}(\tau)\Phi^{\prime\prime}_{q}(\tau)+\frac{1}{\tau}\Phi_{q}(\tau)\Phi^{\prime}_{q}(\tau)\Big{)}{\Phi_{q}^{*}}(\tau^{\prime})^{2}\Big{]}.$
(49)

##  IV Measurements at Future dS Boundary

So far our analysis were general except that we assumed that $n$ is even so we
deal with a single time integral as in Eq. (42). To proceed further, we should
calculate each of $I_{i}(\tau)$ listed above. We start with $I_{3}(\tau)$
which is easier. Let us denote the $(D-2)$ dimensional angular part of the
momentum integral by $V_{D-2}$ as given in Eq. (30). Defining the
dimensionless variables $x\equiv-q\tau^{\prime}$ and $y\equiv-q\tau$ and
switching the orders of the time and momentum integrals, we obtain

$\displaystyle
I_{3}=\frac{V_{D-2}}{(2\pi)^{D-1}}\imaginary[\frac{\pi^{2}H^{D-4}}{16}\int_{0}^{\infty}dyy^{D-2}\big{(}H_{\nu}^{(1)}(y)\big{)}^{2}\int_{y}^{\infty}\frac{dx}{x}\big{(}H_{\nu}^{(2)}(x)\big{)}^{2}\Big{]}\,.$
(50)

Looking at the integral over the $x$ variable we see that it is in the form of
a nested integral. Furthermore, its integrand falls off quickly for large $x$
as the integrand oscillates rapidly with a decaying amplitude. Therefore, one
expects the dominant contribution for the interior integral comes from the
lower bound when $x\rightarrow y$.

To proceed further and to calculate the integrals analytically, we have to
impose some simplification approximations. As a reasonable approximation, we
take $\tau\rightarrow 0$. This corresponds to performing the measurement at
the future boundary of dS. In inflationary models, this corresponds to
performing the vacuum expectation value at the time of end of inflation. In
the limit $\tau\rightarrow 0$, the mode function $\Phi_{q}(\tau)$ in Eq. (11)
simplifies to

$\displaystyle\Phi_{q}(\tau)\rightarrow-\frac{i\Gamma(\nu)}{\pi}\big{(}\frac{\pi}{4H}\big{)}^{\frac{1}{2}}e^{\frac{i\pi}{2}(\nu+\frac{1}{2})}\big{(}-H\tau\big{)}^{\frac{{D}-1}{2}}\big{(}\frac{-2}{k\tau}\big{)}^{\nu}\,.$
(51)

In the context of inflation, this represents the superhorizon limit of
cosmological perturbations when $q\ll aH$ so $y=-q\tau\rightarrow 0$.

Expanding the interior integrand for $x\ll 1$ and taking the lower bound of
integral with $x\rightarrow y\rightarrow 0$, the rest of integral over $y$ can
be taken analytically yielding to

$\displaystyle
I_{3}=\frac{4^{\nu}{\pi}^{\frac{1}{2}}H^{D-4}}{32\nu}\frac{\Gamma(-\nu+\frac{d}{2}-\frac{1}{2})\Gamma(\nu-\frac{d}{2}+1)\Gamma(\frac{d}{2}-\frac{1}{2})}{\Gamma(2\nu-\frac{d}{2}+\frac{3}{2})\Gamma(1-\nu)^{2}}\sin(\pi\nu)\sin\big(\pi\nu-\frac{\pi
D}{2}\big{missing})\,.$ (52)

Following the same strategy, we can calculate $I_{1},I_{2}$ and $I_{4}$
analytically. It turns out that $I_{i}$ are related to each other so all of
them can be expressed in terms of $I_{3}$ as follows:

$\displaystyle I_{1}=\frac{(4\nu^{2}+d-2\nu-1)(-2\nu+d-1)}{4(d-2\nu)}I_{3}\,,$
(53) $\displaystyle I_{2}=-\frac{(d-1-4\nu)(d-1-2\nu)(d-1)}{4(d-2\nu)}I_{3}$
(54)

and

$\displaystyle I_{4}=8\nu^{2}I_{3}\,.$ (55)

With the above analytical values of the in-in integrals at hand, we can
calculate $\langle\rho\rangle_{\Omega},\langle P\rangle_{\Omega}$ and $\langle
T\rangle_{\Omega}$. The expressions for these quantities for a general value
of $\xi$ are too complicated to report here so we consider two spacial limits
for analytical presentations. First, the conformal point where
$\xi=\frac{1}{6}$ and second, the limit of small deviation from the conformal
point with $\delta\xi\equiv\xi-\frac{1}{6}\ll 1$. For general values of $\xi$
where the analytical results are intractable, we present the numerical plots
of $\langle\rho\rangle_{\Omega}$ and $\langle T\rangle_{\Omega}$.

###  IV.1 Conformal point: $\xi=\frac{1}{6}$

At the conformal point with $\xi=\frac{1}{6}$, we obtain

$\displaystyle\Delta\langle\rho\rangle=\frac{\lambda{(3-n)}}{3n}(n-1)!!\Big{(}\frac{-H^{2}}{24\pi^{2}}\Big{)}^{\frac{n}{2}}\,,$
(56)

and

$\displaystyle\Delta\langle
P\rangle=\frac{\lambda(2n-9)}{9n}(n-1)!!\Big{(}\frac{-H^{2}}{24\pi^{2}}\Big{)}^{\frac{n}{2}}\,.$
(57)

From the above formulas for $\Delta\langle\rho\rangle$ and $\Delta\langle
P\rangle$, and noting that when $\lambda=0$ both $\langle\rho\rangle$ and
$\langle P\rangle$ vanish, the equation of state $w=\frac{P}{\rho}$ is
obtained to be

$\displaystyle w=\frac{-2n+9}{3n-9}\,.$ (58)

For example for $n=4$, corresponding to two-loop quantum corrections, we
obtain $w=\frac{1}{3}$ so the quantum corrections in energy momentum tensor
behave like radiation. On the other hand, for large values of $n$, the
equation of state approaches to $w\rightarrow-\frac{2}{3}$. It is intriguing
that the quantum corrections from self-interactions are not in the form of
$w=-1$ which is expected from local Lorentz invariance.

It is instructive to calculate $\Delta\langle T\rangle$ as well, yielding

$\displaystyle\Delta\langle
T\rangle=\frac{\lambda(n-4)}{n}(n-1)!!\Big{(}\frac{-H^{2}}{24\pi^{2}}\Big{)}^{\frac{n}{2}}\,.$
(59)

Interestingly, for the case $n=4$, we see that the quantum corrections in the
trace of energy momentum tensor vanish. This is consistent with the fact that
for $n=4$ the parameter $\lambda$ is dimensionless so the theory is
classically scale invariant and $T=0$ at the classical level. It is intriguing
that the two-loop quantum correction respects this result as well. However, it
is an open question whether or not higher orders loops corrections (i.e.
$\lambda^{2}$ and higher orders corrections) respect this conclusion.

###  IV.2 Small deviation from conformal point

Now suppose we slightly deviate from the conformal point with
$\delta\xi=\xi-\frac{1}{6}\ll 1$. We calculate the quantum corrections to
leading order in $\delta\xi$. By increasing the value of $n$, the analysis
become complicated, so here we present the results for two cases $n=4$ and
$n=6$.

Starting with $n=4$ to linear order in $\big{(}\xi-\frac{1}{6}\big{)}$ we
obtain

$\displaystyle\Delta\langle\rho\rangle\simeq-\frac{\lambda
H^{4}}{2304\pi^{4}}+\frac{\lambda
H^{4}}{32\pi^{4}}\Big{[}-\frac{1}{\epsilon}+\ln\big(\frac{H^{2}}{4\pi\mu^{2}}\big{missing})+\frac{5}{2}-\gamma\Big{]}\big{(}\xi-\frac{1}{6}\big{)}+{\cal
O}\big{(}\big{(}\xi-\frac{1}{6}\big{)}^{2}\big{)}\,,\quad(n=4)$ (60)

and

$\displaystyle\Delta\langle P\rangle\simeq-\frac{\lambda
H^{4}}{6912\pi^{4}}+\frac{\lambda
H^{4}}{96\pi^{4}}\Big{[}-\frac{1}{\epsilon}+\ln\big(\frac{H^{2}}{4\pi\mu^{2}}\big{missing})+\frac{23}{6}-\gamma\Big{]}\big{(}\xi-\frac{1}{6})+{\cal
O}\big{(}\big{(}\xi-\frac{1}{6}\big{)}^{2}\big{)}\,,\quad(n=4)\,,$ (61)

in which $\gamma$ is the Euler number. We have the divergent $\epsilon^{-1}$
term in both $\Delta\langle\rho\rangle$ and $\Delta\langle P\rangle$ which
appears when $\xi\neq\frac{1}{6}$ which should be removed by appropriate
counter terms. Interestingly, in this case we see that the singular terms in
$\Delta\langle\rho\rangle$ and $\Delta\langle P\rangle$ have the equation of
state associated to radiation, $w=\frac{1}{3}$, as suggested in Eq. (58) while
this does not hold for the finite terms. Indeed, we have

$\displaystyle\Delta\langle\rho\rangle-3\Delta\langle
P\rangle\simeq\frac{19\lambda
H^{4}}{96\pi^{4}}\big{(}\xi-\frac{1}{6}\big{)}\quad\quad\quad(n=4)\,,$ (62)

so the divergent $\epsilon^{-1}$ terms are cancelled in the above expression.
In addition, the trace of energy momentum tensor is not zero,

$\displaystyle\Delta\langle T\rangle\simeq-\frac{3\lambda
H^{4}}{16\pi^{4}}\big{(}\xi-\frac{1}{6}\big{)}\quad\quad\quad(n=4)\,.$ (63)

However, we notice that it has no singular part.

Now consider the case $n=6$, corresponding to three-loop bubble diagrams. In
this case the coupling constant $\lambda$ has the dimension of $M^{-2}$ so the
quantum corrections would scale like $\lambda H^{6}$. To linear order in
$\big{(}\xi-\frac{1}{6}\big{)}$ one obtains

$\displaystyle\Delta\langle\rho\rangle\simeq\frac{5\lambda
H^{6}}{27648\pi^{6}}-\frac{15\lambda
H^{6}}{512\pi^{6}}\Big{[}-\frac{2}{3\epsilon}+\ln\big(\frac{H^{2}}{4\pi\mu^{2}}\big{missing})+\frac{71}{54}-\gamma\Big{]}\big{(}\xi-\frac{1}{6}\big{)}+{\cal
O}\big{(}\big{(}\xi-\frac{1}{6}\big{)}^{2}\big{)},\quad(n=6).$ (64)

As expected, we have the singular term $\epsilon^{-1}$. In addition, there
will be singular terms $\epsilon^{-2}$ but it comes at second order
$\big{(}\xi-\frac{1}{6}\big{)}^{2}$. Similarly, for the pressure we obtain

$\displaystyle\Delta\langle P\rangle\simeq\frac{-5\lambda
H^{6}}{82944\pi^{6}}+\frac{5\lambda
H^{6}}{512\pi^{6}}\Big{[}-\frac{2}{3\epsilon}+\ln\big(\frac{H^{2}}{4\pi\mu^{2}}\big{missing})+\frac{127}{54}-\gamma\Big{]}\big{(}\xi-\frac{1}{6}\big{)}+{\cal
O}\big{(}\big{(}\xi-\frac{1}{6}\big{)}^{2}\big{)}\,\quad(n=6).$ (65)

From Eq. (58) for $w$ in conformal point with $n=6$ we obtain
$w=-\frac{1}{3}$, so we expect that the singular terms in
$\Delta\langle\rho\rangle$ and $\Delta\langle P\rangle$ to be related with
this equation of state. Indeed, we have

$\displaystyle\Delta\langle\rho\rangle+3\Delta\langle
P\rangle\simeq\frac{35\lambda
H^{6}}{1152\pi^{6}}\big{(}\xi-\frac{1}{6}\big{)}\quad\quad\quad(n=6)\,,$ (66)

so the singular terms $\epsilon^{-1}$ are cancelled in the above expression.
We have checked that the equation of state $w=-\frac{1}{3}$ also holds for the
most singular terms $\epsilon^{-2}$ which appears at second order
$\big{(}\xi-\frac{1}{6}\big{)}^{2}$.

A conclusion is that the equation of state Eq. (58), which is obtained for the
conformal point, is the equation of state for the singular terms in inverse
powers of $\epsilon$ in $\langle\rho\rangle$ and $\langle P\rangle$ at each
order of $\big{(}\xi-\frac{1}{6}\big{)}$ as well.

###  IV.3 Numerical plots for general value of $\xi$

As the analytical expressions for $\Delta\langle\rho\rangle$ and
$\Delta\langle T\rangle$ for the general values of $\xi$ are very complicated,
here we present their numerical plots for the special case of $n=4$.

Figure 2: The 3D diagrams of $\Delta\langle\rho\rangle$ and $\Delta\langle
T\rangle$ as functions of variables $(\xi,\mu)$ for n=4. We have varied $\mu$
in units of $H$ while $\xi$ is varied in the interval $0<\xi<\frac{3}{16}$. In
the left panel $\Delta\langle\rho\rangle$ is measured in the scale of
$\frac{-\lambda H^{4}}{2304\pi^{4}}$ while in the right panel $\Delta\langle
T\rangle$ is measured in the scale of $-\frac{3\lambda H^{4}}{16\pi^{4}}$. The
green horizontal surface in the left panel represents the surface $z=1$ while
in the right panel it represents the surface $z=\xi-\frac{1}{6}$. In the left
panel near the conformal point $\xi=\frac{1}{6}$, $\Delta\langle\rho\rangle$
approaches the constant value given in Eq. (60) while in the right panel
$\Delta\langle T\rangle$ approaches the formula given in Eq. (63) the near the
conformal point.

In Fig. 2, we have presented the three-dimensional behaviour of
$\Delta\langle\rho\rangle$ and $\Delta\langle T\rangle$ as functions of two
parameters $(\xi,\mu)$. We have varied $\mu$ in units of $H$ while $\xi$ is
varied in the interval $0<\xi<\frac{3}{16}$ in which the index $\nu$ is real.
In the left panel of this figure, we have presented
$\Delta\langle\rho\rangle/(\frac{-\lambda H^{4}}{2304\pi^{4}})$ while in the
right panel we have presented $\Delta\langle T\rangle/(-\frac{3\lambda
H^{4}}{16\pi^{4}})$. In the left panel, the horizontal surface represents the
surface with the value equal to unity while in the right panel the horizontal
surface represents the surface $z=\xi-\frac{1}{6}$, independent of $\mu$. In
the left panel, we see that near the conformal point $\xi=\frac{1}{6}$,
$\Delta\langle\rho\rangle$ approaches the constant value given in Eq. (60). On
the other hand, in the right panel we see that near the conformal point
$\Delta\langle T\rangle$ approaches the formula given in Eq. (63) in which the
quantum corrections in trace vanish at the conformal point.

Figure 3: The diagram shows the behaviour of $\Delta\langle\rho\rangle$ and
$\Delta\langle T\rangle$ as functions of $\mu$ for various fixed values of
$\xi$. As in Fig. 2, $\Delta\langle\rho\rangle$ is measured in the scale of
$\frac{-\lambda H^{4}}{2304\pi^{4}}$ while $\Delta\langle T\rangle$ is
measured in the scale of $\frac{-3\lambda H^{4}}{16\pi^{4}}$. At the conformal
point $\xi=\frac{1}{6}$, $\Delta\langle\rho\rangle$ and $\Delta\langle
T\rangle$ are constant while the behaviours of the curves change for values of
$\xi$ below and above this value.

To have a better visualization, in Figs. 3 and 4 we have presented the two-
dimensional sections of the above plot as functions of $\mu$ ($\xi$) for fixed
values of $\xi$ ($\mu$). In both panels of Fig. 3 we see that
$\xi=\frac{1}{6}$ plays special role in which at this point
$\Delta\langle\rho\rangle=\frac{-\lambda H^{4}}{2304\pi^{4}}$ and
$\Delta\langle T\rangle=0$ independent of the value of $\mu$. This property is
reinforced in Fig. 4 in which all curves merge to each other at the conformal
point $\xi=\frac{1}{6}$.

Figure 4: The diagram shows the behaviour of $\Delta\langle\rho\rangle$ and
$\Delta\langle T\rangle$ as functions of $\mu$ for various fixed values of
$\xi$. As in previous plots $\Delta\langle\rho\rangle$ is measured in the
scale of $\frac{-\lambda H^{4}}{2304\pi^{4}}$ while $\Delta\langle T\rangle$
is measured in the scale of $\frac{-3\lambda H^{4}}{16\pi^{4}}$. At the
conformal point $\xi=\frac{1}{6}$, all curves merge to a fixed value in each
panel.

##  V Summary and Discussions

In this work we have studied the quantum fields with self-interaction
potential $V\propto\lambda\Phi^{n}$ in a dS background and calculated the
vacuum expectation values of the energy density, pressure and the trace of the
energy momentum tensor. We have employed the perturbative in-in formalism to
calculate the corrections in $\langle\rho\rangle_{\Omega}$ etc to first order
in $\lambda$. To simplify the analysis we have considered even values of $n$.
This is because for even values of $n$ the non-zero corrections appear at the
first order of $\lambda$ while for odd values of $n$ the non-zero corrections
appear at the order $\lambda^{2}$. Technically speaking, this corresponds to
having a single time in-in integral for even values of $n$ while for odd
values of $n$, we need to consider nested (double time) in-in integrals.

Our analysis were performed for general values of $\xi$ but to report the
analytical results we have considered two special limits, the conformal point
$\xi=\frac{1}{6}$ and the case with a small deviation from the conformal point
$\delta\xi\ll 1$. At the conformal point we obtain the equation of state Eq.
(58) between $\langle\rho\rangle_{\Omega}$ and $\langle P\rangle_{\Omega}$. In
particular, for $n=4$ we obtain the intriguing result that $w=\frac{1}{3}$
with $\langle T\rangle_{\Omega}=0$. This suggests that the classical conformal
invariance associated to the case $n=4$ is respected under two-loop (i.e.
order $\lambda$) quantum corrections. It is an open question if this symmetry
does hold to higher orders loops corrections. On the other hand, for other
values of $n$, the coupling $\lambda$ is dimensionful so the theory is not
conformally invariant even at the classical level. Therefore, the equation of
state is different that that of radiation. It approaches
$w\rightarrow-\frac{2}{3}$ for large values of $n$. In the case of
$\delta\xi\neq 0$ but being small, we have calculated
$\langle\rho\rangle_{\Omega}$ and $\langle P\rangle_{\Omega}$ to first order
in $\delta\xi$ for cases $n=4$ and $n=6$. Unlike the conformal point, we have
the divergent terms $\epsilon^{-1}$ at the order $\delta\xi$ and
$\epsilon^{-2}$ for $n=6$ at the order $\delta\xi^{2}$. We have found that
these singular terms between $\langle\rho\rangle_{\Omega}$ and $\langle
P\rangle_{\Omega}$ are related by the equation of state Eq. (58) at each order
of $\delta\xi$.

While we presented the analytical results only for the above two special
cases, but we have presented the numerical plots of $\Delta\langle\rho\rangle$
and $\Delta\langle T\rangle$ for general values of $\xi$ for $n=4$ in Figs. 2,
3 and 4. All these figures highlight the special roles played by the conformal
limit $\xi=\frac{1}{6}$.

There are a number of directions in which the current study can be extended.
The first which comes to mind is repeating this investigation for odd values
of $n$. In this case, the in-in analysis become more complicated as one has to
go to second order in in-in integrals, with corrections appearing at order
$\lambda^{2}$. The second question is to extend the current analysis for $n=4$
to second order in $\lambda$ and see if the classical conformal invariance
with $w=\frac{1}{3}$ and $\langle T\rangle_{\Omega}=0$ hold at order
$\lambda^{2}$ or not. Since we were mostly interested in conformal limits, we
have set the mass of the field to be zero. However, one can easily incorporate
the effects of mass in the current analysis as well. In particular, as a
physical example, one may consider a symmetry breaking potential like
$V=\lambda(\Phi^{2}-v^{2})^{2}$ in which the effective mass, the cubic and the
quartic couplings are all related via the coupling constant $\lambda$.
Employing in-in formalism to second order in $\lambda$, one can calculate the
quantum corrections in vacuum zero point energy from both the cubic and
quartic self-interactions. As the theory is not scale invariant even at the
classical level, then one may not expect that the conclusions $w=\frac{1}{3}$
and $\langle T\rangle_{\Omega}=0$ to hold under loop corrections. We would
like to come back to the interesting example of symmetry breaking potential in
future.

Acknowledgments: We are grateful to Richard Woodard and Ali Akbar Abolhasani
for useful discussions and correspondences. The work of H. F. is partially
supported by INSF of Iran under the grant number 4022911. H.S. expresses
gratitude to E. T. Akhmedov, D. V. Bykov, A. Bazarov, and I. Aref’eva for the
constructive discussions they had during conference fields&Strings 2024,
Moscow, and appreciates the warm hospitality of the Steklov Institute.

##  Appendix A Wick Contractions

In this Appendix we briefly review some steps concerning Wick contractions.

Let us consider Eq. (40) in configuration space,

$\Delta\langle{\rho_{3}}(\tau,{{\mathbf{x}}_{\mathbf{0}}})\rangle=\frac{{\lambda{H^{2}}}}{n}\operatorname{Im}\Big{[}\int_{-\infty}^{\tau}d\tau^{\prime}a{(\tau^{\prime})^{D}}\int{{d^{D-1}}}x{\langle
0|\Phi{(\tau,{{\mathbf{x}}_{\mathbf{0}}})^{2}}\Phi{(\tau^{\prime},{\mathbf{x}})^{n}}|0\rangle_{0}}\Big{]}\,.$
(67)

Going to Fourier space, and setting ${\mathbf{x}}_{\mathbf{0}}=0$ because of
the translation invariance, this yields

$\displaystyle\Delta\langle{\rho_{3}}(\tau,{{\mathbf{x}}_{\mathbf{0}}})\rangle$
$\displaystyle=$ $\displaystyle\langle
0|\frac{{{\mu^{(n/2)(4-D)}}}}{n}\lambda{H^{2}}{\mathrm{Im}}\Big{[}\int_{-\infty}^{\tau}d\tau^{\prime}a{(\tau^{\prime})^{D}}\int{\frac{{{d^{D-1}}{\mathbf{p}_{1}}}}{{{{(2\pi)}^{D-1}}}}}\frac{{{d^{D-1}}{\mathbf{p}_{2}}}}{{{{(2\pi)}^{D-1}}}}\frac{{{d^{D-1}}{\mathbf{q}_{1}}}}{{{{(2\pi)}^{D-1}}}}\cdots\frac{{{d^{D-1}}{\mathbf{q}_{n}}}}{{{{(2\pi)}^{D-1}}}}$
$\displaystyle\times$
$\displaystyle\int{{d^{D-1}}}x\left[{{a_{{{\mathbf{p}}_{1}}}}\phi_{{\mathbf{p}_{1}}}(\tau){}+a_{{{\mathbf{p}}_{1}}}^{\dagger}\phi_{{\mathbf{p}_{1}}}^{*}(\tau)}\right]\left[{{a_{{{\mathbf{p}}_{2}}}}{\phi_{{\mathbf{p}_{2}}}}(\tau)+a_{{{\mathbf{p}}_{2}}}^{\dagger}\phi_{{\mathbf{p}_{2}}}^{*}(\tau)}\right]$
$\displaystyle\times$
$\displaystyle\left[{{a_{{{\mathbf{q}}_{1}}}}\phi_{{\mathbf{q}_{1}}}(\tau^{\prime}){e^{i{{\mathbf{q}}_{1}}\cdot{\mathbf{x}}}}+a_{{{\mathbf{q}}_{1}}}^{\dagger}\phi_{{\mathbf{q}_{1}}}^{*}(\tau^{\prime}){e^{-i{{\mathbf{q}}_{1}}\cdot{\mathbf{x}}}}}\right]\left[{{a_{{{\mathbf{q}}_{2}}}}\phi_{{\mathbf{q}_{2}}}(\tau^{\prime}){e^{i{{\mathbf{q}}_{2}}\cdot{\mathbf{x}}}}+a_{{{\mathbf{q}}_{2}}}^{\dagger}\phi_{{\mathbf{q}_{2}}}^{*}(\tau^{\prime}){e^{-i{{\mathbf{q}}_{2}}\cdot{\mathbf{x}}}}}\right]$
$\displaystyle\times$
$\displaystyle\left[{{a_{{{\mathbf{q}}_{3}}}}\phi_{{\mathbf{q}_{3}}}(\tau^{\prime}){e^{i{{\mathbf{q}}_{3}}\cdot{\mathbf{x}}}}+a_{{{\mathbf{q}}_{3}}}^{\dagger}\phi_{{\mathbf{q}_{3}}}^{*}(\tau^{\prime}){e^{-i{{\mathbf{q}}_{3}}\cdot{\mathbf{x}}}}}\right]\left[{{a_{{{\mathbf{q}}_{4}}}}\phi_{{\mathbf{q}_{4}}}(\tau^{\prime}){e^{i{{\mathbf{q}}_{4}}\cdot{\mathbf{x}}}}+a_{{{\mathbf{q}}_{4}}}^{\dagger}\phi_{{\mathbf{q}_{4}}}^{*}(\tau^{\prime}){e^{-i{{\mathbf{q}}_{4}}\cdot{\mathbf{x}}}}}\right]$
$\displaystyle\times$
$\displaystyle\cdots\left[{{a_{{{\mathbf{q}}_{n}}}}\phi_{{\mathbf{q}_{n}}}(\tau^{\prime}){e^{i{{\mathbf{q}}_{n}}\cdot{\mathbf{x}}}}+a_{{{\mathbf{q}}_{n}}}^{\dagger}\phi_{{\mathbf{q}_{n}}}^{*}(\tau^{\prime}){e^{-i{{\mathbf{q}}_{n}}\cdot{\mathbf{x}}}}}\right]|0\rangle\Big{]}\,.$

As a specific example consider $n=6$ where we present the Wick contractions
for one case:

$\begin{gathered}\Delta\left\langle{{\rho_{3}}\left({\tau,{{\mathbf{x}}_{\mathbf{0}}}}\right)}\right\rangle\supset\langle
0\Big{|}\frac{{{\mu^{(3)(4-D)}}}}{6}\lambda{H^{2}}\operatorname{Im}\left[{\int_{-\infty}^{\tau}d{\tau^{\prime}}a{{\left({{\tau^{\prime}}}\right)}^{D}}\int{\frac{{{d^{D-1}}{{\mathbf{p}}_{1}}}}{{{{(2\pi)}^{D-1}}}}}\frac{{{d^{D-1}}{{\mathbf{p}}_{2}}}}{{{{(2\pi)}^{D-1}}}}\frac{{{d^{D-1}}{{\mathbf{q}}_{1}}}}{{{{(2\pi)}^{D-1}}}}\cdots\frac{{{d^{D-1}}{{\mathbf{q}}_{6}}}}{{{{(2\pi)}^{D-1}}}}}\right.\hfill\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\int{{d^{D-1}}}x\left[{{a_{{{\mathbf{p}}_{1}}}}{a_{{{\mathbf{p}}_{2}}}}{a_{{{\mathbf{q}}_{1}}}}{a_{{{\mathbf{q}}_{2}}}}a_{{{\mathbf{q}}_{3}}}^{\dagger}a_{{{\mathbf{q}}_{4}}}^{\dagger}a_{{{\mathbf{q}}_{5}}}^{\dagger}a_{{{\mathbf{q}}_{6}}}^{\dagger}}\right]\,\times\left[{{e^{i({{\mathbf{q}}_{1}}+{{\mathbf{q}}_{2}}-{{\mathbf{q}}_{3}}-{{\mathbf{q}}_{4}}-{{\mathbf{q}}_{5}}-{{\mathbf{q}}_{6}})\cdot{\mathbf{x}}}}}\right]\hfill\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left[{{\phi_{{{\mathbf{p}}_{1}}}}(\tau){\phi_{{{\mathbf{p}}_{2}}}}(\tau){\phi_{{{\mathbf{q}}_{1}}}}\left({{\tau^{\prime}}}\right){\phi_{{{\mathbf{q}}_{2}}}}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{3}}}^{*}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{4}}}^{*}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{5}}}^{*}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{6}}}^{*}\left({{\tau^{\prime}}}\right)}\right]\Bigg{]}\Big{|}0\rangle\,.\\\
\end{gathered}$ (69)

Using Eq. (8) and noting that

$\int{{d^{D-1}}}x\,{e^{i\left({{{\mathbf{q}}_{\mathbf{i}}}\pm{{\mathbf{q}}_{\mathbf{j}}}}\right)\cdot{\mathbf{x}}}}={(2\pi)^{D-1}}{\delta^{(D-1)}}\left({{{\mathbf{q}}_{\mathbf{i}}}\pm{{\mathbf{q}}_{\mathbf{j}}}}\right)\,,$
(70)

then Eq. (69) is cast into the following form,

$\begin{gathered}\Delta\left\langle{{\rho_{3}}\left({\tau,{{\mathbf{x}}_{\mathbf{0}}}}\right)}\right\rangle\supset\langle
0\Big{|}\frac{{{\mu^{3(4-D)}}}}{6}\lambda{H^{2}}\operatorname{Im}\left[{\int_{-\infty}^{\tau}d{\tau^{\prime}}a{{\left({{\tau^{\prime}}}\right)}^{D}}\int{\frac{{{d^{D-1}}{{\mathbf{p}}_{1}}}}{{{{(2\pi)}^{D-1}}}}}\frac{{{d^{D-1}}{{\mathbf{p}}_{2}}}}{{{{(2\pi)}^{D-1}}}}\frac{{{d^{D-1}}{{\mathbf{q}}_{1}}}}{{{{(2\pi)}^{D-1}}}}\cdots\frac{{{d^{D-1}}{{\mathbf{q}}_{6}}}}{{{{(2\pi)}^{D-1}}}}}\right.\hfill\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\
\delta({{\mathbf{q}}_{2}}-{{\mathbf{q}}_{3}})\delta({{\mathbf{q}}_{1}}-{{\mathbf{q}}_{4}})\delta({{\mathbf{p}}_{1}}-{{\mathbf{q}}_{5}})\delta({{\mathbf{p}}_{2}}-{{\mathbf{q}}_{6}})\delta({{\mathbf{q}}_{1}}+{{\mathbf{q}}_{2}}-{{\mathbf{q}}_{3}}-{{\mathbf{q}}_{4}}-{{\mathbf{q}}_{5}}-{{\mathbf{q}}_{6}})\hfill\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\
\left[{{\phi_{{{\mathbf{p}}_{1}}}}(\tau){\phi_{{{\mathbf{p}}_{2}}}}(\tau){\phi_{{{\mathbf{q}}_{1}}}}\left({{\tau^{\prime}}}\right){\phi_{{{\mathbf{q}}_{2}}}}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{3}}}^{*}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{4}}}^{*}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{5}}}^{*}\left({{\tau^{\prime}}}\right)\phi_{{{\mathbf{q}}_{6}}}^{*}\left({{\tau^{\prime}}}\right)}\right]\Bigg{]}\Big{|}0\rangle\,.\end{gathered}$
(71)

After some manipulations, the above equation takes the following form

$\displaystyle\begin{gathered}\left\langle{{\rho_{3}}\left({\tau,{{\mathbf{x}}_{\mathbf{0}}}}\right)}\right\rangle\supset\frac{{{\mu^{3(4-D)}}}}{6}\lambda{H^{2}}\operatorname{Im}\left[{\int_{-\infty}^{\tau}d{\tau^{\prime}}a{{\left({{\tau^{\prime}}}\right)}^{D}}\int{\frac{{{d^{D-1}}{{\mathbf{q}}_{1}}}}{{{{(2\pi)}^{D-1}}}}}\frac{{{d^{D-1}}{{\mathbf{q}}_{2}}}}{{{{(2\pi)}^{D-1}}}}\frac{{{d^{D-1}}{{\mathbf{q}}_{5}}}}{{{{(2\pi)}^{D-1}}}}}\right.\\\
\left[{{{\left|{{\phi_{{{\mathbf{q}}_{1}}}}\left({{\tau^{\prime}}}\right)}\right|}^{2}}{{\left|{{\phi_{{{\mathbf{q}}_{2}}}}\left({{\tau^{\prime}}}\right)}\right|}^{2}}\phi_{{{\mathbf{q}}_{5}}}^{2}(\tau)\phi_{{{\mathbf{q}}_{5}}}^{*2}\left({{\tau^{\prime}}}\right)}\right]\Bigg{]}\,.\end{gathered}$
(74)

It is easy to show that the desired equation for
$\Delta\langle{\rho_{3}}(\tau,{{\mathbf{x}}_{\mathbf{0}}})\rangle$ i.e. Eq.
(41), for $n=4$ is obtained accordingly.

## References

  * (1) B. S. DeWitt, “Quantum Field Theory in Curved Space-Time,” Phys. Rept. 19, 295-357 (1975).
  * (2) N. D. Birrell and P. C. W. Davies, “Quantum Fields in Curved Space,” Cambridge Univ. Press, 1984.
  * (3) S. A. Fulling, “Aspects of Quantum Field Theory in Curved Space-time,” London Math. Soc. Student Texts 17, 1-315 (1989)
  * (4) L. E. Parker and D. Toms, “Quantum Field Theory in Curved Spacetime: Quantized Field and Gravity,” Cambridge University Press, 2009.
  * (5) V. Mukhanov and S. Winitzki, “Introduction to quantum effects in gravity,” Cambridge University Press, 2007.
  * (6) S. Weinberg, “Cosmology,” Oxford University Press, 2008.
  * (7) D. Baumann, “Cosmology,” Cambridge University Press, 2022,
  * (8) H. Kodama and M. Sasaki, Prog. Theor. Phys. Suppl. 78, 1-166 (1984).
  * (9) V. F. Mukhanov, H. A. Feldman and R. H. Brandenberger, Phys. Rept. 215, 203-333 (1992).
  * (10) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)], [arXiv:1807.06209 [astro-ph.CO]].
  * (11) Y. Akrami et al. [Planck], Astron. Astrophys. 641, A10 (2020), [arXiv:1807.06211 [astro-ph.CO]].
  * (12) S. Weinberg, “The Cosmological Constant Problem,” Rev. Mod. Phys. 61, 1-23 (1989).
  * (13) V. Sahni and A. A. Starobinsky, Int. J. Mod. Phys. D 9, 373-444 (2000), [arXiv:astro-ph/9904398 [astro-ph]].
  * (14) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559-606 (2003), [arXiv:astro-ph/0207347 [astro-ph]].
  * (15) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753-1936 (2006), [arXiv:hep-th/0603057 [hep-th]].
  * (16) J. Martin, Comptes Rendus Physique 13, 566-665 (2012) [arXiv:1205.3365 [astro-ph.CO]].
  * (17) S. W. Hawking, Commun. Math. Phys. 43, 199-220 (1975) [erratum: Commun. Math. Phys. 46, 206 (1976)].
  * (18) W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * (19) W. G. Unruh and R. M. Wald, Phys. Rev. D 29, 1047-1056 (1984).
  * (20) T. Jacobson, “Introduction to quantum fields in curved space-time and the Hawking effect,” [arXiv:gr-qc/0308048 [gr-qc]].
  * (21) G. Cozzella, S. A. Fulling, A. G. S. Landulfo and G. E. A. Matsas, Phys. Rev. D 102, no.10, 105016 (2020), [arXiv:2009.13246 [hep-th]].
  * (22) H. Firouzjahi and A. Talebian, Phys. Rev. D 106, no.12, 123505 (2022), [arXiv:2210.15186 [gr-qc]].
  * (23) T. S. Bunch and P. C. W. Davies, Proc. Roy. Soc. Lond. A 360, 117-134 (1978).
  * (24) T. S. Bunch and P. C. W. Davies, J. Phys. A 11, 1315-1328 (1978).
  * (25) S. J. Avis, C. J. Isham and D. Storey, Phys. Rev. D 18, 3565 (1978).
  * (26) L. H. Ford and L. Parker, Phys. Rev. D 16, 245-250 (1977).
  * (27) P. R. Anderson, Phys. Rev. D 32, 1302 (1985).
  * (28) P. R. Anderson and L. Parker, Phys. Rev. D 36, 2963 (1987).
  * (29) I. L. Buchbinder and S. D. Odintsov, Acta Phys. Polon. B 18 , 237-242 (1987).
  * (30) S. D. Odintsov, Phys. Lett. B 213 , 7-10 (1988).
  * (31) S. D. Odintsov, Theor. Math. Phys. 82 , 45-51 (1990).
  * (32) C. Armendariz-Picon and A. Diez-Tejedor, JCAP 11, 030 (2023) [arXiv:2305.16293 [gr-qc]].
  * (33) Y. Zhang, X. Ye and B. Wang, Sci. China Phys. Mech. Astron. 63, no.5, 250411 (2020), [arXiv:1903.10115 [gr-qc]].
  * (34) X. Ye, Y. Zhang and B. Wang, JCAP 09, 020 (2022), [arXiv:2205.04761 [gr-qc]].
  * (35) E. T. Akhmedov and P. A. Anempodistov, Phys. Rev. D 105, no.10, 105019 (2022), [arXiv:2204.01388 [hep-th]].
  * (36) E. T. Akhmedov, K. V. Bazarov, D. V. Diakonov and U. Moschella, Phys. Rev. D 102, no.8, 085003 (2020), [arXiv:2005.13952 [hep-th]].
  * (37) E. T. Akhmedov, K. V. Bazarov, D. V. Diakonov, U. Moschella, F. K. Popov and C. Schubert, Phys. Rev. D 100, no.10, 105011 (2019), [arXiv:1905.09344 [hep-th]].
  * (38) V. K. Onemli and R. P. Woodard, Class. Quant. Grav. 19, 4607 (2002), [arXiv:gr-qc/0204065 [gr-qc]].
  * (39) V. K. Onemli and R. P. Woodard, Phys. Rev. D 70, 107301 (2004), [arXiv:gr-qc/0406098 [gr-qc]].
  * (40) T. Brunier, V. K. Onemli and R. P. Woodard, Class. Quant. Grav. 22, 59-84 (2005), [arXiv:gr-qc/0408080 [gr-qc]].
  * (41) S. P. Miao and R. P. Woodard, Class. Quant. Grav. 23, 1721-1762 (2006), [arXiv:gr-qc/0511140 [gr-qc]].
  * (42) T. Prokopec, N. C. Tsamis and R. P. Woodard, Phys. Rev. D 78, 043523 (2008), [arXiv:0802.3673 [gr-qc]].
  * (43) T. M. Janssen, S. P. Miao, T. Prokopec and R. P. Woodard, Class. Quant. Grav. 25, 245013 (2008), [arXiv:0808.2449 [gr-qc]].
  * (44) E. O. Kahya, V. K. Onemli and R. P. Woodard, Phys. Rev. D 81, 023508 (2010), [arXiv:0904.4811 [gr-qc]].
  * (45) S. P. Miao, N. C. Tsamis and R. P. Woodard, J. Math. Phys. 51, 072503 (2010), [arXiv:1002.4037 [gr-qc]].
  * (46) D. Glavan, S. P. Miao, T. Prokopec and R. P. Woodard, Phys. Rev. D 101, no.10, 106016 (2020), [arXiv:2003.02549 [gr-qc]].
  * (47) D. Glavan, S. P. Miao, T. Prokopec and R. P. Woodard, JHEP 03, 088 (2022), [arXiv:2112.00959 [gr-qc]].
  * (48) R. P. Woodard, Symmetry 15, no.4, 856 (2023), [arXiv:2303.05111 [hep-th]].
  * (49) M. Beneke, P. Hager and A. F. Sanfilippo, [arXiv:2312.06766 [hep-th]].
  * (50) G. ’t Hooft and M. J. G. Veltman, Nucl. Phys. B 44, 189-213 (1972).
  * (51) G. ’t Hooft and M. J. G. Veltman, Ann. Inst. H. Poincare Phys. Theor. A 20, 69-94 (1974).
  * (52) C. G. Bollini and J. J. Giambiagi, Nuovo Cim. B 12, 20-26 (1972).
  * (53) S. Deser and P. van Nieuwenhuizen, Phys. Rev. D 10, 401 (1974).
  * (54) J. S. Dowker and R. Critchley, Phys. Rev. D 13, 3224 (1976).
  * (55) A. O. Barvinsky and G. A. Vilkovisky, Phys. Rept. 119, 1-74 (1985).
  * (56) H. Firouzjahi and H. Sheikhahmadi, Phys. Rev. D 108, no.6, 065002 (2023), [arXiv:2307.00977 [gr-qc]].
  * (57) G. C. Wick, Phys. Rev. 80, 268-272 (1950).
  * (58) S. Weinberg, Phys. Rev. D 72, 043514 (2005), [arXiv:hep-th/0506236 [hep-th]].