# Inefficient acceleration of electrons in the shocked wind of the massive
star $\Theta^{1}$ Ori C within the Trapezium cluster

W. Bednarek
University of Łódź, Faculty of Physics and Applied Informatics, Department of
Astrophysics, ul. Pomorska 149/153, 90-236 Łódź, Poland,
<EMAIL_ADDRESS>

(Accepted . Received ; in original form )

###### Abstract

Shocked winds of massive stars in young stellar clusters have been proposed as
possible sites in which relativistic particles are accelerated. Electrons,
accelerated in such an environment, are expected to efficiently comptonize
optical radiation (from massive stars) and the infra-red radiation (re-
scattered by the dust within the cluster) producing GeV-TeV gamma-rays. We
investigate the time dependent process of acceleration, propagation and
radiation of electrons in the stellar wind of the massive star $\Theta^{1}$
Ori C within the Trapezium cluster. This cluster is located within the nearby
Orion Nebula (M 42). We show that the gamma-ray emission expected from the
Trapezium cluster is consistent with the present observations of the Orion
Molecular Cloud by the Fermi-LAT telescope provided that the efficiency of
energy conversion from the stellar wind to relativistic electrons is very low,
i.e. $\chi<10^{-4}$. For such low efficiencies, the gamma-ray emission from
electrons accelerated in the stellar wind of $\Theta^{1}$ Ori C can be only
barely observed by the future Cherenkov telescopes, e.g. the Cherenkov
Telescope Array (CTA).

###### keywords:

stars: winds, outflows — ISM: bubbles — radiation mechanisms: non-thermal —
gamma-rays: stars

††pagerange: Inefficient acceleration of electrons in the shocked wind of the
massive star $\Theta^{1}$ Ori C within the Trapezium
cluster–References††pubyear: 2015

##  1 Introduction

Winds of massive stars in young clusters are considered since early 80-ties as
likely sites in which cosmic ray particles are able to reach TeV-PeV energies
(e.g. Völk & Forman 1982, Cesarsky & Montmerle 1983, Völk & Forman 1982, Webb,
Axford & Forman 1985). In fact, a few of them, e.g. Cyg OB2 (Aharonian et al.
2005, Abdo et al. 2007, Abeysekara 2021, Cao et al. 2021), Westerlund 1
(Abramowski et al. 2012, Aharonian et al. 2022), Westerlund 2 (Aharonian et
al. 2007, H.E.S.S. Collaboration et al. 2011), have been recently observed at
TeV $\gamma$-rays. These clusters, with typical age estimated on a few mln
years, contain massive stars and their final products of stellar evolution,
i.e. Supernova Remnants (SNRs) and/or pulsars and Pulsar Wind Nebulae (PWNe).
It is not straightforward to indicate which type of objects, either the winds
of massive stars (isolated or cumulative forming super-bubbles), which are
still on the Main Sequence, or the remnants of stellar evolution (SNRs, PWNe)
play the dominant role in acceleration process of cosmic rays. In fact, both
groups of these objects are able to supply comparable amount of mechanical
energy into the cosmic space (of the order of $\sim$1050 erg s-1), ejected in
the form of the stellar winds or the expanding envelopes of supernovae.
However, there are clear differences between proprieties of their potential
acceleration sites. The shocks caused by the stellar winds differ from shocks
around SNRs which also differ from the shocks structures around PWNe. The
winds of massive stars are expected to keep comparable velocities during the
whole Main Sequence stage. On the other hand, the shells around the SNRs are
expected to significantly decelerate on a time scale of a few thousand years.
The structures of the outflows around the PWNe are even much more complicated.
The inner parts of the PWNe contain relativistic plasma from the pulsar. This
part is separated by the relativistic shock from the outer part of the nebula
which contains non-relativistic outflow which finally interacts with the
surrounding material from the supernova explosion and the interstellar matter.
SNRs are expected to accelerate particles to hundreds of TeV or even PeV
energies only at very early stages right after the explosion (see e.g.
Marcowith et al. 2018, Inoue et al. 2021, Brose et al. 2022. These differences
might have consequences for the maximum energies reached by the particles and
their efficiency of acceleration (see e.g. Brose et al 2020, Celli et al.
2019, Aharonian, Yang & de Oña Wilhelmi 2019). The way to distinguish between
the roles of these two types of potential acceleration sites is to search for
the $\gamma$-ray emission from very young stellar clusters in which even the
most massive stars are still observed on the Main Sequence stage. Such
clusters are still free from compact objects. The lack of significant
$\gamma$-ray emission will indicate that the presence of supernovae provide
necessary conditions for the acceleration of particles above TeV energies.

The evolution of an isolated stellar wind bubble within the star forming
region has been studied theoretically already a long time ago (Weaver et al.
1977). Up to now, a simple, one-zone models have been considered for the
radiation processes of accelerated particles in the large scale open cluster
(e.g. Bednarek 2007, Maurin et al. 2016). In these models particles
(electrons, hadrons) are injected into the dense H II regions of the cluster,
with typical size of $\sim$10 pc. These particles are accelerated/re-
accelerated on multiple shocks formed by the whole population of massive stars
within the cluster. In the present work, we consider the high energy radiation
produced by electrons accelerated in the wind shock of the massive star which
dominates the energy output of the stellar cluster.

The interesting candidate for the application of our model is the Trapezium
cluster in the Orion Nebula (M 42) which is one of the closest star formation
region with the radius $\sim$4 pc at the distance of 410 pc (Reid et al.
2009). The Orion Nebula cluster contains a very young cluster of two thousand
stars within the distance of 20 light years. The 4 brightest stars within 1.5
light years form so called the Trapezium cluster. The most luminous star in
the Trapezium cluster, $\theta^{1}$ Ori C, has the luminosity $2\times 10^{5}$
L⊙ (O6Vp type) (Simon Diaz et al. 2006). Its radiation produces dominant
source of ionization of the gas within surrounding Orion Nebula cloud which
mass is $\sim 2.6\times 10^{3}~{}M_{\odot}$ (Pabst et al. 2019).

We investigate the hypothesis in which particles (electrons, hadrons) are
accelerated on a strong shock formed in the stellar wind from the star
$\theta^{1}$ Ori C as a result of its interaction with the surrounding gas of
Orion Nebula. However, due to the strong radiation field from this massive
star, and also infra-red radiation re-processed by the surrounding gas,
electrons should suffer strong energy losses on the Inverse Compton scattering
(ICS) process, producing $\gamma$-rays in the GeV-TeV energy range. We follow
the propagation of electrons in the wind of the massive star taking into
account mentioned above their energy losses and by applying the hydrodynamical
model developed by Weaver et al. 1977). From the comparison of the expected
$\gamma$-ray emission with the observational constraints by the satellite and
ground telescopes, we constrain the efficiency of electron acceleration in the
cavity around massive star $\theta^{1}$ Ori C. The consequences of these
constraints are discussed.

We have chosen the stellar cluster in which a single massive star dominates.
However, often many stars produce comparable wind cavities which can merge
forming a super-bubble around the whole cluster. Possible acceleration of
particles within super-bubbles around young clusters has been considered in
e.g. Bykov (2001), Bykov (2014), or recently by e.g. Vieu et al. (2022) and
Blasi & Morlino (2023).

$\theta^{1}$ Ori C | L⋆ | M⋆ | T⋆ | R⋆ | B$\star$ | Lw | $\dot{M}_{\rm w}$ | vw
---|---|---|---|---|---|---|---|---
unit | L⊙ | $M_{\odot}$ | K | R⊙ | G | erg/s | M⊙/yr | km/s
value | $2\times 10^{5}$ | 33 | $3.9\times 10^{4}$ | $10.6\pm 1.5$ | $1100\pm 100$ | $8\times 10^{35}$ | $4\times 10^{-7}$ | 2500
Table 1: The basic parameters of the massive star $\theta^{1}$ Ori C: Stellar:
the luminosity L⋆ (a), the mass M⋆ (b), the surface temperature T⋆ (a), the
radius R⋆ (a), the surface magnetic field strength (e), the power of the
stellar wind Lw (c,d), the mass loss rate $\dot{M}_{\rm w}$ and the wind
velocity vw (c,d). References: (a) Simon Diaz et al. (2006), (b) Balega et al.
(2014), (c) Howarth et al (1989), (d) Stahl et al. (1996), (e) Donati et al.
(2002).

##  2 A cavity around an early type star within the stellar cluster

The Trapezium cluster contains four luminous stars known as $\theta^{1}$ Ori
A,B,C,D. Around $80\%$ of the cluster ionizing radiation is provided by the
most massive member of the cluster, $\theta^{1}$ Ori C, which is a binary
system of O type star and B type star (Balega et al. 2014 ). The star
$\theta^{1}$ Ori C is immersed in the region of a low density, hot, X-ray
emitting gas which forms a bubble expanding with the velocity of $\sim$13 km
s-1 (see Fig. 3 in Pabst et al. 2019). The star produces a fast wind which
parameters are reported in Table 1. The external pressure of the hot gas
balances the pressure of the stellar wind forming a large scale shock at the
distance $R_{\rm 1}$ from the star (see Fig. 1 in Weaver et al. 1977). This
shock can efficiently accelerate particles since the plasma flows through this
shock with large velocity. On the other hand, the shock on the border between
the shocked interstellar gas and the ambient interstellar gas, at $R_{2}$, is
slow (the velocity of $\sim 13$ km s-1). The slow outer shock is not able to
accelerate electrons to large energies since the acceleration efficiency
depends on the velocity of the plasma flow through the shock (see Sect. 3).
The small scale shocks in the wind from $\theta^{1}$ Ori C can also appear as
a result of the interaction with less massive stars (e.g. $\theta^{1}$ Ori A
and D). These two stars have similar masses producing similar winds. They are
also located at the similar projection distance from $\theta^{1}$ Ori C. If
shocks around these stars are able to accelerate electrons to $\sim$TeV
energies, then they can lose efficiently energy on the ICS of thermal
radiation from the stellar surfaces. Moreover, the star $\theta^{1}$ Ori C is
itself a compact binary system of massive stars. The main component (C1) has
the parameters reported in Table. 1. The companion star (C2) is the B type
star with the mass $11\pm 5$ M⊙ (Balega et al. 2014).

The age of the Trapezium cluster is estimated to be $<$1 mln years, with the
likely value of $\sim$2$\times 10^{5}$ yrs (Pabst et al. 2019). Therefore,
even super-massive stars are expected to be still on the Main Sequence stage.
The Orion Nebula M42 (galactic coordinates: $209^{o},-19.27^{o}$) is one of
the closest, massive star forming regions to the Sun.

Figure 1: The time scale for the energy gains of electrons ($\tau_{\rm acc}$)
in the acceleration process at the shock in the wind around the massive star
$\theta^{1}$ Ori C is compared with the electron time scales for the energy
losses on different process: the advection of electrons from the shock with
the velocity of the wind from the star $\theta^{1}$ Ori C (dot-dashed), the
synchrotron energy losses, $\tau_{\rm syn}$ (dashed), and the IC energy losses
in the stellar radiation, $\tau_{\rm IC/star}$ (dotted), the infra-red
radiation from dust $\tau_{\rm IC/inf}$ (dot-dot-dashed), and the MBR
$\tau_{\rm IC/MBR}$ (dot-dot-dot-dashed) at different moments of the expansion
of the stellar wind nebula counting from its origin: $T=10^{2}$ yrs (on the
left), $10^{3}$ yrs (left-middle), $10^{4}$ yrs (right-middle), and $2\times
10^{5}$ yrs (right). The surface magnetic field of the star is assumed to be
equal to 300 Gs.

