# Characterization of an aerosolized nanoparticle beam beyond the diffraction
limit through strong field ionization

Michael Davino, Tobias Saule, Nora G. Helming, J.A. Powell, Carlos Trallero-
Herrero*
Department of Physics
University of Connecticut
Storrs, CT 06269
<EMAIL_ADDRESS>

###### Abstract

The study of nanomaterials is an active area of research for technological
applications as well as fundamental science. A common method for studying
properties of isolated nanoparticles is by an in-vacuum particle beam produced
via an aerodynamic lens. Despite being common practice, characterization of
such beams has proven difficult as light scattering detection techniques fail
for particles with sizes beyond the diffraction limit. Here we present a new
technique for characterizing such nanoparticle beams using strong field
ionization. By focusing an ultrafast, mJ-level laser into the particle beam, a
nanoparticle within the laser focus is ionized and easily detected by its
ejected electrons. This method gives direct access to the nanoparticle density
at the location of the focus and, by scanning the focus through the transverse
and longitudinal profiles of the particle beam, the 3-dimensional particle
density distribution can be attained. Further, we show that strong field
ionization is effective in detecting spherical nanoparticles as small as 10 nm
in diameter. Additionally, this technique is an effective tool in optimizing
the particle beam for specific applications. As an example we show that the
particle beam density and width can be manipulated by restricting the gas flow
into the aerodynamic lens.

##  1 Introduction

Nanoparticles (NPs) and their associated technologies are active areas of
research across multiple scientific disciplines with important applications in
medicine [1, 2, 3], chemical analysis [4, 5], and solar energy [6, 7]. Apart
from the variety of specific applications, the optical properties of NPs are
an area of interest for fundamental scientific research [4, 8, 9]. While much
of the research and applications involve NPs embedded in substrates [10] or in
dense aqueous suspensions [7], there is a need to study isolated NPs [11]. In
the case of studying isolated aerosolized NPs (as compared to isolated NPs in
a suspension) producing a target sample is particularly difficult. One common
method for generating such a sample is by aerosolization of a NP suspension
and subsequent collimation of the aerosol into an in-vacuum particle beam by
means of an aerodynamic lens [12]. This NP beam serves as the target for
study.
Though NP beams as produced by an aerodynamic lens have been thoroughly
investigated numerically [13, 14, 15], experimental characterization has
proven difficult as the various characterization techniques come with a range
of limitations. In particular, light scattering techniques fail for beams
comprised of NPs beyond the diffraction limit of the light source. Early
characterizations were performed by detecting NP ions via a Faraday cup
electrometer and using a knife edge to physically impede the NP beam, thus
giving access to the transverse spatial distribution of the NPs [12, 16].
While this method can also measure the particle velocities by applying a
transverse electric field to deflect the ionized NP beam, it has the obvious
drawback of requiring the NPs to be ionized. More recently light-sheet imaging
(LSI) has been used to directly image NPs in a beam by light scattered from a
CW laser source [17]. This technique measures both the spatial distribution of
the NPs and the exact number of particles per exposure time, all in a non-
disruptive manner. However, being a scattering based technique, LSI has the
significant drawback that it is limited to larger NPs (> 100 nm in [17]) due
to diffraction limitations. Finally, laser-induced breakdown detection (LIBD)
has also been utilized to characterize NP beams [18]. The LIBD technique
employs a pulsed laser focused into the NP beam and detects the presence of a
NP in the focus by light emitted from its plasma. This detection scheme can be
used to spatially resolve the number of NPs in a beam comparable to the
capabilities of LSI but with the demonstrated capacity to detect particles as
small as 5 nm [19].
Here we present a new technique for acquiring spatially resolved NP density
information for an aerosolized NP beam. By focusing an ultrafast mJ-level
pulsed laser into a NP beam electrons are ejected via strong field ionization
(SFI) and are easily detectable by a channel electron multiplier (CEM).
Distinguishing laser pulses with a NP present in the focus gives direct access
to the NP density of the beam at that position. Further, by scanning the laser
focus through the transverse profile of the NP beam at multiple points along
its axis of propagation, we attain the 3-dimensional NP density distribution
of the particle beam. We show this scheme to be effective for spherical NPs of
various materials and sizes, with NP diameters ranging from 10 nm – 120 nm.
SFI is shown here to be a useful tool in monitoring the output of an
aerodynamic lens setup, which is especially valuable in ensuring optimal and
consistent experimental conditions. As an example, we use the technique to
demonstrate that a NP beam may be manipulated by adjusting the flow of aerosol
input to the aerodynamic lens, effectively allowing for the NP beam to be
tailored to specific applications.

##  2 Experimental Setup

A schematic of the experimental setup is depicted in Figure 1 (a). A well-
established procedure for producing a collimated NP beam is utilized: an
aqueous NP suspension is aerosolized by a syringe pump fed atomizer, dried to
remove excess water, and injected into an aerodynamic lens where a series of
apertures collimates the dry NPs into a beam. The setup used here is similar
to the setup used by Powell et al.[20].

Figure 1: (a) depicts the experimental setup. An atomizer is fed an aqueous NP
suspension by a constant output syringe pump and the solution is aerosolized
by nitrogen gas at a fixed pressure. The aerosol then passes through a series
of dryers consisting of a Nafion™ membrane dryer and two silica bead dryers.
The tubing that couples the output of the dryers to the aerodynamic lens is
equipped with a pinch bolt that allows for choking the flow of gas and
nanoparticles into the lens. The aerodynamic lens is kept at rough vacuum and
focuses the dry nanoparticles into a beam which is sent into the interaction
chamber which is equipped with an ion gauge and CEM. Outside of vacuum the
output of a Titanium:Sapphire laser system is incident on a focusing lens
mounted on a translational stage. The focusing beam reflects off a silver
mirror mounted in a motorized mount and is sent into the interaction chamber.
In vacuum the beam is steered by a series of dielectric mirrors mounted on a
translational stage such that it focuses onto the NP beam a variable distance
L from from the output of the aerodynamic lens. The laser is focused onto the
nanoparticle beam a distance L from the exit of the aerodynamic lens. A HeNe
laser is reflected off the motorized mirror and imaged with a CCD camera to
track the position of the mirror. (b) shows a plot of the raw CEM signal for
two laser pulses, one with a NP present in the interaction region and one
without.

Two atomizers were employed in this work; the first is a design based on a TSI
Inc. aerosol generator (model 3076) which produces an aerosol with an average
droplet size of 0.3 $\mu$m. This atomizer uses impaction to remove large
droplets from the final aerosol and thus a particularly fine mist is attained,
however this means that much of the aerosol recondenses. For stable NP
suspensions (e.g. SiO2) this suspension can be collected and recycled back
through the atomization process, but for less stable NPs the impaction process
can be a catalyst for the NPs to aggregate. In general, we have found that
metallic NPs are particularly susceptible to such aggregation and thus, for
metallic samples a commercial, needle based, atomizer from Fuso Seiki Co.
(Model ST-10SK-0.5X-G2) is used. This atomizer does not use an impactor to
discriminate against large aerosol droplets, thus the average droplet size is
larger, on the order of tens of microns, but allows for the use of the entire
liquid sample without concerns of aggregation. Both atomizers were backed with
30 psi of nitrogen gas. The flow rate of liquid sample into the atomizers is
kept constant by a syringe pump and set to 1 mL/min for the TSI atomizer and
0.1 mL/min for the Fuso atomizer.
The output of either atomizer is coupled directly into a series of dryers used
to remove the solvent (water) from the aerosol and leaving only the NPs
suspended in gas to enter the aerodynamic lens. The first drying stage is a
cylindrical Nafion™ membrane (5 inches diameter, 18 inches length) which
operates by proton transport across the membrane assisted by a counterflow dry
gas (compressed air, < 5% relative humidity). This membrane dryer is followed
by two silica bead diffusion dryers. The relative humidity of the suspended
NPs is measured following the three dryers and is typically < 5% for the TSI
atomizer and < 15% for the Fuso atomizer.
The output of the dryers is coupled into the aerodynamic lens by an orifice
(denoted ‘critical orifice’ in Figure 1) that has an inner diameter of 330
$\mu$m. The flow of NPs through the critical orifice can be limited by a pinch
bolt installed on the hosing between the silica dryers and the orifice.

Material | Size (nm) | Mass (g) | Aqueous Concentration (NP/mL) | Manufacturer
---|---|---|---|---
SiO2 | 60 | 3.0e-16 | 3.3e14 | General Research and Engineering
SiO2 | 120 | 2.4e-15 | 4.2e13 | General Research and Engineering
Ag | 10 | 5.5e-18 | 4.0e12 | Cytodiagnostics
Au | 120 | 1.8e-14 | 2.1e09 | Cytodiagnostics
ZnO | 30 | 7.9e-17 | 2.8e15 | Nanografi
TiO2 | 33 | 8.0e-17 | 4.4e15 | Nanografi
Fe3O4 | 16 | 1.1e-17 | 2.1e16 | Nanografi
Table 1: Detailed NP information for the samples used. The concentration
presented is the commercially available concentration, some samples were
further diluted and are noted as such in the text.

