# Kinetic Simulations of Strongly-Magnetized Parallel Shocks: Deviations from
MHD Jump Conditions

Colby C. Haggerty,1,2 Antoine Bret,3,4 Damiano Caprioli2,5
1Institute for Astronomy, University of Hawai‘i, Honolulu, HI 96822, USA
2Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL
60637, USA
3ETSI Industriales, Universidad de Castilla-La Mancha, 13071 Ciudad Real,
Spain
4Instituto de Investigaciones Energéticas y Aplicaciones Industriales, Campus
Universitario de Ciudad Real, 13071 Ciudad Real, Spain
5Enrico Fermi Institute, The University of Chicago, 5640 S Ellis Ave, Chicago,
IL 60637, USA E-mail<EMAIL_ADDRESS>(CCH)

###### Abstract

Shocks waves are a ubiquitous feature of many astrophysical plasma systems,
and an important process for energy dissipation and transfer. The physics of
these shock waves are frequently treated/modeled as a collisional, fluid MHD
discontinuity, despite the fact that many shocks occur in the collisionless
regime. In light of this, using fully kinetic, 3D simulations of non-
relativistic, parallel propagating collisionless shocks comprised of electron-
positron plasma, we detail the deviation of collisionless shocks form MHD
predictions for varying magnetization/Alfvénic Mach numbers, with particular
focus on systems with Alfénic Mach numbers much smaller than sonic Mach
numbers. We show that the shock compression ratio decreases for sufficiently
large upstream magnetic fields, in agreement with the predictions of Bret &
Narayan (2018). Additionally, we examine the role of magnetic field strength
on the shock front width. This work reinforces a growing body of work that
suggest that modeling many astrophysical systems with only a fluid plasma
description omits potentially important physics.

###### keywords:

editorials, notices – miscellaneous

††pubyear: 2021††pagerange: Kinetic Simulations of Strongly-Magnetized
Parallel Shocks: Deviations from MHD Jump Conditions–Kinetic Simulations of
Strongly-Magnetized Parallel Shocks: Deviations from MHD Jump Conditions

##  1 Introduction

In a collisional fluid, dissipation at the front of a shockwave is provided by
binary collisions. As a consequence, the shock front is a few mean free path
thick (Zel’dovich & Raizer, 2002). Yet, in plasma, shock waves can be
sustained by collective effects, with a front smaller than the mean free path
by several orders of magnitude (Petschek, 1958; Buneman, 1964; Sagdeev, 1966).

A good example is the bow shock formed at the boundary of the Earth’s
magnetosphere and the solar wind, in which the shock front thickness is on the
order of 100 kilometers, while the proton mean free path is comparable to an
AU (Bale et al., 2003; Schwartz et al., 2011).

Such shocks are omnipresent in astrophysical systems and are good potential
candidates to explain the origin of cosmic rays, gamma-ray bursts and fast
radio bursts; as such, they have been under intense scrutiny for decades
(e.g., Blandford & Eichler, 1987; Piran, 2005; Sironi et al., 2021). Because
of their collisionless characteristics the fluid formalism cannot be
straightforwardly applied to study them. Consequently, Particle-In-Cell (PIC)
simulations, based on a model that evolves the equations of collisionless,
kinetic plasmas, has become an important tool employed in the investigation of
shocks in recent years (e.g., Hoshino et al., 1992; Dieckmann et al., 2000;
Spitkovsky, 2008a; Sironi & Spitkovsky, 2011a; Caprioli & Spitkovsky, 2014;
Park et al., 2016). The fluid properties of these kinetic shock simulations,
such as the compression ratio and downstream heating, are often measured by
their corresponding MHD predictions, even though the relevance of MHD for
collisionless shocks is not straightforward (Bret, 2020; Haggerty & Caprioli,
2020).

In this respect, parallel shocks are particular interesting when comparing the
kinetic and MHD picture. In the MHD picture, the strength of a parallel
magnetic field does not affect the hydrodynamics of shock; in this
configuration, the field and the fluid are completely disconnected
(Lichnerowicz (1976), or Kulsrud (2005), p. 141). Therefore, any detectable
effect of a parallel field on a collisionless shock represents a departure
from MHD and must be due to kinetic effects.

In a series of recent theoretical works collisionless (electron-positron)
plasma shock hydrodynamics were reconsidered from a kinetic framework (Bret et
al., 2017; Bret & Narayan, 2018, 2020), in which the density jump of a
collisionless shocks was predicted deviate from the standard MHD, Rankine
Hugoniot jump conditions; namely that an increasing the magnetic field
strength in a parallel shock should decrease the compression ratio. In this
paper, we constrain this prediction with state-of-the-art PIC simulations, and
we examine the effects of increased field strength on the thermal and
nonthermal properties of the shock. The paper is outlined as follows: In Sec.
2 we review the theory and predictions for the departure from MHD, in Sec. 3
we discuss the PIC model and the setup for the shock simulations, in Sec. 4
the fluid properties of the are examined, namely the density, the anisotropic
temperature and the entropy, and are shown to be in good agreement with
theory, in Sec. 5 we detail how the magnetic field strength modified the shock
front’s width , and finally we conclude in Sec. 6.

Figure 1: A diagram of the shock system considered, and setup in the
simulations, with subscripts 1 and 2 denoting upstream and downstream
quantities respectively.

##  2 Theory

For completeness we briefly summarize the model developed for the density jump
in Bret & Narayan (2018), hereafter BN18. We considered a non-relativistic
shock in an electron-positron (pair) plasma. The upstream is assumed
isotropic. The external magnetic field $\mn@boldsymbol{B}_{0}$ is parallel to
the shock normal, as shown in Fig. 1. The upstream density, pressure and
temperature are $n_{1},P_{1}$ and $T_{1}$ respectively. Downstream quantities
are labeled with the subscript “2”, and the downstream are assumed to be
anisotropic with $T_{2x}\neq T_{2y}$. We define the downstream anisotropy
parameter $A_{2}$ as,

$A_{2}=\frac{T_{2y}}{T_{2x}}.$ (1)

The field strength is measured by the $\sigma$ parameter (magnetization),

$\sigma\equiv\frac{B_{0}^{2}/8\pi}{\frac{1}{2}n_{1}mv_{1}^{2}}=\frac{v_{A0}^{2}}{v_{1}^{2}},$
(2)

where $v_{A0}$ is the Alfvén speed and $v_{1}$ the upstream velocity in the
shock frame.

While the standard MHD jump conditions assume isotropy, several authors have
adapted them to account for a downstream pressure anisotropy $A_{2}$
(Karimabadi et al., 1995; Vogl et al., 2001; Erkaev et al., 2000; Gerbig &
Schlickeiser, 2011). Yet, in these studies, $A_{2}$ is a free parameter
leaving the problem under-determined and an additional assumption is required
to pinpoint $A_{2}$ in terms of the upstream quantities (Gary, 1993; Bale et
al., 2009; Maruca et al., 2011).

The additional assumption made in BN18 is that the transit through the shock
front is adiabatic in the perpendicular direction while the entropy increase
goes into the parallel temperature. Regarding the perpendicular temperatures,
the relations derived in (Chew et al., 1956) imply their conservation.
Together with the macroscopic conservation equations, this assumption allows
to fully determine the characteristics of the downstream, including its
temperature anisotropy $A_{2}$.

Depending on the value of $A_{2}$ thus obtained, the downstream can be
firehose unstable if the temperature anisotropy is large relative to the
magnetic field strength. For weak fields, i.e., where $\sigma$ smaller than
some critical value $\sigma_{c}$, the downstream is firehose unstable, while
for strong enough fields, with $\sigma>\sigma_{c}$, the downstream remains
stable and constitutes the end state of the shock. The resulting density jump
($r$) can be derived,