The powerful wind from the dominant star within the Trapezium cluster,
$\theta^{1}$ Ori C, interacts with the surrounding medium forming a large
scale cavity (a bubble) in the Orion Nebula. This cavity is filled with a hot
shocked gas with the temperature of the order of $T\sim 2\times 10^{6}$ K
(Gödel et al. 2008). In the interaction of the fast stellar wind with this
surrounding gas a shock structure appears. Detailed structure of the stellar
wind cavity, at different moments after its formation, is modelled by Castor
et al. (1975) and Weaver et al. (1977). This model gives analytic solutions of
the hydrodynamic equations for the gas content of the cavity at different
stages of its evolution (i.e. as a function of time) and as a function of
basic parameters of the stellar wind (wind velocity and the mass loss rate)
and the density of surrounding medium. Schematic structure of the bubble
around a massive star in stellar cluster is shown in Fig. 1 in Weaver et al.
(1977). Four regions around the star can be distinguished: (1) the fast
stellar wind; (2) the shocked stellar wind; (3) the shocked interstellar gas;
and (4) the ambient interstellar gas. Specific parts are separated by the
shock in the stellar wind (at $R_{1}$), the contact surface (at $R_{\rm c}$),
and the outer radius of the shocked interstellar gas (at $R_{2}$). We assume
that electrons are accelerated at the strong shock at the distance $R_{1}$.
They are advected to the outer part of the bubble with the velocity of the
shocked stellar wind described by Eq. 11 (Weaver et al. 1977). The shocked
stellar wind reaches the boundary between different media at the distance from
the star $R_{\rm c}$ with the velocity $v(R_{\rm c})=3R_{\rm c}/(5t)$, where
$t$ is the time after formation of the star. For the parameters of
$\theta^{1}$ Ori C, this velocity is close to $\sim$10 km s-1. For the
parameters of $\theta^{1}$ Ori C, we derive the formulae for above mentioned
characteristic distances (see Weaver et al. 1977) in order to model the
processes of propagation and radiation of electrons accelerated in the region
of the strong shock in the wind region.

Note however the shortcomings of the Weaver et al. model. For example, it
ignores the effects of the photo-ionization on the dynamics of the wind
bubble. The numerical simulations (e.g. Dwarkardas 2022), and also analytic
calculations (e.g. Chevalier 1997), have shown that due to hydro-dynamic
instabilities the spherical symmetry of the wind can be destroyed. We do not
take these effects into account in our simplified model for the wind bubble
dynamics which is considered to be spherically symmetric.

Electrons, accelerated at the early time in the bubble (i.e for $t<100$ yrs),
are advected with the shocked stellar wind. Electrons with energies above a
few GeV have enough time to cool in the magnetic field and the radiation field
of the massive star and the re-processed infra-red radiation fields (see Fig.
1 on the left, for the electron cooling time scales and their advection time
scale at the example age of the wind bubble equal to 100 yrs). These electrons
cool already within the stellar wind region below the distance $R_{\rm c}$
(calculated for the age of the wind bubble). Therefore, early accelerated
electrons are not transported into the dense, shocked interstellar gas. They
are confined close to the wind shock region. The soft radiation field within
the stellar wind cavity comes mainly from the massive star (optical), from the
surrounding cloud (infra-red), and also the Microwave Background Radiation
(MBR). The density of the optical radiation is assumed to drop with the
distance squared from the star and the density of the infrared radiation is
kept to be constant within the wind bubble. The infra-red emission by dust, in
the swept up shell around the wind bubble, is $6\times 10^{4}$ L⊙. The
temperature of the dust is $\sim$200 K (Pabst et al. 2019). The dust is heated
by the optical emission from $\theta^{1}$ Ori C. These optical photons
propagate freely through the ionized medium of the hot bubble. We calculate
the average density of this infra-red photons assuming the dimension of the
swept up dust layer equal to $\sim$2 pc. The wind bubble with such radius
around $\theta^{1}$ Ori C has been measured by Pabst et al. 2019). This value
is consistent with the theoretical solution for the wind structure around
massive stars (Weaver et al. 1977).

##  3 High energy electrons in the stellar wind cavity

It is assumed that electrons are accelerated at a constant rate (independent
on time, i.e. the amount of energy transferred from the shock to the
relativistic electrons is time independent) at the stellar wind shock from the
moment of the formation of the cavity up to the observation time, i.e. after
$2\times 10^{5}$ yrs. Electrons are advected from the shock region with the
stellar wind. During their propagation in the stellar wind, they lose energy
mainly on the IC and the synchrotron processes. Energy losses of electrons
depend on time due to the time dependent conditions within the cavity (i.e.
the radiation field and the strength of the magnetic field). We take into
account the effects of electron energy losses in the varying medium applying
the step time method.

Let us define the structure of the magnetic field around the massive star
which is needed to determine the acceleration process of electrons at the
stellar wind shock and their synchrotron energy losses. The inner structure of
the magnetic field of the star $\theta^{1}$ Ori C is assumed to be of the
dipole type with an intensity at the surface reported in Table 1 and the
magnetic axis inclined at the angle of $42^{\circ}\pm 6^{\circ}$ with respect
to the rotation axis (Donati et al. 2002), i.e. $B(r)=B_{\star}/r^{3}$, where
$r=R/R_{\star}$ and $B_{\star}=10^{3}B_{3}$ G. The density of the stellar wind
at some distance, $R$, from the star is, $\rho=\dot{M}_{\rm w}/4\pi
R^{2}v_{\rm w}m_{\rm p}\approx 1.8\times 10^{-14}/r^{2}$ g cm-3. Its energy
density depends on the distance from the star as $U_{\rm w}=\rho v_{\rm
w}^{2}/2\approx 550/r^{2}$ erg cm-3. This energy density of the stellar wind
starts to dominate over the energy density of the magnetic field, $\rho_{\rm
B}=B(r)^{2}/(8\pi)\approx 4\times 10^{4}B_{3}^{2}/r^{6}$ erg cm-3, at
distances larger than the Alfven radius which is equal to $r_{\rm A}\approx
2.9B_{3}^{1/2}$ for the parameters of $\theta^{1}$ Ori C. We conclude that the
magnetic field of the star has the dipole structure only below the Alfven
radius, $r_{\rm A}$. At distances above $r_{\rm A}$, the magnetic field
structure in the wind is well described by the dominant radial component,
$B(r)=(B_{\star}/(r_{\rm A})(r_{\rm r/t}/r)^{2})$ where $r_{\rm r/t}=v_{\rm
w}/v_{\rm rot}$ is the radius where the radial magnetic field turns into
toroidal. $v_{\rm rot}$ is the rotational velocity of $\theta^{1}$ Ori C equal
to 53 km/s (Howarth et al. 1997). For $\theta^{1}$ Ori C, we obtain $r_{\rm
r/t}\sim 50$. At large distances from the star, the magnetic field becomes
toroidal due to the rotation of the star, $B(r)\approx B_{\star}/(r_{\rm
A}r_{\rm r/t}r)$. Then, the magnetic field strength at large distance from the
star $\theta^{1}$ Ori C (i.e at the radius of the shock $r_{\rm sh}\gg r_{\rm
r/t}$) obtains toroidal structure. It can be expressed by $B_{\rm sh}\approx
4\times B(r_{\rm sh})\approx 2\times 10^{-5}B_{3}/R_{18}$ G for the toroidal
magnetic field, where the distance from the star is renormalized according to
$R=10^{18}R_{18}$ cm.

Figure 2: Spectral Energy Distribution (SED) of the $\gamma$ ray emission
produced by electrons accelerated in the large scale shock around the star
$\theta^{1}$ Ori C. Electrons are injected with the power law spectrum with
the spectral index equal to -2, extending from 1 GeV up to $E_{\rm max}$.
$E_{\rm max}$ is determined by the balance between their the acceleration time
scale and the dynamical time scale or the synchrotron energy loss time scale
(see Fig. 1). In these calculations we accept that $\chi=0.1$ of the wind
power is transferred to relativistic electrons. Electrons are accelerated
during the period $t_{\rm min}=10$ yrs and $t_{\rm max}=2\times 10^{5}$ yrs
after the formation of the wind cavity. $\gamma$-rays from the comptonization
of the stellar radiation (dot-dashed curves), the infra-red radiation
(dashed), and the MBR (dotted). The surface magnetic field of the star
$\theta^{1}$ Ori C is $B=100$ G (figure on the left), 300 G (middle), and
$10^{3}$ G (right).

Particles are assumed to be accelerated at the large scale shock in the
stellar wind with the power law spectrum up to the maximum energies. In order
to determine these maximum energies of electrons, we compare the time scales
for the energy gains of electrons with their time scales for the energy losses
on different processes. The energy gain time scale from the acceleration
mechanism at the shock can be parametrised by $\tau_{\rm acc}=R_{\rm
L}/\xi/c\approx 10^{2}E_{\rm TeV}/\xi_{-3}B_{\rm G}$ s, where $R_{\rm L}$ is
the Larmor radius of the electrons, $c$ is the velocity of light, the energy
of the electron is $E_{\rm e}=1E_{\rm TeV}$ TeV, the magnetic field strength
at the shock ($B_{\rm G}$) is expressed in Gauss, and $\xi=10^{-3}\xi_{-3}$ is
the acceleration coefficient. The advection time scale of electrons from the
shock with the stellar wind is $\tau_{\rm adv}=3R_{\rm sh}/v_{\rm w}\approx
10^{10}R_{18}$ s. The synchrotron energy loss time scale of electrons is
$\tau_{\rm syn}=3m_{\rm e}\gamma_{\rm e}/(4c\rho_{\rm B}\sigma_{\rm
T}\gamma_{\rm e}^{2})\approx 141/(B_{\rm G}^{2}E_{\rm TeV})$ s, where $m_{\rm
e}$ is the electron rest mass, $\sigma_{\rm T}$ is the Thomson cross section,
$\gamma_{\rm e}=E_{\rm e}/m_{\rm e}c^{2}$ is the Lorentz factor of electrons.
We have also calculated the energy loss time scale of electrons on the Inverse
Compton Scattering of the radiation from the massive star $\theta^{1}$ Ori C,
infrared radiation (reprocessed in the gas) and the Microwave Background
Radiation (MBR) by using the approximate analytic formula from Moderski et al.
(2005). These radiation fields are defined in Sect. 2. All those energy gains
and losses of electrons, at the large scale shock in the wind of $\theta^{1}$
Ori C, are compared for different moments after the formation of the wind
cavity (see Fig. 1). In these calculations, the acceleration coefficient of
electrons in the shock is assumed to be equal to $\xi=0.3(v_{\rm w}/c)^{2}$.
We estimate that electrons can reach energies of the order of several up to a
few tens of TeV (Fig. 1). At the early stage of the nebula, the acceleration
process of electrons is limited mainly by the synchrotron process but at the
latter time by the dynamical time of the shock.

We assume that electrons are accelerated with the power law spectrum,
$dN/(dEdt)=KE^{-\alpha}$ (and the spectral index $\alpha=2$) up to the maximum
energies as constrained in Fig. 1. $K$ is the normalization factor of the
spectrum to the part of the kinetic power of the stellar wind. In fact, the
spectra of particles accelerated at the stellar wind termination shocks could
be much softer in the case of curved shocks around the Sun’s type stars.
However, it has been shown (e.g. see Sect. 5 in Webb et al. 1985) that this is
not the case for the shocks with the parameters typical for the OB type wind
shocks. For the modelling of the $\gamma$-ray emission, we assume that
$\chi=0.1$ of the wind energy is transferred to relativistic electrons with
energies above $E_{\rm e}^{\rm min}=1$ GeV ($L_{\rm e}=0.1L_{\rm w}\approx
8\times 10^{34}$ erg s-1). Electrons are confined within the advection time at
the shock region, losing a part of their energies on the synchrotron and the
IC processes. The electrons in the most energetic part of their spectrum, i.e.
above a few TeV (see Fig. 1), lose energy mainly on the synchrotron process.
Produced synchrotron radiation extends between $\sim$2 eV up to $\sim$20 keV,
for the electrons with energies between $\sim 3-300$ TeV. These electrons can
also produce TeV $\gamma$-rays by scattering optical stellar radiation and the
infra-red radiation (stellar radiation reprocessed by the gas). Electrons with
energies below $\sim$3 TeV interact mainly in the IC process by scattering
optical radiation from the star $\theta^{1}$ Ori C. They produce $\gamma$-rays
with the spectrum extending up to $\sim$TeV energy range (due to the
scattering in the Thomson and the Klein-Nishina regimes). Our model for the
$\gamma$-ray production is applied to the star $\theta^{1}$ Ori C, since the
radiation output and the wind power output of this star dominate the
energetics of the Trapezium cluster.

##  4 $\gamma$-rays from IC scattering

Figure 3: As in Fig. 2 but for the spectral index of electrons equal to -3.