The aerodynamic lens consists of a series of six apertures of decreasing
diameter. The pressure gradient across the lens is established by 1 atm input
into the critical orifice and between 100 and 500 mTorr maintained at the exit
of the lens by an Edwards Vacuum Drystar roughing pump and booster combination
(Models QDP80 and QMB500F respectively). The NPs are collimated as they move
through the aerodynamic lens and the output beam is sent through a 2 mm
diameter skimmer into the interaction chamber which is kept between 8e-5 and
3e-3 Torr. The pressure in the interaction chamber is entirely dependent on
the setting of the pinch bolt, and thus the ion gauge pressure sensor in the
interaction chamber serves as a relative measure of the gas flow (and
subsequently, NPs) into the aerodynamic lens. At a set position of the pinch
bolt the chamber pressure is stable to ± 5e-6 Torr.
The interaction chamber is equipped with a continuous dynode channel electron
multiplier (CEM) set to detect electrons. The output CEM signal is put through
a boxcar integrator with a gate width of 15 $\mu$s and subsequently recorded
by a National Instruments measurement box which is controlled via LabVIEW.
The laser used here is a Titanium:Sapphire system producing 1.5 mJ pulses at a
rate of 1 kHz (1.5 W) with a duration of 30 fs at a central wavelength of 790
nm. The beam is sent through a 500 mm focusing lens mounted on a translational
stage and subsequently reflected by a silver mirror mounted on a 1D motorized
mount which injects the beam into the interaction chamber. In vacuum the beam
is steered by a series of dielectric mirrors mounted on a translational stage
such that it focuses onto the NP beam a variable distance L from from the
output of the aerodynamic lens. The focal spot size was measured to be 45.8 ±
0.5 $\mu$m 1/e2 diameter (peak intensity on the order of 1e15 W/cm2).
The motorized mirror is used to scan the laser focus through the transverse
profile of the NP beam. A Helium-Neon (HeNe) laser is utilized as a reference
to track the position of the motorized mirror by detecting its reflection off
the mirror with a CCD camera. By measuring the distances from the motorized
mirror to both the focus and the camera, the relative displacement of the
focus with respect to the NP beam may be attained.
The NP samples used in this study are commercially available products
purchased as an aqueous suspension and can be diluted to control the NP
concentration. Complete details regarding the NP samples can be found in Table
1.

##  3 Results

A critical prerequisite in characterizing a NP beam by means of SFI is the
ability to identify laser shots for which a NP is present in the interaction
region. Inset Figure 1 (b) shows the raw CEM signal for both the presence and
absence of a NP in the interaction region. There is a clear distinction seen,
but it is important to realize that not every NP ‘hit’ will necessarily
generate as strong of a CEM signal. The number of ejected electrons has been
shown to be strongly dependent on the peak intensity of the laser pulse
impingent on the NP[20]. Further, because the NP is several orders of
magnitude smaller than the focal spot size and thus the local peak intensity
will vary with position of the NP within the focal region[21]. This means that
the integrated CEM signals that are measured here (now referred to as simply
the CEM signal) will constitute a continuous spectrum meaning that the
distinction between the presence and absence of a NP in the focal region is
non-trivial.
To better identify the presence of a NP in the interaction region we have
performed concentration dependent NP studies which are presented in Figure 2.
The laser focus was positioned at L = 5.5 cm and aligned to maximize the
number of clearly distinguishable NP hits in order to focus (approximately) in
the center of the NP beam. Figure 2 (a) is a histogram of the CEM signal for
1e5 laser shots focused into a NP beam consisting of 120 nm SiO2 at three NP
suspension concentrations spanning from 4e11 to 4e13 NP/mL as well as a
solvent (water) only sample. Note that the artifact present in the solvent
only data at a signal value of $\sim$0.5 in Figure 2 (a) is known to be due to
transient noise in the CEM signal line. Because of this, subsequent data
presented in this manuscript rejects laser pulses that records a CEM signal of
this value.

Figure 2: Integrated CEM signal for various suspension concentrations of 120
nm SiO2. (a) is a histogram of 1e5 laser shots for three NP concentrations as
well as a solvent (water) only sample (i.e. 0 NP/mL). The dashed vertical red
line indicates the threshold chosen to distinguish the presence or absence of
a NP in the interaction region, CEM signals larger than the threshold are
considered a NP hit. (b) plots the number of NP hits, normalized per laser
shot, as a function of atomized (aqueous) NP sample concentration. Note the
artifact present in the solvent only data set at a signal value of $\sim$0.5
in (a) is known to be due to transient noise in the CEM signal line. Error
bars represent the standard deviation within the data set which becomes
negligible at higher concentrations.

From Figure 2 (a) it is immediately evident that, when NPs are atomized,
larger CEM signals are observed. As this is the case, we then use the CEM
signal distribution attained from the solvent only and define a threshold that
will be used to identify a NP hit. This threshold is chosen as a CEM signal
value larger than those seen for the solvent only sample and is denoted by a
dashed vertical red line in Figure 2 (a). We then say that any laser shot
which produces a CEM signal larger than this threshold constitutes the
presence of a NP in the interaction region. Figure 2 (b) shows the number of
NP hits as a function of NP concentration produced by utilizing this analysis
technique. The observed relation fits with the expectation that higher
concentrations of the aqueous suspension directly correlates to a higher in-
vacuum hit rate.
We must now note some of the limitations associated with this NP detection
scheme. Examination of Figure 2 (a) shows that, with NPs being atomized, there
are a significant amount of laser shots for which the CEM signal is larger
than the average signal from water only, but smaller than the threshold value.
This is to say that there are a significant number of laser shots where there
is likely a NP present in the interaction region (or more likely, on its
edges) which are not being counted as such. This means that the number of NPs
counted by this scheme represents a lower bound to the actual NP count. This
is compounded with the realization that this scheme, as described, does not
differentiate the presence of a single NP or multiple NPs within the
interaction region for a single laser shot. Instead, we are counting a maximum
of a single NP per laser shot when there may be more. For these reasons the
number of NPs detected in-vacuum will be an undercount to the actual number.
With these limitations in mind, we now aim to quantify the NP beam width. This
is accomplished by scanning the laser focus through the transverse profile of
the NP beam using the motorized mirror which injects the laser into the
interaction chamber (see Figure 1). The reference HeNe allows for a change in
angle of the motorized mirror to be measured and subsequently correlated to a
relative displacement of the laser focus. By measuring the number of NP hits
per laser shot at many points along the transverse profile of the NP beam, the
relative density of the beam is effectively mapped in one dimension. Figure 3
shows a series of such scans for samples of SiO2, Ag, and Au. It is
immediately evident that at L = 5.5 cm each of the four samples produced a NP
beam, each with a width on the order of hundreds of microns. Further, the beam
width is largely preserved through L = 17.7 cm (L = 9 cm for Figure 3 (d)),
but the total number of NPs hit decreases as the laser focus moves further
from the exit of the aerodynamic lens. We also note that Figure 3 (a) and (b)
show that both 60 nm and 120 nm SiO2 retain high-quality beam profiles to L =
17.7 cm, whereas there is a clear degradation in the beam quality for the
metallic samples, 10 nm Ag and 120 nm Au, seen in Figure 3 (c) and (d)
respectively. This is particularly evident in Au, for which a beam profile was
entirely unattainable at L = 17.7 cm. Not presented are beam profiles for the
remaining samples listed in Table 1 (ZnO, TiO2, and Fe3O4) which were only
attained for L = 5.5 cm due to limited sample availability. The beam quality
for the samples not presented are comparable to Ag at L = 5.5 cm (Figure 3
(c)).

Figure 3: 1D scans of NP hits per laser shot for (a) 60 nm SiO2, (b) 120 nm
SiO2, (c) 10 nm Ag, and (d) 120 nm Au taken with a main chamber pressure of
3.5e-4 Torr. Each subfigure presents scans for both the smallest and largest L
value that a scan was performed at. The atomized NP concentration for (a) —
(d) are, respectively: 5.0e13, 6.3e12, 4.0e12, and 2.1e9 NP/mL. CEM signals
were recorded for 5e3, 5e3, 5e4, and 3e5 laser shots for each step of the
motorized mirror in (a) — (d) respectively.

To put a more specific value to the NP beam width for the type of scan
presented in Figure 3 we now define a beam width, w0.8, such that w0.8 is the
range, centered about the center of mass of the NP hit profile, for which the
integral of the profile in that range constitutes 80% of all the NP hits for
that scan. This definition is applied to perform a quantitative analysis of
the NP beam width for 60 nm and 120 nm SiO2 with respect to a variable gas
flow into the aerodynamic lens and is presented in Figure 4 (a). Here the gas
flow is controlled by the pinch bolt (shown in Figure 1) and monitored by the
pressure in the interaction chamber. It is important to note that w0.8 is
unaffected by reasonable choice of threshold.

Figure 4: 60 nm and 120 nm SiO2 NP beams: (a) w0.8, (b) peak density, and (c)
total NP hits as a function of gas flow into the aerodynamics lens as
controlled by the pinch bolt (depicted in Figure 1) and monitored by the
pressure in the interaction chamber. Error bars in (a) represent the step size
of the motorized mirror, whereas those in (b) additionally account for the
uncertainty in the focal spot size diameter.