$r\equiv\frac{n_{2}}{n_{1}}=\left\\{\begin{array}[]{r}\frac{5+\chi_{1}^{2}\left(5-\sigma+\sqrt{\Delta}\right)}{2\left(\chi_{1}^{2}+5\right)},\sigma<\sigma_{c},\\\
\frac{2\chi_{1}^{2}}{\chi_{1}^{2}+3},\sigma>\sigma_{c},\end{array}\right.$ (3)

where,

$\displaystyle\chi_{1}^{2}$ $\displaystyle=$
$\displaystyle\frac{v_{1}^{2}}{P_{1}/n_{1}},$ $\displaystyle\sigma_{c}$
$\displaystyle=$ $\displaystyle
1-\frac{4}{\chi_{1}^{2}+3}-\frac{1}{\chi_{1}^{2}},$ (4) $\displaystyle\Delta$
$\displaystyle=$
$\displaystyle\frac{25}{\chi_{1}^{4}}-\frac{10(\sigma+3)}{\chi_{1}^{2}}+(\sigma-9)(\sigma-1).$

The model explained in BN18 made use of the parameter $\chi_{1}$, which is a
pseudo Mach number. The reason why the model is not adapted to the use of a
proper Mach number can be understood by examining Eq. (3), which shows that
for $\sigma=0$ and $\chi_{1}\gg 1$, the density jump reads $r=4$, while for
$\sigma>\sigma_{c}$ and $\chi_{1}\gg 1$, we find $r=2$. Therefore, the model
presents a transition from a 3D strong shock with $\gamma=5/3$ at $\sigma=0$,
to a 1D strong shock with $\gamma=3$ at $\sigma>\sigma_{c}$ (Bret, 2021). As a
result, within the model of BN18, it is impossible to define a Mach number
using a fixed adiabatic index since the effective adiabatic index of the model
varies between 5/3 and 3.

In the rest of this work, the parameter $\chi_{1}$ is related to the real Mach
number $\mathcal{M}$ through,

$\chi_{1}^{2}=\gamma\frac{v_{1}^{2}}{\gamma
P_{1}/n_{1}}~{}~{}\Rightarrow~{}~{}\chi_{1}=\sqrt{\gamma}\mathcal{M},$ (5)

where $\gamma$ is the adiabatic index of the plasma. Accordingly, the
downstream anisotropy $A_{2}$ reads,

$A_{2}=\left\\{\begin{array}[]{r}\frac{1}{2}\frac{\chi_{1}^{2}(r-2)(r-1)+r(5r-3)}{\chi_{1}^{2}(r-1)+r},\sigma<\sigma_{c},\\\
\frac{4\chi_{1}^{2}}{\chi_{1}^{4}+2\chi_{1}^{2}-3},\sigma>\sigma_{c}.\end{array}\right.$
(6)

Notably, in the strong shock limit $\chi_{1}=\infty$, Eqs. (3,6) reduce to,

$r\equiv
r_{\infty}=\left\\{\begin{array}[]{r}\frac{1}{2}\left(\sqrt{(\sigma-9)(\sigma-1)}+5-\sigma\right),\sigma<1,\\\
2,\sigma>1,\end{array}\right.$ (7)

and

$A_{2}\equiv
A_{2\infty}=\left\\{\begin{array}[]{r}\frac{1}{4}\left(\sqrt{(\sigma-9)(\sigma-1)}+1-\sigma\right),\sigma<1\\\
0,\sigma>1.\end{array}\right.$ (8)

with, as expected, $r_{\infty}(\sigma=0)=4$ and $A_{2\infty}(\sigma=0)=1$.

##  3 Simulations

In order to test these predictions we perform a survey of non-relativistic 3D
particle-in-cell (PIC) simulations of electron-positron (pair-plasma) shocks
with TRISTAN-MP (Buneman, 1993; Spitkovsky, 2005). The choice to model a pair-
plasma was made to more directly compare with the theory developed in BN18
which was derived for an electron-positron plasma. This choice was motivated
by two reasons: the temperature of the two species would have the same value
downstream of the shock and the instabilities, relevant for this problem, are
unchanged in the electron-positron limit (Gary & Karimabadi, 2009;
Schlickeiser, 2010). In the ion-electron limit, the temperature of the two
species is not necessarily equal (Guo et al., 2017, 2018; Tran & Sironi,
2020), however, examining this problem in this framework would require a
description of the ion/electron energy partition. While this is an important
step in understanding the physics of non-relativistic collisionless shocks, it
is beyond the scope of this work

The setup of the shock simulations are similar to those described in
Spitkovsky (2008b); Sironi & Spitkovsky (2011b); Sironi et al. (2013), where
the shock is formed by initializing the plasma with a supersonic flow towards
a reflecting wall (in the negative $x$ direction) on the left hand side of the
simulation. The plasma is reflected and forms a shock that is mediated by some
counter streaming instabilities. The plasma is continuously injected by an
injector receding away from the shock which extends the simulation domain. The
boundary conditions along $y$ and $z$ are periodic.

Values in the simulations are measured in normalized code units: speeds are
normalized to the speed of light, times to the inverse plasma frequency
$\omega_{pe}^{-1}=\sqrt{m/4\pi ne^{2}}$, where $e$ is the fundamental charge
and $n$ is the number density of the inflowing plasma, and lengths,
consistently, to the electron inertial length $d_{e}=c/\omega_{pe}$. The
maximum allowed simulation length (i.e., where the receding injector stops) is
$L_{x}=1200c/\omega_{pe}$, and the transverse width and depth are
$L_{y}=L_{z}=16c/\omega_{pe}$ for all simulations. The simulations are
performed with $2$ electron/positron pairs per cell, $4$ grid cells per
$c/\omega_{pe}$, and with the time step chosen such that $\Delta t=.45\Delta
x/c$. A convergence simulation was run for the $M_{A}=M_{s}$ case, with double
the transverse box size, twice the spatial resolution and $\times 4$ particles
per cell, and no significant difference in the results found.

The inflowing electrons and positrons are initialized as a drifting Maxwellian
with equal temperature, $\Delta\gamma\equiv k_{B}T/mc^{2}=10^{-4}$. Because of
the reflecting wall, the downstream plasma is at rest in the simulation frame.
In this frame, the upstream plasma is in-flowing at $v_{sh}/c=0.1$, which
corresponds to a Mach number of $\mathcal{M}=v_{sh}/\sqrt{\gamma P/n}\approx
7.75$. However, the shock predictions and theory outlined above are determined
in the rest frame of the shock. The two frames are ultimately related by the
density jump $r$, i.e., $v_{1}=v_{sh}/(1-1/r)$, as discussed in many previous
simulations of shocks formed by reflecting walls (Caprioli & Spitkovsky, 2014;
Haggerty & Caprioli, 2019, 2020).

The simulations are initialized with a uniform magnetic field anti-parallel to
the upstream flow (i.e., $+x$). The initial magnetic field strength varies
between simulations and is characterized by the shock energy density in the
downstream frame $\sigma^{\prime}=B_{1}^{2}/4\pi n_{1}mv_{sh}^{2}$ and is
related to the value described by Eq. 2 by the density jump,
$\sigma=\sigma^{\prime}(1-1/r)^{2}$. The survey consists of 9 simulations
performed over a range of $\sigma^{\prime}$ values from $10^{-3}$ to $25$,
allowing for the consideration of both super-Alfvénic and sub-Alfvénic shocks.
The simulations are evolved for 1000’s of $\omega_{pe}^{-1}$ but are stopped
before particles can escape the simulation.

##  4 Hydrodynamics

Various hydrodynamic shock quantities can be determined from the survey of
simulations, and in the following section they are presented and compared with
the theory of BN18.

###  4.1 Density Jump

Fig. 2 shows the shock density profile integrated in the transverse
directions, for $\sigma$ up to 25. Data are plotted at time
$6000\omega_{pe}^{-1}$ to insure that the shock is fully formed.

Figure 2: Panel (a): Shock density profile integrated in the transverse
directions, for $\chi_{1}=10$ and various $\sigma$, for this and all following
figures, the data points are color-coded by value of $\sigma$ and correspond
to the values in the color bar. Panel (b): The average downstream density
normalized to the upstream density for various $\sigma$. Data are plotted at
time $6000\omega_{pe}^{-1}$ and the dashed line pertains to Eq. (7).

The jump obtained for values of $\sigma\ll 1$ is 4 consistent with the MHD
prediction for strong shocks. Although the critical magnetization is close to
unity ($\sigma_{c}\sim 1$), a clear departure from the standard MHD prediction
is observed for $\sigma$ small as $0.1$. In this regime, where
$0.1\lesssim\sigma\lesssim 1$, the Mach number based on the Alfvén speed would
range between $1$ and $10$, and so the decreasing compression ratios seems
consistent with the decreasing Mach number.

The lower plot in figure 2 displays the density jump for a range of $\sigma$
values, contrasted against Eq. (7) represented by the black dashed line.
Simulations show the downstream density enhancement decreasing from 4 to 2,
consistent with prediction presented above. The theoretical predictions
attribute the departure from the MHD prediction to the stabilization of
downstream temperature anisotropy, which can be directly examined form the
summations and is discussed in Sec. 4.4

Figure 3: Panel (a): The change in the kinetic entropy, as defined by Eq. 9,
across the shock in the simulations with different values of $\sigma$ (color
coded). Panel (b): The downstream parallel distribution function ($\int
fdv_{y}dv_{z}$) from simulations with different values of $\sigma$. The
upstream distribution shown by the black dashed line at time
$6000\omega_{pe}^{-1}$. Simulations with smaller values of $\sigma$ show
larger entropy increases with more normal downstream distributions, while
simulations with larger $\sigma$ values have a smaller entropy increase, with
distributions that are double peeked.

###  4.2 Entropy

For simulations with relatively small values of $\sigma\ll\sigma_{c}$, it is
clear that a shock in the typical sense is forming, i.e., the upstream and
downstream regions are separated by a sharp discontinuity, and the downstream
plasma is compressed consistent with the MHD jump conditions; however, for
simulations with $\sigma\gtrsim\sigma_{c}$ it is less clear if a bona fide
shock has developed, or we are witnessing interpenetrating fluids. Ultimately,
shocks must increase entropy, and as such the issue can be addressed by the
entropy difference between the upstream and the downstream. While there remain
unsettled issues in the interpretation of entropy and reversibility in a
collisionless system, or even effective collisions introduced by poor counting
statistics of the PIC model (e.g. Montgomery & Nielson, 1970; Shalaby et al.,
2017; Kadanoff, 2017; Yang et al., 2017; Liang et al., 2019; Du et al., 2020,
and references therein), we limit considerations to the definition of kinetic
entropy density per particle:

$s=\frac{S}{n}=-\frac{\int f\ln fd^{3}v}{\int fd^{3}v}$ (9)

where $S$ is the kinetic entropy density, $f$ is the distribution function.
This calculation was explicitly performed in Sec. 4 of BN18 and considered
again in Bret (2021) where the proof is made that for large values of
$\chi_{1}$, the entropy jump for both the large and small $\sigma$ regimes is
definitely positive. Since the marginal firehose stability is observed in the
present simulations at small $\sigma$ (see Fig. 4), while the perpendicular
temperature conservation is also observed far larger $\sigma$’s (see Sec 4.4),
the calculations of BN18 and Bret (2021) should apply, ensuring an entropy
increase at the transition. The measured entropy jump ($\Delta s$) for all of
the simulations as a function of $\sigma$ is shown in panel (a) of Fig. 3,
with the black dashed line showing the theoretical prediction. The value is
calculated by taking the difference of integrated discretized distribution
functions averaged over a region of $100\ d_{e}$ in the downstream and
upstream. There is good quantitative agreement for simulations with
$\sigma<\sigma_{c}$, and as $\sigma$ is increased and the downstream remains
firehose stable across the shock, there is qualitative agreement with theory,
i.e., the difference in entropy reduces. It is worth noting that the entropy
increase for the $\sigma\approx 6$ simulation agrees well with theory,
however, in the simulation with an even larger magnetization, $\sigma\approx
14$, the entropy goes to nearly zero ($\Delta s\approx 0$). This difference
can be understood by considering the parallel downstream distribution
functions for these simulations, shown in panel (b) of Fig. 3. In the entropy
derivation of BN18, the distribution was assumed to have the form of a bi-
maxwellian, while it is clear from the dark red line in the lower panel, that
the distribution has two clear peaks associated with counter-streaming beams.
However, it is likely that —if the simulation were run for longer times in a
larger box— the double-peaked distribution function would tend towards the bi-
maxwellian state, as the distribution function has already gaussianized
relative to the initial distribution. From this it seems reasonable that the
entropy will increase across the discontinuity, thus constituting a
collisionless shock even in large-$\sigma$ systems.

Figure 4: Profiles of anisotropic parameters for strong shocks
($\chi_{1}=10$), taken at $6000\omega_{pe}^{-1}$, with varying values of
$\sigma$ (color-coded). Panel (a): Ratios of the perpendicular to parallel
electron temperature ($A=T_{yy}/T_{xx}$), averaged over the transverse
directions and shown relative to the shock front. Panel (b): Profiles of both
sides of Eq. 10 (left side, $A$, green line; right side, firehose, maroon)
superimposed on the density normalized to the downstream value (black line)
for $\sigma^{\prime}=0.1$ ($\sigma\approx 0.056$). The system is firehose
unstable in regions where $1-1/\beta_{\parallel}>A_{2}$, corresponding to the
yellow shaded region in the upstream, in contrast to the marginal stability in
the downstream. Panel (c): The averaged downstream temperature anisotropy
(right axis, red line with squares) and firehose parameter (left axis, black
line with triangles) shown as a function of $\sigma$. The red dashed line
correspond to prediction for $A_{2}$ in the strong shock limit, as defined in
Eq. 6.

###  4.3 Firehose stability

The density jump measured in the simulations is in good agreement with the
theoretical predictions from BN18 and to further verify the consistency
between simulations an BN18, we also test the assumptions that form the basis
of the model. The main ingredient of the BN18 model is that for weak magnetic
field strengths, the plasma turns firehose unstable at the shock front
crossing. As a result, it migrates to the stability threshold so that the
downstream is marginally firehose stable. For strong fields, we expect the
downstream to never becomes unstable, the condition for which is (Gary &
Karimabadi, 2009; Schlickeiser, 2010),

$A_{2}>1-\frac{1}{\beta_{\parallel 2}},$ (10)

with,

$\beta_{\parallel 2}=\frac{n_{2}k_{B}T_{\parallel
2}}{B_{0}^{2}/8\pi}=\frac{n_{2}k_{B}T_{x2}}{B_{0}^{2}/8\pi}=\frac{4}{\chi_{2}^{2}\sigma}.$
(11)

Fig. 4(a) shows the anisotropy parameter $A(x)=T_{y}/T_{x}$ with the x-axis
re-centered on the shock front for various values of $\sigma$ (color coded).
The anisotropy far upstream for all simulations, namely $x\gtrsim
700c/\omega_{pe}$, is 1 but it decreases in the shock foot111In order to
correctly render the shock front and the downstream, we do not show this far
upstream where $A$ reaches unity. In contrast with what was assumed in BN18,
the structure of the shock is not a discontinuity, which will be discussed
further in Sec. 5. The anisotropy parameter downstream behaves as expected; at
low $\sigma$, the field cannot sustain a stable anisotropy and the downstream
is nearly isotropic and $A_{2}\sim 1$. As the field strength increases and
$\sigma>\sigma_{c}$, $A_{2}$ progressively goes to 0. This is further
demonstrated in Fig. 4(c), in which the red line (corresponding to the right
y-axis) shows the average downstream value of $A_{2}$ as a function of
$\sigma$. In Fig. 4(c) the red dashed line shows the prediction for $A_{2}$ as
a function of $\sigma$ (Eq. 8), and the two are found to be in good agreement.

Another ingredient of BN18 was that for weak fields, the plasma turns firehose
unstable when crossing the shock front before it migrates to marginal
stability. In this respect, Fig. 4(b) shows both sides of Eq. 10 superimposed
on the shock density profile for $\sigma=0.056$ ($\sigma^{\prime}=0.1$). While
the downstream shows marginal stability with $A_{2}\sim 1-1/\beta_{\parallel
2}$, we find $A_{2}<1-1/\beta_{\parallel 2}$ for
$100c/\omega_{pe}<x<250c/\omega_{pe}$. As expected, the plasma turns firehose
unstable as it passes to the downstream, before it settles to marginal
stability.

Finally, we can determine the downstream firehose condition for the different
$\sigma$ simulations; Fig. 4(c) shows the averaged downstream values from the
two sides of Eq. 10 (with $A_{2}$ denoted by the red line with squares
measured by the right axis and the black line with triangles measured by the
left axis). For small values of $\sigma$ Both terms are close to 1. As
$\sigma$ increases to $\sigma\gtrsim 0.1$, the LHS of Eq. 10 drops and even
becomes negative implying that the downstream is firehose stable, while
simultaneously the temperature ratio ($A_{2})$, also decreases from 1 to 0
showing that a stable field aligned temperature anistropy has developed.

###  4.4 Structure of the Downstream Temperature Anisotropy

Figure 5: Panel (a) and (b): Profiles of the parallel ($T_{xx}$, (a)) and
perpendicular ($T_{zz},(b)$) temperatures, averaged along the transverse
directions, showing the development of a strong temperature anisotropy for
increasing values of $\sigma$, color-coded. Panel (c): Components of the
average downstream temperatures, in units normalized to
$mv^{2}_{1}=mv_{A}^{2}/\sigma$, as a function of $\sigma$, with the two black
dashed lines corresponding to $1$ and $1/3$ respectively, with $T_{xx}$ as
pink, $T_{yy}$ as dark green, $T_{zz}$ as purple and the trace of the
temperature tensor, $\mathrm{Tr}[T]=\sum_{i}T_{i}$ as grey.

The predictions laid out in BN18 ultimately derived from the temperature
anisotropy generated downstream of the shock. The anisotropy is predicted to
grow for strong magnetic fields, as the inflowing plasma is expected to remain
magnetized, such that the perpendicular temperature should not change
appreciably across the shock. This assumption is examined in Fig. 5, which
shows the parallel ($T_{xx}$, panel a) and out-of-plane ($T_{zz}$, panel b)
temperature profiles for the different shock simulations. Note that the
temperatures in this figure are normalized to twice the kinetic energy of the
inflowing plasma, rather than to the rest mass energy, and the x-axis is
centered around the shock. For the simulations with smaller values of $\sigma$
($\lesssim 0.3$, blue to green) both the parallel and perpendicular
temperatures increase to roughly equal values downstream, resulting in an
isotropic downstream temperature, consistent with Fig. 4(a). For increasing
magnetization ($\sigma\gtrsim 0.3$, yellow to red), The parallel temperature
increases dramatically, while the perpendicular temperature become nearly
constant across the shock, verifying the assumption of BN18.

Further details about the downstream anisotropy and be found from Fig. 5(c),
which analyses how the thermal energy is split among the degrees of freedom.
In the plot, the average downstream temperatures are shown as a function of
$\sigma$, with the parallel temperature represented by the pink line with
squares, perpendicular with a green line with triangles (y)/purple line with
circles (z) and the average with a gray line with ‘x’s (denoted by
$\mathrm{Tr}[T]=\sum T_{i}$.) 222The choice of normalization makes it clear
that the total downstream temperature is equal to the second moment of two
counter streaming beam populations moving at $\pm v_{sh}$. This fact is likely
a consequence of the simulation design, as the downstream population must be
at rest in simulation frame, and so all of the energy in this frame must be
thermal. It seems natural that the downstream temperature must be equal to the
composite energy of the counter streaming cold beams. For $\sigma<0.1$, the
different components of the temperature are equal. As the only stable
configuration has $A_{2}\sim 1$, the downstream is isotropic with its thermal
energy evenly distributed along the 3 degrees of freedom. As a consequence,
each one carries 1/3 of the total, and all together they sum up to
$mv^{2}_{sh}$ which is the total amount of energy initially available in the 2
cold counter streaming plasmas. For $\sigma>0.1$, the downstream becomes
increasingly anisotropic with increasing $\sigma$ according to the mechanism
described in BN18. The parallel temperature $T_{x}$ rises while $T_{y}$ and
$T_{z}$ decrease in such a way that $\mathrm{Tr}[T]$ is kept constant and
equal to $mv^{2}_{sh}$ (momentum conservation).

##  5 Shock Front Width

From the density profiles shown in Fig. 2(a), it is clear that the width of
the shock has some dependence on $\sigma$. To further examine this effect, we
determined the width of the shock fronts for each simulation at
$6000\,\omega_{pe}^{-1}$. This was achieved by fitting the shifted, average
density profiles with a ${\rm tanh}$ function of the form,

$\frac{n(x)}{n_{1}}=\frac{1+r}{2}+\frac{r-1}{2}\tanh(-x/\lambda).$ (12)

The value of $\lambda$ for each simulation is determined with a non-linear
least squares fit and the results are shown in Fig. 6.

Although the width of the shock varies non-monotonically with $\sigma$, two
distinct patterns appear to emerge when considering simulations with
$\sigma\ll 0.1$ or $\sigma\gg 0.1$. In the small $\sigma$ case, the shock
width seems to reach a constant size for decreasing $\sigma$. This is
consistent with the shock formation being mediated by the two stream
instability for the present parameters (Bret et al., 2008), which has a growth
rate independent of $\sigma$.

Figure 6: The width of the shock as a function of $\sigma$, determined from
the best fit of a ${\rm tanh}$ function shown in Eq. (12), color-coded by
$\sigma$.

For $\sigma\gg 0.1$, the shock width increases with $\sigma$. This increase is
potentially due to the increasing number of pitch angle scattering required to
return particles traveling up-streaming back towards the downstream. Each
deflection is expected to alter a particle’s pitch angle by $\Delta
mu\sim\delta B/B_{0}\sim 1/\sqrt{\sigma}$. For a larger value of $\sigma$,
more deflections will be required to send particles back potentially leading
to a larger shock width. This process is shown in an animation of self-
consistent particle tracing for the $\sigma^{\prime}=2$ simulation (orange
line in this papers figures) in the supplementary material.

##  6 Conclusion

In this work, using self-consistent kinetic PIC simulations, we verify the
departure of parallel collisionless shocks from fluid MHD predictions; by
varying the magnetization on collisionless pair plasma shocks we show that a
parallel propagating shock wave is suppressed for sufficiently strong magnetic
field strengths, consistent with kinetic theoretical predictions of BN18. For
a fixed (parallel) shock inclination and sonic Mach number
($v_{1}/v_{th}=10$), we test the predictions of BN18 of how the shock
properties depend on the strength of the upstream magnetic field relative to
the upstream shock speed, a quantity parametrized by $\sigma\equiv
B_{0}^{2}/4\pi n_{1}mv_{1}^{2}$. The simulations are found to be in good
agreement with the theoretical predictions, namely that as $\sigma$ is
increased to $\gtrsim 1$, the shocked plasma remains well magnetized,
suppressing the firehose instability and leading to a large field aligned
temperature anisotropy downstream. The entropy jump between the downstream and
upstream is considered and found to be approximate agreement with predictions.

Additionally we examine how the kinetic properties of the shock are modified
by the magnetization. The effective shock width is determined for a fixed time
for each of the simulations, and the width is found to be a minimum for an
upstream plasma beta of unity ($\beta\sim 1$). For decreasing and increasing
values of $\sigma$ the shock width broadens, but saturating for small values
of $\sigma$ as the instability mediating the shock reformation becomes two
stream like and insensitive to the magnetic field strength.

In the theory and simulations of this work, only a pair plasma was considered;
this choice simplified the problem by ensuring that the downstream
temperatures would be the same for both species. In shocks comprised of
electron-proton plasmas, it is likely that the temperature of each species
will be different (Tran & Sironi, 2020). This added complexity should be
further explored to insure the applicability of this work to space and
astrophysical shocks and will be considered in future investigations.

These results contribute to the building consensus collisionless shocks
fundamentally require a kinetic description to accurately capture both the
thermal and non-thermal plasma physics.

## Data Availability

The datasets used in this manuscript were derived using the open source
particle-in-cell software Tristan-MP v2 (https://ntoles.github.io/tristan-
wiki/). The the average output for each simulation requires $\sim 14$ GB of
space ($\gtrsim 200$ GB for all of the simulations and convergence studies)
and as such either the data underlying this article or the simulation
initialization files will be shared on reasonable request to the corresponding
author.

## Acknowledgements

We thank the anonymous referee for thorough comments which helped improve the
manuscript. This work has been achieved under projects ENE2016-75703-R from
the Spanish Ministerio de Economía y Competitividad and SBPLY/17/180501/000264
from the Junta de Comunidades de Castilla-La Mancha. As well as the FDSS NSF
AGS-1936393 grant to the Institute for Astronomy of the University of Hawaii.
Simulations were performed on computational resources provided by the
University of Chicago Research Computing Center and XSEDE TACC (TG-AST180008).

## References

  * Bale et al. (2003) Bale S. D., Mozer F. S., Horbury T. S., 2003, Phys. Rev. Lett., 91, 265004
  * Bale et al. (2009) Bale S. D., Kasper J. C., Howes G. G., Quataert E., Salem C., Sundkvist D., 2009, Phys. Rev. Lett., 103, 211101
  * Blandford & Eichler (1987) Blandford R., Eichler D., 1987, Phys. Rep., 154, 1
  * Bret (2020) Bret A., 2020, ApJ, 900, 111
  * Bret (2021) Bret A., 2021, arXiv e-prints, p. arXiv:2107.10324
  * Bret & Narayan (2018) Bret A., Narayan R., 2018, Journal of Plasma Physics, 84, 905840604
  * Bret & Narayan (2020) Bret A., Narayan R., 2020, Laser and Particle Beams, 38, 114
  * Bret et al. (2008) Bret A., Gremillet L., Bénisti D., Lefebvre E., 2008, Phys. Rev. Lett., 100, 205008
  * Bret et al. (2017) Bret A., Pe’er A., Sironi L., Sa̧dowski A., Narayan R., 2017, Journal of Plasma Physics, 83, 715830201
  * Buneman (1964) Buneman O., 1964, Physics of Fluids, 7, S3
  * Buneman (1993) Buneman O., 1993, in H. Matsumoto Y. Omura eds, Computer Space Plasma Physics. Tokyo: Terra Scientific, p. 67
  * Caprioli & Spitkovsky (2014) Caprioli D., Spitkovsky A., 2014, Astrophys. J., 783, 91
  * Chew et al. (1956) Chew G. F., Goldberger M. L., Low F. E., 1956, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 236, 112
  * Dieckmann et al. (2000) Dieckmann M. E., Chapman S. C., McClements K. G., Dendy R. O., Drury L. O., 2000, Astronomy & Astrophysics, 356, 377
  * Du et al. (2020) Du S., Zank G. P., Li X., Guo F., 2020, Phys. Rev. E, 101, 033208
  * Erkaev et al. (2000) Erkaev N. V., Vogl D. F., Biernat H. K., 2000, Journal of Plasma Physics, 64, 561
  * Gary (1993) Gary S., 1993, Theory of Space Plasma Microinstabilities. Cambridge Atmospheric and Space Science Series, Cambridge University Press
  * Gary & Karimabadi (2009) Gary S. P., Karimabadi H., 2009, Physics of Plasmas, 16, 042104
  * Gerbig & Schlickeiser (2011) Gerbig D., Schlickeiser R., 2011, The Astrophysical Journal, 733, 32
  * Guo et al. (2017) Guo X., Sironi L., Narayan R., 2017, The Astrophysical Journal, 851, 134
  * Guo et al. (2018) Guo X., Sironi L., Narayan R., 2018, The Astrophysical Journal, 858, 95
  * Haggerty & Caprioli (2019) Haggerty C. C., Caprioli D., 2019, ApJ, 887, 165
  * Haggerty & Caprioli (2020) Haggerty C. C., Caprioli D., 2020, ApJ, 905, 1
  * Hoshino et al. (1992) Hoshino M., Arons J., Gallant Y. A., Langdon A. B., 1992, ApJ, 390, 454
  * Kadanoff (2017) Kadanoff L., 2017, J. Stat. Phys, 167, 1079
  * Karimabadi et al. (1995) Karimabadi H., Krauss-Varban D., Omidi N., 1995, Geophysical Research Letters, 22, 2689
  * Kulsrud (2005) Kulsrud R. M., 2005, Plasma Physics for Astrophysics. Princeton Univ. Press, Princeton, NJ
  * Liang et al. (2019) Liang H., et al., 2019, Physics of Plasmas, 26, 082903
  * Lichnerowicz (1976) Lichnerowicz A., 1976, Journal of Mathematical Physics, 17, 2135
  * Maruca et al. (2011) Maruca B. A., Kasper J. C., Bale S. D., 2011, Phys. Rev. Lett., 107, 201101
  * Montgomery & Nielson (1970) Montgomery D., Nielson C. W., 1970, Physics of Fluids, 13, 1405
  * Park et al. (2016) Park H.-S., et al., 2016, Journal of Physics Conference Series, 688, 012084
  * Petschek (1958) Petschek H. E., 1958, Rev. Mod. Phys., 30, 966
  * Piran (2005) Piran T., 2005, Rev. Mod. Phys., 76, 1143
  * Sagdeev (1966) Sagdeev R. Z., 1966, Reviews of Plasma Physics, 4, 23
  * Schlickeiser (2010) Schlickeiser R., 2010, The Open Plasma Physics Journal, 3, 1
  * Schwartz et al. (2011) Schwartz S. J., Henley E., Mitchell J., Krasnoselskikh V., 2011, Phys. Rev. Lett., 107, 215002
  * Shalaby et al. (2017) Shalaby M., Broderick A. E., Chang P., Pfrommer C., Lamberts A., Puchwein E., 2017, ApJ, 841, 52
  * Sironi & Spitkovsky (2011a) Sironi L., Spitkovsky A., 2011a, The Astrophysical Journal, 726, 75
  * Sironi & Spitkovsky (2011b) Sironi L., Spitkovsky A., 2011b, ApJ, 741, 39
  * Sironi et al. (2013) Sironi L., Spitkovsky A., Arons J., 2013, The Astrophysical Journal, 771, 54
  * Sironi et al. (2021) Sironi L., Plotnikov I., Nättilä J., Beloborodov A. M., 2021, Phys. Rev. Lett., 127, 035101
  * Spitkovsky (2005) Spitkovsky A., 2005, in Bulik T., Rudak B., Madejski G., eds, American Institute of Physics Conference Series Vol. 801, Astrophysical Sources of High Energy Particles and Radiation. pp 345–350 (arXiv:astro-ph/0603211), doi:10.1063/1.2141897
  * Spitkovsky (2008a) Spitkovsky A., 2008a, Astrophys. J. Lett., 673, L39
  * Spitkovsky (2008b) Spitkovsky A., 2008b, Astrophys. J. Lett., 682, L5
  * Tran & Sironi (2020) Tran A., Sironi L., 2020, ApJ, 900, L36
  * Vogl et al. (2001) Vogl D. F., Biernat H. K., Erkaev N. V., Farrugia C. J., Mühlbachler S., 2001, Nonlinear Processes in Geophysics, 8, 167
  * Yang et al. (2017) Yang Y., et al., 2017, Physics of Plasmas, 24, 072306
  * Zel’dovich & Raizer (2002) Zel’dovich I., Raizer Y., 2002, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena. Dover Books on Physics, Dover Publications


# On the quasinormal modes of the de Sitter spacetime

A. López-Ortega
Centro de Investigación en Ciencia Aplicada y Tecnología Avanzada.
Unidad Legaria. Instituto Politécnico Nacional.
Calzada Legaria # 694. Colonia Irrigación. Delegación Miguel Hidalgo.
México, D. F., México.
C. P. 11500<EMAIL_ADDRESS>
###### Abstract

Modifying a method by Horowitz and Hubeny for asymptotically anti-de Sitter
black holes, we establish the classical stability of the quasinormal modes of
the de Sitter spacetime. Furthermore using a straightforward method we
calculate the de Sitter quasinormal frequencies of the gravitational
perturbations and discuss some properties of the radial functions of these
quasinormal modes.

KEYWORDS: de Sitter spacetime; Quasinormal modes; Classical stability

##  1 Introduction

The de Sitter spacetime is an extensively studied solution of the
gravitational field equations. It is a maximally symmetric spacetime and due
to its simplicity it is the appropriate background to explore and test several
ideas and methods [1], [2].

If a spacetime with an event or cosmological horizon is perturbed, it is well
known that the perturbation oscillates with a set of decaying modes, usually
called quasinormal modes (QNM), whose complex frequencies are determined by
the physical parameters of the background. They have been calculated for
several spacetimes (see the reviews [3]-[6]) and they are useful in the study
of several physical phenomena [3]-[6].

The quasinormal frequencies (QNF) of the de Sitter spacetime have been
determined in several references [7]–[13]. We know the QNF of classical fields
propagating in lower [7]–[9] and higher [10]–[13] dimensional de Sitter
spacetimes. Moreover, the recent reviews [5], [6] examine some of these
results. For $D$-dimensional ($D\geq 4$) de Sitter spacetimes, in Ref. [10] it
is asserted that for $D$ even, the QNF of the gravitational perturbations are
not well defined, and we notice that this affirmation is reproduced in the
review [5] (see also Appendix A of [11] for a similar conclusion for the
massless Klein-Gordon field in the four-dimensional de Sitter spacetime).
Nevertheless in Ref. [12] it is asserted that for even and odd dimensional de
Sitter spacetimes the QNF of massless fields are well defined (a similar
affirmation appears in the review [6]).

The values of the de Sitter QNF have found several applications, for example,
they have been used in the calculation of the Green functions that are
necessary to discuss some aspects of static patch holography [14], in the
computation of the one-loop thermal partition function for the scalar field
[15], and in the determination of the entropy spectrum for the de Sitter
horizon [16], [17]. Thus we must state clearly the existence and properties of
the de Sitter QNM.

Here using a different method we prove the classical stability of the de
Sitter QNM. This method is motivated by a procedure used to show the stability
of the QNM in asymptotically anti-de Sitter spacetimes [18]. Moreover, to
solve the mentioned issue on the existence of the de Sitter QNF in even
dimensions and employing a straightforward procedure we determine the de
Sitter QNF of the gravitational perturbations. Finally we discuss some
properties of the radial functions for the gravitational de Sitter QNM.

We organize this paper as follows. In Sect. 2 we establish the classical
stability of the de Sitter QNM. To achieve this goal we modify a method by
Horowitz and Hubeny that is suitable for asymptotically anti-de Sitter black
holes [18]. To prove the classical stability of the QNM for the Klein-Gordon
and Dirac fields we find that it is necessary to elaborate on the method by
Horowitz and Hubeny. In Sect. 3 we calculate the QNF of the gravitational
perturbations propagating in the $D$-dimensional de Sitter spacetime. We
compute these with a simpler method than the procedure used in the previous
references and discuss some properties of the radial functions for the de
Sitter QNM. In Sect. 4 we enumerate our results. Finally, in Appendix we
analyze the absence of QNF for the Klein-Gordon field with special values of
the mass and the azimuthal number that propagates in four-dimensional de
Sitter spacetime.

##  2 Stability of the de Sitter QNM

In what follows we consider the $D$-dimensional ($D\geq 4$) spherically
symmetric de Sitter spacetime [1], [2]

$\textrm{d}s^{2}=(1-r^{2})\textrm{d}t^{2}-(1-r^{2})^{-1}\textrm{d}r^{2}-r^{2}\textrm{d}\Omega^{2}_{D-2},$
(1)

where $\textrm{d}\Omega^{2}_{D-2}$ is the round metric on the $(D-2)$-sphere.
As is well known, this background has a cosmological horizon. In the units
chosen in the metric (1) the horizon is located at $r=1$. We recall that in a
spherically symmetric background the equations of motion for the perturbations
simplify to Schrödinger type equations [3]–[6], [19], [20]

$\frac{\textrm{d}^{2}U}{\textrm{d}r_{*}^{2}}+\omega^{2}U-VU=0,$ (2)

where $r_{*}$ is the tortoise coordinate or is related to it, $V$ is an
effective potential that depends on the perturbation, and $U$ is related to
the radial function of the perturbation. For the discussion that follows, in
Table 1 we enumerate the effective potentials of several fields propagating in
de Sitter spacetime [10]–[13].

Table 1: Effective potentials in $D$-dimensional de Sitter spacetime. We denote the mass of the field by $\mu$, $q$ is a nonnegative integer, $f=1-r^{2}$, and $\lambda=\pm(q+(D-2)/2)$ [10]–[13]. Field | Effective potential
---|---
Klein-Gordon | $V_{KG}=f\left(\frac{(D-2)(D-4)}{4r^{2}}f-(D-2)+\frac{q(q+D-3)}{r^{2}}+\mu^{2}\right)$
Electromagnetic scalar | $V_{EMS}=f\left(\frac{q(q+D-3)}{r^{2}}+\frac{(D-2)(D-4)}{4r^{2}}-\frac{(D-6)(D-4)}{4}\right)$
Electromagnetic vector | $V_{EMV}=f\left(\frac{q(q+D-3)-1}{r^{2}}+\frac{(D-2)^{2}-2(D-2)+4}{4r^{2}}-\frac{(D-2)(D-4)}{4}\right)$
Gravitational tensor | $V_{GPT}=f\left(\frac{(2q+D-4)(2q+D-2)}{4r^{2}}-\frac{D(D-2)}{4}\right)$
Gravitational vector | $V_{GPV}=f\left(\frac{(2q+D-4)(2q+D-2)}{4r^{2}}-\frac{(D-4)(D-2)}{4}\right)$
Gravitational scalar | $V_{GPS}=f\left(\frac{(2q+D-4)(2q+D-2)}{4r^{2}}-\frac{(D-4)(D-6)}{4}\right)$
Dirac | $V_{\pm}=\frac{f(\lambda^{2}+\mu^{2}r^{2})^{3}}{r^{2}(\lambda^{2}+\mu^{2}r^{2}+\lambda\mu f/(2\omega))^{2}}\pm\frac{f^{3/2}(\lambda^{2}+\mu^{2}r^{2})^{5/2}}{\lambda^{2}+\mu^{2}r^{2}+\lambda\mu f/(2\omega)}$
|
$\times\left[\frac{\frac{1}{2f}\frac{\textrm{d}f}{\textrm{d}r}+\frac{3\mu^{2}r}{\lambda^{2}+\mu^{2}r^{2}}}{r(\lambda^{2}+\mu^{2}r^{2}+\lambda\mu
f/(2\omega))}-\frac{\lambda^{2}+3\mu^{2}r^{2}+\frac{\lambda\mu
r}{2\omega}\frac{\textrm{d}f}{\textrm{d}r}+\frac{\lambda\mu
f}{2\omega}}{r^{2}(\lambda^{2}+\mu^{2}r^{2}+\lambda\mu
f/(2\omega))^{2}}\right]$

Here we define the de Sitter QNM as the solutions to Eq. (2) that satisfy:

(a) They are purely outgoing near the de Sitter horizon.

(b) They go to zero at $r=0$.

Notice that the boundary condition at the origin (b) is slightly different
from that of [10], [12]. In those references it is imposed the boundary
condition “the radial function is regular at $r=0$”, in particular the radial
function can be equal to a constant different from zero. Here we are imposing
the boundary conditions (a) and (b) on the function $U$ that appears in the
Schrödinger type equation (2).

The motivation for the boundary conditions (a) and (b) is as follows. Near the
cosmological horizon $(r=1)$ the effective potentials of Table 1 behave as

$\lim_{r\to 1}V=0,$ (3)

hence near the cosmological horizon the solution of Eq. (2) is

$U=C_{1}\textrm{e}^{-i\omega r_{*}}+C_{2}\textrm{e}^{i\omega r_{*}},$ (4)

where $C_{1}$ and $C_{2}$ are constants. Thus near the horizon the solution
(4) represents ingoing and outgoing waves. For the de Sitter QNM we must take
the purely outgoing wave.

Near $r=0$, almost all the effective potentials of Table 1 behave in the form
(see below for some exceptions)

$\lim_{r\to 0}V\to+\infty,$ (5)

and therefore for the potentials that satisfy this condition it is reasonable
to impose the boundary condition $U(0)=0$.

One exception to the behavior near $r=0$ given in the formula (5) is that of
the effective potential $V_{+}$ for the Dirac field with $\lambda=1$ moving in
four-dimensional de Sitter spacetime, since in this example we get

$\lim_{r\to
0}V_{+}=\frac{\frac{3\mu^{2}}{2}-\frac{1}{2}+\frac{7\mu^{3}+3\mu}{4\omega}}{\left(1+\frac{\mu}{2\omega}\right)^{3}},$
(6)

and for $\mu=0$ we obtain

$\lim_{r\to 0}V_{+}=-\frac{1}{2}.$ (7)

For the effective potential $V_{+}$ both limits (6) and (7) are finite at
$r=0$. In Fig. 1, for the massless Dirac field with $\lambda=1$, we plot the
effective potential $V_{+}$ in four-dimensional de Sitter spacetime. We see
that for these values of $\mu$, $\lambda$, and $D$ the effective potential
$V_{+}$ is negative and that near $r=0$ and $r=1$ its behavior is consistent
with the analytical limits (3) and (7).

Figure 1: For the four-dimensional de Sitter spacetime we draw the effective
potential $V_{+}$ of the massless Dirac field with $\lambda=1$. Figure 2: For
the four-dimensional de Sitter spacetime we plot the effective potentials
$V_{KG}$ of the Klein Gordon field with $q=0$, $\mu^{2}=3$ (dashed line), and
$\mu^{2}=1$ (continuous line).

Furthermore in four-dimensional de Sitter spacetime, for the massive Klein-
Gordon field with $q=0$ we find that its effective potential becomes (see
Table 1)

$V_{KG}=(\mu^{2}-2)(1-r^{2}),$ (8)

and we point out that this potential is finite at $r=0$. Also for $\mu^{2}<2$
the effective potential (8) is negative (see Fig. 2). As far as we have
analyzed there is no another effective potential of Table 1 that behaves in
this form for the allowed values of its physical parameters.111For the Klein-
Gordon field the effective potential (8) is similar to that for the same field
propagating in two-dimensional de Sitter spacetime [21]. The difference is
that the factor multiplying to $(1-r^{2})$ is $(\mu^{2}-2)$ in four dimensions
and $\mu^{2}$ in two dimensions. Notice that in two-dimensional de Sitter
spacetime the effective potential of the massive Klein-Gordon field is always
positive, in contrast to the behavior of the effective potential (8) (see Fig.
2).

For the effective potentials that behave in this way we must discuss whether
the boundary condition $U(0)=0$ can be imposed, and we believe that for these
examples we need a separate analysis. In what follows we do not analyze these
potentials, because we like to show some results about the de Sitter QNM that
satisfy the boundary conditions (a) and (b) enumerated above.

Establishing a bound on the imaginary parts of the QNF, Horowitz and Hubeny
prove that the QNM of the Schwarzschild anti-de Sitter black hole are stable
[18]. Although the exact values of the de Sitter QNF are known, and using
these we can discuss the classical stability of the QNM, it is convenient to
have an independent proof of their stability and prove it with a different
method. Hence we modify the method by Horowitz and Hubeny [18] to show the
classical stability of the de Sitter QNM.

Taking the time dependence as $\textrm{exp}(-i\omega t)$ and following
Horowitz and Hubeny we propose that the function $U$ of Eq. (2) satisfies the
de Sitter QNM boundary condition near the horizon, that is, this function
takes the form

$U=\textrm{e}^{i\omega r_{*}}\Phi,$ (9)

where the function $\Phi$ is a solution of the differential equation

$f\frac{\textrm{d}^{2}\Phi}{\textrm{d}r^{2}}+\left(\frac{\textrm{d}f}{\textrm{d}r}+2i\omega\right)\frac{\textrm{d}\Phi}{\textrm{d}r}-\frac{V}{f}\Phi=0,$
(10)

and it behaves as $\Phi(0)=0$ and near the horizon $\Phi$ goes to a constant.

As in Ref. [18] we multiply Eq. (10) by the complex conjugate of $\Phi$.222In
what follows $\Phi^{*}$ denotes the complex conjugate of $\Phi$. Integrating
the result we obtain

$\int_{0}^{1}\textrm{d}r\left[\Phi^{*}f\frac{\textrm{d}^{2}\Phi}{\textrm{d}r^{2}}+\left(\frac{\textrm{d}f}{\textrm{d}r}+2i\omega\right)\Phi^{*}\frac{\textrm{d}\Phi}{\textrm{d}r}-\frac{V}{f}\left|\Phi\right|^{2}\right]=0.$
(11)

Next we integrate by parts the first term of the expression (11) and observe
that the surface terms vanish. Thus we get that the previous formula becomes

$\int_{0}^{1}\textrm{d}r\left[\left|\frac{\textrm{d}\Phi}{\textrm{d}r}\right|^{2}f-2i\omega\Phi^{*}\frac{\textrm{d}\Phi}{\textrm{d}r}+\frac{V}{f}\left|\Phi\right|^{2}\right]=0.$
(12)

We take the imaginary part of the previous expression [18] to get

$\int_{0}^{1}\textrm{d}r\,\Phi^{*}\frac{\textrm{d}\Phi}{\textrm{d}r}=-\frac{\omega^{*}\left|\Phi(1)\right|^{2}}{\omega-\omega^{*}},$
(13)

and from the formulas (12) and (13) we find

$\int_{0}^{1}\textrm{d}r\left[f\left|\frac{\textrm{d}\Phi}{\textrm{d}r}\right|^{2}+\frac{V}{f}\left|\Phi\right|^{2}\right]=2i\omega\int_{0}^{1}\textrm{d}r\,\Phi^{*}\frac{\textrm{d}\Phi}{\textrm{d}r}=-\frac{\left|\omega\right|^{2}\left|\Phi(1)\right|^{2}}{{\mathbb{I}}{\mathrm{m}}(\omega)}.$
(14)

Hence if $V\geq 0$ then ${\mathbb{I}}{\mathrm{m}}(\omega)<0$ and the de Sitter
QNM are stable. Thus for the fields of Table 1 we study whether their
effective potentials are nonnegative to prove the classical stability of their
de Sitter QNM. For the gravitational and electromagnetic perturbations it is
possible to show that for the allowed values of the parameters, their
effective potentials are always positive (see Table 1). Thus for the
gravitational and electromagnetic fields we assert that their de Sitter QNM
are stable.

For the Klein-Gordon and Dirac fields propagating in de Sitter spacetime their
effective potentials are not strictly positive. For example, for the massless
Klein-Gordon field with $q=0$, when $D>4$, if $1>r>\sqrt{\frac{D-4}{D}}$ then
$V_{KG}<0$. Moreover, for the massless Dirac field, if $\lambda>1$,
$1>r^{2}>1-\frac{1}{\lambda^{2}}$ then $V_{+}<0$.333In four-dimensional de
Sitter spacetime, for the Klein-Gordon with $q=0$ and $\mu^{2}<2$, we find
$V_{KG}<0$ in $r\in(0,1)$ (see Fig. 2) and for the massless Dirac field with
$\lambda=1$ we get that $V_{+}<0$ for $r\in(0,1)$ (see Fig. 1). Nevertheless
note that these effective potentials do not diverge at $r=0$.  Notice that for
the massless Klein-Gordon field only for the value $q=0$ the effective
potential is not strictly positive, whereas for the massless Dirac field the
effective potential $V_{+}$ is not strictly positive for the allowed values of
$\lambda$. (See Figs. 3 and 4 for some examples.)

Figure 3: For the five-dimensional de Sitter spacetime we draw the effective
potential $V_{KG}$ of the Klein-Gordon field with $q=0$ and $\mu=0$. Figure 4:
For the five-dimensional de Sitter spacetime we plot the effective potential
$V_{+}$ of the massless Dirac field with $\lambda=3/2$.

From the formula (14), for these values of the parameters for the Klein-Gordon
and Dirac fields we can not assert that the corresponding de Sitter QNM are
stable. Hence we need to modify the previous method. We follow to Ishibashi
and Kodama [22] to find

$\int_{0}^{1}\textrm{d}r\left[f\left|\frac{\textrm{d}\Phi}{\textrm{d}r}\right|^{2}+\frac{V}{f}\left|\Phi\right|^{2}\right]=\int_{0}^{\infty}\textrm{d}r_{*}\left[\left|\frac{\textrm{d}\Phi}{\textrm{d}r_{*}}\right|^{2}+V\left|\Phi\right|^{2}\right]=\int_{0}^{\infty}\textrm{d}r_{*}\Phi^{*}\left[-\frac{\textrm{d}^{2}}{\textrm{d}r_{*}^{2}}+V\right]\Phi.$
(15)

We notice that in the third step the boundary terms cancel. In this formula
and in what follows $r_{*}$ is the tortoise coordinate of the de Sitter
spacetime.

Defining the differential operator [22]

$\tilde{D}=\frac{\textrm{d}}{\textrm{d}r_{*}}+S,$ (16)

where $S$ is a function of $r$ such that $S(1)=0$, we get from the expression
(15)

$\displaystyle\int_{0}^{\infty}\textrm{d}r_{*}\Phi^{*}\left[-\frac{\textrm{d}^{2}}{\textrm{d}r_{*}^{2}}+V\right]\Phi$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\textrm{d}r_{*}\left[(\tilde{D}\Phi)^{*}(\tilde{D}\Phi)+\Phi^{*}\left(V+\frac{\textrm{d}S}{\textrm{d}r_{*}}-S^{2}\right)\Phi\right]$
(17) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\textrm{d}r_{*}\left[(\tilde{D}\Phi)^{*}(\tilde{D}\Phi)+\Phi^{*}\tilde{V}\Phi\right],$

with $\tilde{V}$ given by

$\tilde{V}=V+\frac{\textrm{d}S}{\textrm{d}r_{*}}-S^{2}.$ (18)

Thus for a field such that its effective potential $V$ is not strictly
positive, if we find a new potential $\tilde{V}$ such that $\tilde{V}\geq 0$,
then we show the stability of the corresponding de Sitter QNM.

For example, for the Klein-Gordon field we choose [22]

$S=-\frac{(D-2)f}{2r},$ (19)

to obtain

$\tilde{V}=f\left(\frac{q(D+q-3)}{r^{2}}+\mu^{2}\right).$ (20)

Since for the allowed values of the physical parameters in de Sitter spacetime
this new potential satisfies $\tilde{V}\geq 0$, we find that for the Klein-
Gordon field its de Sitter QNM are stable. Notice that this result is valid
for the effective potentials of the massive and massless Klein-Gordon fields.

For the effective potential $V_{+}$ of the massless Dirac field we choose

$S=-\frac{\lambda f^{1/2}}{r},$ (21)

to get that the new potential is equal to444Writing the effective potentials
of the massless Dirac field as [19]
$V_{\pm}=W^{2}\pm\frac{\textrm{d}W}{\textrm{d}r_{*}},$ where $W=\lambda
f^{1/2}/r$, we obtain that the function $S$ is equal to $S=-W$.

$\tilde{V}_{+}=0.$ (22)

Thus the formula (17) simplifies to

$\int_{0}^{\infty}\textrm{d}r_{*}\left|(\tilde{D}\Phi)\right|^{2}>0,$ (23)

and therefore the de Sitter QNM of the massless Dirac field are stable.

We note that the functions (19) and (21) satisfy $S(1)=0$. Furthermore the
results (20) and (22) are valid for an arbitrary metric function $f$ that
satisfies $f(1)=0$, $f(0)$ is a finite quantity, and $f$ is positive in
$r\in(0,1)$. In particular these results are valid for the de Sitter
spacetime.

It is convenient to notice that for the Klein-Gordon, Dirac, electromagnetic,
and gravitational perturbations our results on the classical stability of
their de Sitter QNM are valid for $D\geq 4$. In contrast, it is known that for
asymptotically de Sitter spherically symmetric black holes, their classical
stability depends on the spacetime dimension in a non trivial way [22]-[26].
For example, only for $D=4,5,6$ we know analytical proofs of the classical
stability of the $D$-dimensional Schwarzschild de Sitter spacetime [22], [23].
For $D=5,\dots,11$ there are numerical calculations that state their classical
stability [24]. We do not know analytical or numerical calculations
determining the classical stability of Schwarzschild de Sitter black holes for
$D\geq 12$.

Furthermore, based on numerical calculations, in [25] it is showed that for
$D=7,\dots,11$ and for sufficiently large charge and cosmological constant,
the Reissner-Nordström de Sitter black holes are unstable under coupled
electromagnetic and gravitational scalar type perturbations. Only for $D=4,5$
we are aware of analytical proofs of the classical stability for Reissner-
Nordström de Sitter black holes [22], [23]. For additional arguments about the
instability of these black holes see Ref. [26]. Hence the classical stability
of the de Sitter QNM is stated in a more general way than the stability of the
$D$-dimensional Schwarzschild de Sitter and Reissner-Nordström de Sitter black
holes.

##  3 Gravitational QNM of the de Sitter spacetime

As we comment in the previous sections, the QNM of the $D$-dimensional de
Sitter spacetime ($D\geq 4$) are determined in [7]-[13] and in these
references, on their existence we find at least two different results. To
settle this problem, based on Refs. [18], [27], [28] and with a different
method to that of Refs. [10], [12], in what follows we calculate the QNF of
the three types of gravitational perturbations propagating in $D$-dimensional
de Sitter spacetimes. Notice that the method we use in this work is simpler
than the procedure of Refs. [10], [12], and we hope that it will be useful to
solve the mentioned issue. We choose the gravitational perturbations due to
its intrinsic physical relevance and since these perturbations are analyzed in
[10], [12].

First we note that in $D$-dimensional de Sitter spacetime the Schrödinger type
equations (2) for the three types of gravitational perturbations simplify to
the common form [10]

$(1-r^{2})\frac{{\rm d}^{2}U}{{\rm
d}r^{2}}-2r\frac{\textrm{d}U}{\textrm{d}r}+\left(\frac{\omega^{2}}{1-r^{2}}-\frac{B(B+1)}{r^{2}}+A\right)U=0,$
(24)

where

$\displaystyle B=\frac{2q+D-4}{2},\qquad A=\left\\{\begin{array}[]{l
l}\frac{(D-2)D}{4}&\,\,\,\,\textrm{tensor type},\\\ \\\
\frac{(D-4)(D-2)}{4}&\,\,\,\,\textrm{vector type},\\\ \\\
\frac{(D-6)(D-4)}{4}&\,\,\,\,\textrm{scalar type}.\end{array}\right.$ (30)

To solve Eq. (24) we make the change of variable $y=r^{2}$ [10], [12] to
obtain

$y(1-y)\frac{\textrm{d}^{2}U}{\textrm{d}y^{2}}+\left(\frac{1}{2}-\frac{3y}{2}\right)\frac{\textrm{d}U}{\textrm{d}y}+\left(\frac{\omega^{2}}{4(1-y)}-\frac{B(B+1)}{4y}+\frac{A}{4}\right)U=0.$
(31)

Proposing that the function $U$ takes the form

$U=(1-y)^{E}y^{C}U_{2},$ (32)

we find that if the quantities $E$ and $C$ are solutions of the equations

$E^{2}+\frac{\omega^{2}}{4}=0,\qquad\quad
C^{2}-\frac{C}{2}-\frac{B(B+1)}{4}=0,$ (33)

then the function $U_{2}$ must be a solution of the differential equation

$\displaystyle y(1-y)\frac{\textrm{d}^{2}U_{2}}{\textrm{d}y^{2}}$
$\displaystyle+$
$\displaystyle\left(2C+\frac{1}{2}-\left(2E+2C+\frac{3}{2}\right)y\right)\frac{\textrm{d}U_{2}}{\textrm{d}y}$
(34) $\displaystyle-$
$\displaystyle\left(C^{2}+2EC+\frac{C}{2}+E^{2}+\frac{E}{2}-\frac{A}{4}\right)U_{2}=0,$

which is a hypergeometric type differential equation [29], [30]

$y(1-y)\frac{{\rm d}^{2}U_{2}}{{\rm d}y^{2}}+(c-(a+b+1)y)\frac{{\rm
d}U_{2}}{{\rm d}y}-abU_{2}=0.$ (35)

Motivated by the Horowitz-Hubeny method [18], we choose $E=-i\omega/2$ to
satisfy the boundary condition (a) of the de Sitter QNM near the horizon and
$C=(1+B)/2$ to fulfill the boundary condition (b) at $r=0$. Making these
choices for $E$ and $C$ we get that the parameters $a$, $b$, and $c$ of Eq.
(35) are equal to

$\displaystyle a$
$\displaystyle=\frac{B}{2}+\frac{3}{4}-\frac{1}{2}\left(\frac{1}{4}+A\right)^{1/2}-\frac{i\omega}{2},\quad\quad
c=B+\frac{3}{2},$ $\displaystyle b$
$\displaystyle=\frac{B}{2}+\frac{3}{4}+\frac{1}{2}\left(\frac{1}{4}+A\right)^{1/2}-\frac{i\omega}{2}.$
(36)

As in Refs. [27], [28] to obtain the appropriate behavior at the boundaries we
impose that the function $U_{2}$ be a polynomial. Therefore we must satisfy
the conditions [29], [30]

$a=-p,\qquad\quad\textrm{or}\qquad\quad b=-p,\quad\qquad p=0,1,2,\dots$ (37)

From these equations we get that the QNF of the $D$-dimensional de Sitter
spacetime are equal to

$\omega=-i(q+D-1-\eta+2p),\qquad\qquad\omega=-i(q+\eta+2p),$ (38)

where

$\displaystyle\eta=\left\\{\begin{array}[]{l l}0,&\,\,\,\textrm{tensor
type},\\\ 1,&\,\,\,\textrm{vector type},\\\ 2,&\,\,\,\textrm{scalar
type}.\end{array}\right.$ (42)

We point out that the QNF (38) are identical to those of Refs. [10], [12] (see
for example the formula (37) of [12]). We expect that this method allows us to
calculate the QNF of other fields moving in de Sitter spacetime.

Notice that the method used here to calculate the QNF (38) is valid for even
and odd $D$, and hence these QNF are valid for all $D$. In contrast to the
result that we find in [10] asserting that the QNF (38) are valid only for $D$
odd. In Ref. [12] it is asserted that the QNF (38) are valid for even and odd
$D$, but in the procedure of this reference it is necessary to study different
solutions for $D$ even and $D$ odd. Furthermore in the method of Ref. [12] the
employed solutions and the analysis of the problem are more complicated. We
expect that our previous results establish the existence of well defined QNF
of the gravitational perturbations in even and odd de Sitter spacetimes.

Furthermore we note that our conclusion on the existence of well defined de
Sitter QNF in even and odd dimensions it is valid for the massless Klein-
Gordon field, since in the Schrödinger type equations the effective potential
for this field is equal to that of the tensor type gravitational perturbations
[10]. In contrast to gravitational perturbations for which $q\geq 2$, for the
Klein-Gordon field the angular eigenvalues $q$ take the values
$q=0,1,2,\dots$, but this fact does not modify our calculation method or
change the result for the de Sitter QNF (38). Hence the de Sitter QNF of the
massless Klein-Gordon field are well defined in even and odd dimensions.

Making the change of variable

$y=\frac{1+x}{2},$ (43)

we find that for the chosen values of the quantities $C$ and $E$, Eq. (34)
becomes

$\displaystyle(1-x^{2})\frac{\textrm{d}U_{2}}{\textrm{d}x^{2}}$
$\displaystyle+\left(B+\frac{1}{2}+i\omega-\left(\frac{3}{2}+B-i\omega+1\right)x\right)\frac{\textrm{d}U_{2}}{\textrm{d}x}$
$\displaystyle-\left(\frac{B(B+1)}{4}+\frac{2(B+1)(1-i\omega)}{4}-\frac{A+\omega^{2}+i\omega}{4}\right)U_{2}=0,$
(44)

where $x\in(-1,1)$. This equation is a differential equation of Jacobi type
[30], [31]

$(1-x^{2})\frac{\textrm{d}^{2}U_{n}}{\textrm{d}x^{2}}+(\beta-\alpha-(\alpha+\beta+2)x)\frac{\textrm{d}U_{n}}{\textrm{d}x}+n(n+\alpha+\beta+1)U_{n}=0,$
(45)

with parameters $\alpha$, $\beta$, and $n$ equal to

$\alpha=-i\omega,\quad\qquad\beta=B+\frac{1}{2},\quad\qquad
n=\frac{i\omega-B-3/2\pm(1/4+A)^{1/2}}{2}.$ (46)

For $n$ a nonnegative integer, from the last formula of (46) we get that the
de Sitter QNF (38). Furthermore from the QNF for the three types of
gravitational perturbations we obtain that the parameter $\alpha$ satisfies
$\alpha\leq-2$. These values that we find for $\alpha$ have a noticeable
consequence.

Let $U_{n_{1}}$ and $U_{n_{2}}$ two solutions of Eq. (45) for different
eigenvalues. Following the usual procedure to show the orthogonality of the
eigenfunctions [30], we arrive to the formula

$\int_{-1}^{1}(1-x)^{\alpha}(1+x)^{\beta}U_{n_{1}}U_{n_{2}}\textrm{d}x=\frac{(1-x)^{\alpha+1}(1+x)^{\beta+1}\left[U_{n_{2}}\frac{\textrm{d}U_{n_{1}}}{\textrm{d}x}-U_{n_{1}}\frac{\textrm{d}U_{n_{2}}}{\textrm{d}x}\right]^{1}_{-1}}{n_{2}(n_{2}+\alpha+\beta+1)-n_{1}(n_{1}+\alpha+\beta+1)},$
(47)

which is divergent at $x=1$ for $\alpha<-1$, as for the de Sitter radial
functions that we find above.

Hence, the radial functions of the de Sitter QNM are not orthogonal solutions
of Eq. (45) in the interval $(-1,1)$ with weight $(1-x)^{\alpha}(1+x)^{\beta}$
and they can not be the usual Jacobi polynomials which are orthogonal with the
previous weight. In Ref. [14] it is said that the de Sitter QNM radial
functions involve Jacobi polynomials, without additional comments or
discussion. As we say above the radial functions of the de Sitter QNM do not
involve the usual Jacobi polynomials [30], [31]. As far as we know this fact
is not noted in previous papers.

##  4 Summary

Based on the work by Horowitz and Hubeny, for the Klein-Gordon, Dirac,
electromagnetic, and gravitational perturbations we calculate a bound for the
imaginary parts of their de Sitter QNF and using this bound we prove the
stability of their de Sitter QNM. Without a doubt, we believe that this
alternative proof shows the stability of the de Sitter QNM. Notice that for
the Klein-Gordon and Dirac fields we employ a slightly different method that
for the electromagnetic and gravitational perturbations.

Furthermore, we remark that for the de Sitter QNM (satisfying the boundary
conditions (a) and (b) of Sect. 2) of the Klein-Gordon, Dirac,
electromagnetic, and gravitational perturbations, our result on their
classical stability is valid for all the spacetime dimensions satisfying
$D\geq 4$.

Moreover, as a noticeable example, we obtain that for $q=0$ and $\mu^{2}=2$
the Klein-Gordon field propagates freely in the static patch of the four-
dimensional de Sitter spacetime (1) (see Sect. 2 and Appendix). Recently for
odd-dimensional de Sitter spacetimes in global coordinates, it is found that
several fields propagate freely when they satisfy some conditions [32]. Our
result shows that for $q=0$ and $\mu^{2}=2$ the Klein-Gordon field behaves in
a similar way in the static patch of the four-dimensional de Sitter spacetime.

Using a simple method we calculate the de Sitter QNF of the gravitational
perturbations. Our results coincide with those of Ref. [12], but the method of
Sect. 3 is simpler than the procedure of Refs. [10], [12]. We expect that
these calculations establish, in odd and even spacetime dimensions, the
existence of well defined de Sitter QNF for massless fields.

We also find that for the gravitational perturbations, the radial functions of
their de Sitter QNM are solutions of the Jacobi differential equation, but
these radial functions are not orthogonal in the interval $(-1,1)$ with the
usual weight $(1-x)^{\alpha}(1+x)^{\beta}$, since for the de Sitter QNF (38)
we get that the parameter $\alpha$ satisfies $\alpha\leq-2$. We think that
this fact deserves further research.

##  5 Acknowledgments

This work was supported by CONACYT México, SNI México, EDI-IPN, COFAA-IPN, and
Research Projects SIP-20120773 and SIP-20121648.

##  Appendix A On the de Sitter QNF of the Klein-Gordon field with $q=0$,
$\mu^{2}=2$, and $D=4$

For the Klein-Gordon field with $q=0$ and propagating in four-dimensional de
Sitter spacetime, from the formula (8) we observe that if $\mu^{2}=2$ the
effective potential of the Schrödinger type equation goes to zero. Therefore
the solutions of Eq. (2) take the form (4) for all $r\in(0,1)$, and hence we
can not satisfy the boundary conditions of the de Sitter QNM. Thus for these
values of the parameters we find that the formula (51) of Ref. [12] does not
give valid QNF of the de Sitter spacetime.

To see where the procedure of Ref. [12] fails for these values of the physical
parameters, we note that for $q=0$ the radial function regular at $r=0$ is
(see the formulas (46) and (48) of [12])

$R=C_{1}(1-y)^{i\omega/2}{}_{2}F_{1}(a,b;c;y),$ (48)

with

$a=\frac{i\omega}{2}+1,\qquad\qquad
b=\frac{i\omega}{2}+\frac{1}{2},\qquad\qquad c=\frac{3}{2}.$ (49)

We point out that the quantities $a$, $b$, and $c$ of the formulas (49)
satisfy

$a-b+1=c,$ (50)

and therefore we write the radial function (48) as (see the formulas (15.4.14)
and (8.6.9) of [29])

$\displaystyle R$
$\displaystyle=C_{1}\Gamma(3/2)y^{-1/4}(1-y)^{-1/2}P^{-1/2}_{-i\omega/2-1/2}\left(\tfrac{1+y}{1-y}\right)$
$\displaystyle=-\frac{C_{1}}{2i\omega
y^{1/2}}\left((1+y^{1/2})^{-i\omega}(1-y)^{i\omega/2}-(1+y^{1/2})^{i\omega}(1-y)^{-i\omega/2}\right),$
(51)

where $\Gamma$ denotes the gamma function and $P^{a}_{b}$ the associated
Legendre functions.

From the last expression we find that near the horizon the radial function
(48) contains ingoing and outgoing waves. Notice that we can not cancel the
ingoing part of the radial function without cancel the outgoing part. Thus in
four-dimensional de Sitter spacetime, for the Klein-Gordon field with these
values of $q$ and $\mu$, the method of Ref. [12] does not work, since we can
not satisfy the boundary conditions of the de Sitter QNM.

It is convenient to note that the absence of the de Sitter QNF of the Klein-
Gordon field for $q=0$, $\mu^{2}=2$, and $D=4$ is an exceptional situation,
and it is caused by the fact that for these values of $q$, $\mu$, and $D$ the
effective potential $V$ of the Klein-Gordon field is equal to zero. As far as
we can see, only for these values of the physical parameters we expect that
the de Sitter QNM do not exist.

## References

  * [1] M. Spradlin, A. Strominger and A. Volovich, arXiv:hep-th/0110007.
  * [2] Y. b. Kim, C. Y. Oh and N. Park, Jour. of the Korean Phys. Soc. 42, 573 (2003) [arXiv:hep-th/0212326].
  * [3] K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999) [arXiv:gr-qc/9909058].
  * [4] H. P. Nollert, Class. Quantum Grav. 16, R159 (1999).
  * [5] E. Berti, V. Cardoso and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009) [arXiv:0905.2975 [gr-qc]].
  * [6] R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011) [arXiv:1102.4014 [gr-qc]].
  * [7] D. P. Du, B. Wang and R. K. Su, Phys. Rev. D 70, 064024 (2004) [arXiv:hep-th/0404047].
  * [8] A. López-Ortega, Gen. Rel. Grav. 38, 743 (2006) [arXiv:gr-qc/0605022].
  * [9] A. Zelnikov, JHEP 0807, 010 (2008) [arXiv:0805.4031 [hep-th]].
  * [10] J. Natario and R. Schiappa, Adv. Theor. Math. Phys. 8 (2004) 1001 [arXiv:hep-th/0411267].
  * [11] T. R. Choudhury and T. Padmanabhan, Phys. Rev. D 69, 064033 (2004) [arXiv:gr-qc/0311064].
  * [12] A. López-Ortega, Gen. Rel. Grav. 38 (2006) 1565 [arXiv:gr-qc/0605027].
  * [13] A. López-Ortega, Gen. Rel. Grav. 39 (2007) 1011 [arXiv:0704.2468 [gr-qc]].
  * [14] D. Anninos, S. A. Hartnoll and D. M. Hofman, Class. Quant. Grav. 29, 075002 (2012) [arXiv:1109.4942 [hep-th]].
  * [15] F. Denef, S. A. Hartnoll and S. Sachdev, Class. Quant. Grav. 27, 125001 (2010) [arXiv:0908.2657 [hep-th]].
  * [16] A. Lopez-Ortega, Phys. Lett. B 682, 85 (2009) [arXiv:0910.5779 [gr-qc]].
  * [17] A. Lopez-Ortega, Int. J. Mod. Phys. D 20, 2525 (2011) [arXiv:1112.6211 [gr-qc]].
  * [18] G. T. Horowitz, V. E. Hubeny, Phys. Rev. D 62, 024027 (2000). [hep-th/9909056].
  * [19] S. Chandrasekhar, The Mathematical Theory of Black Holes, (Oxford University Press, Oxford, 1983).
  * [20] A. Lopez-Ortega, Lat. Am. J. Phys. Educ. 3, 578 (2009) [arXiv:0906.2754 [gr-qc]].
  * [21] A. Lopez-Ortega and I. Vega-Acevedo, Gen. Rel. Grav. 43, 2631 (2011) [arXiv:1105.2802 [gr-qc]].
  * [22] A. Ishibashi and H. Kodama, Prog. Theor. Phys. 110, 901 (2003) [arXiv:hep-th/0305185].
  * [23] A. Ishibashi and H. Kodama, Prog. Theor. Phys. Suppl. 189, 165 (2011) [arXiv:1103.6148 [hep-th]].
  * [24] R. A. Konoplya and A. Zhidenko, Nucl. Phys. B 777, 182 (2007) [hep-th/0703231 [HEP-TH]].
  * [25] R. A. Konoplya and A. Zhidenko, Phys. Rev. Lett. 103, 161101 (2009) [arXiv:0809.2822 [hep-th]].
  * [26] V. Cardoso, M. Lemos and M. Marques, Phys. Rev. D 80, 127502 (2009) [arXiv:1001.0019 [gr-qc]].
  * [27] S. Musiri and G. Siopsis, Phys. Lett. B 563, 102 (2003) [arXiv:hep-th/0301081].
  * [28] S. Musiri and G. Siopsis, Phys. Lett. B 576, 309 (2003) [arXiv:hep-th/0308196].
  * [29] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Graphs, and Mathematical Table, (Dover Publications, New York, 1965).
  * [30] Z. X. Wang and D. R. Guo, Special Functions, (World Scientific Publishing, Singapore, 1989).
  * [31] G. Szegö, Orthogonal Polynomials, (American Mathematical Society, Providence, Rhode Island, 1975).
  * [32] P. Lagogiannis, A. Maloney and Y. Wang, arXiv:1106.2846 [hep-th].