As we have shown above, the energies of electrons injected at the shock in the
stellar wind evolve in time due to their time dependent energy losses in
varying magnetic field and the radiation field within the wind cavity. The
calculation of the high energy radiation from electrons, which were
accelerated in the past, have been done numerically since the physical
conditions for these electrons (and so their injected spectrum but not their
injected power) depend on time. It is assumed that the bubble has the age of
$t_{\rm max}=2\times 10^{5}$ yrs. We follow the evolution of the wind bubble
from the moment of its formation up to, $t_{\rm max}$, with the step time
method in the log scale with the step $\Delta$ log(t). At the time ’t’ in the
past (counted from the formation of the wind bubble, with the step in time
$\Delta~{}\log(t)=0.001$ (yr)), we inject relativistic electrons with the
power law spectrum up to the maximum energies which are determined from the
comparison of their acceleration time scale with their energy loss time scales
or their advection time scale. The spectrum of electrons at the time, t, is
fragmented onto a sequence of bunches with the size in energy
$\Delta\log(E_{\rm(n)})=\Delta\log(E_{\rm(n-1)}$) + 0.025 (TeV), starting from
$\Delta\log(E_{\rm 1})$ = -3, up to the maximum energies allowed by their
energy losses. The number of electrons in the bunch with energy, E, is
calculated from $\Delta N=(dN/dEdt)E(\Delta\log(E))t(\Delta\log(t))(\ln
10)^{2}$. Different radiation processes such as, the synchrotron and the IC
processes in the stellar radiation, infra-red radiation, and also the
Microwave Background Radiation (at late periods after formation of the wind
bubble), are taken into account. The fate of these specific bunches of
electrons are followed from the time ’t’ up to ’$t_{\rm max}$’, including
their energy losses and the location within the bubble at different moments
’t’. The evolution of the energy of electrons in the bunch is followed in time
in linear scale with the time step, $\Delta t=10^{-4}\tau_{\rm loss}$, where
$\tau_{\rm loss}$ is the energy dependent time scale for the electron energy
losses on considered radiation processes (i.e. $1/\tau_{\rm loss}=1/\tau_{\rm
IC}+1/\tau_{\rm syn}$). In this way, we arrive with a specific electron bunch
to the present moment $t_{\rm max}$. Gamma-ray emission is calculated as the
sum of emission from all the bunches of electrons which arrived to the time,
$t_{\rm max}$ with energies above 1 GeV . So, we do not calculate the spectrum
of electrons at the time $t_{\rm max}$ but only radiation from the whole
population of the electron bunches which fate is considered in time from the
moment of their acceleration, at the time ’t’, up to the present moment at
’$t_{\rm max}$’.

We calculate the $\gamma$-ray spectra for a few different values of the
surface magnetic field of the massive star (see Fig. 2). The $\gamma$-ray
spectra are dominated by the IC scattering of the stellar radiation with
significant contribution of the IC scattering of the infra-red radiation at
the largest energy range. The IC scattering of the MBR plays minor role. The
spectra show steepening at around $\sim$100 GeV due to the transition between
the IC scattering in the Thomson and the Klein-Nishina regimes and also time
dependent cooling. The $\gamma$-ray spectra, produced by the electrons
injected with a steeper spectrum, show strong emission at GeV energies but low
level of emission at TeV energies (see Fig. 3).

We also investigate the contribution to the observed at present $\gamma$-ray
spectrum from the electrons which were accelerated at different ages of the
wind cavity. The contribution to the electron spectrum is shown on the right
of Fig. 4. The $\gamma$-ray spectra, produced by these electrons, are shown on
the left of Fig. 4. It is clear that electrons accelerated early in the wind
bubble lose energy very efficiently. They are not able to contribute to the
presently observed spectrum at TeV energy range. However, electrons,
accelerated at the bubble with the age $t=10^{4}$ yrs, can still contribute
significantly to the GeV energy range (see the upper Fig. 4). On the other
hand, $\gamma$-rays with energies above $\sim$100 GeV can only be produced by
electrons accelerated within the last half of the lifetime of the stellar wind
cavity (see the bottom Fig. 4). Therefore, our estimate should be considered
as an upper limit.

Figure 4: Spectral Energy Distribution (SED, figures on the left) of the
$\gamma$-ray emission (parameters as in Fig. 2) but for the electrons
accelerated at different early periods after formation stellar wind cavity
around $\theta^{1}$ Ori C (upper figure): between $t=10^{5}$ yrs and
$t=2\times 10^{5}$ yrs (dot-dashed curve), between $t=10^{4}$ yrs and
$t=10^{5}$ yrs (dashed), between $t=10^{3}$ yrs and $t=10^{4}$ yrs (dotted),
and between $t=10^{2}$ yrs and $t=10^{3}$ yrs (dot-dot-dashed). The
$\gamma$-ray emission from the whole period of electron acceleration, between
$t_{\rm min}=10$ yrs and $T_{\rm max}=2\times 10^{5}$ yrs is marked by the
solid curve. The emission from electrons accelerated at late periods is shown
in the bottom figure, between $t=1.85\times 10^{5}$ yrs and $t=2\times 10^{5}$
yrs (dot-dot-dashed curve), between $t=1.65\times 10^{5}$ yrs and
$t=1.85\times 10^{5}$ yrs (solid), between $t=1.45\times 10^{5}$ yrs and
$t=1.65\times 10^{5}$ yrs (dotted), between $t=1.25\times 10^{5}$ yrs and
$t=1.45\times 10^{5}$ yrs (dashed), and between $t=10^{5}$ yrs and
$t=1.25\times 10^{5}$ yrs (dot-dashed). The corresponding spectra of
electrons, for which these $\gamma$-ray spectra have been calculated, are
shown on the right figures.

For the stellar wind cavity around the star $\theta^{1}$ Ori C, with the
radius of $R_{\rm neb}=2$ pc (Pabst et al. 2019), and the distance of the
Trapezium cluster $\sim 410$ pc, the angular size of the TeV $\gamma$-ray
emission is equal to $\sim 0.28$ degree. This size of the wind cavity is
comparable to the diffusion distance of the electrons with energies 1 TeV in
the magnetic field at the shock region (estimated on $B=3\times 10^{-6}$ G,
see Sect. 3), assuming the simple Bohm prescription for the diffusion
coefficient, which should be considered as an optimistic assumption on the
diffusion model of electrons within the wind bubble. However note, that the
magnetic field is expected to have toroidal structure at the wind bubble
boundary. Therefore, we conclude that in fact electrons accelerated at the
stellar wind shock are likely confined within the shocked wind around
$\theta^{1}$ Ori C.

In principle, these $\gamma$-rays might be absorbed in the radiation from the
star $\theta^{1}$ Ori C. We estimate the optical depth for the $\gamma$-rays
as a function of distance, $D=d\times R_{\star}$, from the star, on
$\tau_{\gamma-\gamma}\approx D\times n_{\rm ph}\times\sigma_{\gamma-\gamma}$,
where $n_{\rm ph}\approx 1.5\times 10^{15}/d^{2}$ ph./cm3 is the density of
stellar photons as a function of distance from the star for the surface
temperature and the radius of the star given in Table 1.
$\sigma_{\gamma-\gamma}\approx\sigma_{\rm T}/3$ is the cross section for
$\gamma-\gamma$ collision. Then, $\tau_{\gamma-\gamma}\approx 187/d$. The
optical depth for the $\gamma$-rays is below one for distances from the star
$D>1.4\times 10^{14}$ cm. Therefore, we conclude that $\gamma$-rays, produced
at the present location of the shock in the stellar wind region equal to
$r_{1}\approx 10^{18}$ cm (see Eq. 12 in Weaver et al. 1977), do not suffer
strong absorption.

##  5 Constraints on the acceleration efficiency of electrons

The Orion Molecular Complex (OMC) is one of the largest, nearby concentration
of gas in the solar vicinity. It also contains many interesting stellar
objects including the Trapezium cluster. Therefore, it is one of the favourite
targets for the $\gamma$-ray observatories. It was detected already by the
COS-B (Caraveo et al. 1980) and the EGRET (Digel et al. 1995) Observatories.
Recently, OMC is also observed in the GeV $\gamma$-ray energy range by the
Fermi-LAT detector (Ackermann et al. 2012, Yang et al. 2014, Neronov et al.
2012, Peng et al. 2019, Baghmanyan et al. 2020). The $\gamma$-ray emission
from OMC is extended forming two separate parts OMC A and OMC B. This
$\gamma$-ray emission is directionally consistent with the mass concentrations
within OMC. Therefore, it is likely due to the interaction of cosmic ray
hadrons with the matter of the complex (e.g. see Young et al. 2014). The upper
limit on the TeV $\gamma$-ray emission from the OMC has been recently reported
by the HAWC Collaboration (Albert et al. 2021). The level of the $\gamma$-ray
emission, predicted in terms of our model from the Trapezium cluster (for the
electron acceleration efficiency $\chi<10^{-3}$), is in fact clearly below the
upper limit on the extended emission from the OMC A (Albert et al. 2021).
Recently, the upper limits on the GeV $\gamma$-ray emission from the direction
of the Trapezium cluster has been also obtained by Maurin et al. (2016), based
on the analysis of the Fermi-LAT data. Maurin et al. (2016) considered a
hadronic model as responsible for the $\gamma$-ray emission from the stellar
wind cavity around $\theta^{1}$ Ori C. In their model, hadrons, accelerated at
the stellar wind shock, produce $\gamma$-rays via decay of pions produced in
their collisions with the matter of the surrounding gas. They constrained the
efficiency of hadron acceleration in the wind of $\theta^{1}$ Ori C below a
few percent.

In Fig. 5, we compare the $\gamma$-ray spectra, expected in terms of our
leptonic model from the stellar wind cavity around $\theta^{1}$ Ori C star,
with the upper limits on the Fermi-LAT emission observed towards the direction
of the Trapezium cluster (Maurin et al. 2016). Based on this comparison, we
conclude that in order to be consistent with this upper limit, the energy
conversion efficiency from the stellar wind to energetic electrons cannot be
larger than $\chi\sim 10^{-4}$. The efficiency of the electron acceleration is
clearly stronger constrained than the efficiency of acceleration of hadrons
obtained by Maurin et al. (2016). This is due to the fact that electrons
suffer much stronger energy losses than hadrons. They are not able to
accumulate efficiently within the stellar wind nebula around $\theta^{1}$ Ori
C star. Future observations with the Cherenkov Telescope Array (CTA) will be
able only marginally constrain the TeV $\gamma$-ray emission from the wind
cavity around $\theta^{1}$ Ori C star provided that they have leptonic origin
(see Fig. 5). Based on the comparison of the Fermi-LAT upper limits with the
results of our modelling, we conclude that electrons can not be efficiently
accelerated within the wind cavity around the massive star $\theta^{1}$ Ori C.

Figure 5: Spectral Energy Distribution (SED) of the $\gamma$ ray emission
produced by electrons accelerated in the large scale shock around the star
$\theta^{1}$ Ori C (solid curve) is scalled by the electron acceleration
efficiency $\chi=10^{-4}$. It is compared with the Fermi-LAT upper limits on
the GeV $\gamma$-ray flux from the Trapezium cluster (dot-dashed line, Maurin
et al. 2016) and also the upper limit on the TeV $\gamma$-ray emission from
the Orion Molecular Cloud obtained by the HAWC Collaboration (dashed curve,
see Table 3 in Albert et al. 2021). We also show the sensitivity of the future
Cherenkov Telescope Array (50 hr sensitivity) of the CTA South by the thin
dotted curve (Maier et al. 2017).

##  6 The role of shocks around other stars within Trapezium cluster

The brightest, visible stars in the Trapezium cluster (mentioned above) are in
fact close binary systems of massive stars. Therefore, particles (electrons
and hadrons), accelerated on the shocks between the winds of these companion
stars, could also produce $\gamma$-ray emission as discussed by e.g. White &
Chen (1992), Eichler & Usov (1993), Benaglia & Romero (2003), Torres et al.
(2004), Bednarek (2005), Reimer et al. (2006), Pittard & Dougherty (2006).