The NP beam width analysis presented in Figure 4 (a) is now utilized to
approximate a particularly interesting value, the in-vacuum NP volume density
(Figure 4 (b)). By using the thresholding method we have attained the number
of NP hits per laser shot and thus, by knowing the interaction volume, the NP
volume density can be calculated. Once again, note that the limitations of the
NP detection scheme imply that this approximation will be a lower bound to the
actual NP density. In Figure 4 (b) we present the peak NP density as a
function of interaction chamber pressure. Specifically, we estimate the
density for each pressure at the focal position corresponding to the maximum
number of NP hits for the given 1D scan. The interaction region is
approximated as a cylinder with diameter equal to the 1/e2 spot size of the
laser focus (45.8 ± 0.5 $\mu$m) and length equal to the w0.8. Note that we use
w0.8 to be the length of the interaction region as the Rayleigh range of the
laser focus is millimeters long, an order of magnitude larger than the
effective width of the NP beam. Finally, Figure 4 (c) shows the total number
of NP hits for the scan performed at each pressure.
The analysis described above culminates in a detailed picture regarding the
behavior of the aerodynamic lens with respect to the input gas flow, where
each of the three subfigures that comprise Figure 4 provides key information.
Figure 4 (a) depicts a lens with a ‘preferred’ gas flow input, that is, the NP
beam is focused mostly tightly at the gas flow rate corresponding to a
pressure in the interaction chamber between 3e-4 and 4e-4 Torr for both 60 nm
and 120 nm SiO2. Consistent with this observation, Figure 4 (b) then shows
that the in-vacuum NP density is also maximized here. Finally, despite the NP
beam being ‘optimized’ at these pressures, the total number of NPs comprising
the beam is seen in Figure 4 (c) as generally increasing with the gas flow
throughput (as measured by chamber pressure). Though none of these
observations are particularly surprising, Figure 4, as a whole, demonstrates
two important points. First, the methods we have described thus far are self-
consistent, validating the described SFI based technique as a reliable
characterization method for an in-vacuum NP beam. Further, it demonstrates the
technique as a powerful and practical tool when using an aerodynamic lens.
Depending on the application an aerodynamic lens is being used for, a user may
prefer a tight NP beam for a reliable target, whereas another application may
require that the total number of NPs in-vacuum be maximized and a small target
is of little consequence. The techniques presented here allow the user freedom
to explore their relevant parameter space and tailor their NP beam to the
specific application.
Finally, a similar study is now applied to a variable L (Figure 5) which
allows for a 3-dimensional characterization of the NP beam to be achieved. As
seen in Figure 5 (a) w0.8 nearly doubles from L = 5.5 cm to L = 9 cm and then
remains largely constant through the range allowed by the experimental
apparatus. Figure 5 (b) presents the complementary NP beam densities to the
beam widths presented in Figure 5 (a). It is immediately evident that the NP
density decreases the further from the exit of the aerodynamic lens, dropping
by a factor of $\sim$2 for both 60 nm and 120 nm SiO2 from L = 5.5 cm to L =
18 cm. This is the case despite the NP beam width being largely constant with
increasing L. Further, though the atomized NP samples vary by an order of
magnitude for 60 nm and 120 nm SiO2 (5e13 NP/mL and 6.3e12 NP/mL
respectively), the maximum in-vacuum densities differ by less than a factor of
2.

Figure 5: NP beam characterization along its axis of propagation for 60 nm and
120 nm SiO2 atomized at suspension concentrations of 5e13 and 6.3e12 NP/mL
respectively. (a) presents w0.8 for five L values between 5.5 cm and 18 cm and
each 1D scan uses the same parameters as those used for the SiO2 scans
presented in Figure 2. (b) shows the approximate in-vacuum NP volume density
for these same five L values. This density is attained by modelling the
interaction region as a cylinder of diameter equal to the 1/e2 spot size of
the laser focus and length equal to the w0.8 of the NP beam. Error bars are
computed by the same methods as in Figure 4 (a) and (b) respectively.

Before giving our concluding remarks, we should mention that the SFI
characterization technique is largely comparable to LIBD with the key
difference being the detection of electrons (SFI) versus plasma photons
(LIBD). Both detection schemes utilize a pulsed laser to determine the
presence of a NP within the laser focus for each laser shot and are effective
for smaller NP sizes ($\geq$ 10 nm). This gives both techniques access to the
spatial dependence of the NP density with resolution equal to the laser focal
volume and are subject to similar limitations. It is important to highlight
some differences between SFI used here and LIBD as performed by Barreda et
al.[18]. First, we have made use of a fs laser which allows for focal peak
intensities more than 1000 times larger than those produced by the ns laser
used in [18]. This ensures that not only are the NPs in the focal region
ionized, but they are ionized to an extent which allows electrons to be
detected via CEM. Further, this extra peak power permits flexibility for an
optical setup where the focusing lens is outside of the vacuum chamber. In our
opinion, the biggest strength of SFI vs LIBD is that the latter relies on the
emission of plasma photons from the NP which is a material dependent process
and may not always be realizable, whereas ionization is possible with
sufficient laser intensity, something very achievable with modest laser
systems. This big advantage could make SFI a more general prospect compared to
LIBD.

##  4 Conclusion

We have demonstrated the use of strong field ionization for the detection of
isolated NPs in vacuum and the utility of this technique for 3-dimensional
characterization of NP beams produced by an aerodynamic lens. By using a CEM
in conjunction with a fs pulsed laser focused into the NP beam we determine if
there is a NP present in the laser focus for each laser pulse by detecting its
ejected electrons. With this information and with a known laser focal volume
the NP density in the focus is determined. Scanning the laser focus through
the NP beam allows for simple and effective mapping of the 3-dimensional NP
density distribution. As discussed in the text, the thresholding technique
undercounts the number of events in the beam and thus the calculated NP beam
density represents a lower bound to the real NP density.
Finally, we have applied this scheme to optimize the NP beam. The parameter to
optimize in this instance is the aerosol flow into the aerodynamic lens. It is
shown that the aerosol input influences both the collimation (width) of the NP
beam and the number of NPs that comprise the beam, both of which effect the NP
volume density and is crucial to many applications.

##  5 Acknowledgments

The authors would like to recognize the Air Force Office of Scientific
Research (AFOSR) for financial support of this work under grant no.
FA9550-21-1-0387. T.S. was partially funded by the U.S. Department of Energy
(DOE) under grant no. DE-SC0019098.

##  6 Author contributions statement

M.D., T.S., N.G.H., J.A.P., and C.T. conceived the technique. M.D. and J.A.P.
constructed the experimental apparatus. M.D. and T.S. conducted the
experiments. M.D. analysed the results and prepared the manuscript. All
authors contributed in discussions and reviewed the manuscript.

##  7 Availability of data and materials

The datasets used/analyzed in this manuscript are available from the
corresponding author on reasonable request.

##  8 Additional information

The authors declare no competing financial interests.

## References

  * [1] Daniel Bobo, Kye J Robinson, Jiaul Islam, Kristofer J Thurecht, and Simon R Corrie.  Nanoparticle-based medicines: a review of fda-approved materials and clinical trials to date.  Pharmaceutical research, 33(10):2373–2387, 2016.
  * [2] Stephanie D Steichen, Mary Caldorera-Moore, and Nicholas A Peppas.  A review of current nanoparticle and targeting moieties for the delivery of cancer therapeutics.  European journal of pharmaceutical sciences, 48(3):416–427, 2013\.
  * [3] Naresh Kumar Rajendran, Sathish Sundar Dhilip Kumar, Nicolette Nadene Houreld, and Heidi Abrahamse.  A review on nanoparticle based treatment for wound healing.  Journal of Drug Delivery Science and Technology, 44:421–430, 2018\.
  * [4] Jyoti Boken, Parul Khurana, Sheenam Thatai, Dinesh Kumar, and Surendra Prasad.  Plasmonic nanoparticles and their analytical applications: A review.  Applied Spectroscopy Reviews, 52(9):774–820, 2017.
  * [5] Jan Krajczewski, Karol Kołątaj, and Andrzej Kudelski.  Plasmonic nanoparticles in chemical analysis.  RSC advances, 7(28):17559–17576, 2017.
  * [6] Keisuke Nakayama, Katsuaki Tanabe, and Harry A Atwater.  Plasmonic nanoparticle enhanced light absorption in gaas solar cells.  Applied Physics Letters, 93(12):121904, 2008.
  * [7] Yuki Kameya and Katsunori Hanamura.  Enhancement of solar radiation absorption using nanoparticle suspension.  Solar Energy, 85(2):299–307, 2011.
  * [8] Marcello F Ciappina, José A Pérez-Hernández, Alexandra S Landsman, William A Okell, Sergey Zherebtsov, Benjamin Förg, Johannes Schötz, Lennart Seiffert, Thomas Fennel, Tahir Shaaran, et al.  Attosecond physics at the nanoscale.  Reports on Progress in Physics, 80(5):054401, 2017.
  * [9] Yu-xi Zhang and Yu-hua Wang.  Nonlinear optical properties of metal nanoparticles: a review.  RSC advances, 7(71):45129–45144, 2017.
  * [10] Vasyl G Kravets, Andrei V Kabashin, William L Barnes, and Alexander N Grigorenko.  Plasmonic surface lattice resonances: a review of properties and applications.  Chemical reviews, 118(12):5912–5951, 2018.
  * [11] Emilie Ringe, Bhavya Sharma, Anne-Isabelle Henry, Laurence D Marks, and Richard P Van Duyne.  Single nanoparticle plasmonics.  Physical Chemistry Chemical Physics, 15(12):4110–4129, 2013.
  * [12] Peng Liu, Paul J Ziemann, David B Kittelson, and Peter H McMurry.  Generating particle beams of controlled dimensions and divergence: Ii. experimental evaluation of particle motion in aerodynamic lenses and nozzle expansions.  Aerosol Science and Technology, 22(3):314–324, 1995.
  * [13] Peng Liu, Paul J Ziemann, David B Kittelson, and Peter H McMurry.  Generating particle beams of controlled dimensions and divergence: I. theory of particle motion in aerodynamic lenses and nozzle expansions.  Aerosol Science and Technology, 22(3):293–313, 1995.
  * [14] Xuefeng Zhang, Kenneth A Smith, Douglas R Worsnop, Jose L Jimenez, John T Jayne, Charles E Kolb, James Morris, and Paul Davidovits.  Numerical characterization of particle beam collimation: Part ii integrated aerodynamic-lens–nozzle system.  Aerosol Science and Technology, 38(6):619–638, 2004.
  * [15] Xiaoliang Wang, Ashok Gidwani, Steven L Girshick, and Peter H McMurry.  Aerodynamic focusing of nanoparticles: Ii. numerical simulation of particle motion through aerodynamic lenses.  Aerosol science and technology, 39(7):624–636, 2005.
  * [16] Xiaoliang Wang and Peter H McMurry.  An experimental study of nanoparticle focusing with aerodynamic lenses.  International Journal of Mass Spectrometry, 258(1-3):30–36, 2006\.
  * [17] Lena Worbs, Jannik Lübke, Nils Roth, Amit K Samanta, Daniel A Horke, and Jochen Küpper.  Light-sheet imaging for the recording of transverse absolute density distributions of gas-phase particle-beams from nanoparticle injectors.  Optics express, 27(25):36580–36586, 2019.
  * [18] F-A Barreda, Christophe Nicolas, J-B Sirven, F-X Ouf, J-L Lacour, Emmanuel Robert, S Benkoula, J Yon, Catalin Miron, and Olivier Sublemontier.  In-situ characterization of nanoparticle beams focused with an aerodynamic lens by laser-induced breakdown detection.  Scientific reports, 5(1):1–7, 2015.
  * [19] JR Lead and KJ Wilkinson.  Environmental Colloids and Particles: Behaviour, Structure and Characterisation.  John Wiley & Sons Hoboken, NJ, USA, 2006.
  * [20] Jeffrey A Powell, Adam M Summers, Qingcao Liu, Seyyed Javad Robatjazi, Philipp Rupp, Johannes Stierle, Carlos Trallero-Herrero, Matthias F Kling, and Artem Rudenko.  Interplay of pulse duration, peak intensity, and particle size in laser-driven electron emission from silica nanospheres.  Optics express, 27(19):27124–27135, 2019.
  * [21] J Powell, SJ Robatjazi, A Vajdu, V Makhija, Johannes Stierle, X Li, Y Malakar, WL Pearson, C Sorensen, Matthias F Kling, et al.  Intensity-dependent photoelectron spectroscopy of gas-phase nanoparticles without focal volume averaging.  In CLEO: QELS_Fundamental Science, pages FTh4M–3. Optical Society of America, 2016.