In fact, $\theta^{1}$ Ori C is a binary system of an O type massive star and a
B type star with the semi-major axis equal to $a=4.1\times 10^{13}$ cm (for
the distance of $D=415.0\pm 6$ pc), and the eccentricity of 0.6 (Rzaev et al.
2021). The distance between the stars at the periastron is $D=a(1-e)\approx
1.6\times 10^{13}$ cm and at the apoastron is $D=a(1+e)\approx 6.6\times
10^{13}$ cm. The wind parameters of the B type star in this binary system are
unfortunately unknown. Therefore, we will apply the typical values for B type
stars, i.e. the wind velocity of the order of $v_{\rm w}=3\times 10^{8}$ cm
s-1, and their mass loss rates of the order of $\dot{M}_{\rm w}=10^{-9}$ M⊙
s-1. As a result of collisions of these winds, a double shock structure
appears within the $\theta^{1}$ Ori C binary system at the distance from the
more massive star equal to, $X_{\rm p,a}=R_{p,a}/(1+\sqrt{\eta})\approx
0.053R_{p,a}$, where the parameter $\eta$ for the binary system $\theta^{1}$
Ori C is estimated on, $\dot{M}_{\rm w}^{\rm O}v_{\rm w}^{\rm O}/\dot{M}_{\rm
w}^{\rm B}v_{\rm w}^{\rm B}\approx 320$, where $\dot{M}_{\rm w}^{\rm O}$ and
$\dot{M}_{\rm w}^{\rm B}$ are the mass loss rates of the O type and B type
stars and $v_{\rm w}^{\rm O}$ and $v_{\rm w}^{\rm B}$ are their wind
velocities, respectively. We estimate the solid angle overtaken by the wind of
the Be type star as seen from the location of the more massive star in
$\theta^{1}$ Ori C. It is equal to $\Delta\Omega=\pi\theta^{2}/(4\pi)\sim
7\times 10^{-4}$, where $\theta=1/(1+\sqrt{\eta})$. Therefore, we conclude
that only a very small part of the power of the wind from the more massive
star in $\theta^{1}$ Ori C, of the order of $\Delta\Omega L_{\rm w}\approx
6\times 10^{31}$ erg s-1, can be eventually converted to relativistic
electrons in the region of collisions of these winds with the binary system
$\theta^{1}$ Ori C. This is much less than the power available for the
electrons accelerated at the large scale shock around the Trapezium cluster
considered above. Therefore, we can neglect the role of the acceleration of
electrons within the binary system $\theta^{1}$ Ori C itself.

Other main members of the cluster, i.e. the stars $\theta^{1}$ Ori A and
$\theta^{1}$ Ori D, could in principle also form obstacles for the dominant
wind from the star $\theta^{1}$ Ori C. We assume that B type stars in
$\theta^{1}$ Ori A and $\theta^{1}$ Ori D have typical parameters expected for
B type stars. The massive star in multiple binary system $\theta^{1}$ Ori A is
a B0.5V star with the mass of $15.3$ M⊙. It is at the projected distance of
12.7 mas from $\theta^{1}$ Ori C. This angular distance corresponds to the
physical separation of $7.6\times 10^{16}$ cm, assuming the distance to
Trapezium cluster of $\sim$415 pc. The massive companion in $\theta^{1}$ Ori A
star has the radius $4.5$ R⊙, the luminosity $1.6\times 10^{4}$ L⊙ and the
surface temperature $3\times 10^{4}$ K (Nieva & Przybilla 2014). In collisions
of the winds from $\theta^{1}$ Ori C and $\theta^{1}$ Ori A, the shock
structure defined by the parameter $\eta$ equal to $\sim$320 should form. So
then, as in the case of the shock around the B type star in the binary system
$\theta^{1}$ Ori C, this shock can extract only a small part of the wind
energy from the stellar wind of $\theta^{1}$ Ori C. The star $\theta^{1}$ Ori
D has quite similar parameters to $\theta^{1}$ Ori A, i.e the mass of $18\pm
6$ M⊙ at the projected distance of 13.4 mas from $\theta^{1}$ Ori C, the
radius $5.6$ R⊙, the luminosity $2.95\times 10^{4}$ L⊙, and the surface
temperature $3.2\times 10^{4}$ K (Simón-Díaz et al. 2006). Therefore, the
effect of the star $\theta^{1}$ Ori D can be also safely neglected in respect
to the radiation processes occurring at the large scale shock around Trapezium
cluster.

In the case of strongly magnetized, massive stars of the $\theta^{1}$ Ori C
type, also particles accelerated in other mechanisms might contribute to the
expected $\gamma$-ray emission. For example, particles can be also accelerated
in the reconnection region of the stellar magnetic field forced by the strong
winds from the luminous star (e.g. Usov & Melrose 1992). Electrons,
accelerated in the induced electric field, can produce synchrotron non-thermal
radiation (Trigilio et al. 2004, Leto et al. 2017) and also up-scatter dense
stellar radiation up to $\sim$TeV $\gamma$-ray energies (e.g. Bednarek 2021).
We conclude that possible contribution from other processes to the
$\gamma$-ray emission from the stellar wind bubble around $\theta^{1}$ Ori C
can only farther limit estimated above energy conversion efficiency from the
stellar wind to relativistic electrons.

##  7 Conclusion

We elaborated a model for the acceleration of electrons in the wind cavity
around the most massive star in the Trapezium cluster (i.e. $\theta^{1}$ Ori
C) within the Orion Nebula (M 42). We choose a very young cluster (the age
$\sim$2$\times 10^{5}$ yrs) in which even the most massive stars have not
exploded yet as a supernovae. In this way, we are able to investigate the high
energy radiation which originates in the winds of massive stars. This
radiation is not contaminated by the presence of shell type supernova remnants
or Pulsar Wind Nebulae and the compact objects (in the binary systems) formed
during supernova explosions. It is assumed that electrons are accelerated at a
constant rate during the whole lifetime of the cavity around the most massive
star in Trapezium cluster, i.e. $\theta^{1}$ Ori C. For reasonable
acceleration efficiency, electrons can reach TeV energies. However, electrons
suffer significant energy losses, due to the synchrotron and the Inverse
Compton processes, which strongly modify their spectra observed at present
within the stellar wind cavity. Applying the time dependent hydrodynamic model
for the evolution of the wind cavity (Weaver et al. 1977), we calculate the
$\gamma$-ray emission from electrons at the present moment within the cavity
around the star $\theta^{1}$ Ori C. This $\gamma$-ray emission mainly arises
due to the Inverse Compton scattering of the dominant optical radiation from
the $\theta^{1}$ Ori C and the strong infra-red radiation which appears as a
result of the reprocessing of the stellar optical radiation by the gas
surrounding the massive star.

We confront this predicted GeV-TeV $\gamma$-ray emission with the available
constraints on the $\gamma$-ray emission from this region of the sky in the
GeV and TeV energy range. This $\gamma$-ray emission cannot overcome the level
of the observed GeV $\gamma$-ray emission from the part of the Giant Molecular
Complex called Orion A (see e.g. Young et al. 2014). This emission is likely
due to the interaction of the galactic cosmic rays with the matter of the huge
Orion Nebula. Predicted by us $\gamma$-ray emission from electrons cannot also
overcome the upper limits on the TeV $\gamma$-ray emission from the Orion
Cluster (obtained by HAWC Collaboration, see Albert et al. 2021). Recently,
the upper limit on the GeV emission from the direction of the Trapezium
cluster has been also reported by Maurin et al. (2016). Based on the
comparison of the results of our calculations with these Fermi-LAT upper
limits, we conclude that the electron acceleration efficiency, i.e. the
conversion factor of the wind power of the star $\theta^{1}$ Ori C to the
relativistic electrons, cannot be larger than $\sim$10-4. This indicates that
electrons cannot be efficiently accelerated in the wind cavity around the
massive star $\theta^{1}$ Ori C. Predicted level of the TeV $\gamma$-ray
emission from the wind cavity around $\theta^{1}$ Ori C (Fig. 5) can stay
marginally within the sensitivity limits of the future Cherenkov observatories
(e.g. the Cherenkov Telescope Array).

The constraint on the efficiency of electron acceleration in the wind of
$\theta^{1}$ Ori C is clearly lower than the constraint on the acceleration
efficiency of hadrons. This constraint (Maurin et al. 2016), on the level of a
few percent, is derived based on a similar scenario but in terms of the
hadronic model. The difference in the constraints on the efficiency of hadron
and electron acceleration is due to the more efficient energy loses of
electrons versus hadrons in the environment around the star $\Theta^{1}$ Ori
C.

## Acknowledgments

I would like to thank the Referee for many valuable comments. This work is
supported by the grant through the Polish National Research Centre No.
2019/33/B/ST9/01904.

## Data Availability

The simulated data underlying this article will be shared on reasonable
request to the corresponding author.

## References

  * [1] Abdo, A.A. et al. 2007 ApJ 658, L33
  * [2] Abdo, A.A. et al. 2010 Science 329, 817
  * [3] Abeysekara, A.U., Albert, A., Alfaro, R. et al. 2021 Nature Astronomy 5, 465
  * [4] Abramowski, A. et al. 2012 A&A 537, A114
  * [5] Acciari, V.A. et al. 2022 Nature Astronomy, 6, 689
  * [6] Ackermann, M. et al. 2012 ApJ 756, 4
  * [7] Aharonian, F., Yang, R., de Oña Wilhelmi, E. 2019 Nature Astronomy 3, 561
  * [8] Aharonian, F. et al. 2005 A&A 431, 197
  * [9] Aharonian, F. et al. 2007 A&A 467, 1075
  * [10] Aharonian, F. et al. 2022 A&A 666, A124
  * [11] Albert, A. et al. 2021 ApJ 914, 106
  * [12] Baghmanyan et al. 2020 ApJL 901, L4
  * [13] Bednarek, W. 2005 MNRAS 363, L46
  * [14] Bednarek, W. 2021 MNRAS 507, 3292
  * [15] Benaglia, P., Romero, G.E. 2003 A&A 399, 1121
  * [16] Blasi, P., Morlino, G. 2023, MNRAS, in press (arXiv:2306.03762), https://doi.org/10.48550/arXiv.2306.03762
  * [17] Brose, R., Pohl, M., Sushch, I., Petruk, O., Kuzyo, T. 2020 A&A 634, 59
  * [18] Brose, R., Sushch, I., Mackey, J. 2022 MNRAS 516, 492
  * [19] Bykov, A.M. 2001 SSRv 99, 317
  * [20] Bykov, A.M. 2014 A&ARv 22, 77
  * [21] Cao, Z., Aharonian, F. A., An, Q., et al. 2021 Nature 594, 33
  * [22] Caraveo, P. et al. 1980 A&A 91, L3
  * [23] Castor, J., McCray, R., Weaver, R. 1975 ApJ 200, L107
  * [24] Celli, S., Morlino, G., Gabici, S., Aharonian, F. A. 2019 MNRAS 490, 4317
  * [25] Cesarsky, C.J., Montmerle, T. 1983 SSRv 36, 173
  * [26] Chevalier, R.A. 1997 ApJ 488, 263
  * [27] Digel, S.W., Hunter, S.D., Mukherjee, R. 1995 ApJ 441, 270
  * [28] Donati, J.-F., Babel, P.J., Harries, P.T.J., Howarth, P.I.D., Petit P., Semel M. 2002 MNRAS 333, 55
  * [29] Dwarkadas, V.V. 2022 Galaxies, 10, 37
  * [30] Eichler, D., Usov, V. 1993 ApJ 40, 271
  * [31] Güdel, M., Briggs, K.R., Montmerle, T., Audard, M., Rebull, L., Skinner, S.L. 2008 Science 319, 309
  * [32] H.E.S.S. Collaboration, Abramowski, A., Acero, F., et al. 2011 A&A 525, A46
  * [33] Howarth, I.D., Prinja, R.K. 1989 ApJS 69, 527
  * [34] Howarth, I.D., Siebert, K.W., Hussain, G.A.J., et al. 1997, MNRAS, 284, 265
  * [35] Inoue, T., Marcowith, A., Giacinti, G., Jan van Marle, A., Nishino, S. 2021 ApJ 2022, 7
  * [36] Leto, P. et al. 2017 MNRAS 457, 2820
  * [37] Maier, G. et al. (CTA Consortium) 2017, in Proc. 35th ICRC (Busan, Korea), PoS(ICRC2017), 846
  * [38] Marchili, N. et al. 2018 A&A 615, A82
  * [39] Marcowith, A., Dwarkadas, V.V., Renaud, M., Tatischeff, V., Giacinti, G. 2018 MNRAS 479, 4470
  * [40] Maurin, G. et al. 2016 A&A 591, A71
  * [41] Moderski, R., Sikora, M., Coppi, P.S., Aharonian, F. 2005 MNRAS 363, 954
  * [42] Neronov, A., Semikoz, D., Taylor, A. 2012 PRL 108, 051105
  * [43] Nieva, M.-F., Przybilla, N. 2014 A&A, 566, A7
  * [44] Pabst, C., Higgins, R., Goicoechea, J.R. et al. 2019 Nature 565, 618
  * [45] Peng, F.-K. et al. 2019A&A 621, A70
  * [46] Pittard, J.M., Dougherty, S.M. 2006 MNRAS 372, 801
  * [47] Reid, M.J. et al. 2009 ApJ 700, 137
  * [48] Reimer, A., Pohl, M.K., Reimer, O. 2006 ApJ 644, 1118
  * [49] Rzaev, A.Kh, Shimansky, V.V., Kolbin, A.I. 2021 MNRAS 504, 3787
  * [50] Simon-Diaz, S., Herrero, A., Esteban, C., Najarro, F., 2006, A&A 448, 351
  * [51] Stahl, O. et al. 1996 A&A 312, 539
  * [52] Torres, D.F., Domingo-Santamaria, E., Romero, G.E. 2004 ApJ 601, L75
  * [53] Trigilio, C. Leto, P., Umana, G., Leone, F., Buemi, C. 2004 A&A 418, 593
  * [54] Usov, V.V., Melrose, D.B. 1992 ApJ 395, 575
  * [55] Yang, R., de Oña Wlhelmi, E., Aharonian, F. 2014 A&A 566, A142
  * [56] Vieu, T., Gabici, S., Tatischeff, V., Ravikularaman, S. 2022 MNRAS 512, 1275
  * [57] Völk, H.J., Forman, M. 1982 ApJ 253, 188
  * [58] Weaver, R., McCray, R., Castor, J., Shapiro, P., Moore, R. 1977 ApJ 218, 377
  * [59] Webb, G.M., Axford, W.I., Forman, M. A. 1985 ApJ 298, 684
  * [60] White, R.L., Chen, W. 1992 ApJ 387, L81