The refractive index of the relic gravitons and the nHz band
Massimo Giovannini 111e-mail address<EMAIL_ADDRESS>
Department of Physics, CERN, 1211 Geneva 23, Switzerland
INFN, Section of Milan-Bicocca, 20126 Milan, Italy

Abstract

If the refractive index of the relic gravitons increases during a conventional
stage of inflationary evolution the spectral energy density is blue at
intermediate frequencies above the fHz and then flattens out after a knee that
is typically smaller than the mHz. We investigate here the conditions leading
to a sufficiently large spectral energy density in the nHz range where some
peculiar signatures observed with the pulsar timing arrays have been recently
attributed to cosmic gravitons. If these potential evidences are combined with
the most recent bounds provided by wide-band interferometers in the audio
range (i.e. between few Hz and the kHz) the allowed regions of the parameter
space are compatible with both determinations and also with all the other
constraints associated with the background of relic gravitons produced during
inflation. The present analysis suggests that the pulsar timing arrays are
sensitive to the evolution of the refractive index during early stages of the
inflationary evolution. This physical interpretation of the preliminary
empirical evidence is distinguishable from other perspectives since the high-
frequency normalization, the blue spectral index and the tensor to scalar
ratio cannot be independently assigned but are all related to the frequency of
the knee that is ultimately determined by the competition between the rate of
evolution of the refractive index and the slow-roll corrections.

The pumping action of the space-time curvature is the primary source of relic
gravitons [2, 3] but when a conventional stage of inflationary expansion is
replaced by a radiation-dominated epoch, the spectral energy density in
critical units at the present conformal time $\tau_{0}$ (denoted hereunder by
$\Omega_{gw}(\nu,\tau_{0})$) is quasi-flat for comoving frequencies $\nu$
larger than $100$ aHz [4]. The transition across the epoch of matter-radiation
equality leads to an infrared branch where
$\Omega_{gw}(\nu,\tau_{0})\propto\nu^{-2}$ between the aHz and $100$ aHz [5].
The conventional inflationary evolution sems to exclude the possibility of a
drastic deviation from scale-invariance in the direction of blue spectral
indices for typical frequencies larger than fHz where, in a conservative
perspective, $h_{0}^{2}\Omega_{gw}(\nu,\tau_{0})={\mathcal{O}}(10^{-16.5})$.
On the contrary the corrections typically go into the direction of a
decreasing spectral slope as it happens, for instance, in the single-field
case where, thanks to the consistency relations, the tensor spectral index
$n_{T}$ is notoriously related to the tensor to scalar ratio $r_{T}$ as
$n_{T}\simeq-r_{T}/8$. Since $r_{T}$ is currently assessed from the analysis
of the temperature and polarization anisotropies of the Cosmic Microwave
Background (CMB) it is difficult to argue that $n_{T}$ could be positive. For
typical frequencies larger than the mHz the spectral energy density increases
provided the post-inflationary evolution is stiffer than radiation. For
frequencies between the fHz and the mHz the spectral index can be blue in
generalized gravity theories and when the dominant energy condition is
violated (see, for instance, [6] for a recent review focussed on these
issues).

Gravitational waves also acquire an effective index of refraction when they
travel in curved space-times [7]. In this context the blue spectral slopes
arise from a more mundane effect associated with the variation of the
refractive index even if the background geometry evolves according to a
conventional stage of expansion possibly supplemented by a standard
decelerated epoch [8]. The general idea behind this suggestion is that the
effective action of single-field inflationary models involves all the
different terms that include four derivatives and are suppressed by the
negative powers of a large mass scale [9]. If parity breaking terms are
included in the effective action [10, 11], the relic graviton background may
be polarized but this possibility has been already discussed elsewhere [12].
There are however non-generic models of inflation where the higher-order
corrections assume a particular form since the inflaton has some particular
symmetry (like a shift symmetry $\varphi\to\varphi+c$) or because the rate of
inflaton roll remains constant (and possibly larger than $1$), as it happens
in certain fast-roll scenarios [13] (see also, for instance, [14, 15]). Other
examples include the situation where the higher-order curvature corrections
are given in terms of the Gauss-Bonnet combination weighted by some inflaton
dependent-coupling [16, 17, 18]. In [8] (see also [20, 19]) it has been argued
that in all these situations the effective action of the relic gravitons is be
modified and ultimately assumes the following general form:

$S_{g}=\frac{1}{8\ell_{P}^{2}}\int
d^{4}x\biggl{[}A(\tau)\,\,\partial_{\tau}h_{ij}\,\partial_{\tau}h^{ij}-B(\tau)\,\partial_{k}h_{ij}\partial^{k}h^{ij}\biggr{]},$
(1)

where we assume that the background is conformally flat and
$\partial_{i}h_{j}^{\,\,\,\,i}=h_{i}^{\,\,\,\,i}=0$. While both terms
$A(\tau)$ and $B(\tau)$ depend on the conformal time coordinate $\tau$ we can
always factor $A(\tau)$ and parametrize and introduce an effective refractive
index $n(\tau)$ associated with the interactions with the background geometry,
as originally argued in Ref. [7, 8]):

$S_{g}=\frac{1}{8\ell_{P}^{2}}\int
d^{4}x\,\,a^{2}(\tau)\,\,\biggl{[}\partial_{\tau}h_{ij}\,\partial_{\tau}h^{ij}-\frac{\partial_{k}h_{ij}\partial^{k}h^{ij}}{n^{2}(\tau)}\biggr{]},$
(2)

where $a(\tau)$ is the scale factor in the conformal time coordinate; we are
assuming here a conformally flat background metric
$\overline{g}_{\mu\nu}=a^{2}(\tau)\,\eta_{\mu\nu}$. After Eq. (2) has been
proposed in Ref. [8] apparently different parametrizations have been later
introduced (see [19, 20] and discussion theierin). The difference between
these strategies consists in modifying the first term (rather than the second)
inside the squared bracket of Eq. (2). This choice is immaterial since the two
parametrizations of the effect are related by a rescaling of the four-
dimensional metric through a conformal factor that involves the refractive
index itself. Equation (2) simplifies by introducing a new time coordinate
conventionally referred to as the $\eta$-time:

$S_{g}=\frac{1}{8\ell_{P}^{2}}\int\,\,d^{4}x\,\,b^{2}(\eta)\,\,\biggl{[}\partial_{\eta}h_{ij}\,\,\partial_{\eta}h^{ij}-\partial_{k}h_{ij}\,\,\partial^{k}h^{ij}\biggr{]},\qquad\qquad
b(\eta)=\frac{a(\eta)}{\sqrt{n(\eta)}},\qquad\ell_{P}=\sqrt{8\pi G}.$ (3)

Since the $\eta$-time is defined by $n(\eta)\,d\eta\,=\,d\tau$ and the scale
factor is everywhere continuous with its first derivative222Note that the
continuity of the extrinsic curvature implies the continuity of the first
(conformal) time derivative of the scale factor., the evolution of the
refractive index is specified unambiguously by assigning $n(a)$. Equation (3)
generalizes the standard Ford-Parker action [21] to the case of a dynamical
refractive index. Even if the phase velocity of the relic gravitons is not
required to be sub-luminal we shall consider the situation where $n(a)\geq 1$.
In what follows for the sake of illustration we shall be assuming that $n(a)$
changes appreciably during inflation and it goes to $1$ during the standard
decelerated stage of expansion:

$n(a)=n_{\ast}\frac{(a/a_{\ast})^{\alpha}\,\,e^{-\xi(a/a_{1})}}{(a/a_{*})^{\alpha}+1}+1,\qquad\qquad
n_{\ast}=n_{i}(a_{\ast}/a_{i})^{\alpha}=n_{i}e^{\alpha\,N_{\ast}},$ (4)