On unitarity in singlet inflation with a non-minimal coupling to gravity

Oleg Lebedev${}^{\,a}$, Yann Mambrini${}^{\,b}$ and Jong-Hyun Yoon${}^{\,b}$

aDepartment of Physics and Helsinki Institute of Physics,
Gustaf Hällströmin katu 2a, FI-00014 Helsinki, Finland
bUniversité Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

Abstract

We study inflationary models based on a non-minimal coupling of a singlet
scalar to gravity, focussing on the preheating dynamics and the unitarity
issues in this regime. If the scalar does not have significant couplings to
other fields, particle production after inflation is far less efficient than
that in Higgs inflation. As a result, unitarity violation at large non-minimal
couplings requires a different treatment. We find that collective effects in
inflaton scattering processes during preheating make an important impact on
the unitarity constraint. Within effective field theory, the consequent upper
bound on the non-minimal coupling is of order a few hundreds.

###### Contents

  1. 1 Introduction
  2. 2 Inflaton background dynamics
    1. 2.1 Basics of inflation in the Einstein frame
    2. 2.2 Postinflationary background dynamics in the Jordan frame
      1. 2.2.1 Spike-like feature in the inflaton evolution
      2. 2.2.2 Non-adiabaticity
  3. 3 Postinflationary dynamics on the lattice and unitarity
    1. 3.1 Unitarity violation: simple estimate
    2. 3.2 Impact of collective effects on the unitarity bound
    3. 3.3 Effect of a significant Higgs-inflaton coupling
  4. 4 Conclusion

##  1 Introduction

Inflation is one of the cornerstones of modern cosmology [1, 2, 3, 4]. An
attractive class of inflationary models is based on a non-minimal scalar
coupling to curvature $\xi$ [5, 6, 7],

$\Delta{\cal L}={1\over 2}\xi R\,\phi^{2}\;,$ (1)

where $R$ is the scalar curvature and $\phi$ is the inflaton. A well-known
example is ‘‘Higgs inflation’’ [8]. An important feature of these models is
that they provide an excellent fit to the current PLANCK data [9].
Furthermore, the presence of the coupling $\xi$ is expected on general
grounds, e.g. it is generated radiatively, making the framework theoretically
well motivated.

One of the constraints on the size of the non-minimal coupling to gravity is
imposed by unitarity considerations [10, 11]. The inflationary models are
based on the effective field theory (EFT) description and thus are valid up to
a certain energy scale. After inflation, when the inflaton background has
become negligible, the corresponding EFT cut-off is of order [10, 11]

$\Lambda\sim{M_{\rm Pl}\over\xi}\;,$ (2)

while during inflation, the presence of the large scalar background raises it
to

$\Lambda_{\rm infl}\sim{M_{\rm Pl}}$ (3)

in the Einstein frame [12]. If $\xi\gg 1$, the effective description may be
problematic after inflation bringing consistency of such models in question.

In particular, after Higgs inflation, a ‘‘violent’’ reheating stage takes
place, which is driven mainly by the Higgs couplings to the gauge bosons. Much
of the energy density of the system is converted into energetic gauge bosons
with momenta close to the Planck scale [13]. Clearly, this is problematic
since the effective field theory cut-off during reheating is far below the
Planck scale. This phenomenon was explored further in [14].

In our current work, we study the unitarity issues in singlet scalar inflation
driven by $\xi\gg 1$, where the singlet does not have significant couplings to
other fields, unlike the Higgs. In this case, postinflationary particle
production is far less efficient than that in Higgs inflation and the
corresponding unitarity considerations do not apply. Nevertheless, it is clear
that very large $\xi$ would be problematic: due to inflaton quanta production,
the energy density after inflation decreases in time slower than the inflaton
background does, thus violating the bound (2) eventually. We quantify this
statement using lattice simulations. More importantly, we evaluate the impact
of inflaton collective effects on scattering amplitudes during preheating and
find that these make the unitarity bound substantially stronger. Finally, we
study the effect of a significant non-derivative inflaton coupling to another
field. Our main result is that the unitarity bound requires $\xi$ to be below
a few hundreds, which is important for inflationary model building.

We confine ourselves to unitarity studies in Higgs-like inflationary models in
the metric formulation. Analogous unitarity bounds become much looser in
Palatini formalism [15], although this set-up is not as minimalistic. The
original Higgs inflation model can also be unitarized, i.e. made UV-complete,
for the price of Higgs inflation being a ‘‘mirage’’ in the sense that the
inflationary dynamics gets dominated by a different field [16, 17, 18, 19].
Other relevant studies of unitarity issues during and after Higgs inflation
can be found in [20, 21, 22, 23]. In particular, [22] focuses on explicit
calculation of the scattering amplitudes, while [23] discusses the relevant
field redefinitions which make the computations transparent. Numerical studies
of the preheating dynamics in Higgs-like inflation models have been performed
in [24, 25, 26, 27].

##  2 Inflaton background dynamics

We start with the inflaton action in the Jordan frame,

${\cal S}=\int d^{4}x\sqrt{-g}\left({1\over 2}M_{\rm Pl}^{2}R+{1\over 2}\xi
R\,\phi^{2}-{1\over
2}g^{\mu\nu}\partial_{\mu}\phi\,\partial_{\nu}\phi-V\right)\;,$ (4)

where the metric convention is $-+++$ and $\xi>0$.111The $\xi$ convention
differs from that of our previous work [28]. The potential is

$\displaystyle V={1\over 4}\lambda_{\phi}\phi^{4}\;,$ (5)

where we have neglected the inflaton mass term irrelevant at large field
values. In what follows, we will use the Planck units

$M_{\rm Pl}=1\;.$ (6)

and consider a large $\xi$ regime,

$\xi\gg 1\;.$ (7)

###  2.1 Basics of inflation in the Einstein frame

The transition to the Einstein frame, where the curvature-dependent term is
${1\over 2}R$, is accomplished by the metric transformation [29]

$g_{\mu\nu}^{\rm E}=\Omega\,g_{\mu\nu}\;,$ (8)

with

$\Omega=1+\xi\phi^{2}\;.$ (9)

This eliminates the scalar coupling to curvature for the price of a non-
canonical $\phi$ kinetic term. The canonically normalized inflaton $\chi$
satisfies ${d\chi\over
d\phi}=\sqrt{1+\xi(1+6\xi)\phi^{2}\over(1+\xi\phi^{2})^{2}}\;,$ which is
solved by [30]

$\,\chi(\phi)=\sqrt{1+6\xi\over\xi}\,\sinh^{-1}\left(\sqrt{(1+6\xi)\xi}\,\phi\right)-\sqrt{6}\,\sinh^{-1}\left({\sqrt{6}\xi\phi\over\sqrt{1+\xi\phi^{2}}}\right)\;.$
(10)

Inflation takes place when $\xi\phi^{2}\gg 1$, in which case
$\chi\simeq{\sqrt{3/2}}\,\ln\xi\phi^{2}$ and the Einstein frame potential
$V^{\rm E}=V/\Omega^{2}$ has the form

$V_{E}(\chi)\simeq{\lambda_{\phi}\over 4\xi^{2}}\,\left(1-e^{-\sqrt{2\over
3}|\chi|}\right)^{2}\;.$ (11)

It is exponentially close to a flat potential and provides an excellent fit to
the PLANCK data [8, 9]. The energy density during inflation can be
approximated by $V^{\rm E}=3H^{2}\simeq{\lambda_{\phi}\over 4\xi^{2}}\;.$
Computing the $\epsilon$ parameter from the above potential and applying the
COBE normalization [8], one gets

${\lambda_{\phi}\over 4\xi^{2}}=4\times 10^{-7}\,{1\over N^{2}}\;,$ (12)

where $N$ is the number of $e$-folds of inflation. For $N=60$,
${\lambda_{\phi}\over 4\xi^{2}}\simeq 10^{-10}$. This constraint implies that
there is only one independent variable in our analysis, e.g. $\xi$.

Inflation ends when $\chi\sim 1$. After that, the potential becomes quadratic
and, subsequently, quartic [30]:

$V_{E}(\chi)\simeq\left\\{\begin{tabular}[]{ccc}${\lambda_{\phi}\over
6\xi^{2}}\chi^{2}$&for&${1\over 2\xi}\ll|\chi|\ll 1~{},$\\\
${\lambda_{\phi}\over 4}\chi^{4}$&for&$|\chi|\ll{1\over
2\xi}~{}.$\end{tabular}\right.$ (13)

In our work, we focus primarily on these regimes.

###  2.2 Postinflationary background dynamics in the Jordan frame

To study the inflaton background dynamics after inflation, it is often
convenient to use the Jordan frame. This also applies to lattice simulations
of postinflationary dynamics [31] since the equations of motion are easier to
handle in the Jordan frame. In this section, we analyze how the system behaves
in terms of two field variables: the homogeneous inflaton field $\phi(t)$ and
the scale factor $a(t)$.

In the Friedmann metric

$ds^{2}=-dt^{2}+a(t)^{2}\,dx^{i}\,dx^{i}\;,$ (14)

the scalar curvature has the form

$R=6\left[{\ddot{a}\over a}+\left({\dot{a}\over a}\right)^{2}\right]\;,$ (15)

where the dot denotes the time derivative. As usual, the Hubble rate is
defined by $H=\dot{a}/a$. The individual components of the Ricci tensor
$R_{\mu\nu}$ can be found, for instance, in [31]. The system satisfies the
Einstein equation $R_{\mu\nu}-{1\over 2}g_{\mu\nu}R=\,T_{\mu\nu}\,,$ where
$T_{\mu\nu}$ is the scalar energy-momentum tensor. Due to the non-minimal
coupling to gravity, it contains a curvature-dependent contribution.

Since there are two field variables in the system, it is sufficient to use two
independent equations of motion (EOM). The scalar EOM for the homogeneous
background $\phi$ and the Einstein equation $R_{00}-{1\over 2}g_{00}R=T_{00}$
can be written as

$\displaystyle\ddot{\phi}+3H\,\dot{\phi}-{1\over a^{2}}\nabla^{2}\phi-\xi
R\phi+{\partial V\over\partial\phi}=0\;,$ (16) $\displaystyle
3H^{2}=\rho(\phi)\;.$ (17)

These equations can be recast in a more convenient form. Tracing over the
indices of the Einstein equation, $R=-T_{\mu}^{\mu}\;,$ and computing
$T_{00}(\phi)=\rho(\phi)$ for the inflaton, one finds [31, 28]

$\displaystyle
R=\frac{1}{1+(6\xi+1)\xi\,\phi^{2}}\,\left[(1+6\xi)\,\partial^{\mu}\phi\partial_{\mu}\phi+4V+6\xi\,\phi
V^{\prime}_{\phi}\right]\;,$ (18)
$\displaystyle\rho(\phi)=\frac{1}{2}\dot{\phi}^{2}+V(\phi)-3\xi
H^{2}\phi^{2}-6\xi H\phi\dot{\phi}\;,$ (19)

accounting for the difference between our definition of $\xi$ and that of
[28]. Note the unusual presence of the Hubble-dependent terms in the scalar
energy density, which appear due to the kinetic scalar-graviton mixing in the
Jordan frame (see the discussion in [28]).

Plugging the expression for $R$ in the scalar EOM, we get for $6\xi\gg 1$,

$\ddot{\phi}+3H\,\dot{\phi}+{6\xi^{2}\dot{\phi}^{2}+\lambda_{\phi}\phi^{2}\over
1+6\xi^{2}\phi^{2}}\;\phi=0\;,$ (20)

such that ${6\xi^{2}\dot{\phi}^{2}+\lambda_{\phi}\phi^{2}\over
1+6\xi^{2}\phi^{2}}$ can be viewed as the effective mass squared. We observe
that there is no tachyonic instability and the effective mass squared at large
$\phi$ is $\lambda_{\phi}/(6\xi^{2})$. To solve the scalar EOM, we need an
expression for $H$ in terms of $\phi$. Denoting the usual $\xi=0$ Hubble rate
by $H_{0}$,

$3H_{0}^{2}={1\over 2}\dot{\phi}^{2}+V\;,$ (21)

one solves the quadratic equation (17) to get

$H=\sqrt{{H_{0}^{2}\over 1+\xi\phi^{2}}+\left({\xi\phi\dot{\phi}\over
1+\xi\phi^{2}}\right)^{2}}-{\xi\phi\dot{\phi}\over 1+\xi\phi^{2}}\;.$ (22)

Plugging this into (20), we obtain an equation just for $\phi$, which can be
solved numerically. Note that at large $\phi$, $H=H_{0}/(\sqrt{\xi}\phi)$
which is not constant during inflation, unlike the Hubble rate in the Einstein
frame.

Figure 1:  The inflaton and Hubble rate evolution in the Jordan frame for
$\xi=10$. The scale factor is normalized to 1 when $R=6H^{2}$ ($\ddot{a}=0$)
holds true numerically.

The inflationary solution is $\ddot{\phi}\simeq 0$ such that $\phi(t)$ is
linear, with a small velocity
$\dot{\phi}\sim\sqrt{\lambda_{\phi}}\,\xi^{-3/2}$. A numerical example of the
field and Hubble rate evolution in the Jordan frame is shown in Fig. 1.222 The
$\lambda_{\phi}$ dependence can be eliminated from the scalar EOM by
introducing a rescaled time variable $t^{\prime}=\sqrt{\lambda_{\phi}}\,t$. At
large field values, both $\phi$ and $H$ are linear in time. Inflation ends
when $\phi$ becomes of order

$\phi\sim{1\over\sqrt{\xi}}\;,$ (23)

such that $1+\xi\phi^{2}$ can no longer be approximated by $\xi\phi^{2}$. At
this point, the velocity contribution to the effective inflaton mass in (20)
becomes important, $\xi^{2}\dot{\phi}^{2}\sim\lambda_{\phi}\phi^{2}$.
Similarly, the velocity terms in (22) start playing a significant role. The
energy density at the end of inflation is of order
$\rho\sim\lambda_{\phi}/\xi^{2}$. Finally, when $\phi$ reaches zero, it starts
oscillating inducing oscillations in the Hubble rate as well.

####  2.2.1 Spike-like feature in the inflaton evolution

Ema et al. [13] have observed that the inflaton field derivative exhibits a
spike-like feature relevant to particle production. Such a spike is caused by
the shape of the inflaton potential, as viewed in the Einstein frame, because
it combines a quadratic potential at larger field values with a quartic
potential at smaller inflaton values. An earlier study of preheating in this
system can be found in [32].

Let us examine the appearance of the spike using our equations of motion. To
understand the essence of the mechanism, one may neglect the Universe
expansion. We then have

$\ddot{\phi}+\omega^{2}(\phi)\,\phi=0~{}~{},~{}~{}\omega^{2}={6\xi^{2}\dot{\phi}^{2}+\lambda_{\phi}\phi^{2}\over
1+6\xi^{2}\phi^{2}}\;.$ (24)

The quantity $\omega$ has a different behaviour at large and small field
values. Consider first the larger $\phi$ values.

(a) $\xi^{2}\phi^{2}\gg 1$.

In this regime,
$\omega^{2}\simeq\lambda_{\phi}/(6\xi^{2})+(\dot{\phi}/\phi)^{2}$. Let us use
the Ansatz $\phi=\phi_{0}\,\cos\omega t$, where $\phi_{0}$ and $\omega$ are
treated as approximately constant. As long as the field value is large enough,
i.e. $\cos^{2}\omega t\gg\sin^{2}\omega t$, the oscillation frequency is
indeed almost constant and

$\omega^{2}\simeq{\lambda_{\phi}\over 6\xi^{2}}\;.$ (25)

This expression therefore gives the frequency of large-amplitude oscillations.

Figure 2:  The inflaton field and its velocity evolution in the Jordan frame
at $\xi=10^{4}$. The scale factor is normalized to 1 when $R=6H^{2}$
($\ddot{a}=0$) holds true numerically.

(b) $\xi^{2}\phi^{2}\ll 1$.

Consider the small field limit $\phi\rightarrow 0$. In this case,

$\omega^{2}\simeq 6\xi^{2}\dot{\phi}^{2}$ (26)

and $\omega$ cannot a priori be treated as a slowly varying function. The EOM
becomes

$\ddot{\phi}+6\xi^{2}\dot{\phi}^{2}\,\phi=0\;.$ (27)

This is easily solved for $\dot{\phi}$,

$\dot{\phi}={\rm const}\times e^{-3\xi^{2}\phi^{2}}\;,$ (28)

where the constant is determined by the boundary conditions. In the small
field regime, the exponent is close to one, so $\dot{\phi}\simeq{\rm const}$.
The size of the constant can be evaluated using approximate energy
conservation. The energy density at the end of inflation,
$\lambda_{\phi}/\xi^{2}$, gets converted into the kinetic energy at $\phi=0$,
so $\dot{\phi}^{2}\sim\lambda_{\phi}/\xi^{2}$. Hence

$\dot{\phi}\sim\sqrt{\lambda_{\phi}}/\xi$ (29)

and

$\omega\sim\sqrt{\lambda_{\phi}}\,.$ (30)

This corresponds to a much sharper variation of $\phi$ compared to the large
field case, by a factor of $\xi$. As $|\phi|$ increases away from the zero
crossing, the velocity becomes exponentially suppressed according to (28),
which generates a spike in $\dot{\phi}$.

The Hubble friction does not change this behaviour fundamentally, so the full
system exhibits similar features. The main point is that the frequency
$\omega$ characterizing the variation of $\phi$ in time changes abruptly
around inflaton zero crossings $\phi=0$, which indicates non-adiabaticity. A
numerical example illustrating the presence of spike-like features in shown in
Fig. 2. Due to the logarithmic scale, there also appear ‘‘negative’’ spikes
corresponding to $\dot{\phi}\rightarrow 0$.

####  2.2.2 Non-adiabaticity

Semiclassical particle production is characterized by the non-adiabaticity
parameter $\dot{\omega}/\omega^{2}$. The oscillation frequency $\omega$
changes abruptly at the border of the regions considered above, so let us
focus on $\xi|\phi|\lesssim 1$. In this case,

$\omega\simeq{\sqrt{6}\xi|\dot{\phi}|\over\sqrt{1+6\xi^{2}\phi^{2}}}\;.$ (31)

This is because $\lambda_{\phi}\phi^{2}\ll\xi^{2}\dot{\phi}^{2}$ in this
region. Indeed, approximate energy conservation requires
$\lambda_{\phi}/\xi^{2}\sim\dot{\phi}^{2}+\lambda_{\phi}\phi^{4}$, while the
potential term is negligible at $\phi\sim 1/\xi$.

Figure 3:  The non-adiabaticity parameter for $\xi=10^{4}$. The scale factor
is normalized to 1 when $R=6H^{2}$ ($\ddot{a}=0$) holds true numerically.

One then finds for $\dot{\phi}>0$,

${\dot{\omega}\over\omega^{2}}\simeq-2\,{\sqrt{6}\xi\phi\over\sqrt{1+6\xi^{2}\phi^{2}}}\;.$
(32)

The second factor in this expression is bounded by one, so
$|\dot{\omega}/\omega^{2}|>1$ at $\xi\phi\sim 1$ with the maximal value of 2.
This behaviour is seen in Fig. 3.

These considerations can be extended to modes with non-zero momenta by
replacing

$\omega^{2}\rightarrow{k^{2}\over a^{2}}+\omega^{2}\;.$ (33)

As long as $k^{2}/a^{2}\lesssim\omega^{2}$, the above conclusions apply. In
the non-adiabatic region, $\omega\lesssim\sqrt{\lambda_{\phi}}$. Therefore,
the momentum cut-off for the produced particles is of order

${k\over a}\lesssim\sqrt{\lambda_{\phi}}\;.$ (34)

##  3 Postinflationary dynamics on the lattice and unitarity

The dynamics of the inflaton field after inflation are complicated by
backreaction and rescattering effects. These lead to fast decay of the zero
mode and growth of inhomogeneities, both of which have an impact on unitarity
violation considerations. To take such effects into account, we resort to
lattice simulations using the tool CosmoLattice [33, 34].

Figure 4:  Evolution of the energy density scale $\rho^{1/4}$ in Planck units.
Left: $\xi=100$, $\lambda_{\phi}=4.4\times 10^{-6}$. Right: $\xi=2000$,
$\lambda_{\phi}=1.8\times 10^{-3}$.

One of the features of CosmoLattice is a rearrangement of the second-order
equations into the first-order ones. However, in the presence of the non-
minimal coupling, this becomes problematic since the pressure and energy
density terms in the second-order Friedmann equation have an explicit
dependence on the scale factor. Instead, we use a customized version of
CosmoLattice as in [31], which solves for the scale factor with the help of
Eq. 18 and an explicit Runge-Kutta algorithm. We employ a lattice with 320
sites along each axis, which is sufficient to capture the characteristic power
spectra. We observe good convergence of the results and consistency with other
studies such as [31], when applicable.

The initial conditions for the simulation are chosen as follows: at the end of
inflation, $\phi\sim 1/\sqrt{\xi}$ in Planck units, we neglect $\dot{\phi}$
and solve Eq. 20 to find $\phi,\dot{\phi}$ at later times.333We find that the
results are insensitive to the initial value of $\dot{\phi}$ as long as it is
small. At the same time, we evolve the field fluctuations with a Python code
to set up the boundary conditions for the CosmoLattice simulations. For the
latter, we use $\phi=0.01$ as the initial value ($a=1$) and adopt the
normalization (12) with $N=60$. This implies, for example, that at $a=1$,
$\rho^{1/4}=0.0023\,$ for $\xi=100$ and $\rho^{1/4}=0.00247\,$ for $\xi=500$.
At this stage, we set vacuum fluctuations as the initial conditions for the
inflaton momentum modes and evolve these with the help of lattice simulations.
We find that using inflationary fluctuations instead does not make a visible
difference since inflation affects superhorizon modes $k/a\lesssim H$, while
the hard modes are much more important for our purposes (see also [31]).

Given the shape of the inflaton potential, one expects the energy density to
scale as non-relativistic matter initially, while at later times, the scaling
becomes radiation-like. In practice, the field fluctuations play an important
role such that the truly non-relativistic regime shortens significantly and
most of the time at $\phi>1/\xi$ the system behaves as semirelativistic, i.e.
exhibits a scaling law between $a^{-3}$ and $a^{-4}$. For $\phi\lesssim
1/\xi$, the system becomes radiation-like, so

$\displaystyle\rho\propto
a^{-(3+x)}~{}~{}~{},~{}~{}~{}1/\xi\lesssim\phi\lesssim 1/\sqrt{\xi}~{}~{},$
$\displaystyle\rho\propto a^{-4}~{}~{}~{},~{}~{}~{}\phi\lesssim 1/\xi~{}~{},$
(35)

where $0<x<1$. Numerical examples of the energy density evolution are shown in
Fig. 4. The apparent strong oscillatory behaviour is due to lattice
limitations, yet the above scaling applies to $\rho$ averaged over a few
oscillation periods.

Figure 5:  Evolution of the inflaton momentum spectrum (red to blue). Left:
$\xi=100$, $a_{\rm end}=300$. Right: $\xi=500$, $a_{\rm end}=70$. $n_{k}$ is
the occupation number and the comoving momentum $k$ is in units of
$\sqrt{\lambda_{\phi}}$, where $\lambda_{\phi}$ is fixed by the COBE
normalization.

The oscillating inflaton background evolves non-adiabatically (Sec. 2.2.2) and
resonantly excites non-zero momentum modes via an analog of Eq. 20.444The EOM
for the different momentum modes do not decouple and cannot be written in a
simple way, e.g. due to momentum-dependent contributions to the curvature.
Instead, one solves the system on the lattice in coordinate space as detailed
in [31, 28] and then decomposes the solution in the Fourier modes. The process
is highly non-linear and can be studied reliably only on the lattice.
Expanding $\phi(x)$ in 3-momentum modes $\phi_{k}$ as in [35], it is
convenient to define a rescaled variable $Y_{k}\equiv a\,\phi_{k}$ as a
function of conformal time $\tau$: $d\tau=dt/a$. Then, the oscillator number
$n_{k}$ and the corresponding frequency $\omega_{k}$ are defined by

$\displaystyle
n_{k}\equiv\frac{1}{2}\left(\omega_{k}|Y_{k}|^{2}+\frac{1}{\omega_{k}}|\dot{Y}_{k}|^{2}\right)\;,$
(36) $\displaystyle\omega^{2}_{k}\equiv
k^{2}+a^{2}\left\langle\frac{\partial^{2}V(\phi)}{\partial\phi^{2}}\right\rangle\;,$
(37)

where $\langle...\rangle$ denotes spacial averaging. $n_{k}$ changes in time
signifying particle production and rescattering. Evolution of the (comoving)
momentum spectrum for two representative examples $\xi=100$ and $\xi=500$ is
shown in Fig. 5. One observes the momentum cut-off of order
$\sqrt{\lambda_{\phi}}$, as expected. The typical occupation numbers are very
large, $n_{k}\gg 1$, manifesting the importance of collective effects. At late
times, one also observes an isolated peak in the spectrum. It is due to the
momentum modes excited in the $\phi^{4}$ potential via the Lamé equation [35].
The magnitude of this momentum band is given approximately by

$k\sim\sqrt{\lambda_{\phi}}\,\phi_{*}\;,$ (38)

where $\phi_{*}\sim 1/\xi$ corresponds to the onset of the quartic potential
regime. This only gives a rough estimate since the velocity of the field at
$\phi_{*}$ is non-zero, unlike in the analysis of [35], and $\phi_{*}$ is not
precisely defined. Finally, the bump is expected to be smoothed out on a
larger time scale due to rescattering.

An important ingredient in our study is the evolution of the inflaton zero
mode. The inflaton background affects the unitarity bound, hence it is
essential to determine its decay time. Fig. 6 displays the evolution of the
background $|\langle\phi\rangle|$ and the fluctuation amplitude
$\sqrt{\langle\phi^{2}\rangle-\langle\phi\rangle^{2}}$ for $\xi=100,\,2000$.
We observe that at large $\xi$, the fluctuations become larger than the
background when the zero mode is of the size

$\phi_{*}\sim{1\over\xi}\;.$ (39)

This inflaton value signifies the transition to the quartic potential regime.
From this point on, the system becomes fully dominated by fluctuations and the
background can be neglected for most purposes. (For small $\xi\sim 5$, these
conclusions do not apply and the fluctuations are much suppressed.)

Figure 6:  Evolution of the inflaton zero mode and the inflaton variance.
Left: $\xi=100$, $\lambda_{\phi}=4.4\times 10^{-6}$. Right: $\xi=2000$,
$\lambda_{\phi}=1.8\times 10^{-3}$.

###  3.1 Unitarity violation: simple estimate

When the inflaton background value is negligible, a ballpark estimate of the
unitarity bound on $\xi$ can be obtained by analyzing the energy density of
the system. Since $\xi\langle\phi^{2}\rangle\ll 1$, the difference between the
Jordan and Einstein frames is insignificant in terms of the energy scales,
hence one can use either frame to estimate the unitarity bound.

The canonically normalized inflaton variable $\chi$ in the $Einstein$ frame is
defined by the differential equation

${d\chi\over
d\phi}=\sqrt{1+\xi_{\phi}(1+6\xi_{\phi})\phi^{2}\over(1+\xi_{\phi}\phi^{2})^{2}}\;.$
(40)

The solution is given by (10), although it is often more convenient to work
directly with (40) in specific limits.

Since the kinetic and the curvature terms are canonical, all the non-trivial
physics resides in the scalar potential
$V=\lambda_{\phi}\phi^{4}/\left[4(1+\xi\phi^{2})^{2}\right]$ formulated in
terms of $\chi$. Given the inflaton background value at the particular
evolution stage, one expands the potential in terms of fluctuations over the
background. Perturbative unitarity within the effective field theory
description then imposes a constraint on the energy scale of the fluctuations.
During inflation, the expansion is of the form $(\delta\chi)^{n}$ so that the
unitarity violation scale is $M_{\rm Pl}$ [12]. The same bound also applies at
the end of inflation when $\chi\sim\xi\phi^{2}\sim 1$. This can be shown by
examining (40) in the vicinity of $\xi\phi^{2}\sim 1$. However, at very small
$\phi\ll 1/\xi$ and large $\xi\gg 1$,

$\phi\simeq\chi\,(1-\xi^{2}\chi^{2})\;.$ (41)

The expansion around the vacuum then has the form $(\xi\chi)^{n}$, which
implies that the perturbative scattering amplitudes based on $V(\chi)$ blow up
at the energy scale $1/\xi$ .555The size of the coupling $\lambda_{\phi}$ does
not matter for a sufficiently large $n$. The characteristic energy scale for
the scattering processes can be taken to be $\rho^{1/4}$ at the time when the
background becomes negligible, so the unitarity limit corresponds to

$\rho^{1/4}(\phi_{*})\sim{1\over\xi}\;.$ (42)

The resulting bound on $\xi$ can be estimated by
$\left({\lambda_{\phi}/4\xi^{2}}\right)^{1/4}\times{a_{*}^{-3/4}}\sim{1/\xi}$,
where $a_{*}$ is the scale factor corresponding to the decay of the zero mode
and we have assumed non-relativistic scaling of the energy density before
$a_{*}$. Since typically $a_{*}$ is between 5 and 10, this gives the critical
value of $\xi$ around 1000. A more careful numerical analysis yields

$\xi_{\rm max}^{0}\simeq 2000\;.$ (43)

Here the superscript $0$ serves to emphasise that this result is based on
simple dimensional analysis.

###  3.2 Impact of collective effects on the unitarity bound

After inflation, the inflaton field exhibits complicated dynamics where
collective effects are important. The system can be treated as a collection
inflaton quanta with given occupation numbers $n_{k}$, which can be very
large. On the other hand, the typical momentum can be much below the naive
estimate $\rho^{1/4}$. In what follows, we estimate the impact of collective
effects on the unitarity bound within the effective field theory approach.

To be specific, let us study the Higgs quanta production by the inflaton
field. Since the system must be reheated at some stage, consider a small
trilinear coupling in the Jordan frame,

$\Delta V={\sigma_{\phi h}}\,\phi\,H^{\dagger}H\;,$ (44)

which leads to inflaton decay at late times. In the Einstein frame, the
corresponding interaction is $\Delta V/\left[(1+\xi\phi^{2})^{2}\right]$
written in terms of $\chi\ll 1$. We may choose ${\sigma_{\phi h}}$ to be so
small that it does not affect the preheating dynamics. When the inflaton
background becomes negligible, one can use (41) or, more precisely, a high
order expansion of (10) around zero, which gives $\phi$ as a series in
$(\xi\chi)^{n}$ plus subleading terms. As a result, we obtain a set of
interactions

${\sigma_{\phi h}}\,H^{\dagger}H\;{\chi^{n}\over\Lambda^{n-1}}~{},$ (45)

where

$\Lambda\equiv{1\over\xi}\;.$ (46)

Such interactions induce, in particular, $n\rightarrow 2$ Higgs production
processes. Their efficiency depends on the mode occupation numbers or, in
other words, momentum distribution function. For highly occupied modes, the
corresponding amplitude can become large violating perturbative unitarity.

To make simple estimates, let us approximate the inflaton momentum
distribution function by a step function,

$f(p)=f\times\theta(p_{\rm max}-|p|)\;,$ (47)

where $f$ is a constant and $p_{\rm max}$ corresponds to the ‘‘maximal’’
3-momentum. The inflaton number density and the energy density are defined by
$n=\int{d^{3}p\over(2\pi)^{3}}\,f(p)$ and
$\rho=\int{d^{3}p\over(2\pi)^{3}}\,E_{p}\,f(p)\;,$ respectively. In practice,
the inflaton quanta can be treated as relativistic, $E_{p}\simeq|p|$ since the
induced inflaton mass $\sqrt{\lambda_{\phi}\langle\phi^{2}\rangle}$ is smaller
than the typical momenta. Then,

$n={1\over 6\pi^{2}}\;p_{\rm max}^{3}\,f~{}~{},~{}~{}\rho={1\over
8\pi^{2}}\;p_{\rm max}^{4}\,f\;.$ (48)

We observe that for $f\gg 1$, the characteristic momentum $p_{\rm max}$ is far
below $\rho^{1/4}$.

The probability of the Higgs pair production with fixed momenta via (45) is
given by

$\rho(n\rightarrow 2)={1\over
n!}\int\left(\prod_{i=1}^{n}{d^{3}k_{i}\over(2\pi)^{3}2E_{i}}\,f(k_{i})\,\right)(2\pi)^{4}\delta^{(4)}\left(\sum_{i}K_{i}-\sum_{j}P_{j}\right)\;|{\cal
M}(n\rightarrow 2)|^{2}\;,$ (49)

where $K_{i}$ and $P_{j}$ are the initial and final 4-momenta, respectively;
${\cal M}(n\rightarrow 2)$ is the standard QFT transition amplitude, and $n!$
accounts for identical particles in the initial state. Let us normalize the
coefficient of operator (45) to $1/n!$ to avoid large combinatorial factors in
the amplitude. Then, performing the integrals and using the Stirling
approximation for $n!$, one finds that for sufficiently large $n$, the
transition probability grows as

$\rho(n\rightarrow 2)\propto\left(c_{n}\,{p_{\rm
max}\sqrt{f}\over\Lambda}\right)^{2n}\;,$ (50)

with $c_{n}=\sqrt{e/(8\pi^{2}n)}$. The prefactor in (50) is unimportant for
large enough $n$ so that the precise value of $\sigma_{\phi h}$ does not play
any role. If the factor in the parentheses is significantly larger than one,
the probability grows uncontrollably violating unitarity. We note that the
precise value of $c_{n}$ depends on the coefficient of operator (45): for
instance, if it is order 1 instead of $1/n!$, the $\sqrt{n}$ factor in $c_{n}$
would appear in the numerator instead of the denominator. Given this
ambiguity, in what follows we will simply assume $c_{n}\sim{\cal
O}(1)$.666Although the coefficients of higher dimensional operators can be
computed, their signs alternate. Given that our system is a superposition of
states with fixed particle numbers (“squeezed” state), the collective impact
of these operators is difficult to estimate, hence we simply focus on a single
term assuming $c_{n}\sim 1$.

We conclude that the figure of merit for the effective field theory expansion
is roughly

$\kappa\sim{p_{\rm max}\sqrt{f}\over\Lambda}$ (51)

for moderately large $n$, and unitarity requires

$\kappa\lesssim 1\;.$ (52)

Hence, highly occupied states can cause problems even for low characteristic
energy of the inflaton quanta. To give an example, we find that for $\xi=500$
at $a\sim 6$ corresponding to the decay of the inflaton zero mode,777We obtain
these numbers from the lattice output by parametrizing the physical energy and
particle densities in terms of $p_{\rm max}$ and $f$. These quantities are
obtained in the Jordan frame, however, at this stage $\phi\sim\chi$ and the
difference between the frames is insignificant for our purposes.

$p_{\rm max}\sim 10^{-2}\sqrt{\lambda_{\phi}}~{}~{},~{}~{}f\sim 10^{4}\;,$
(53)

which makes $\kappa\sim 5$. The inflaton effective mass at this stage is of
order
$\sqrt{\lambda_{\phi}\langle\phi^{2}\rangle}\sim\sqrt{\lambda_{\phi}}/\xi$,
which is below $p_{\rm max}$ justifying our relativistic approximation. For
larger $\xi$, $\kappa$ grows further: the main factors are $p_{\rm
max}\propto\sqrt{\lambda_{\phi}}\propto\xi$ due to the inflationary constraint
and $1/\Lambda\propto\xi$, making $\kappa\propto\xi^{2}$. Even though
$\sqrt{f}$ decreases tending to ${\cal O}(1)$ at very large $\xi$ and the
decay time for the zero mode increases slowly with $\xi$ (Fig. 6), these
factors do not overcome the $\xi^{2}$ growth. We conclude that unitarity is
violated at

$\xi_{\rm max}\sim{\rm few}\,\times\,100\;$ (54)

according to the above criterion, keeping in mind the uncertainty associated
with $c_{n}$.

We see that collective effects lead to a stronger unitarity bound. This
tendency can be understood as follows. If we impose a bound on the energy
density of the system, this constrains $p_{\rm max}\,f^{1/4}$ according to
(48). On the other hand, the combination that appears in scattering processes
is $p_{\rm max}\,f^{1/2}$. Thus, we have

$p_{\rm max}\,f^{1/2}\gg p_{\rm max}\,f^{1/4}\,,$ (55)

and the scattering bound is stronger by the factor $f^{1/4}\gg 1$. For small
occupation numbers, on the other hand, the results are similar. Below we
illustrate the above statement with a simple example.

Example. Let us apply an analogous analysis to Higgs production by an inflaton
background in the $\phi^{2}$ potential. The process can be viewed as
annihilation of the non-relativistic inflaton quanta. Suppose $V={1\over
2}m^{2}\phi^{2}$ and the interaction term is $\phi^{4}h^{2}/\Lambda^{2}$. The
momentum distribution function for the inflaton can be written as
$f(p)=(2\pi)^{3}\delta^{(3)}(p)\,n$, where $n$ is the number density. Although
this system is non-relativistic, $n$ and $\rho$ can be expressed in terms of
the parameters of our box-shaped distribution function. Since $\rho=mn$, one
finds $p_{\rm max}=4m/3$ which plays the role of the characteristic energy
scale. The average occupation number is then found via
$f=8\pi^{2}\,\rho/p_{\rm max}^{4}$, which yields $f\sim 15\,\phi^{2}/m^{2}\gg
1$. In the case at hand, we can also take into account the effect of $c_{n}$:
with the above normalization, the analog of $c_{n}$ is about $1/7$ such that
$c_{n}\;{p_{\rm max}\sqrt{f}}/\Lambda$ becomes unity at $\phi\sim\Lambda$,
which signals unitarity violation. Therefore, we obtain the expected result.
The cutoff $\Lambda$ is far above the typical energy of the inflaton quantum
$m$ and the large occupation number plays a crucial role. Note that requiring
$\rho^{1/4}$ of the inflaton field to be below the cut-off would impose a very
different constraint $\phi<\Lambda^{2}/m$.

Similar considerations apply to inflaton self-interaction
$\chi^{l+4}/\Lambda^{l}$. The corresponding $n\rightarrow 2$ reaction rate is
enhanced by a further factor $f^{2}$ from the final state, which is however
inconsequential for a sufficiently large $l$. Thus, one again obtains the
unitarity bound (54). Note also that, if the bound is violated, the system
becomes sensitive to arbitrarily small couplings to new fields, which signals
pathology, at least within effective field theory.

It is important to understand the limitations of our analysis. We have relied
on perturbative description within effective field theory, yet violation of
perturbative unitarity does not always imply fundamental problems [36, 37].
Also, we have dealt with $n$-particle states, whereas the system in reality is
a complicated superposition of such states. However, we believe that our
analysis captures the main physical features and results in a reasonable
estimate of $\xi_{\rm max}$, above which the effective field theory
description becomes inadequate.

###  3.3 Effect of a significant Higgs-inflaton coupling

Figure 7:  Evolution of the inflaton zero mode and the variances in the
presence of the inflaton-Higgs coupling $\lambda_{\phi h}=0.5$ at $\xi=1000$
with 4 Higgs d.o.f.

A significant inflaton coupling to other fields, e.g. the Higgs boson [38,
39], leads to more efficient particle production and quicker background decay.
Fig. 7 shows an example for $\xi=1000$ with $\lambda_{\phi h}=0.5$, where the
coupling is defined by

$\Delta V={1\over 2}{\lambda_{\phi h}}\,\phi^{2}\,H^{\dagger}H\;.$ (56)

The zero mode decays by $a\sim 3$, which makes its lifetime shorter by a
factor of 2-3 compared to the zero coupling case. This occurs due to strong
backreaction effects of the produced Higgses on the inflaton background. By
the time $a\sim 4$, the Higgs field variance becomes larger than that of the
inflaton. At this stage, the system reaches quasiequilibrium such that each
degree of freedom carries the same fraction of the total energy. Since the
Higgs field has 4 d.o.f. at high energies, it dominates. We find that the
Higgs self-coupling between 0 and $10^{-2}$ at this scale does not affect the
results in any significant way. We also note that the presence of a non-
minimal Higgs coupling to gravity $H^{\dagger}HR$ would accelerate the decay
of the inflaton zero mode, yet its effect can largely be captured by
increasing the Higgs-inflaton coupling.

Since the background decays very quickly, the simple unitarity bound based on
$\rho^{1/4}(\phi_{*})\sim{1/\xi}$ becomes stronger. For the above parameters,
it requires $\xi<1000$. Another effect of a strong inflaton-Higgs coupling is
that the spectrum shifts to the UV which increases $\kappa$. For instance,
$\xi=1000$ yields $\kappa\sim 20$ at the background decay time. The resulting
$\xi_{\rm max}$ decreases by about a factor of 2. Although this makes the
unitarity constraint stronger, the improvement does not change the estimate
(54) fundamentally, given the uncertainties involved.

A significant Higgs-inflaton coupling is interesting in the context of
inflaton dark matter [40] and inflaton thermalization [41]. As shown recently
in [42], the minimal inflaton dark matter model is only viable if
$m_{\phi}\simeq m_{h}/2$. This conclusion was based on the assumption that the
non-minimal coupling $\xi$ cannot exceed about 300 without violating
unitarity. Our current work supports this assumption and reinforces the
conclusions of [42]. We note also that the direct dark matter detection limits
have become significantly stronger [43], which puts many variants of Higgs
portal dark matter under extra pressure and confines inflaton dark matter to a
narrow resonance region.

Particle production becomes even more efficient if the inflaton has derivative
couplings, as in original Higgs inflation [8]. This case was studied in [13]
with the conclusion that, for $\xi\gg 1$, much of the inflaton energy can be
converted into particles already within a single inflaton oscillation.
Clearly, this is problematic for unitarity since Higgs inflation requires
$\xi\sim 5\times 10^{4}$.

Finally, let us note that this framework offers interesting reheating, dark
matter and leptogenesis phenomenology [44, 45, 46]. It does not require an
explicit coupling of dark matter to the inflaton. Since the inflationary scale
is high, dark matter is abundantly produced by gravitational effects during
inflation and/or preheating [47].

##  4 Conclusion

Models of inflation driven by a non-minimal scalar coupling to curvature are
among the most attractive and experimentally viable models of the Early
Universe. We have studied postinflationary dynamics in a singlet scalar model
of this type. In particular, employing lattice simulations, we have focused on
perturbative unitarity constraints after inflation with a large non-minimal
coupling $\xi$. The main ingredients in our study are the decay time of the
inflaton zero mode and collective effects associated with large occupation
numbers of the inflaton momentum modes. We find that for values of $\xi$ above
a few hundreds, such collective effects lead to large scattering amplitudes
within effective field theory, thereby violating perturbative unitarity. This
is the case even if the inflaton does not have any substantial couplings to
other fields. We thus obtain an upper bound on $\xi$, at least, in the
effective field theory description. If the inflaton has a significant coupling
to other fields, particle production becomes more efficient and the unitarity
bound becomes stronger.

Acknowledgements. We are grateful to the authors of [31] for providing us with
a customized version of CosmoLattice. This work was performed using HPC
resources from the ‘‘Mésocentre’’ computing center of CentraleSupélec, École
Normale Supérieure Paris-Saclay and Université Paris-Saclay supported by CNRS
and Région Île-de-France (https://mesocentre.universite-paris-scalay.fr/).
This project has received support from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 860881-HIDDeN, and the IN2P3 Master Project UCMN.

## References

  * [1] A. A. Starobinsky, Phys. Lett. B 91 (1980) 99-102.
  * [2] A. H. Guth, Phys. Rev. D 23 (1981) 347-356.
  * [3] A. D. Linde, Phys. Lett. B 108 (1982) 389-393; Phys. Lett. B 129 (1983), 177-181.
  * [4] V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33, 532-535 (1981).
  * [5] N. A. Chernikov and E. A. Tagirov, Ann. Inst. H. Poincare Phys. Theor. A 9, 109 (1968).
  * [6] R. Fakir and W. G. Unruh, Phys. Rev. D 41, 1783-1791 (1990).
  * [7] I. L. Buchbinder, S. D. Odintsov and I. L. Shapiro, ‘‘Effective action in quantum gravity,’’ Bristol, UK: IOP (1992).
  * [8] F. L. Bezrukov and M. Shaposhnikov, Phys. Lett. B 659, 703-706 (2008).
  * [9] P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016).
  * [10] C. P. Burgess, H. M. Lee and M. Trott, JHEP 09 (2009), 103.
  * [11] J. L. F. Barbon and J. R. Espinosa, Phys. Rev. D 79 (2009), 081302.
  * [12] F. Bezrukov, A. Magnin, M. Shaposhnikov and S. Sibiryakov, JHEP 01, 016 (2011).
  * [13] Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, JCAP 02, 045 (2017).
  * [14] Y. Ema, R. Jinno, K. Nakayama and J. van de Vis, Phys. Rev. D 103, no.10, 103536 (2021).
  * [15] F. Bauer and D. A. Demir, Phys. Lett. B 698, 425-429 (2011).
  * [16] G. F. Giudice and H. M. Lee, Phys. Lett. B 694, 294-300 (2011).
  * [17] O. Lebedev and H. M. Lee, Eur. Phys. J. C 71, 1821 (2011).
  * [18] J. L. F. Barbon, J. A. Casas, J. Elias-Miro and J. R. Espinosa, JHEP 09, 027 (2015).
  * [19] M. He, A. A. Starobinsky and J. Yokoyama, JCAP 05, 064 (2018).
  * [20] J. Fumagalli and M. Postma, JHEP 05, 049 (2016).
  * [21] Y. Hamada, K. Kawana and A. Scherlis, JCAP 03, 062 (2021).
  * [22] A. Ito, W. Khater and S. Rasanen, JHEP 06, 164 (2022).
  * [23] G. K. Karananas, M. Shaposhnikov and S. Zell, JHEP 06, 132 (2022).
  * [24] J. Repond and J. Rubio, JCAP 07, 043 (2016).
  * [25] J. van de Vis, R. Nguyen, E. I. Sfakianakis, J. T. Giblin and D. I. Kaiser, Phys. Rev. D 102, no.4, 043528 (2020).
  * [26] F. Dux, A. Florio, J. Klarić, A. Shkerin and I. Timiryasov, JCAP 09, 015 (2022).
  * [27] C. Joana, Phys. Rev. D 106, no.2, 023504 (2022).
  * [28] O. Lebedev, T. Solomko and J. H. Yoon, JCAP 02, 035 (2023).
  * [29] D. S. Salopek, J. R. Bond and J. M. Bardeen, Phys. Rev. D 40, 1753 (1989).
  * [30] J. Garcia-Bellido, D. G. Figueroa and J. Rubio, Phys. Rev. D 79, 063531 (2009).
  * [31] D. G. Figueroa, A. Florio, T. Opferkuch and B. A. Stefanek, [arXiv:2112.08388 [astro-ph.CO]].
  * [32] S. Tsujikawa, K. i. Maeda and T. Torii, Phys. Rev. D 61, 103501 (2000).
  * [33] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, JCAP 04, 035 (2021).
  * [34] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, [arXiv:2102.01031 [astro-ph.CO]].
  * [35] P. B. Greene, L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 6175-6192 (1997).
  * [36] G. Dvali, G. F. Giudice, C. Gomez and A. Kehagias, JHEP 08, 108 (2011).
  * [37] G. Dvali and L. Eisemann, Phys. Rev. D 106, no.12, 125019 (2022).
  * [38] O. Lebedev, Prog. Part. Nucl. Phys. 120, 103881 (2021).
  * [39] C. Gross, O. Lebedev and M. Zatta, Phys. Lett. B 753, 178-181 (2016).
  * [40] R. N. Lerner and J. McDonald, Phys. Rev. D 80, 123507 (2009).
  * [41] O. Lebedev, T. Nerdi, T. Solomko and J. H. Yoon, Phys. Rev. D 106, no.4, 043537 (2022).
  * [42] O. Lebedev and J. H. Yoon, Phys. Lett. B 821, 136614 (2021).
  * [43] J. Aalbers et al. [LZ], [arXiv:2207.03764 [hep-ex]].
  * [44] S. Clery, Y. Mambrini, K. A. Olive, A. Shkerin and S. Verner, Phys. Rev. D 105, no.9, 095042 (2022).
  * [45] R. T. Co, Y. Mambrini and K. A. Olive, Phys. Rev. D 106, no.7, 075006 (2022).
  * [46] B. Barman, S. Cléry, R. T. Co, Y. Mambrini and K. A. Olive, JHEP 12, 072 (2022).
  * [47] O. Lebedev, JCAP 02, 032 (2023).