where $a_{i}$ and $a_{1}$ mark, respectively, the beginning and the end of the
inflationary epoch; $a_{*}$ defines the boundary of the refractive stage. What
matters is not the specific analytical form (that only affects the transition
regions) but the three successive physical regimes described by Eq. (4). In
particular we have that $n(a)\to 1$ when $a\gg a_{1}$ and the sharpness of
this transition is controlled by $\xi\geq 1$. When $a_{*}<a<a_{1}$ $n(a)$ is
constant (but still larger than $1$), i.e. $n(a)\simeq n_{\ast}>1$. Finally
for $a<a_{\ast}$ we have the truly refractive stage where $n(a)\simeq
n_{\ast}(a/a_{\ast})^{\alpha}$. Even if in Ref. [8] more general possibilities
have been examined (including a post-inflationary evolution of $n(a)$), for
the present ends it is enough to consider the minimal case illustrated by Eq.
(4).

The different dynamical stages of the model are illustrated in the cartoon of
Fig. 1 where the common logarithm of the effective horizon $\dot{b}/b$ is
reported. In contrast with the standard notations, the overdot denotes here a
derivation with respect to the $\eta$ coordinate (and not a derivation with
respect to the cosmic time as usually implied). From the shape of Fig. 1 we
can deduce that the spectral energy density will have three different
branches. The modes $k<k_{eq}$ correspond to wavelengths leaving the Hubble
radius for $b<b_{*}$ and reentering during the matter-dominated stage. We
remind that, as usual,
$k_{eq}=0.0732\,h_{0}^{2}\Omega_{M0}\,\mathrm{Mpc}^{-1}$ where $\Omega_{M0}$
is the present fraction in dusty matter333The comoving frequency associated to
$k_{eq}$ is given by $\nu_{eq}=k_{eq}/(2\pi)=1.597\times
10^{-17}\,\mathrm{Hz}$ for $h_{0}^{2}\Omega_{M0}=0.1411$ and
$h_{0}^{2}\Omega_{R0}=4.15\times 10^{-5}$..

Figure 1: The common logarithm effective horizon is approximately illustrated
in terms of the scale factor. The different wavenumbers hit twice the
effective horizon and the two crossings determine the slope of the spectrum in
a given dynamical stage.

The modes $k_{eq}<k<k_{\ast}$ (where $k_{\ast}=1/\eta_{\ast}$) correspond to
the wavelengths leaving the Hubble radius during the refractive stage and
reentering when the background is dominated by radiation. In this range of
wavenumbers the spectral energy density is either blue or violet, depending on
the value of $\alpha$ which is always positive as long as the refractive index
is larger than $1$. Finally in the region $k_{\ast}<k<k_{1}$ the spectrum is
quasi-flat since the corresponding wavelengths cross the Hubble radius during
the inflationary phase (when the refractive index is not dynamical anymore)
and reenter in the radiation stage of expansion. The comoving frequencies
$\nu_{max}=k_{max}/(2\,\pi)$ and $\nu_{\ast}=k_{*}/(2\,\pi)$ cannot be
assigned at wish depending on the convenience of the potential signal but they
are determined from the various stages of the model. In particular the maximal
frequency is typically ${\mathcal{O}}(200)$ MHz:

$\nu_{max}=0.26\,\biggl{(}\frac{\epsilon}{0.003}\biggr{)}^{1/4}\biggl{(}\frac{{\mathcal{A}}_{\mathcal{R}}}{2.41\times
10^{-9}}\biggr{)}^{1/4}\,\biggl{(}\frac{h_{0}^{2}\Omega_{R0}}{4.15\times
10^{-5}}\biggr{)}^{1/4}\,\mathrm{GHz},$ (5)

where $\epsilon$ denotes the slow-roll parameter and
${\mathcal{A}}_{{\mathcal{R}}}$ is the amplitude of the scalar power spectrum;
according to the current determinations we have that the tensor to scalar
ratio should be smaller than $0.064$ [23, 24]; in a conservative perspective
we shall be requiring $\epsilon\leq 0.06/16\simeq 0.003$ as assumed in Eq.
(5). We remind that the typical value of the pivot scale used in this
discussion is $k_{p}=0.002\,\mathrm{Mpc}^{-1}$ corresponding to a typical
comoving frequency $\nu_{p}=3.092\,\,\mathrm{aHz}$. Given $\nu_{max}$ the
value of $\nu_{\ast}$ defines the knee of the spectrum:

$\nu_{\ast}=\biggl{(}1+\frac{\alpha}{1-\epsilon}\biggr{)}e^{N_{\ast}(\alpha+1)-N_{t}}\,\,\nu_{max},$
(6)

where $N_{t}$ denotes the total number of $e$-folds while, as in Eq. (4),
$N_{\ast}=\ln{(a_{\ast}/a_{i})}$. In Eq. (6) as well as in the forthcoming
discussion we shall always be assuming that $n_{i}\to 1$; different choices
are possible (provided $n_{i}\geq 1$) but their effect does not modify the
conclusions since the value of $n_{i}$ can always be traded for a shorter
refractive phase. The relevant point to appreciate is that $\nu_{\ast}$ will
control the typical frequency of the knee of the spectrum and it does depend
on $N_{*}$, $N_{t}$ and $\alpha$.

The Hamiltonian associated with Eq. (3) is directly expressible in the
$\eta$-time and it is:

$H_{g}(\eta)=\int
d^{3}x\biggl{[}\frac{8\ell_{P}^{2}}{b^{2}}\pi_{i\,j}\,\pi^{\,i\,j}+\frac{b^{2}}{8\ell_{P}^{2}}\partial_{k}h_{i\,j}\,\partial^{k}h^{\,i\,j}\biggr{]},\qquad\qquad\pi_{i\,j}=\frac{b^{2}}{8\ell_{P}^{2}}\partial_{\eta}h_{i\,j}.$
(7)

After promoting the field operators and their conjugate momenta to the status
of quantum operators, from Eq. (7) the governing equations for
$\widehat{h}_{i\,j}$ and $\widehat{\pi}_{i\,j}$ are:

$\partial_{\eta}\widehat{h}_{i\,j}=\frac{8\ell_{P}^{2}}{b^{2}}\,\widehat{\pi}_{i\,j},\qquad\qquad\partial_{\eta}\widehat{\pi}_{i\,j}=\frac{b^{2}}{8\ell_{P}^{2}}\,\nabla^{2}\widehat{h}_{i\,j}.$
(8)

In Fourier space Eq. (8) implies a specific solution of the corresponding mode
functions $F_{k}(\eta)$ and $G_{k}(\eta)$; more specifically we have

$\displaystyle\widehat{h}_{ij}(\vec{x},\eta)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}\,\ell_{P}}{(2\pi)^{3/2}}\sum_{\lambda}\int\,d^{3}k\,\,e^{(\lambda)}_{ij}(\vec{k})\,\biggl{[}F_{k,\lambda}(\eta)\,\widehat{a}_{\vec{k},\,\lambda}e^{-i\vec{k}\cdot\vec{x}}+F^{\ast}_{k,\lambda}(\eta)\,\widehat{a}_{\vec{k},\,\lambda}^{\dagger}e^{i\vec{k}\cdot\vec{x}}\biggr{]},$
(9) $\displaystyle\widehat{\pi}_{ij}(\vec{x},\eta)$ $\displaystyle=$
$\displaystyle\frac{1}{4\,\sqrt{2}\,\ell_{P}\,(2\pi)^{3/2}}\sum_{\lambda}\int\,d^{3}k\,\,e^{(\lambda)}_{ij}(\vec{k})\,\biggl{[}G_{k,\lambda}(\eta)\,\widehat{a}_{\vec{k},\,\lambda}\,e^{-i\vec{k}\cdot\vec{x}}+G^{\ast}_{k,\lambda}(\eta)\,\widehat{a}_{\vec{k},\,\lambda}^{\dagger}\,e^{i\vec{k}\cdot\vec{x}}\biggr{]}.$
(10)

In Eqs. (9)–(10) $\lambda=\oplus,\,\otimes$ runs over the tensor polarizations
but, as in the conventional situation, the evolution of the mode functions is
the same for each of the two values of $\lambda$:

$\dot{G}_{k,\,\lambda}=-k^{2}\,b^{2}\,F_{k,\,\lambda},\qquad\dot{F}_{k,\,\lambda}=\frac{G_{k,\,\lambda}}{b^{2}},$
(11)

where we recall, as already mentioned in connection with Fig. 1 that the
overdot now denotes a derivation with respect to $\eta$ (i.e.
$\dot{F}_{k,\,\lambda}=\partial_{\eta}\,F_{k,\,\lambda}$). From Eqs. (9)–(10)
the canonical commutation relations at equal $\eta$-times become

$[\widehat{h}_{i\,j}(\vec{x},\eta),\,\widehat{\pi}_{m\,n}(\vec{y},\eta)]=i\,{\mathcal{S}}_{i\,j\,m\,n}(\hat{k})\,\delta^{(3)}(\vec{x}-\vec{y}),$
(12)

In Eq. (12) ${\mathcal{S}}_{i\,j\,m\,n}(\hat{k})$ follows from the sum over
the tensor polarizations and the Wronskian of the mode functions is fixed:

$\displaystyle{\mathcal{S}}_{i\,j\,m\,n}(\hat{k})=\frac{1}{4}\biggl{[}p_{i\,m}(\hat{k})\,p_{j\,n}(\hat{k})+p_{j\,m}(\hat{k})\,p_{i\,n}(\hat{k})-p_{i\,j}(\hat{k})\,p_{m\,n}(\hat{k})\biggr{]},$
(13) $\displaystyle
F_{k,\,\lambda}(\eta)\,G_{k,\,\lambda}^{\ast}(\eta)-F_{k,\,\lambda}^{\ast}(\eta)\,G_{k,\,\lambda}(\eta)\,=\,i/b^{2}(\eta),$
(14)

where $p_{i\,j}(\hat{k})=(\delta_{i\,j}-\hat{k}_{i}\,\hat{k}_{j})$ denotes the
standard projector and $\hat{k}^{i}=k^{i}/|\vec{k}|$. Note that if the
normalization of the Wronskian differs from Eq. (14) the commutation relation
of Eq. (12) is not canonical anymore.

The mode functions are normalized during the refractive phase where, according
to Eq. (4), the $\eta$-time and the conformal time are related as
$(-\eta/\eta_{\ast})=(-\tau/\tau_{\ast})^{\alpha/(1-\epsilon)+1}$ and
$\eta_{\ast}=\tau_{\ast}(1-\epsilon)/[n_{\ast}(1+\alpha-\epsilon)]$; both
relations follow from the definition of the $\eta$-time (i.e.
$n(\eta)d\eta=d\tau$) and also from $n(a)=n_{*}(a/a_{\ast})^{\alpha}$, as
implied by Eq. (4) when $a<a_{\ast}$. In the $\eta$-time the explicit
expression of $b(\eta)$ is:

$b(\eta)=b_{\ast}\biggl{(}-\frac{\eta}{\eta_{\ast}}\biggr{)}^{-\sigma},\qquad\sigma=\frac{2-\alpha}{2(1+\alpha-\epsilon)},\qquad\mathrm{for}\qquad\eta<-\eta_{\ast},$
(15)

where $b_{\ast}=a_{\ast}/\sqrt{n_{\ast}}$. In the refractive regime the
solution of Eq. (11) is therefore given by:

$F_{k}(\eta)=\frac{{\mathcal{N}}}{\sqrt{2k}\,b(\eta)}\,\sqrt{-k\eta}\,H_{\mu}^{(1)}(-k\,\eta),\qquad\qquad
G_{k}(\eta)=-\frac{{\mathcal{N}}}{b(\eta)}\,\sqrt{\frac{k}{2}}\,\sqrt{-k\eta}\,H_{\mu-1}^{(1)}(-k\eta),$
(16)

where $\mu=\sigma+1/2$ and $H_{\mu}^{(1)}(-k\eta)$ is the Hankel function of
first kind [25]; note that $|{\mathcal{N}}|=\sqrt{\pi/2}$ and, with this
choice, the Wronskian normalization condition of Eq. (14) is satisfied. Thanks
to Eq. (16) we can compute the tensor power spectrum that eventually
determines the low-frequency normalization of the spectral energy density:

$\displaystyle{\mathcal{P}}_{T}(k,\eta)$ $\displaystyle=$
$\displaystyle\frac{4\ell_{P}^{2}\,k^{3}}{\pi^{2}}\bigl{|}F_{k}(\eta)\bigr{|}^{2}\to\biggl{(}\frac{H_{1}}{M_{P}}\biggr{)}^{2}\,\,{\mathcal{C}}(n_{T},N_{*},N_{t},\epsilon)\,\,\biggl{(}\frac{k}{a_{1}H_{1}}\biggr{)}^{n_{T}},$
(17) $\displaystyle\qquad{\mathcal{C}}(n_{T},N_{*},N_{t},\epsilon)$
$\displaystyle=$
$\displaystyle\frac{2^{6-n_{T}}}{\pi^{2}}\biggl{|}1+\frac{\alpha}{1-\epsilon}\biggr{|}^{2-n_{T}}\,\Gamma^{2}\biggl{(}\frac{3-n_{T}}{2}\biggr{)}e^{\alpha\,N_{\ast}(3-n_{T})-n_{T}(N_{\ast}-N_{t})},$
(18)

The second expression at the right hand side of Eq. (17) corresponds to the
limit $|k\eta|\ll 1$ where ${\mathcal{P}}_{T}(k,\eta)$ becomes constant in
time and the tensor to scalar ratio is:

$r_{T}(\nu)=\frac{\epsilon}{\pi}\,{\mathcal{C}}(n_{T},N_{*},N_{t},\epsilon)\biggl{(}\frac{\nu}{\nu_{max}}\biggr{)}^{n_{T}},\qquad\qquad\lim_{\alpha\to
0}{\mathcal{C}}(n_{T},N_{*},N_{t},\epsilon)=16\pi.$ (19)

For $\alpha\to 0$ we have that $r_{T}\to 16\epsilon$ and the standard
consistency relation is recovered. In the same limit the tensor spectral index
goes to $-2\epsilon$:

$n_{T}=\frac{3\alpha-2\epsilon}{1+\alpha-\epsilon}=\frac{3\alpha}{1+\alpha}+\frac{\epsilon(\alpha-2)}{(1+\alpha)^{2}}+{\mathcal{O}}(\epsilon^{2}).$
(20)

Equation (20) defines the tensor spectral index in the intermediate frequency
range where $\alpha$ is larger than $\epsilon$ and this is why the exact
result can always be expanded in the limit $\epsilon\ll 1$. The slope of the
tensor power spectrum evaluated for the wavelengths larger than the Hubble
radius coincides, after reentry, with the slope of spectral energy density
before the knee:

$h_{0}^{2}\,\,\Omega_{\mathrm{gw}}(\nu,\tau_{0})={\mathcal{T}}^{2}(\nu,\nu_{eq})\,\,r_{\mathrm{T}}(\nu_{p})\,\,\biggl{(}\frac{\nu}{\nu_{p}}\biggr{)}^{n_{\mathrm{T}}},\qquad\nu<\nu_{\ast},$
(21)

where $\nu_{p}=3.092\,\,\mathrm{aHz}$ and ${\mathcal{T}}(\nu,\nu_{eq})$ is:

${\mathcal{T}}(\nu,\nu_{eq})=6.453\times
10^{-8}\sqrt{1+c_{2}\biggl{(}\frac{\nu_{eq}}{\nu}\biggr{)}+b_{2}\biggl{(}\frac{\nu_{eq}}{\nu}\biggr{)}^{2}},\qquad
c_{2}=0.5238,\qquad b_{2}=0.3537,$ (22)

where the constants $c_{2}$ and $b_{2}$ follow from the transfer function for
the spectral energy density [22] and are obtained by integrating the mode
functions across the radiation-matter transition and by computing
$\Omega_{\mathrm{gw}}(\nu,\tau)$ for different frequencies. Above $\nu_{\ast}$
the spectral energy density is instead quasi-flat, as anticipated in Fig. 1
and it is approximately given by:

$h_{0}^{2}\Omega_{\mathrm{gw}}(\nu,\tau_{0})={\mathcal{T}}^{2}(\nu,\nu_{eq})\,\,r_{\mathrm{T}}(\nu_{p})\,\,\biggl{(}\frac{\nu_{\ast}}{\nu_{p}}\biggr{)}^{n_{\mathrm{T}}}\,\,\biggl{(}\frac{\nu}{\nu_{\ast}}\biggr{)}^{-2\epsilon}\,\qquad\nu_{\ast}<\nu<\nu_{max},$
(23)

where the contribution of ${\mathcal{T}}^{2}(\nu,\nu_{eq})$ is, in practice,
frequency-independent for $\nu>\nu_{\ast}$ and it can be approximated as
$4.165\times 10^{-15}$. The results of Eqs. (21)–(22) and (23) agree with the
previous results of Ref. [8] and are a natural candidate for exploring a
potential signal in the nHz range: the blue spectrum in the intermediate range
and the flat slope at high-frequencies have a clear dynamical rationale which
is encoded in Fig. 1. It is however clear that the spectral energy density
obtained in this concrete scenario does not support the viewpoint that
$n_{T}$, $r_{T}$ and $\nu_{\ast}$ could be independently assigned without
worrying about their mutual interplay.

There are two classes of recent constraints that are relevant for the present
discussions. The $12.5$ yrs data set of the NANOGrav pulsar timing array (PTA)
experiment reported evidence of a stochastic signal in the nHz range. The
features of this purported signal would imply, in the present notations,
that444See also Ref. [27] for the $11$ yrs and Ref. [28] for the $9$ yrs data
set where some useful upper limits on the relic gravitons have been
reported.[26]:

$10^{-8.6}<\,h_{0}^{2}\Omega_{gw}(\nu,\tau_{0})<\,10^{-9.8},\qquad\qquad
3\,\mathrm{nHz}<\nu<100\,\mathrm{nHz}.$ (24)

The Kagra, LIGO and Virgo collaborations, in their attempt to constraint the
stochastic backgrounds of relic gravitons also reported a constraint [29]
implying, in the case of a flat spectral energy density,

$\Omega_{gw}(\nu,\tau_{0})<5.8\times 10^{-9},\qquad\qquad
20\,\,\mathrm{Hz}<\nu_{KLV}<76.6\,\,\mathrm{Hz},$ (25)

where $\nu_{KLV}$ denotes the Kagra-LIGO-Virgo frequency. This limit improves
on a series of bounds previously deduced by the wide-band interferometers (see
[6] for a recent review). In particular in Ref. [30] the analog of Eq. (25)
implied $\Omega_{gw}(\nu,\tau_{0})<6\times 10^{-8}$ for a comparable frequency
interval and always in the case of a flat spectral energy density. While the
bound of Eq. (25) could be immediately used also in our case since at high-
frequency the spectral energy density is nearly scale-invariant, it is useful
to elaborate on the result of Ref. [29]. The collaboration actually reports a
threefold bound which could be parametrized as

$\Omega_{gw}(\nu,\tau_{0})=\overline{\Omega}(\delta)\biggl{(}\frac{\nu}{\nu_{ref}}\biggr{)}^{\delta},\qquad\qquad\nu_{ref}=25\,\,\mathrm{Hz}.$
(26)

In terms of Eq. (26) the results of Ref. [29] read
$\overline{\Omega}(0)<5.8\times 10^{-9}$ (valid in the case $\delta=0$),
$\overline{\Omega}(2/3)<3.4\times 10^{-9}$ (when $\delta=2/3$) and
$\overline{\Omega}(3)<3.9\times 10^{-10}$ (when $\delta=3$). When the value of
$\delta$ increases the bound becomes more restrictive once the reference
frequency has been kept fixed. The three results are unified in the following
interpolating formula

$\log{\overline{\Omega}}(\delta)<-\,8.236-\,0.335\,\delta-0.018\,\delta^{2}.$
(27)

While the quadratic fit is slightly more accurate the linear one, the essence
of the arguments does not change in the two cases since the different points
fall approximately on the same straight line (i.e. $-8.223-0.393\,\delta$).
Since in the present case the bound (27) should be applied at high-frequencies
we will have $\delta=-2\epsilon/(1-\epsilon)$ with $\epsilon<0.01$. We can
immediately see from Eq. (27) that, to leading order in $\epsilon$, Eq. (27)
implies that
$\log{\overline{\Omega}}(\epsilon)<-\,8.236-\,0.335\,\epsilon-0.393\epsilon^{2}$
giving, at most, a negligible correction in the exponent. The pivotal
parameters that determine the spectrum are $\alpha$, $N_{\ast}$ and $N_{t}$.
If $N_{*}$ is of the order of $N_{t}$ the transition to normalcy occurs at the
end of inflation555 It is not difficult to show that in this case, as
previously discussed, it is impossible to get a large signal in the nHz range
without jeopardizing the big-bang nuclosynthesis constraint (see [8, 19]).
This point is quite well known also from other scenarios where the blue
spectra are not induced by the variation of the refractive index [6]. It is
clear that this conclusion holds, in the present context, because $\nu_{max}$
is cannot be arbitrarily fixed. When $N_{*}<N_{t}$ the transition to normalcy
takes place well before the onset of the radiation-dominated epoch (i.e. when
the background is still inflating deep inside the quasi-de Sitter stage of
expansion). The scale-dependent constraint on the tensor to scalar ratio is
conventionally applied at the pivot scale $k_{p}$ (see the discussion after
Eq. (5)) corresponding to the comoving frequency $\nu_{p}={\mathcal{O}}(3)$
aHz. In a conservative perspective we therefore require the limit [23, 24]
$r_{T}(\nu_{p})\leq 0.06$. The bounds coming from big-bang nucleosynthesis
[31, 32, 33] imply instead a constraint on the integral of the spectral energy
density of the gravitons:

$h_{0}^{2}\int_{\nu_{bbn}}^{\nu_{max}}\Omega_{\mathrm{gw}}(\nu,\tau_{0})d\ln{\nu}=5.61\times
10^{-6}\Delta N_{\nu}\biggl{(}\frac{h_{0}^{2}\Omega_{\gamma 0}}{2.47\times
10^{-5}}\biggr{)},$ (28)

where $\Omega_{\gamma 0}$ is the (present) critical fraction of CMB photons.
The lower limit of integration in Eq. (28) is given by the frequency
corresponding to the Hubble rate at the nucleosynthesis epoch:

$\nu_{bbn}=2.252\times
10^{-11}\biggl{(}\frac{N_{eff}}{10.75}\biggr{)}^{1/4}\biggl{(}\frac{T_{\mathrm{bbn}}}{\,\,\mathrm{MeV}}\biggr{)}\biggl{(}\frac{h_{0}^{2}\Omega_{\mathrm{R}0}}{4.15\times
10^{-5}}\biggr{)}^{1/4}\,\,\mathrm{Hz}\simeq 0.01\,\mathrm{nHz},$ (29)

where $N_{eff}$ denotes the effective number of relativistic degrees of
freedom entering the total energy density of the plasma and $T_{\mathrm{bbn}}$
is the temperature of big-bang nucleosynthesis. The limit of Eq. (28) sets an
indirect constraint on the extra-relativistic species (and, among others, on
the relic gravitons). Since Eq. (28) is relevant in the context of neutrino
physics (when applied to massless fermionic species), the limit is often
expressed for practical reasons in terms of $\Delta N_{\nu}$ representing the
contribution of supplementary neutrino species. The actual bounds on $\Delta
N_{\nu}$ range from $\Delta N_{\nu}\leq 0.2$ to $\Delta N_{\nu}\leq 1$; the
integrated spectral density in Eq. (28) is thus between $10^{-6}$ and
$10^{-5}$.

Figure 2: The allowed region of the parameter space is illustrated together
with the explicit form of the spectral energy density when the corresponding
parameters are selected within the shaded region where all the relevant are
satisfied. The different contours appearing in the plot at the left correspond
to the values of $h_{0}^{2}\,\Omega_{gw}(\nu_{KLV},\tau_{0})$ where
$\nu_{KLV}$ approximately denotes most sensitive frequency domain of the wide-
band interferometers (i.e. Kara-LIGO-Virgo) that recently provided interesting
results that can be interpreted as upper limits on the spectral energy density
of the relic gravitons. The shaded box highlights the region suggested by the
PTA. We finally recall that $\log{}$ denotes the common logarithm (as opposed
to the natural logarithm) employed, for instance, in the definition of the
number of $e$-folds (see Eq. (6) and comments thereafter).

The limits stemming from CMB analyses (i.e. $r_{T}(\nu_{p})<0.06$), the
requirement of Eq. (24), the Kagra-LIGO-Virgo limit (KLV) of Eq. (25) (see
also Eqs. (26)–(27)) and the BBN constraints of Eq. (28) have all been
enforced in Fig. 2. The shaded area (appearing in the plot at the left of Fig.
2 and resembling a circular segment) illustrates the region where are all the
various constraints are concurrently satisfied when $N_{t}$ is fixed to a
fiducial value. On top of the mentioned bounds, for consistency it is
preferable to require, in a model-independent perspective, that
$\nu_{\ast}>\nu_{bbn}$: if this condition is not included the allowed region
in Fig. 2 gets even larger. The shaded box appearing in Fig. 2 illustrates the
region defined by Eq. (24) where the tentative signal of the PTA should be
approximately located. In the right plot of Fig. 2 the spectral energy density
is reported for different values of $\alpha$ and $N_{\ast}$ all falling within
the shaded region that appears in the left plot. When the parameters are
selected within the circular segment the corresponding spectral energy density
falls within the PTA box while the limits coming from wide-band
interferometers are satisfied. The results of Fig. 2 suggest that a large
signal in the nHz range is obtained when the variation of the refractive index
occurs sufficiently early during the inflationary stage and anyway not beyond
the first $20$ $e$-folds. In units of the inflationary Hubble rate, the rate
of variation of the refractive index must fall in the range
$0.25<\alpha<0.55$. For this portion of the parameter space the frequency of
the knee is in the nHz range.

It is now interesting to fix $\alpha$ to some fiducial value and to deduce the
constraints in the plane $(N_{\ast},\,N_{t})$: this strategy is illustrated in
Fig. 3. It is clear that for $\alpha\simeq 2/7$ the allowed region and the
maximal signal of the model (see the right plot of Fig. 3) occur for
$N_{t}={\mathcal{O}}(60)$ and $N_{\ast}={\mathcal{O}}(20)$. Furthermore for
$\alpha=2/7\simeq{\mathcal{O}}(0.28)$, $N_{t}={\mathcal{O}}(60)$ and
$N_{\ast}={\mathcal{O}}(20)$ the spectral energy density passes through the
PTA box and undershoots the $KLV$ bound. The rationale for the apparently
peculiar choice of $\alpha$ appearing in Fig. 3 stems, in short, from the
following consideration. The PTA results are often reported in terms of a
chirp amplitude scaling as $\nu^{-2/3}$ for a typical reference frequency
${\mathcal{O}}(\mathrm{yr}^{-1})$. If, at the present time,
$h_{c}(\nu,\tau_{0})$ scales as $\nu^{\beta}$ the spectral energy density has
the following dependence on the comoving frequency666See Ref. [6] for a
discussion of the chirp amplitude and of its relation with the spectral energy
density. In general terms the chirp amplitude is simply related to the tensor
power spectrum as $2h_{c}^{2}(k,\eta)={\mathcal{P}}_{T}(k,\eta)$. :

$h_{c}(\nu,\tau_{0})=Q\biggl{(}\frac{\nu}{\nu_{ref}}\biggr{)}^{\beta},\qquad\Rightarrow\qquad\Omega_{gw}(\nu,\tau_{0})=\frac{2\,\pi^{2}\,\nu^{2}}{3H_{0}^{2}}Q^{2}\biggl{(}\frac{\nu}{\nu_{ref}}\biggr{)}^{2\alpha}\propto\biggl{(}\frac{\nu}{\nu_{ref}}\biggr{)}^{2+2\beta},$
(30)

where, as already mentioned, $\nu_{ref}={\mathcal{O}}(\mathrm{yr}^{-1})$ and
$Q$ just denotes a constant amplitude. If $\beta=-2/3$ we have that
$\Omega_{gw}\propto\nu^{2/3}$. The value of $\alpha$ corresponding to
$\beta=-2/3$ is obtained by setting $n_{T}\simeq 2/3$ in Eq. (20).
Consequently we have that $\alpha=(2+4\epsilon)/7$ which can be approximated
as $\alpha=2/7+{\mathcal{O}}(\epsilon)$ since $\epsilon<10^{-3}$. In this case
the spectral energy density is illustrated in the right plot of Fig. 3 for a
set of values selected within the allowed region of the parameter space. If we
would really pretend that the PTA measures a real signal with chirp amplitude
proportional to $\nu^{-2/3}$ the region of the parameter space illustrated in
Fig. 3 should be the most promising one.

Figure 3: In the plot at the left we compute the allowed region in the
parameter space when $\alpha$ is fixed to $2/7$. This choice implies that the
chirp amplitude scales as $\nu^{-2/3}$ while the spectral energy density goes
as $\nu^{2/3}$. In the plot at the right we illustrate the spectral energy
density for a set of parameters selected within the shaded region appearing in
the exclusion plot at the right. As in Fig. 2 common logarithms are employed
on both axes of the right plot.

If the results of the PTA are interpreted as an upper limit, the allowed
region of the parameter space gets wider as we can argue by comparing the left
plot of Fig. 4 with the corresponding result reported in Fig. 2, always in the
plot at the left. The logic of the upper limit has been encouraged by the
observational collaborations themselves in the 11-yr release [27] where the
PTA data implied, within our notations, that
$h_{0}^{2}\Omega_{gw}(\nu_{ref},\tau_{0})<3.4(1)\times 10^{-10}$. If we apply
this viewpoint to the present case (and consequently avoid the requirements
associated with Eq. (24)) we obtain the exclusion plot reported at the left of
Fig. 4 where the shaded area corresponds, as usual, to the allowed region.
From Fig. 4 the area where $\alpha<0.25$ is now viable and it corresponds to a
value of $N_{\ast}$ that is comparatively larger than in Fig. 2. It is however
clear that, in this region, the purported “signal” attributed to relic
gravitons cannot be reproduced. Indeed in the right plot of Fig. 4 the
spectral energy density of the relic gravitons always undershoots the PTA box
and the upper limits of wide-band detectors.

All in all the evolution of refractive index of the gravitons during the early
stages of the inflationary evolution leads to a blue (i.e. slightly
increasing) slope of the spectral energy density at intermediate frequencies
above the fHz. While the specific values of the slopes are determined by the
competition of the slow-roll parameter and of the rate of variation of the
refractive index, the general idea suggests that increasing frequency spectra
can also be obtained in the framework of conventional inflationary scenarios.
Some time ago [6] (see also [19, 20]) it was suggested that the allowed
regions of the parameter space could only be probed by some of the planned
detectors operating in the audio band or in the mHz region. Broadly speaking
the previous results assumed
$h_{0}^{2}\Omega_{gw}(\nu_{pulsar},\tau_{0})<10^{-8}$ for
$\nu_{pulsar}={\mathcal{O}}(10)$ nHz. The generous bounds previously applied
(and motivated by the early stage of the PTA) are now being improved. In this
investigation, following some recent claims, we focussed the attention on the
nHz domain and demonstrated that the recent findings, if confirmed, are
compatible with the variation of the refractive index of the relic gravitons
during the early stages of the inflationary evolution. The results reported
here have been purposely derived in the minimal framework where the high-
frequency plateau ranges between the nHz and the MHz region. If the slow-roll
condition is violated towards the end of inflation the spectral energy density
is not constant at high-frequencies and the flat plateau is replaced by a
decreasing branch. In this case the spectral energy density may even develop a
maximum between the nHz and the kHz. To address a large signal in the nHz
domain the present analysis suggests that the realistic spectra should have at
least three branches: the low-frequency branch in the aHz region originally
discussed in [5], an intermediate region where the spectral index increases
and the high-frequency domain (extending up to the MHz) where the spectrum is
nearly scale-invariant. We could easily have spectra with four different
branches [34]. This behaviour arises since the post-inflationary evolution is
not dominated by radiation but it contains either a stiff epoch or a delayed
reheating.

Figure 4: The shaded area is obtained by imposing all the previous constraints
except for Eq. (24) which is instead replaced by the limit
$h_{0}^{2}\Omega_{gw}(\nu_{ref},\tau_{0})<3.4(1)\times 10^{-10}$ for
$\nu_{ref}={\mathcal{O}}(10^{-8})$ Hz. In the plot at the right a specific
choice of the parameters. This result shows that if the PTA results are
interpreted as upper limits (as originally done when introducing the variation
of the refractive index [6]) the range of variation of $N_{*}$ and $\alpha$
gets larger even in the minimal model examined here.

The inflationary evolution of the refractive has been discussed in a minimal
scenario where the spectral energy density of cosmic gravitons is blue at
intermediate frequencies and then flattens out at high-frequency. The obtained
spectra satisfy all the constraints associated with the stochastic backgrounds
of relic gravitons and may also explain the recent findings of the PTA that
however need a convincing confirmation. While the EPTA (European Pulsar Timing
Array) [35, 36] and the PPTA (Parkes Pulsar Timing Array) [37, 38] should
eventually offer an independent (and convincing) evidence for the claim, it
would also be quite interesting to see if the present bounds coming from wide-
band interferometers could be improved by, at least, a couple of orders of
magnitude. The improved limits would allow to decide if the spectral energy-
density at high-frequencies is really flat or should instead decrease. The
theoretical perspective conveyed here suggests that the blue spectral index,
the tensor to scalar ratio and the knee of the spectra cannot be independently
assigned but their mutual relations are essential to determine the allowed
regions of the parameter space. For this reason we expect that the potential
signatures of the refractive index will be distinguishable from other
potential sources of power at intermediate frequencies. If the forthcoming
measurements will not support the evidences currently claimed by the PTA
measurements we will have to conclude that the variation of the refractive
index may be even longer and potentially lead to a large signal in the MHz
domain. In this case, however, also other spikes are potentially expected and
the improved bounds of the PTA arrays (together with the upper limits of the
wide-band detectors) will clarify the origin of the signals in the high-
frequency and ultra-high-frequency domains.

The author wishes to thank T. Basaglia, A. Gentil-Beccot, S. Rohr and J. Vigen
of the CERN Scientific Information Service for their valuable help.

## References

  * [1]
  * [2] L. P. Grishchuk, Sov. Phys. JETP 40, 409 (1975) [Zh. Eksp. Teor. Fiz. 67, 825 (1974)].
  * [3] L. P. Grishchuk, Annals N. Y. Acad. Sci. 302, 439 (1977).
  * [4] A. A. Starobinsky, JETP Lett. 30, 682 (1979) [Pisma Zh. Eksp. Teor. Fiz. 30, 719 (1979)].
  * [5] V. A. Rubakov, M. V. Sazhin and A. V. Veryaskin, Phys. Lett. B 115, 189 (1982).
  * [6] M. Giovannini, Prog. Part. Nucl. Phys. 112, 103774 (2020).
  * [7] P. Szekeres, Annals Phys. 64, 599 (1971); P. C. Peters, Phys. Rev. D 9, 2207 (1974).
  * [8] M. Giovannini, Class. Quant. Grav. 33, 125002 (2016) [arXiv:1507.03456 [astro-ph.CO]].
  * [9] S. Weinberg, Phys. Rev. D 77, 123541 (2008).
  * [10] S.-Y. Pi and R. Jackiw, Phys. Rev. D 68, 104012 (2003).
  * [11] A. Lue, L. Wang, and M. Kamionkowski, Phys. Rev. Lett. 83, 1506 (1999).
  * [12] M. Giovannini, Phys. Rev. D 99, 083501 (2019).
  * [13] H. Motohashi and A. A. Starobinsky, JCAP 11, 025 (2019).
  * [14] M. Guerrero, D. Rubiera-Garcia and D. Saez-Chillon Gomez, Phys. Rev. D 102, 123528 (2020).
  * [15] A. Mohammadi, T. Golanbari, S. Nasri and K. Saaidi, Phys. Rev. D 101, 123537 (2020).
  * [16] M. Gasperini and M. Giovannini, Phys. Lett. B 287, 56 (1992).
  * [17] I. Antoniadis, J. Rizos and K. Tamvakis, Nucl. Phys. B 415, 497 (1994).
  * [18] Z. K. Guo and D. J. Schwarz, Phys. Rev. D 80, 063523 (2009).
  * [19] M. Giovannini, Eur. Phys. J. C 78, 442 (2018).
  * [20] M. Giovannini, Phys. Rev. D 98, 103509 (2018).
  * [21] L. H. Ford and L. Parker, Phys. Rev. D 16, 1601 (1977).
  * [22] M. Giovannini, Phys. Lett. B 668, 44 (2008); Class. Quant. Grav. 26, 045004 (2009).
  * [23] Y. Akrami et al. [Planck Collaboration], Astron. Astrophys. 641, A10 (2020).
  * [24] N. Aghanim et al. [Planck Collaboration], Astron. Astrophys. 641, A6 (2020).
  * [25] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1972).
  * [26] Z. Arzoumanian et al., Astrophys. J. Lett. 905, L34 (2020).
  * [27] Z. Arzoumanian et al., Astrophys. J. 859, 47 (2018).
  * [28] Z. Arzoumanian et al. [NANOGrav], Astrophys. J. 821, 13 (2016)
  * [29] R. Abbott et al. [KAGRA, Virgo and LIGO Scientific], Phys. Rev. D 104, 022004 (2021).
  * [30] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. D 100, 061101 (2019).
  * [31] V. F. Schwartzmann, JETP Lett. 9, 184 (1969).
  * [32] M. Giovannini, H. Kurki-Suonio and E. Sihvola, Phys. Rev. D 66, 043504 (2002).
  * [33] R. H. Cyburt, B. D. Fields, K. A. Olive, and E. Skillman, Astropart. Phys. 23, 313 (2005).
  * [34] M. Giovannini, Phys. Lett. B 789, 502 (2019).
  * [35] L. Lentati et al., Mon. Not. Roy. Astron. Soc. 453, 2576 (2015).
  * [36] G. Desvignes et al., Mon. Not. Roy. Astron. Soc. 458, 3341 (2016).
  * [37] D. R. B. Yardley, W. A. Coles, G. B. Hobbs, J. P. W. Verbiest, R. N. Manchester et al. Mon. Not. Roy. Astron. Soc. 414, 1777 (2011).
  * [38] B. B. P. Perera et al., Mon. Not. Roy. Astron. Soc. 490, 4666 (2019).