[table]capposition=top capbtabboxtable[][]

# Growth of power spectrum due to decrease of sound speed during inflation

Rongrong Zhai, Hongwei Yu<EMAIL_ADDRESS>and Puxun Wu
<EMAIL_ADDRESS>Department of Physics and Synergetic Innovation Center
for Quantum Effects and Applications, Hunan Normal University, Changsha, Hunan
410081, China

###### Abstract

We study the amplification of the curvature perturbations due to a small sound
speed and find that its origin is different completely from that due to the
ultraslow-roll inflation. This is because when the sound speed is very small
the enhancement of the power spectrum comes from the fact that the curvature
perturbations at the scales smaller than the cosmic microwave background (CMB)
scale becomes scale-variant, rather than growing that leads to the
amplification of the curvature perturbations during the ultraslow-roll
inflation. At large scales the power spectrum of the curvature perturbations
remains to be scale invariant, which is consistent with the CMB observations,
and then it will have a transient $k^{2}$ growth and finally approach a
$k^{4}$ growth as the scale becomes smaller and smaller. Thus the power
spectrum can be enhanced to generate a sizable amount of primordial black
holes. Furthermore, when the high order correction in the dispersion relation
of the curvature perturbations is considered the growth of the power spectrum
of the curvature perturbations has the same origin as that in the case without
this correction.

##  I Introduction

Black holes, as one of the most mysterious objects in nature, can be formed
from the supernova explosion of massive stars. It can also be generated in the
very early Universe from the matter collapse resulting from large density
perturbations and such a black hole is dubbed a primordial black hole (PBH)
Hawking1971 ; Carr1974 ; Carr1975 . These density perturbations originate from
the curvature fluctuations during inflation Starobinsky1980 ; Guth1980 ;
Linde1982 ; Albrecht1982 , which are stretched to outside the Hubble horizon
by the exponential cosmic expansion with the amplitudes frozen at certain
values. When these superhorizon fluctuations reenter the Hubble horizon during
the radiation- or matter-dominated era, they will lead to the formation of
large scale cosmic structures and at the same time may also result in the PBH
production.

To generate a sizable amount of PBHs, it is required that the amplitude
$\mathcal{P_{R}}$ of the power spectrum of primordial curvature perturbations
reaches the $\mathcal{O}(10^{-2})$ order, as has been demonstrated by
extensive studies carried out in, for example, Refs. S. Chongchitnan ; K.
Kohri ; Byrnes2019 ; E. Bugaev . However, the cosmic microwave background
(CMB) observations have limited the primordial curvature perturbations to be
about $10^{-5}$ with the amplitude of the corresponding power spectrum being
about $10^{-9}$ at the CMB scale Aghanim2020 . Thus, to be consistent with the
CMB observations and at the same time to generate a sizable amount of PBHs,
the primordial curvature perturbations can only be enhanced at scales smaller
than the CMB scale during inflation. Since
$\mathcal{P_{R}}\propto\frac{1}{c_{s}\epsilon}$, where
$\epsilon=-\frac{\dot{H}}{H^{2}}$ is the slow-roll parameter, with $H$ being
the Hubble parameter and $\dot{H}=\mathrm{d}H/\mathrm{d}t$, and $c_{s}$ is the
sound speed of the curvature perturbations, the enhancement of the curvature
perturbations during inflation can be achieved by reducing the rolling speed
of the inflaton to obtain a phase of ultraslow-roll inflation Germani2017 ;
Motohashi2017 ; Ezquiaga2017 ; H. Di2018 ; Ballesteros2018 ; Dalianis2019 ;
Gao2018 ; Drees2021 ; C.Fu2020 ; Xu2020 ; Lin2020 ; Dalianis2021 ; Yi2021 ;
Gao2021 ; Yi2021b ; TGao2021 ; Solbi2021 ; Gao2021b ; Solbi2021b ; Zheng2021 ;
Teimoori2021a ; Cai2021 ; Wang2021 ; Fuchengjie2019 ; fuchengjie2020 ;
Dalianis(2020) ; Teimoori2021 ; Heydari2022 ; Heydari2022b ; Bhaumik2020 ;
Karam2022  or the sound speed of the curvature perturbations Ballesteros2022 ;
Gorji2022 ; Ballesteros2019 ; Romano2016 ; Romano3 ; Ozsoy2018 .

Slowing down the inflaton can be naturally realized by flattening the
inflationary potential as well as by increasing the friction Fuchengjie2019 ;
Germani2011_1 ; Germani2011_2 ; Tsujikawa2012  or introducing a noncanonical
kinetic term Lin2020 . When the cosmic evolution changes from the slow-roll
inflation to the ultraslow-roll one, the slow-roll parameter $\eta$, which is
defined as $\eta=\frac{\dot{\epsilon}}{\epsilon H}$, will change from $\sim 0$
to $-6$. Solving the equation of motion for the curvature perturbations, one
can find that the solution for modes outside the Hubble horizon during the
slow-roll inflation consists of a constant term and a decaying one, which
leads to a nearly scale-invariant spectrum. Assuming that the Universe
transitions suddenly at a time from the slow-roll inflation to the ultraslow-
roll one, we can match the solution of the curvature perturbations and its
first derivative at this time by the Israel junction conditions Israel1966 ;
Deruelle1995 . After this transition, the term, which is decaying with time
before the transition, becomes “growing”. At first, the solution is dominated
by the constant term and the power spectrum of the curvature perturbations
remains to be scale invariant [$\mathcal{P_{R}}(k)\propto k^{0}$]. If the
ultraslow-roll inflation lasts sufficiently long, the growing term will become
dominant, which results in a $k^{4}$-dependent power spectrum Garrilho2019 ;
Byrnes2019 ; Liu2020 . Because the constant and growing terms cancel each
other, the power spectrum has a dip preceding the $k^{4}$ dependence. If an
$\eta=-1$ middle phase is added between the slow- and ultraslow-roll
inflations, it has been found that the power spectrum has a $k^{5}(\log
k)^{2}$ growth Garrilho2019 ; Cole2022 .

Similar to the case of a tiny $\epsilon$, a smaller value of the sound speed
can also lead to an amplification of the primordial curvature perturbations.
For a scalar field with a standard kinetic term, the sound speed equals the
light speed ($c_{s}=1$). A varying $c_{s}$ can be realized in many
inflationary models including the inflation with a noncanonical kinetic term
Armendariz-Picon1999 ; Garriga1999 ; Ballesteros2019 ; A. Y. Kamenshchik , the
Dirac-Born-Infeld inflation Kamenshchik2022 , the multifield inflationary
model Romano1 ; Romano2 , the inflationary models in modified gravity theories
CGermaniprl ; S.Tsujikawa , and so on. Using the effective field theory of
inflation, the generation of PBHs due to a small $c_{s}$ has been studied in
Ballesteros2022 ; Gorji2022 ; Ballesteros2019 ; Romano3 . However, whether the
growth of the power spectrum due to decreasing sound speed is similar to that
due to ultraslow-roll inflation and whether this growth results from the
appearance of a growing term in the solution of the curvature perturbations
remains unclear. Recently, it was pointed out that when the sound speed is
very small the high order $k^{4}$ term in the dispersion relation will become
dominant Ballesteros2022 ; Gorji2022 . Although this high order term leads to
a new solution of the curvature perturbations, the power spectrum is enhanced
by small $c_{s}$ and it has a $k^{4}$ growth, which is the same as that in the
case of the ultraslow-roll inflation. Can both $k^{4}$ growths of the power
spectrum have the same origin? The purpose of the present paper is to answer
the above mentioned questions by studying in detail the growth of the power
spectrum of the curvature perturbations due to the decrease of the sound
speed.

The paper is organized as follows: In Sec. II, we investigate the growth of
the power spectrum due to a sudden decrease of the sound speed. In Sec. III,
the effect of the high order $k^{4}$ term in the dispersion relation is
analyzed. Finally, we give our conclusions in Sec. IV. Throughout this paper,
we set $c=\hbar=M_{\mathrm{Pl}}=1$.

##  II the growth of power spectrum from a small sound speed

To study the growth of the power spectrum of the curvature perturbations due
to a sudden decrease of the sound speed, we need to consider the Sasaki-
Mukhanov equation

$\displaystyle
v^{\prime\prime}_{k}+\left(c_{s}^{2}k^{2}-\frac{z^{\prime\prime}}{z}\right)v_{k}=0$
(1)

in Fourier space, where a dash denotes a differentiation with respect to the
conformal time $\tau$, $v_{k}=z\mathcal{R}_{k}$, $\mathcal{R}$ is the
curvature perturbation, and $z$ is defined to be

$\displaystyle z^{2}\equiv\frac{2a^{2}\epsilon}{c_{s}^{2}}$ (2)

with $a$ being the cosmic scale factor. From Eq. (2), one can obtain that

$\displaystyle\frac{z^{\prime\prime}}{z}=(aH)^{2}\left(2-\epsilon+\frac{3}{2}\eta-3s+s^{2}+s\epsilon-s\eta+\frac{1}{4}\eta^{2}-\frac{1}{2}\eta\epsilon\right)\;,$
(3)

where $s=\frac{\dot{c}_{s}}{c_{s}H}$. Equation (3) can be reexpressed to be

$\displaystyle\frac{z^{\prime\prime}}{z}=\frac{\nu^{2}-1/4}{(-\tau)^{2}}\;,$
(4)

where $aH=-\frac{1}{\tau}$ has been used, and
$\nu\simeq\frac{3}{2}+\frac{1}{3}\epsilon+\frac{1}{2}\eta-s$. If $c_{s}$ is a
constant, the general solution of Eq. (1) has the form

$\displaystyle
v_{k}(\tau)=A\sqrt{-\tau}H_{\nu}^{(1)}(-c_{s}k\tau)+B\sqrt{-\tau}H_{\nu}^{(2)}(-c_{s}k\tau)\;.$
(5)

Here $A$ and $B$ are two constants. Thus, the solution of the curvature
perturbations can be obtained through the relation
$\mathcal{R}_{k}=\frac{v_{k}}{z}$.

Now we consider that during the slow-roll inflation the sound speed $c_{s}$
will change suddenly from $1$ to a very small value denoted by $1/A_{s}$ at
conformal time $\tau_{1}$, where $A_{s}$ is a constant and $A_{s}\gg 1$. Since
$\nu\simeq-3/2$ for the slow-roll inflation with a constant sound speed, we
find that when $c_{s}=1$ the solution of the curvature perturbations has the
form

$\displaystyle\mathcal{R}^{(1)}_{k}(\tau)=i\frac{H}{2\sqrt{\epsilon
k^{3}}}e^{-ik\tau}(1+ik\tau)$ (6)

after choosing the adiabatic Bunch-Davis vacuum condition.

For $c_{s}=1/A_{s}$, from Eq. (5), we obtain the general solution of the
curvature perturbations

$\displaystyle\mathcal{R}^{(2)}_{k}(\tau)=\frac{H}{2\sqrt{\pi A_{s}\epsilon
k^{3}}}\left[A_{2}e^{-\frac{ik\tau}{A_{s}}}(A_{s}+ik\tau)-B_{2}e^{\frac{ik\tau}{A_{s}}}(A_{s}-ik\tau)\right]\;,$
(7)

where $A_{2}$ and $B_{2}$ are two constants. Matching $\mathcal{R}^{(1)}_{k}$
and $\mathcal{R}^{(2)}_{k}$ at $\tau=\tau_{1}$ by using the conditions
$\mathcal{R}^{(1)}_{k}(\tau_{1})=\mathcal{R}^{(2)}_{k}(\tau_{1})$ and
$\mathcal{R}^{\prime(1)}_{k}(\tau_{1})=\mathcal{R}^{\prime(2)}_{k}(\tau_{1})$,
one can obtain that

$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{A_{s}\pi}}{2k}(1+A_{s})(-1+ik+A_{s})e^{ik\big{(}1-\frac{1}{A_{s}}\big{)}},$
(8) $\displaystyle B_{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{A_{s}\pi}}{2k}(-1+A_{s})(1-ik+A_{s})e^{ik\big{(}1+\frac{1}{A_{s}}\big{)}},$
(9)

where we have set $\tau_{1}=-1$, which simplifies the expressions but does not
change any physical result. Substituting $A_{2}$ and $B_{2}$ into Eq. (7)
gives the expression of the curvature perturbations during the phase of
$c_{s}=1/A_{s}$. Expanding this expression in the infrared limit
($c_{s}k\tau\rightarrow 0$), we arrive at

$\displaystyle\mathcal{R}^{(2)}_{k}(\tau)\simeq\frac{iHA_{s}e^{ik}}{2\sqrt{\epsilon
k^{3}}}C_{2}-\dfrac{iHA_{s}e^{ik}}{2\sqrt{\epsilon
k^{5}}}D_{2}-\frac{iHe^{ik}}{4A_{s}\sqrt{\epsilon
k}}(-\tau)^{2}D_{2}+\cdots\;.$ (10)

where

$\displaystyle C_{2}$ $\displaystyle=$ $\displaystyle
A_{s}\cos\left(\frac{k}{A_{s}}\right)-i\sin\left(\frac{k}{A_{s}}\right)$
$\displaystyle D_{2}$ $\displaystyle=$
$\displaystyle\left(A_{s}^{2}-1\right)\sin\left(\frac{k}{A_{s}}\right).$ (11)

For modes that are superhorizon at $\tau_{1}$, we can further expand Eq. (II)
at $k\tau_{1}=-k\rightarrow 0$ and then obtain

$\displaystyle\mathcal{R}^{(2)}_{k}(\tau)\simeq\frac{iHe^{ik}}{2\sqrt{\epsilon
k^{3}}}+\dfrac{He^{ik}}{2\sqrt{\epsilon
k}}-\frac{iH(2A_{s}^{2}+1)e^{ik}k^{1/2}}{12A_{s}^{2}\sqrt{\epsilon}}-\frac{iH(A_{s}^{2}-1)e^{ik}k^{1/2}}{4A_{s}^{2}\sqrt{\epsilon}}(-\tau)^{2}+\cdots\;.$
(12)

It is easy to see that the leading part is constant independent of $\tau$,
which contains three different $k$-dependent terms and the subleading part
decays with time since $|\tau|$ decreases during inflation. These characters
are different from that in the case of the transition from the slow-roll
inflation to the ultraslow-roll one, where there is a growing part and the
constant part only has a $k^{-3/2}$-dependent term Ballesteros2022 .

Equation (12) gives that the power spectrum of the curvature perturbations has
the form

$\displaystyle\mathcal{P}_{\mathcal{R}^{(2)}_{k}}=\frac{k^{3}}{2\pi^{2}}|\mathcal{R}^{(2)}_{k}|^{2}\simeq\frac{H^{2}}{8\pi^{2}\epsilon}+\frac{(A_{s}^{2}-1)H^{2}}{24\pi^{2}A_{s}^{2}\epsilon}k^{2}+\frac{(2A_{s}^{2}+1)^{2}H^{2}}{2\left(12\pi
A_{s}^{2}\right)^{2}\epsilon}k^{4}\;.$ (13)

Apparently, the power spectrum of the curvature perturbations consists of a
$k$-independent term and two $k$-dependent ones. At the CMB scale, the first
term dominates, which leads to a scale-invariant spectrum consistent with the
CMB observations. Going to the scales that are smaller than the CMB scale, the
second term begins to play a role. The power spectrum becomes scale dependent
and has a short era with a $k^{2}$ growth. Then, the third term finally
dominates and the power spectrum displays a $k^{4}$ growth. These results are
shown clearly in Fig 1. This figure indicates that there is no dip in the
power spectrum, although the dip appears in the power spectrum in the
ultraslow-roll inflation. This is because the second term in Eq. (13) does not
cancel the first term. From Fig. 1, one can see that the power spectrum
oscillates at the small scales. The reason is that although these small scales
are superhorizon at $\tau$, they are sub-horizon at $\tau_{1}$. For subhorizon
scales at $\tau_{1}$, Eq. (II) shows clearly that $C_{2}$ and $D_{2}$ are
oscillating functions, which results in an oscillating power spectrum.

Figure 1:  The evolution of power spectrum $\mathcal{P_{R}}$ as functions of
wave number $k$ with $A_{s}=50$. The dashed red line indicates the $k^{4}$
slope.

##  III The effect of higher order correction in dispersion relation

When $c_{s}$ is very small, the high order correction in the dispersion
relation of the curvature perturbations needs to be considered Ballesteros2022
; Gorji2022 ; Ballesteros2019 . When the $k^{4}$ corrected term in the
dispersion relation is included, the numerical calculation shows that the
power spectrum has a $k^{4}$ growth Ballesteros2022 ; Gorji2022 . In this
section, we discuss the origin of this growth. We find that $v_{k}$ satisfies
the following equation

$\displaystyle v_{k}^{\prime\prime}+\left[c_{s}^{2}k^{2}+\alpha
k^{4}(-\tau)^{2}-\frac{\nu^{2}-1/4}{(-\tau)^{2}}\right]v_{k}=0\;.$ (14)

where $\alpha$ is a constant. When $c_{s}=1$, this $k^{4}$ corrected term can
be neglected and $v_{k}$ has a general solution shown in Eq. (5). Thus,
choosing the adiabatic Bunch-Davis vacuum condition the solution of
$\mathcal{R}_{k}^{(1)}$ is given in Eq. (6).

For a very small $c_{s}$, which implies a very large $A_{s}$ after setting
$c_{s}=\frac{1}{A_{s}}$, the $k^{4}$ corrected term in Eq. (14) will dominate
over $c_{s}^{2}k^{2}$ and the latter can be neglected. Solving Eq. (14) and
using $\mathcal{R}_{k}=\frac{v_{k}}{z}$ and $\nu=-3/2$, we find that the
curvature perturbation has the general solution

$\displaystyle\mathcal{R}^{(2)}_{k}(\tau)=-\frac{\alpha^{1/8}H}{2A_{s}\sqrt{\epsilon}}k^{1/2}(-\tau)^{3/2}\left[A_{3}H_{3/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\tau^{2}\right)+B_{3}H_{3/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\tau^{2}\right)\right]\;,$
(15)

where $A_{3}$ and $B_{3}$ are two constants. Using the matching condition:
$\mathcal{R}^{(1)}_{k}(\tau_{1})=\mathcal{R}^{(2)}_{k}(\tau_{1})$ and
$\mathcal{R}^{\prime(1)}_{k}(\tau_{1})=\mathcal{R}^{\prime(2)}_{k}(\tau_{1})$
at time $\tau_{1}$ which is set to be $-1$, one can obtain that

$\displaystyle A_{3}$
$\displaystyle=\frac{A_{s}e^{ik}~{}C}{\alpha^{5/8}k^{4}~{}E}\;,$ (16)
$\displaystyle B_{3}$
$\displaystyle=\frac{A_{s}e^{ik}~{}D}{\alpha^{5/8}k^{4}~{}E}\;,$ (17)

where

$\displaystyle C=$ $\displaystyle
i\sqrt{\alpha}(i+k)k^{2}H_{-1/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)-(3-3ik-2k^{2})H_{3/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)$
$\displaystyle+\sqrt{\alpha}(1-ik)k^{2}H_{7/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right),$
$\displaystyle D=$
$\displaystyle\sqrt{\alpha}(1-ik)k^{2}H_{-1/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)+(3-3ik-2k^{2})H_{3/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)$
$\displaystyle-\sqrt{\alpha}(1-ik)k^{2}H_{7/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right),$
$\displaystyle E=$ $\displaystyle
H_{3/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)\left[H_{-1/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)-H_{7/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)\right]$
$\displaystyle+H_{3/4}^{(1)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)\left[H_{7/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)-H_{-1/4}^{(2)}\left(\frac{1}{2}\sqrt{\alpha}k^{2}\right)\right].$
(18)

In the infrared region, we have

$\displaystyle\mathcal{R}^{(2)}_{k}(\tau)\simeq\frac{He^{ik}}{2\sqrt{\epsilon
k^{3}}}-\dfrac{iHe^{ik}}{2\sqrt{\epsilon
k}}-\frac{He^{ik}k^{1/2}}{6\sqrt{\epsilon}}+\frac{He^{ik}k^{1/2}}{6\sqrt{\epsilon}}(-\tau)^{3}+\cdots\;.$
(19)

This solution consists of three constant terms and the decaying terms, but
there is no growing term. This is similar to the case studied in the previous
section. Then, we obtain the power spectrum during the small sound speed
phase,

$\displaystyle\mathcal{P}_{\mathcal{R}^{(2)}_{k}}=\frac{k^{3}}{2\pi^{2}}|\mathcal{R}^{(2)}_{k}|^{2}\simeq\frac{H^{2}}{8\pi^{2}\epsilon}+\frac{H^{2}}{24\pi^{2}\epsilon}k^{2}+\frac{H^{2}}{72\pi^{2}\epsilon}k^{4}\;.$
(20)

This is the same as that given in Eq. (13) since $A_{s}\gg 1$. So, the high
order correction in the dispersion relation has no contribution to the growth
of the power spectrum of the curvature perturbations.

##  IV conclusion

A generation of abundant PBHs in the early Universe requires that the
amplitude of the power spectrum of the primordial curvature perturbations is
enhanced by about 7 orders during inflation, which can be realized by the
ultraslow-rolling of inflaton, a very small sound speed of the curvature
perturbations, and even some other mechanisms. It has been found that in the
ultraslow-roll inflation the amplification of the power spectrum can be
attributed to the appearance of a growing solution of the curvature
perturbations. In this paper, we find that the origin of the amplification of
the curvature perturbations due to a small sound speed is different completely
from that in the case of the ultraslow-roll inflation. This is because when
the sound speed is very small the enhancement of the power spectrum comes from
the fact that the curvature perturbations at the scales smaller than the CMB
scale becomes scale variant, rather than growing. The power spectrum of
curvature perturbations remains to be scale invariant at large scales, and
then it has a short time $k^{2}$ growth and finally approaches a $k^{4}$
growth as the scale becomes smaller and smaller. Thus, the curvature
perturbations during inflation with a decrease of the sound speed can be
consistent with the CMB observations and at the same time enhanced to generate
a sizable amount of PBHs. Furthermore, we find that when the high order
correction in the dispersion relation of the curvature perturbations is
considered, the growth of the power spectrum of the curvature perturbations
has the same origin as that without this correction.

###### Acknowledgements.

We appreciate very much the insightful comments and helpful suggestions by the
anonymous referee. This work is supported by the National Key Research and
Development Program of China Grant No. 2020YFC2201502, the National Natural
Science Foundation of China under Grants No. 12075084, No. 11775077 and No.
11805063, and the Science and Technology Innovation Plan of Hunan Province
under Grant No. 2017XK2019.

## References

  * (1) S. Hawking, Mon. Not. R. Astron. Soc. 152, 75 (1971).
  * (2) B. J. Carr and S. W. Hawking Mon. Not. R. Astron. Soc.168, 399 (1974).
  * (3) B. J. Carr,  Astrophys. J. 201, 1 (1975).
  * (4) A. H. Guth, Phys. Rev. D 23, 347 (1980).
  * (5) A. D. Linde, Phys. Lett. 108B, 389 (1982).
  * (6) A. A. Starobinsky, Phys. Lett. 91B, 99 (1980).
  * (7) A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982).
  * (8) S. Chongchitnan and G. Efstathiou, J. Cosmol. Astropart. Phys. 01 (2007) 011
  * (9) K. Kohri, D. H Lyth and A. Melchiorri, J. Cosmol. Astropart. Phys. 04 (2008) 038.
  * (10) C. T. Byrnes, P. S. Cole and S. P. Patil, J. Cosmol. Astropart. Phys. 06, (2019) 028
  * (11) E. Bugaev and P. Klimai, Phys. Rev. D 78, 063515 (2008).
  * (12) N. Aghanim et al. (Planck Collaboration), Astron. Astrophys. 641, A6 (2020).
  * (13) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6 (2017).
  * (14) H. Motohashi and W. Hu, Phys. Rev. D 96, 063503 (2017).
  * (15) J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales, Phys. Lett. 039B, 11 (2017).
  * (16) H. Di and Y. Gong, J. Cosmol. Astropart. Phys. 07 (2018) 007.
  * (17) G. Ballesteros and M. Taoso, Phys. Rev. D 97, 023501 (2018).
  * (18) I. Dalianis, A. Kehagias, and G. Tringas, J. Cosmol. Astropart. Phys. 01 (2019) 037.
  * (19) T.-J. Gao and Z.-K. Guo, Phys. Rev. D 98, 063526 (2018).
  * (20) M. Drees and Y. Xu, Eur. Phys. J. C 81, 182 (2021).
  * (21) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 102, 043527 (2020).
  * (22) W. Xu, J. Liu, T. Gao, and Z. Guo, Phys. Rev. D 101, 023505 (2020).
  * (23) J. Lin, Q. Gao, Y. Gong, Y. Lu, and C. Zhang, Phys. Rev. D 101, 103515 (2020).
  * (24) I. Dalianis and K. Kritos, Phys. Rev. D 103, 023505 (2021).
  * (25) Z. Yi, Y. Gong, B. Wang, and Z. Zhu, Phys. Rev. D 103, 063535 (2021).
  * (26) Q. Gao, Y. Gong, and Z. Yi, Nucl. Phys. B969, 115480 (2021).
  * (27) Z. Yi, Q. Gao, Y. Gong, and Z. Zhu Phys. Rev. D 103, 063534 (2021).
  * (28) T. Gao and X. Yang, Eur. Phys. J. C 81, 494 (2021).
  * (29) M. Solbi and K. Karami, J. Cosmol. Astropart. Phys. 08 (2021) 056.
  * (30) Q. Gao, Sci. China Phys. Mech. Astron. 64, 280411 (2021).
  * (31) M. Solbi and K. Karami, Eur. Phys. J. C 81, 884 (2021).
  * (32) R. Zheng, J. Shi, and T. Qiu, arXiv: 2106.04303.
  * (33) Z. Teimoori, K. Rezazadeh, M. A. Rasheed, and K. Karami, J. Cosmol. Astropart. Phys. 10 (2021) 018.
  * (34) R. Cai, C. Chen, and C. Fu, Phys. Rev. D 104, 083537 (2021).
  * (35) Q. Wang, Y. Liu, B. Su, and N. Li, Phys. Rev. D 104, 083546 (2021).
  * (36) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 100, 063532 (2019).
  * (37) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 101, 023529 (2020).
  * (38) I. Dalianis, S. Karydas, and E. Papantonopoulos, J. Cosmol. Astropart. Phys. 06 (2020) 040.
  * (39) Z. Teimoori, K. Rezazadeh, and K. Karami, Astrophys. J. 915, 118 (2021).
  * (40) S. Heydari, and K. Karami,  Eur. Phys. J. C 82, 83 (2022).
  * (41) S. Heydari, and K. Karami, J. Cosmol. Astropart. Phys. 03 (2022) 033.
  * (42) N. Bhaumik and R. K. Jain, J. Cosmol. Astropart. Phys. 01 (2020) 037 .
  * (43) A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen, and H. Veermäe, arXiv: 2205.13540.
  * (44) G. Ballesteros, S. Céspedes, and L. Santoni, J. High Energy Phys. (2022) 074.
  * (45) M. A. Gorji, H. Mothhashi, and S. Mukohyama, J. Cosmol. Astropart. Phys. 02 (2022) 030.
  * (46) G. Ballesteros, J. B. Jiménez, and M. Pieroni, J. Cosmol. Astropart. Phys. 06 (2019) 016.
  * (47) A. E. Romano, arXiv: 2006.07321.
  * (48) A. E. Romano, S. Mooij, and M. Sasaki, Phys. Lett. B 761, 119 (2016).
  * (49) O. Özsoy, S. Parameswaran, G. Tasinato, and I. Zavala, J. Cosmol. Astropart. Phys. 07 (2018) 005.
  * (50) C. Germani and A. Kehagias, Phys. Rev. Lett. 106, 161302 (2011).
  * (51) C. Germani and Y. Watanabe, J. Cosmol. Astropart. Phys. 07 (2011) 031.
  * (52) S. Tsujikawa, Phys. Rev. D 85, 083518 (2012).
  * (53) W. Israel, Nuovo Cim. B 44, 1 (1966).
  * (54) N. Deruelle and V. F. Mukhanov, Phys. Rev. D 52, 5549 (1995).
  * (55) P. Carrilho, K. A. Malik and D. J. Mulryne, Phys. Rev. D 100, 103529 (2019).
  * (56) J. Liu, Z. Guo, and R. Cai, Phys. Rev. D 101, 083535 (2020).
  * (57) P. S. Cole, A. D. Gow, C. T. Byrnes, and S. P. Patil, arXiv: 2204.07573.
  * (58) C. Armendariz-Picon, T. Damour and V. F. Mukhanov, Phys. Lett. B 458, 209 (1999).
  * (59) J. Garriga and V.F. Mukhanov, Phys. Lett. 458B, 219 (1999).
  * (60) A. Y. Kamenshchik, A. Tronconi, T. Vardanyan, and G. Venturi, Phys. Lett. 791B, 201(2019).
  * (61) A. Y. Kamenshchik, A. Tronconi, and G. Venturi, J. Cosmol. Astropart. Phys. 01 (2022) 051.
  * (62) A. E. Romano, S. A. Vallejo Peña, and K. Turzyński, arXiv:2006.00969.
  * (63) A. E. Romano and S. A. Vallejo Peña, Phys. Lett. 784B, 367 (2018).
  * (64) C. Germani and A. Kehagias, Phys. Rev. Lett. 105, 011302 (2010).
  * (65) S. Tsujikawa, Phys. Rev. D 85, 083518 (2012).


# Exploring the dynamics of finite-energy Airy beams: A trajectory analysis
perspective

Ángel S. Sanz* and Rosario Martínez-Herrero  Department of Optics, Faculty of
Physical Sciences, Universidad Complutense de Madrid,
Pza. Ciencias 1, 28040 Madrid, Spain<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In practice, Airy beams can only be reproduced in an approximate manner, with
a limited spatial extension and hence a finite energy content. To this end,
different procedures have been reported in the literature, based on a
convenient tuning of the transmission properties of aperture functions. In
order to investigate the effects generated by the truncation and hence the
propagation properties displayed by the designed beams, here we resort to a
new perspective based on a trajectory methodology, complementary to the
density plots more commonly used to study the intensity distribution
propagation. We consider three different aperture functions, which are
convoluted with an ideal Airy beam. As it is shown, the corresponding
trajectories reveals a deeper physical insight about the propagation dynamics
exhibited by the beams analyzed due to their direct connection with the local
phase variations undergone by the beams, which is in contrast with the global
information provided by the usual standard tools. Furthermore, we introduce a
new parameter, namely, the escape rate, which allow us to perform piecewise
analyses of the intensity distribution without producing any change on it,
e.g., determining unambiguously how much energy flux contributes to the
leading maximum at each stage of the propagation, or for how long self-
accelerating transverse propagation survives. The analysis presented in this
work thus provides an insight into the behavior of finite-energy Airy beams,
and therefore is expected to contribute to the design and applications
exploiting this singular type of beams.

††journal: opticajournal††articletype: Research Article

## Introduction

In 1979, Berry and Balazs found [1] that the free-propagation Schrödinger
equation, a rather simple and even trivial equation, admitted a rather
counter-intuitive nondispersive solution. Provided the initial ansatz is
described by an Airy function, these authors proved that the time-evolution of
the corresponding probability density remains shape-invariant during the whole
propagation. At the same time, without the presence of any external potential,
these wave packets undergo a self-accelerating motion pushing them forwards.
When this behavior is analyzed within the framework of the equivalence
principle, as it was done by Greenberger [2], it can be shown that
Schrödinger’s equation can be recast as the time-independent equation
describing a free falling mass, of which the eigensolutions correspond
precisely to Airy functions. Later on, Unnikrishnan and Rau [3] proved the
uniqueness of these solutions to Schrödinger’s equation as the only free-space
self-accelerating nondispersive one.

The beams conjectured by Berry and Balazs have not become a reality until
their experimental realization, formerly produced in 2007 by Siviloglou et al.
[4], taking advantage of the isomorphism between the free-space Schrödinger
equation for a non-relativist particle with a mass $m$, and the paraxial
Helmholtz equation for monochromatic light propagating in a linear medium with
a refractive index $n$. Ever since, Airy beams have become one of the relevant
components in the puzzle of structured light [5, 6], a very active and
promising research field strongly developed over the last decade due to the
interesting fundamental and technological applications brought in by
conveniently designed light beams. In this regard, different control
mechanisms for their spatial manipulation [7] and applications [8] have been
discussed in the literature.

From the point of view of their design, Airy beams introduce an important
technical question, namely, their spatial extension in real laboratory
realizations is limited, unlike the infinite tail that characterizes Airy
functions. In other words, although the energy content of an ideal Airy beam
is infinite (this is another interesting peculiarity of this singular
solutions), in practice their energy content is finite due to the unavoidable
spatial truncation under experimental conditions. A way to overcome this
drawback was proposed by Siviloglou and Christodoulides [9], which consisted
in using an exponential aperture, such that it has no influence on the beam at
large positive values of $x$, while its tail of secondary maxima is gradually
annihilated for negative $x$-values. This method allows a propagation of the
beam for rather long distances, keeping its distinctive features, but it does
not prevent eventually the diffraction and, consequently, the loss of those
features [4].

The introduction of such aperture functions gives rise to beams physically
realizable in practice. Another method to produce finite-energy beams consists
of considering partially coherent Airy beams, as shown in [10]. Here, however,
motivated by the above procedure, we analyze the effect of different aperture
functions in order to determine which one shows a better performance. However,
instead of directly introducing a spatial truncation, we proceed in the
Fourier space, which translates into convoluting the aperture function with an
ideal Airy beam. A similar procedure has been recently considered
experimentally by Xu et al. [11], achieving off-axis autofocusing and
transverse self-acceleration. In our case, to investigate the properties of
these aperture function, we consider a trajectory-based methodology, which
allows us to follow the transverse energy flux at a local level, in a
hydrodynamic manner. These trajectories are directly connected to the local
spatial variations of the phase undergone by the beam, describing the typical
parabolic evolution in the ideal case. Furthermore, as it is also shown, they
make possible to determine how much a given feature degrades as the beam
propagates forward, like the main maximum of a finite-energy Airy beam, in a
very precise manner. In this work, we thus merge the information obtained from
global behaviors, studied by means of the standard propagation of the
intensity distribution, with local phase effects, which are examined by virtue
of the computation and subsequent analysis of the corresponding flux
trajectories.

## Theory

### General aspects

Consider the paraxial monochromatic scalar field

$\Psi({\bf r}_{\perp},z)=\psi({\bf r}_{\perp},z)e^{ikz},$ (1)

where ${\bf r}_{\perp}=(x,y)$ and $k=2\pi n/\lambda_{0}$ (with
$\lambda_{0}=\omega/c$). The complex-valued amplitude $\psi({\bf
r}_{\perp},z)$ satisfies the paraxial Helmholtz equation,

$i\ \frac{\partial\psi({\bf r}_{\perp},z)}{\partial z}=-\frac{1}{2k}\
\nabla_{\perp}^{2}\psi({\bf r}_{\perp},z),$ (2)

where $\nabla_{\perp}^{2}=\partial^{2}/\partial x^{2}+\partial^{2}/\partial
y^{2}$ is the transverse Laplacian. If this equation is multiplied on both
sides by $\psi^{(}{\bf r}_{\perp},z)$, and then we subtract to the resulting
equation its complex conjugate version, we obtain

$\frac{\partial|\psi|^{2}}{\partial
z}=\frac{1}{2ik}\left[\psi^{*}\left(\nabla_{\perp}^{2}\psi\right)-\left(\nabla_{\perp}^{2}\psi^{*}\right)\psi\right].$
(3)

Rearranging terms on the right-hand side of this equation, it can be further
simplified in the form of the following transport equation

$\frac{\partial|\psi|^{2}}{\partial
z}=-\nabla_{\perp}\left[|\psi|^{2}\frac{1}{k}\ {\rm
Re}\left(\frac{-i\nabla_{\perp}\psi}{\psi}\right)\right],$ (4)

which describes the transversal changes undergone by the intensity
distribution, $|\psi|^{2}$, as it propagates along the $z$-direction. Such
changes are determined by the term on the right-hand side, where the flux
density [12] is associated here with a transverse electromagnetic flux vector,

$\displaystyle{\bf j}_{\perp}({\bf r}_{\perp},z)$ $\displaystyle=$
$\displaystyle\frac{1}{2ik}\left\\{\psi^{*}({\bf
r}_{\perp},z)\left[\nabla_{\perp}\psi({\bf
r}_{\perp},z)\right]-\left[\nabla_{\perp}\psi^{*}({\bf
r}_{\perp},z)\right]\psi({\bf r}_{\perp},z)\right\\}$ (5) $\displaystyle=$
$\displaystyle|\psi({\bf r}_{\perp},z)|^{2}{\bf v}_{\perp}({\bf
r}_{\perp},z).$

Actually, according to the usual definition of a flux density, the second line
of Eq. (5) is specified in terms of the intensity distribution and the
transport field

${\bf v}_{\perp}({\bf r}_{\perp},z)=\frac{1}{k}\ {\rm
Re}\left(\frac{-i\nabla_{\perp}\psi}{\psi}\right).$ (6)

This field quantifies the local rate of change of the flux, analogous to a
local velocity field describing the flow of energy on a given point ${\bf
r}_{\perp}$ at a given value of the $z$ distance. Hence, from now on, we shall
refer to ${\bf v}_{\perp}$ as the effective transverse velocity field.

As it can be noticed, if the complex field amplitude is recast in the
following polar form

$\psi({\bf r}_{\perp},z)=A({\bf r}_{\perp},z)e^{iS({\bf r}_{\perp},z)},$ (7)

where both $A({\bf r}_{\perp},z)$ and $S({\bf r}_{\perp},z)$ are real-valued
fields, the above defined velocity field reads as

${\bf v}_{\perp}({\bf r}_{\perp},z)=\frac{1}{k}\ \nabla_{\perp}S({\bf
r}_{\perp},z).$ (8)

This expression unambiguously shows that any change in the transverse flux is
directly associated with the local (transverse) spatial variations undergone
by the phase of the field phase during its propagation, which eventually leads
to measurable changes in the intensity distribution. Accordingly, to better
understand these physical consequences, one can compute ensembles of
trajectories obtained from the integration of the equation of motion

$\frac{d{\bf r}_{\perp}(z)}{dz}={\bf v}_{\perp}({\bf r}_{\perp},z)=\frac{{\bf
j}_{\perp}({\bf r}_{\perp},z)}{|\psi({\bf r}_{\perp},z)|^{2}},$ (9)

which will provide us with first-hand local information directly associated
with the phase of the field (or, more strictly speaking, its spatial phase
variations) rather than from the profile displayed by the intensity
distribution (from which any phase effect is indirectly inferred a
posteriori). In this regard, these trajectories have already been used to
describe light propagation through optical fibers [13], to account for the
differences between Gaussian beams and Maxwell-Gauss beams [14], or, more
recently, to describe the dynamics of ideal and finite-energy Airy beams [15]
or to investigate light focusing [16]. Furthermore, it has also been used to
establish a connection between quantum and optical analogs given the
isomorphism between the paraxial Helmholtz equation and the time-dependent
Schrödinger equation [17].

In analogy to their quantum-mechanical counterparts [17], the trajectories
obtained after integrating Eq. (9) have three distinctive properties with
practical interest concerning the problem that we are dealing with here:

  * (i)

The trajectories cannot get across the same point ${\bf r}_{\perp}$ at the
same value of $z$. This is a direct consequence that arises from the single
valuedness of the solutions of the paraxial Helmholtz equation (2), except for
a constant phase factor. From the solution in polar form (7), we note that
this phase factor corresponds to an integer multiple of $2\pi$, which does not
affect the trajectory dynamics, as it can be inferred from the Eq. (8).

  * (ii)

Consider a closed curve $\mathcal{C}({\bf r}_{\perp},0)$ at $z=0$ defined by a
set of initial conditions ${\bf r}_{\perp}(0)$ for Eq. (9). At any other value
of $z$, $\mathcal{C}({\bf r}_{\perp},0)$ will uniquely evolve into another
closed curve $\mathcal{C}({\bf r}_{\perp},z)$ consisting of the positions
${\bf r}_{\perp}(z)$ of the trajectories started at ${\bf r}_{\perp}(0)$.

  * (iii)

Following properties (i) and (ii), the partial energy content within the
surface $\Sigma({\bf r}_{\perp},z)$ determined by the curve $\mathcal{C}({\bf
r}_{\perp},z)$, defined as

$E_{\Sigma({\bf r}_{\perp},z)}=\int_{\Sigma({\bf r}_{\perp},z)}|\psi({\bf
r}_{\perp},z)|^{2}d{\bf
r}_{\perp}\propto\lim_{N\to\infty}\sum_{i=1}^{N}\delta\left[\\{{\bf
r}_{\perp}\in\Sigma({\bf r}_{\perp},z)\\}-{\bf r}_{\perp}(z)\right],$ (10)

remains constant regardless of the value of $z$ and the shape acquired by
$\mathcal{C}({\bf r}_{\perp},z)$ along its propagation [18]. Note, on the
right-hand side of the second equality, that $E_{\Sigma({\bf r}_{\perp},z)}$
is directly related to the set of all trajectories laying within
$\mathcal{C}({\bf r}_{\perp},z)$ and with initially distributed within
$\mathcal{C}({\bf r}_{\perp},0)$ according to $|\psi({\bf r}_{\perp},0)|^{2}$.

Specifically, as it is shown below, here we use these properties to
quantitatively analyze some aspects associated with the leading peak of the
input finite-energy Airy beams.

### Finite-energy Airy beams

Let us consider the particular case of the one-dimensional Helmholtz equation,
although this will be done, for simplicity, in terms of the reduced
coordinates $\tilde{x}=x/x_{0}$ ($x_{0}$ denotes an arbitrary typical value
for the transverse coordinate) and $\tilde{z}=z/kx_{0}^{2}=\lambda_{0}z/2\pi
nx_{0}^{2}$ [17],

$i\ \frac{\partial\psi(x,z)}{\partial
z}=-\frac{1}{2}\frac{\partial^{2}\psi(x,z)}{\partial x^{2}}$ (11)

(to simplify notiation, from now on tildes will be removed). As Berry and
Balazs showed [1], for wave packets with an initial amplitude described by an
Airy function, propagating according to Eq. (11), their amplitude remains
shape invariant and undergoes a transverse displacement that increases rapidly
with $z^{2}$:

$\psi(x,z)=e^{i(x-z^{2}/6)z/2}Ai(x-z^{2}/4).$ (12)

Alternatively, the propagation of these beams can be described at a local
level in terms of trajectories associated with the transverse energy flux
[14],

$j(x,z)={\rm Im}\left[\psi^{*}(x,z)\frac{\partial\psi(x,z)}{\partial
x}\right].$ (13)

Specifically, the trajectories arise from the integration along $z$ of the
equation of “motion”

$\frac{dx(z)}{dz}=v(x,z)=\frac{j(x,z)}{I(x,z)},$ (14)

where $v(x,z)$ is the local field describing the energy flow along the
transverse direction ($x$) at a given value of the longitudinal coordinate
($z$), and $I(x,z)=|\psi(x,z)|^{2}=|Ai(x-z^{2}/4)|^{2}$ is the beam intensity
distribution.

Because the field (14) is directly related to the beam local phase variations
along the transverse direction, it becomes a suitable tool to explore any
property associated with the phase and, more specifically, with effects due to
possible variations induced on it (see, for instance, Refs. [15, 16]). Note
that for ideal Airy beams, $v(x,z)$ depends linearly on the $z$-coordinate and
hence the associated trajectories are fully analytical:

$x(z)=x(0)+\frac{z^{2}}{4}.$ (15)

These flux trajectories are parabolas, in agreement with the renowned
parabolic transverse acceleration undergone by the beam during its propagation
along $z$.

In the case of non-ideal Airy beams, that is, realistic beams with a finite
energy content, because the infinite tail of the ideal Airy beam is not
physically realizable in practice, the corresponding truncation introduces
extra phase factors that eventually lead to the gradual dispersion of the
beam. In other words, the beam will gradually lose both its shape invariance
and its transverse self-accelerated propagation beyond a certain spatial
range. In this regard, the above trajectory-based approach constitutes a
suitable tool to detect and explore any change directly related to the phase
at a local level.

The truncation can be introduced by directly acting on the tail of the beam,
as it was formerly proposed by Siviloglou and Christodoulides [9], who
considered an exponentially decaying factor acting on an input (ideal) Airy
beam, i.e.,

$\psi(x,0)=Ai(x)e^{\gamma x}.$ (16)

The propagated beam in this case can be analytically determined and reads as

$\psi(x,z)=e^{i(x-z^{2}/6)z/2+\gamma(y-z^{2}/4)-i\gamma^{2}z/2}Ai(y),$ (17)

where the complex argument $y=x-z^{2}/4+i\gamma z$ in the Airy function makes
the latter to become also complex, thus leading to the appearance of another
phase factor, apart from those arising from the multiplying exponential
prefactor. More specifically, the expression for the velocity field becomes
[15]

$v(x,z)=\frac{z}{2}+\frac{\partial}{\partial x}\left\\{{\rm
arg}\left[Ai(y)\right]\right\\},$ (18)

which leads to a non-analytical equation of motion, where we clearly observe
that the extra term (second term on the right-hand side) is responsible for
the deviation of the trajectories (flux) from the typical parabolic shape
associated with an ideal Airy beam.

In the above case, the field can be represented as a product of the ideal Airy
beam and an exponential. Hence, the corresponding transverse spectrum,

$\tilde{\psi}(k_{x})=e^{ik_{x}^{3}/3-\gamma(k_{x}-i\gamma/2)^{2}+\gamma^{3}/12},$
(19)

is the direct product of the transverse spectrum corresponding to the ideal
Airy beam,

$\tilde{\psi}(k_{x})=e^{ik_{x}^{3}/3},$ (20)

and a complex Gaussian function. Taking into account this result and that,
according to the convolution theorem, the Fourier transform of the convolution
of two functions is the direct product of the Fourier transforms of such
functions, we proceed now as follows to investigate the effects induced by
known aperture functions on the flux dynamics. Consider that the finite-energy
Airy beam results from the convolution function

$\psi(x,z)=\int f(x^{\prime})Ai(x-x^{\prime},z)dx^{\prime},$ (21)

where $w(x)$ represents a given aperture function. From a formal point of
view, here we are going to assume that $w(x)$ satisfies the following two
properties: ($i$) it is $L^{2}$ in order to warrant the finiteness of the
energy content of $\psi(x)$, and ($ii$) under some control parameter it
approaches a Dirac $\delta$-function, which allows us to recover the original
(ideal) infinite energy Airy beam.

Thus, according to the convolution theorem, the Fourier transform of (21)
reads as

$\tilde{\psi}(k_{x})=\tilde{f}(k_{x})e^{ik_{x}^{3}/3},$ (22)

from which we already infer that any additional phase effect is going to be
associated with the presence of the new term $\tilde{f}(k_{x})$, as the Airy
prefactor only governs the curvature of the flux. More importantly, in the
current case, note that the transverse power spectrum of the ideal Airy beam
is totally flat, so the parabolic curvature of the flux only comes from the
cubic phase factor of the Fourier transform. The appearance of
$\tilde{f}(k_{x})$ in (22) is not only going to affect the phase, as seen in
(18), but it may also induce a modulation of the spectral components, thus
removing the flatness of the power spectrum. Both factor are therefore
responsible for the eventual diversion of the propagation of the finite-energy
beam from the one displayed by an ideal Airy beam.

Taking into account the functional form of a propagated Airy beam, Eq. (11),
the amplitude (21) can be recast as

$\psi(x,z)=e^{i(x-z^{2}/6)z/2}\mathcal{A}i(x,z),$ (23)

where

$\mathcal{A}i(x,z)\equiv\int f(x^{\prime})Ai(x-x^{\prime}-z^{2}/4)\
e^{-izx^{\prime}/2}dx^{\prime}.$ (24)

The amplitude (23) leads to a generalized equation of motion that now reads as

$v(x,z)=\frac{z}{2}+\frac{\partial}{\partial x}\left\\{{\rm
arg}\left[\mathcal{A}i(x,z)\right]\right\\},$ (25)

containing the action of an extra term on the right-hand side. This extra term
is directly connected with the spectrum of the aperture functions, which
induces a new local phase configuration apart from the truncation of the beam,
thus leading to the distortion of the beam with respect to its ideal
counterpart.

So far the trajectories obtained after integrating Eq. (25) provide us with an
idea of the newly induced phase effects on the flux dynamics. However, we can
also extract valuable quantitative information regarding the finite-energy
beams if we consider particular ensembles of such trajectories and the above
properties satisfied by the trajectories arising from Eq. (25). Here, in
particular, we consider two quantities to investigate in a more quantitative
manner the dynamics associated with the leading maximum of the input finite-
energy Airy beans generated:

  * (i)

The escape rate, $N(z)$, defined as the partial energy content that remains
confined within two parabolic trajectories delimiting the main maximum of an
ideal Airy beam, given that the main maximum for the three input beams here
considered basically overlaps with that of the ideal Airy beam. This quantity
is computed according to Eq. (10) and the $\mathcal{C}(x,0)$ is determined by
the above mentioned limiting trajectories. In order to produce a fair sampling
in this single-event dynamics, the initial conditions are randomly chosen
according to the weight $|\psi(x,0)|^{2}$ (see Materials and methods for
additional technical details). Depending on the functional form of the
aperture function, the escape rate will fall faster or slower, thus indicating
either a faster or a slower deviation with respect to the ideal Airy beam.

  * (ii)

The average position, $x_{\rm av}(z)$, determined at each $z$-value as the
direct average taking into account all trajectories in the ensemble. Again,
like the escape rate, this quantity supplies local information about how fast
the main peak disappears, and hence how fast the finite-energy beam looses the
properties that characterize an ideal Airy beam.

### Aperture functions

Specifically, to explore the feasibility of aperture functions with different
functional forms and then to determine which one shows a better performance,
we have considered the following aperture functions:

$\displaystyle w_{G}(x)$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{\pi
w_{G}^{2}}}\ e^{-(x-x_{c,G})^{2}/w_{G}^{2}},$ (26) $\displaystyle w_{L}(x)$
$\displaystyle=$ $\displaystyle\frac{1}{\pi
w_{L}}\frac{1}{1+(x-x_{c,L})^{2}/w_{L}^{2}},$ (27) $\displaystyle w_{s}(x)$
$\displaystyle=$ $\displaystyle\frac{1}{\pi w_{s}}\ {\rm
sinc}\left[(x-x_{c,s})/w_{s}\right].$ (28)

In particular, for a comparison on the same footing conditions, we have
considered $x_{c}=0$ in the three cases, while $w$ has been chosen in a way
that the three functions have the same FWHM. In this latter regard, we have
fixed the width of the Gaussian aperture function to be $w_{G}=26.5$ $\mu$m
($0.5$ in reduced units, in terms of $x_{0}$), which renders a ${\rm
FWHM}=44.13$ $\mu$m and the following width values for the two other aperture
functions:

$\displaystyle w_{L}$ $\displaystyle=$ $\displaystyle w_{G}\sqrt{\ln 2}\approx
22.06\ \mu{\rm m},$ (29) $\displaystyle w_{s}$ $\displaystyle\approx$
$\displaystyle\frac{w_{G}\sqrt{\ln 2}}{u_{+}}\approx 11.64\ \mu{\rm m}.$ (30)

where

$u_{+}=\frac{\pi}{2}-4+\sqrt{16+4\pi-\pi^{2}}\approx 1.895$ (31)

is the (positive) solution for ${\rm sinc}\ u$ arising from the Bhāskara I’s
sine approximation formula [19],

${\rm sinc}\ u\approx\frac{16(\pi-u)}{5\pi^{2}-4u(\pi-u)},$ (32)

for the equation ${\rm sinc}\ u=0.5$, with $u=(x-x_{c,s})/w_{s}$.

Figure 1:  (a) Spatial dependence of the three aperture functions considered
here, all with the same FWHM, $44.13$ $\mu$m. The Gaussian aperture function
is represented with black dotted line, the Lorentzian aperture function with
red dashed line, and the sinc-type aperture function with blue solid line. For
a better comparison, the three functions are normalized to the unity. (b)
Transverse power spectrum, $|\tilde{\psi}(k_{x})|^{2}$, of the input
convoluted wave functions arising after applying the convolution with the
aperture functions displayed in part (a) [type and color lines are in
correspondence with the respective figures in (a)]. Again, for a better
comparison, the three spectra have been normalized to the unity at their
center.

The three aperture functions are shown in Fig. 1a, all with the same FWHM,
$44.13$ $\mu$m, which is nearly half the FWHM of the leading maximum of the
ideal (infinite energy) Airy beam considered, $96.13$ $\mu$m. For a better
comparison, the three functions are plotted without the constant prefactor, so
that their maxima coincide with 1, while they asymptotically approach 0 at
large $x$. We notice that, although the three cases present nearly the same
trend around $x_{c}$, the spatial extent of the sinc-type function is larger
than in the other two cases. Nonetheless, if we inspect the transverse power
spectrum, $|\tilde{\psi}(k_{x})|^{2}$, corresponding to the initial ansätze
constructed with these aperture functions (see Fig. 1b), we notice that the
three cases essentially span the same range of momenta, even though the
profile of the respective spectra looks very different (again, for a better
comparison, the three spectra have been normalized to the unity at their
center). Given that the power spectrum of an ideal Airy beam is flat, with an
infinite extension, we note that, while the Gaussian and Lorentzian aperture
functions are going to introduce a gradually decreasing of the amplitude
components involved in the Airy beam, the sinc-type aperture function, on the
contrary, is going to essentially preserve the same amplitude for the spectral
range spanned. In the current case, as it can be seen in Fig. 1b, this
spectral range is $\Delta k_{x}=2/w_{s}\approx 171.76$ mm-1, with all
components, low and large spatial frequencies, essentially contributing the
same. This is in sharp contrast with the other two aperture functions, where
the low spatial frequencies will play a more prominent role than the large
ones (more remarkable in the Lorentzian case than in the Gaussian one).

## Results

### Numerical dynamical simulations

To check the feasibility of the aperture functions here proposed, we have
conducted a series of numerical simulations due to the non-analiticity of the
propagation process in the three cases. As it is mentioned above, the
parameter values considered are analogous to those used in the original
experiments [4], which enables a more direct comparison with such experiments.
Accordingly, in Figs. 2, 3, and 4 we show both an overall view of the
propagation process in terms of trajectories superimposed on the corresponding
density plots of the intensity distribution (see part a in these figures), and
a series of snapshots of the intensity distribution at given values of $z$,
which allow us to elucidate in a qualitative manner how fast the properties of
ideal Airy beams are gradually lost (see part b in the figures).

Figure 2:  (a) Flux trajectories (white solid lines) associated with an Airy
beam convoluted with a Gaussian aperture function. The initial conditions
cover homogeneous the range spanned by the first foremost 10 maxima (including
the leading one). The thick black line denotes a trajectory with its initial
condition at the same position where the main maximum lies. (b) Profiles of
the intensity distribution at: $z=0$ (black solid line), $z=10$ cm (red solid
line), $z=20$ cm (blue solid line), and $z=30$ cm (green solid line). The
curves have been normalized to unity taking into account the maximum value of
the input beam (black solid line). In the inset, the same curves but each
normalized to unity with respect to their corresponding maxima.

In Fig. 2a we observe the propagation for a beam obtained with the Gaussian
aperture function. As it is seen from the density plot, the beam behaves
similarly to an ideal Airy beam for a relatively long distance, with a well-
defined main leading maximum surviving up to nearly the largest $z$-value
considered ($z=30$ cm). The snapshots displayed in Fig. 2b show clear evidence
of such a behavior, particularly in the inset in this figure, where we have
renormalized the three distribution for a better comparison. Nonetheless,
unlike an ideal Airy beam, there is a progressive decrease in the intensity as
well as a loss of the characteristics oscillations that form the tail of the
beam. To understand this behavior, we resort to the trajectory description
(see Fig. 2a). Accordingly, we note that, as it was mentioned in [15], the
limited extension of the tail releases the beam from the infinite energy
influx coming from the back. Consequently, not only the beam propagates
forward, but it can also get diffracted backwards, which is what the
trajectories (white solid lines) clearly show. Actually, we can observe that,
if we choose an initial condition at the position of the maximum of the main
peak of the initial beam, the corresponding trajectory (black solid line)
start deviating from a parabolic-type curve for distances beyond $z=15$ cm,
approximately. In other words, typical Airy beam properties are warranted for
about $15$ cm; afterwards, the beam starts getting distorted.

Figure 3:  (a) Flux trajectories (white solid lines) associated with an Airy
beam convoluted with a Lorentzian aperture function. The initial conditions
cover homogeneous the range spanned by the first foremost 10 maxima (including
the leading one). The thick black line denotes a trajectory with its initial
condition at the same position where the main maximum lies. (b) Profiles of
the intensity distribution at: $z=0$ (black solid line), $z=10$ cm (red solid
line), $z=20$ cm (blue solid line), and $z=30$ cm (green solid line). The
curves have been normalized to unity taking into account the maximum value of
the input beam (black solid line). In the inset, the same curves but each
normalized to unity with respect to their corresponding maxima.

The results for the Lorentzian aperture function are shown in Fig. 3. In this
case, we observe a faster degradation, at shorter distances (around $10$ cm).
In Fig. 3a the density plot reveals a large transfer of flux from the main
peak to the rear secondary maxima, although with the seemingly net effect of
getting stacked around $x=0$ mm. If we inspect Fig. 3b, effectively we clearly
notice an accumulation of intensity around $x=0$ mm (see also inset),
decreasing towards larger values of $x$. From the trajectory description, we
find an overall picture that resembles that of the Gaussian aperture, but with
the subtle difference that now the trajectory starting at the position of the
main maximum remains nearly around the same position, without following the
foremost ensemble of trajectories that still show a parabolic-type shape.

Figure 4:  (a) Flux trajectories (white solid lines) associated with an Airy
beam convoluted with a sinc-type aperture function. The initial conditions
cover homogeneous the range spanned by the first foremost 10 maxima (including
the leading one). The thick black line denotes a trajectory with its initial
condition at the same position where the main maximum lies. (b) Profiles of
the intensity distribution at: $z=0$ (black solid line), $z=10$ cm (red solid
line), $z=20$ cm (blue solid line), and $z=30$ cm (green solid line). The
curves have been normalized to unity taking into account the maximum value of
the beam at $z=10$ cm (red solid line). In the inset, the same curves but each
normalized to unity with respect to their corresponding maxima.

The two previous cases can be understood in the light of the strong filtering
produced by the aperture functions considered in each case, which remove
dramatically the high frequencies from the spectrum (see Fig. 1a). In sharp
contrast, the case of the sinc-type aperture function, represented in Fig. 4,
shows a very different picture. First, except from some intensity bursts
propagating backwards, the density plot shows that all maxima, both the main
peak and the secondary maxima, exhibit the typical parabolic shape expected
for an Airy beam. The results in Fig. 4b clearly indicate that, effectively,
the beam behaves like an ideal Airy beam for longer distances. On the other
hand, the behavior of the trajectories also reveals that there is a larger
survival of the parabolic shape for longer distances. Unlike the two previous
cases, note that here the trajectories show a back and forth displacement,
which keeps confined along parabolic-like channels for longer distances, until
the backward energy influx reduces and they can escape (backwards). This
behavior makes that a large portion of the foremost trajectories will keep
showing their parabolic shape for rather long $z$-distances. Indeed, note that
the trajectory started at the position of maximum of the main peak is nearly
parabolic for the full propagation.

### Dynamics of the main maximum

Apart from a local description of the propagation process, the trajectories
also render interesting information that cannot be directly obtained from the
density plots. In particular, to further understand physically the performance
of the three aperture functions, we can investigate the local dynamics
associated with a piece of the beam without perturbing its propagation. As it
is inferred from [20], the intensity content within two any initial conditions
will remain constant along $z$ provided we propagate those conditions along
the corresponding trajectories. That is, if we consider two initial conditions
delimiting a given intensity maximum (for instance, setting such initial
conditions near the adjacent nodes, but not on top of them), we can determine
at any subsequent $z$-value how the corresponding intensity distributes and
within which boundaries.

Figure 5:  (a) Escape rate, $N(z)$, as a function of $z$ for the three
aperture functions here considered. (b) Average position, $x_{\rm av}(z)$,
over a sampling of 1,000 flux trajectories randomly distributed, at $z=0$,
within the area covered by the leading maximum for the three aperture
functions. The distribution of initial conditions have been taken according to
a biGaussian distribution (See Materials and methods section). To compare
with, the $z$-dependence of the leading maximum of an ideal (infinite energy)
Airy beam is also shown, represented with green dashed line.

Accordingly, we have considered pairs of initial conditions between which we
collect most of the intensity associated with the main peaks in the three
cases (the ranges are given in Table 1 in Materials and methods). Then, in
order to understand how fast the corresponding beams deviate from an ideal
Airy beam, we have computed two quantities. On the one hand, we have computed
the escape rate, $N(z)$, which measures, along $z$, the intensity confined
within two boundaries given by two parabolic trajectories. It is clear that
the better the preservation of the ideal Airy beam properties the larger the
value of $N(z)$ for larger and larger $z$-values. This quantity is plotted in
Fig. 5a. As it can be seen, the decay of $N(z)$ is relatively fast for the
Gaussian and the Lorentzian aperture functions, although the falloff for the
latter is still slightly larger. While the energy (proportional to the number
of trajectories) decays to about a half at $z\approx 18.8$ cm for the Gaussian
aperture function, in the Lorentizan case the same takes place at $z\approx
13.5$ cm. This corresponds to the degradation of the intensity distribution
profile observed in Figs. 2a and 3a, respectively. On the contrary, for the
sinc-type aperture function, $N(z)$ remains nearly constant until $z\approx
19.8$ cm and, extrapolating data, a decrease to about $50\%$ is only expected
for larger distances ($z\approx 33.5$ cm), which is in correspondence with the
behavior observed in Fig. 4a.

On the other hand, taking into account that those initial conditions have been
randomly distributed according to the intensity distribution, their average,
$x_{av}$, can provide an objective measure of the self-accelerated propagation
of only this beam. Again, the better the preservation of the ideal behavior,
the closer this average will be to a parabola. This quantity is represented in
Fig. 5b (see solid lines), where we have also included the trajectories
associated with the initial condition taken at the position of the maximum of
the respective main peaks (dashed lines). To compare with, we have also
included the trajectory of an ideal Airy beam with its initial position at the
$x$-value where the main peak of such a beam reaches its maximum (see green
dotted line). As we can see, the trajectories corresponding to the maximum in
the Gaussian and Lorentzian cases start deviating from the average value at
short and medium distances, respectively, in agreement with what we already
saw in Figs. 2 and 3. More specifically, such a deviation starts being
remarkable at about $20$ cm for the Gaussian aperture function, and about $10$
cm for the Lorentzian one. These distances can be used as distinctive
objective criterion to determine when the performance of the corresponding
aperture function becomes poor. In therms of the average, $x_{av}$, this loss
of efficiency is indicated by the transition from a nearly parabolic behavior,
at the initial stages of the propagation, to a seemingly linear behavior at
larger distances, which is produced by the gradual drain of trajectories
(flux) backwards.

In the case of the sinc-type aperture function, again we observe a very
different behavior, since both the average and the trajectory of the maximum
basically propagate parallel one another. More interestingly, these two
quantities overlap in a good approximation the trajectory of the maximum for
an ideal Airy beam; only at large $z$-values we start observing a discrepancy,
which is a signature of the loss of the ideal behavior.

## Discussion

We have introduced a trajectory-based methodology to investigate the effect of
limiting the energy content of Airy beams by means of aperture functions with
different shapes, and hence to elucidate which functional form better
preserves the properties of an ideal Airy beam at the input. Thus, we have
combined conventional propagation simulations, to determine the intensity
distribution at difference distances from the input plane, with the above
mentioned non-conventional approach based on an on-the-fly computation of the
transverse energy flux, which provides us with information at a more local
level. The transverse energy flux, described in terms of sets of trajectories,
enables a description and analysis of the propagation process on a single-
event basis, making more intuitive the properties of Airy beams (shape
invariance and self-accelerated propagation) at the same time that it is also
possible to perform some additional quantitative analyses of specific
features.

From our analysis, we conclude that, among the three types of aperture
functions investigated, sinc-type ones show a better performance, because they
basically introduce a cut in the beam spectrum, without modifying the
amplitude of its components, unlike other types of aperture functions, such as
Gaussian or Lorentzian profiles. The trajectories show a very interesting
picture in this regard. On the one hand, we note that the trajectories
associated with the Gaussian or the Lorentzian aperture functions make evident
a monotonous energy dissipation backwards. This continuous flow produces a
rather fast lost of fringe visibility of the intensity distribution, as it is
manifested by the corresponding escape rates, and hence the also fast loss of
the typical properties of an Airy beam at relatively close distances from the
input plane. In particular, it has been noticed that the self-accelerated
displacement that should undergo the average position of the energy
accumulated in the main maximum is gradually lost.

On the other, in the case of the sinc-type aperture function, the fact that
all spatial spectral components play have nearly the same amplitude provokes
that the presence of a wiggling displacement of the trajectories and hence a
back-and-forth energy flow between adjacent intensity maxima. This oscillatory
behavior enables the preservation of the typical Airy beam features for much
longer distances, even though the gradual degradation observed in the pattern
along them. Furthermore, unlike the two previous cases, here it has been shown
that the average position of the energy accumulated in the main maximum
follows very nicely the behavior expected for an ideal Airy beam.

In sum, the present analysis is aimed at providing new insights into the
behavior of finite-energy Airy beams, and therefore also expected to
contribute to the design and development of novel applications exploiting
finite-energy Airy beams. Moreover, it is also worth noting that, by virtue of
the tight analogy between the paraxial Helmholtz equation and the Schrödinger
equation, the results shown and discussed here are also directly transferable
to the generation and exploitation of Airy beams with nonzero mass particles,
as it is the case of electron Airy beams, first produced in 2013 by Voloch-
Bloch et al. [21].

## Materials and methods

The simulations conducted in this work to obtain the beam propagation and the
corresponding trajectories, have been carried out following the same procedure
formerly used in [15], based on the split-step for beam propagation and a
spectral recast of the equation of motion for the trajectories [22]. For the
trajectories shown in Figs. 2, 3, and 4, we have considered an evenly-spaced
distribution of 50 initial conditions, taken from a range that covers the
first 10 intensity maxima of the input beam (the leading maximum and the nine
subsequent secondary maxima). Concerning the random distribution required to
compute the quantities shown in Fig. 5, first we have fitted the main maximum
of the input intensity distributions (see black line curves in Figs. 2, 3, and
4) to the biGaussian distribution,

$f(x)=\left\\{\begin{array}[]{lcc}\exp[-(x-x_{0})^{2}/2\sigma_{l}^{2}],&&x<x_{0}\\\
\exp[-(x-x_{0})^{2}/2\sigma_{r}^{2}],&&x\geq x_{0}\\\ \end{array}\right.,$
(33)

which captures fairly well the asymmetry of such a maximum. The values for the
parameters $x_{0}$, $\sigma_{l}$, and $\sigma_{r}$ obtained from the
corresponding fitting procedures for the three cases are given in Table 1.

Table 1:  Parameters used in the fitting of the main maximum to a biGaussian distribution. aperture function |  $x_{c}$ (mm) |  $w_{1}$ (mm) |  $w_{2}$ (mm) | i.c. range (mm)
---|---|---|---|---
Gaussian | $-0.05688$ | $0.02647$ | $0.04844$ | $(-0.123,0.115)$
Lorentzian | $-0.05363$ | $0.02887$ | $0.05412$ | $(-0.128,0.161)$
sinc-type | $-0.06016$ | $0.02341$ | $0.04536$ | $(-0.124,0.112)$

Regarding the specific random choice, first, we select randomly from a $[0,1)$
uniform distribution a number. If this number is smaller than the ratio
$\sigma_{l}/(\sigma_{l}+\sigma_{r})$, then we associate the initial condition
to the left side of the biGaussian distribution (i.e., $x<x_{0}$); otherwise,
the initial condition is associated with the right side. A number is then
picked up randomly following a Gaussian distribution centered at $x_{0}$ and
with a width $\sigma_{l}$, if the previous choice laid on the left side of the
distribution, or $\sigma_{r}$, if the choice led to the right side. It might
happen that the number obtained from this procedure will lay on the opposite
side where it is supposed to be; in such a case, the algorithm re-start the
choice of the side and repeats the process.

## Data Availability

Data underlying the results presented in this paper are not publicly
available, but may be obtained from the authors upon reasonable request.

## Funding

Grants PID2019-104268GB-C21, PID2021-127781NB-I00, and network
RED2022-134391-T “Foundations and Applications of Structured Light” funded by
the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033.

## Disclosures

The authors declare no conflicts of interest.

## References

  * [1] M. V. Berry and N. L. Balazs, “Nonspreading wave packets,” Am. J. Phys. 47, 264–267 (1979).
  * [2] D. M. Greenberger, “Comment on “nonspreading wave packets”,” Am. J. Phys. 48, 256–256 (1980).
  * [3] K. Unnikrishnan and A. R. P. Rau, “Uniqueness of the airy packet in quantum mechanics,” Am. J. Phys. 64, 1034–1035 (1996).
  * [4] G. A. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of accelerating airy beams,” Phys. Rev. Lett. 99, 213901 (2007).
  * [5] A. Forbes, M. de Oliveira, and M. R. Dennis, “Structured light,” Nat. Photonics 15, 253–262 (2021).
  * [6] J. Wang and Y. Liang, “Generation and detection of structured light: A review,” Front. Phys. 9, 688284 (2021).
  * [7] Y. Zhang, H. Zhong, M. R. Belić, and Y. Zhang, “Guided self-accelerating airy beams—a mini-review,” Appl. Sci. 7 (2017).
  * [8] N. K. Efremidis, Z. Chen, M. Segev, and D. N. Christodoulides, “Airy beams and accelerating waves: an overview of recent advances,” Optica 6, 686–701 (2019).
  * [9] G. A. Siviloglou and D. N. Christodoulides, “Accelerating finite energy airy beams,” Opt. Lett. 32, 979–981 (2007).
  * [10] R. Martínez-Herrero and A. S. Sanz, “Partially coherent airy beams: A cross-spectral-density approach,” Phys. Rev. A 106, 053512 (2022).
  * [11] D. Xu, Y. Liu, Z. Mo, J. Jiang, J. Shi, Z. Liang, Y. Wu, J. Zhao, H. Yang, H. Huang, H. Liu, L. Shui, and D. Deng, “Shaping autofocusing airy beams through the modification of fourier spectrum,” Opt. Express 30, 232–242 (2022).
  * [12] H. Margenau and G. M. Murphy, _The Mathematics of Physics and Chemistry_ (D. Van Nostrand, Princeton, NJ, 1943).
  * [13] A. S. Sanz, J. Campos-Martínez, and S. Miret-Artés, “Transmission properties in waveguides: An optical streamline analysis,” J. Opt. Am. Soc. A 29, 695–701 (2012).
  * [14] A. S. Sanz, M. Davidović, and M. Božić, “Bohmian-based approach to gauss-maxwell beams,” Appl. Sci. 10, 1808(1–25) (2020).
  * [15] A. S. Sanz, “Flux trajectory analysis of airy-type beams,” J. Opt. Soc. Am. A 39, C79–C85 (2022).
  * [16] A. S. Sanz, L. L. Sanchez-Soto, and A. Aiello, “A quantum trajectory analysis of singular wave functions,” arXiv:2301.13207v1 (2023).
  * [17] A. S. Sanz and S. Miret-Artés, _A Trajectory Description of Quantum Processes. I. Fundamentals_ , vol. 850 of _Lecture Notes in Physics_ (Springer, Berlin, 2012).
  * [18] A. S. Sanz and S. Miret-Artés, “On the unique mapping relationship between initial and final quantum states,” Ann. Phys. 339, 11–21 (2013).
  * [19] A. García-Sánchez and A. S. Sanz, “Analysis of the gradual transition from the near to the far field in single-slit diffraction,” Phys. Scr. 97, 055507 (2022).
  * [20] A. S. Sanz and S. Miret-Artés, “Setting-up tunneling conditions by means of bohmian mechanics,” J. Phys. A: Math. Theor. 44, 485301(1–17) (2011).
  * [21] N. Voloch-Bloch, Y. Lereah, Y. Lilach, A. Gover, and A. Arie, “Generation of electron airy beams,” Nature 494, 331–335 (2013).
  * [22] A. S. Sanz and S. Miret-Artés, _A Trajectory Description of Quantum Processes. II. Applications_ , vol. 831 of _Lecture Notes in Physics_ (Springer, Berlin, 2014).


# Anatomically informed Bayesian spatial priors for fMRI analysis

###### Abstract

Existing Bayesian spatial priors for functional magnetic resonance imaging
(fMRI) data correspond to stationary isotropic smoothing filters that may
oversmooth at anatomical boundaries. We propose two anatomically informed
Bayesian spatial models for fMRI data with local smoothing in each voxel based
on a tensor field estimated from a $T_{1}$-weighted anatomical image. We show
that our anatomically informed Bayesian spatial models results in posterior
probability maps that follow the anatomical structure.

Index Terms—  Bayesian statistics, functional MRI, activation mapping,
adaptive smoothing

##  1 Introduction

The analysis of functional magnetic resonance imaging (fMRI) data has
generally relied on frequentist statistics to perform inferences about brain
activity. After fitting a general linear model (GLM) to the individual voxel
time series, t-values are calculated for contrasts of interest at every voxel,
resulting in t-maps that can be thresholded at an appropriate significance
level. Isotropic Gaussian smoothing is the most common way of preprocessing
fMRI data, but several adaptive smoothing approaches for detecting brain
activity have been proposed [1, 2, 3, 4, 5, 6, 7].

During the 2000s, a Bayesian framework for fMRI analysis was developed by
Penny et al. [8, 9, 10, 11, 12], which provided increased flexibility compared
to the classical frequentist approach, by allowing the estimation of
individual smoothness parameters for each regressor and autoregressive (AR)
noise coefficient. Sidén et al. [13, 14] extended this work to an efficient
Markov Chain Monte Carlo (MCMC) implementation which works in 2D as well as
3D. However, we are not aware of any work that performs anatomically adaptive
Bayesian spatial modeling. In this work we propose two approaches for
performing 2D Bayesian spatial modeling in an anatomically-adaptive way.

##  2 Methods

###  2.1 Bayesian spatial GLM

The most common way to model fMRI data is as a voxel-wise general linear model
with serial correlations in the residuals, which are modeled as an AR($p$)
process. This combined model is referred to as GLM-AR($p$). In the case of
single subject data, with $T$ volumes, $N$ voxels, $K$ regressors and an AR
model of order $p$, this model is written as

$\underset{[T\times N]}{\mathbf{Y}}=\underset{[T\times
K]}{\mathbf{X}}\underset{[K\times N]}{\mathbf{W}}+\underset{[T\times
N]}{\mathbf{E}},$ $\underset{[T\times N]}{\mathbf{E}}=\underset{[T\times
P]}{\widetilde{\mathbf{E}}}\underset{[P\times
N]}{\mathbf{R}}+\underset{[T\times N]}{\mathbf{Z}},$

where $\mathbf{Y}$ is the observation matrix, $\mathbf{X}$ is the design
matrix, $\mathbf{W}$ is the regressor matrix and $\mathbf{E}$ is the residual
matrix. The residual matrix $\mathbf{E}$ is itself expressed as the product of
a lagged prediction error matrix $\widetilde{\mathbf{E}}$ with an AR
coefficient matrix $\mathbf{R}$ plus a matrix of i.i.d. zero mean Gaussian
errors $\mathbf{Z}$, where $\mathbf{Z}_{\cdot,n}\stackrel{{\scriptstyle
iid}}{{\sim}}N(\mathbf{0},\lambda_{n}^{-1}\mathbf{I}_{T})$, and $\lambda_{n}$
is the noise precision at voxel $n$. In this Bayesian framework, smoothness of
the regression coefficients enters the model through a spatial precision
matrix $\mathbf{D}_{w}$ in their prior distribution, according to

$\mathbf{W}_{k,\cdot}^{\prime}\sim\mathcal{N}\left(\mathbf{0},\alpha_{k}^{-1}\mathbf{D}_{w}^{-1}\right),$

where $\mathbf{W}^{\prime}_{k,\cdot}$ is the transpose of the $k$-th row of
the regression coefficient matrix $\mathbf{W}$, representing the $k$-th
regression coefficient at every voxel, and $\alpha_{k}$ is a smoothness
hyperparameter for the $k$-th regressor (since the following also applies to
the AR coefficients, we refer to their precision matrices generically as
$\mathbf{D}$). The precision matrix $\mathbf{D}$ encodes the conditional
dependencies between every pair of voxels. The smoothness assumption
constrains these dependencies to exist only between a voxel and its immediate
neighbors, which has the advantage of making this matrix very sparse.

###  2.2 Uniform graph Laplacian precision matrix

Graphs constitute a natural way of describing relationships between sets of
elements. Therefore, the precision matrices employed by Penny et al. and Sidén
et al. take the form of graph Laplacian matrices. In their formulation the
precision matrix $\mathbf{D}$ is an unweighted graph Laplacian (UGL), where
$d_{i,i}$ equals the number of pixels neighboring pixel $i$ and ${d_{i,j}=-1}$
if pixels $i,j$ are cardinal neighbors. For inner pixels this corresponds to
the Laplacian operator, which in 2D takes the form

$\mathbf{H}_{\text{UGL}}=\begin{pmatrix}[r]0&-1&0\\\ -1&4&-1\\\
0&-1&0\end{pmatrix}.$

Generically, graph Laplacian matrices can be constructed with [15]

$\mathbf{L}=\mathbf{B}-\mathbf{A},$

where $\mathbf{A}$ is an _adjacency matrix_ whose $a_{i,j}$ element, called a
_weight_ , represents in this context the strength of the conditional
dependence between pixels $i$ and $j$, and $\mathbf{B}$ is a diagonal _degree
matrix_ with $b_{i,i}=\sum_{j}a_{i,j}$, that is, the sum of all the incoming
weights to pixel $i$. Since $\mathbf{B}$ can be constructed from $\mathbf{A}$,
the adjacency matrix is sufficient for generating the graph Laplacian. Under
this formulation the same UGL prior can be constructed by considering its
dependency structure (neighborhood), which for inner pixels takes the shape

$\mathbf{N}_{\text{UGL}}=\begin{pmatrix}[r]0&1&0\\\ 1&0&1\\\
0&1&0\end{pmatrix},$

that is, $a_{i,j}=1$ if pixels $i,j$ are cardinal neighbors. It can be seen
that this results in the same precision matrix $\mathbf{D}\equiv\mathbf{L}$.

Due to the uniform neighborhood shape and constant weights between neighbors
this type of prior is incapable of encoding anatomical information. Such a
prior will not respect anatomical tissue boundaries, mixing, for example,
signals from white and gray brain matter. Conversely, in order for the prior
to encode relevant anatomical information it is necessary for the dependencies
between neighboring pixels to be specified independently at each location.

###  2.3 The structure tensor

In order to establish pixel relationships that respect anatomical boundaries,
we must first estimate the position and orientation of these boundaries. Our
main tool for achieving this is the structure tensor, a tensor estimated at
every pixel (using quadrature filters along different directions [16, 17])
whose eigenvalues and eigenvectors reveal the degree to which spatial
structure is present and the orientation of any such structure, here
understood as lines or edges. A detailed explanation of the estimation and
interpretation of structure tensors falls outside the scope of this paper, see
[18, 16, 19] for a thorough treatment of this topic. For our purposes it is
sufficient to point out that the eigenvector of the structure tensor
corresponding to the largest eigenvalue is perpendicular to the main structure
orientation at the given pixel, while the eigenvalue is proportional to the
degree of certainty of this orientation. The second eigenvalue indicates the
amount of structure in the orientation perpendicular to the first eigenvector.
The structure tensor is, therefore, a flexible model capable or representing
structure in zero (noise, uniform data), one (line, edge), or two (crossing
lines, crossing edges) orientations (for the $2$D case).

As fMRI data lacks contrast, the structure tensor is calculated from a
registered $T_{1}$-weighted volume. Figure 1 shows the eigenvectors of the
structure tensor for a typical brain image. The tensor is small and isotropic
in very uniform regions, while next to lines or edges it becomes anisotropic
and oriented perpendicularly to the line or edge.

Fig. 1: Local orientation represented using the structure tensor. Red vectors
indicate the main local orientation, while green vectors indicate the second
(perpendicular) orientation. If no vector is present, there is no orientation
information available (e.g. due to uniform intensity).

###  2.4 Spatial model with four orientations

In our first proposed solution we assign to every pixel one of four oriented
neighborhood structures according to the local structural orientation in that
pixel. Such a model offers limited angular resolution, since only four
different orientations can be represented. However, it simplifies the
formulation of the possible neighborhoods and provides increased sparsity in
the precision matrix $\mathbf{D}$. We refer to this model as 4DIR. We start by
considering four possible orientations: horizontal, vertical, and two
diagonals. The corresponding orientation vectors are

$\mathbf{d}_{x}=\begin{pmatrix}1\\\
0\end{pmatrix},\quad\mathbf{d}_{y}=\begin{pmatrix}0\\\ 1\end{pmatrix},$
$\mathbf{d}_{xy}=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
1\end{pmatrix},\quad\mathbf{d}_{-xy}=\frac{1}{\sqrt{2}}\begin{pmatrix}-1\\\
1\end{pmatrix}.$

For computational efficiency we want to avoid having to calculate the
eigenvectors of the structure tensor. We can find which of the four
orientations is closest to the main tensor orientation at each point by first
defining tensors corresponding to each of the four orientations

$\mathbf{T}_{x}=\mathbf{d}_{x}\mathbf{d}_{x}^{T},\quad\mathbf{T}_{y}=\mathbf{d}_{y}\mathbf{d}_{y}^{T},$
$\mathbf{T}_{xy}=\mathbf{d}_{xy}\mathbf{d}_{xy}^{T},\quad\mathbf{T}_{-xy}=\mathbf{d}_{-xy}\mathbf{d}_{-xy}^{T},$

and then projecting, through an inner product, the structure tensor at each
point onto the four orientation tensors. The maximum projection value will
correspond to the orientation closest to that of the structure tensors.

Having found at each pixel which of the four orientations is closest to that
of the structure tensor, we define neighborhoods (filters) along each of the
four orientations

$\mathbf{N}_{x}=\begin{pmatrix}[r]0&0&0\\\ 1&0&1\\\
0&0&0\end{pmatrix},\quad\mathbf{N}_{y}=\begin{pmatrix}[r]0&1&0\\\ 0&0&0\\\
0&1&0\end{pmatrix},$ $\mathbf{N}_{xy}=\begin{pmatrix}[r]1&0&0\\\ 0&0&0\\\
0&0&1\end{pmatrix},\quad\mathbf{N}_{-xy}=\begin{pmatrix}[r]0&0&1\\\ 0&0&0\\\
1&0&0\end{pmatrix}.$

As the structure tensor is aligned across the orientation of lines and edges
(see Figure 1), which we are trying to preserve, we associate to each pixel a
neighborhood perpendicular to the orientation of the structure at that point
(e.g. if $\mathbf{d}_{x}$ is the orientation for some pixel, then the
appropriate neighborhood would be $\mathbf{N}_{y}$).

Having assigned a specific neighborhood to each pixel, a new adjacency matrix
$\mathbf{A}_{\text{4DIR}}$ can be constructed and used to define a graph
Laplacian $\mathbf{D}_{\text{4DIR}}$. However, due to the specific process for
calculating the orientation at each position, it cannot be guaranteed that the
relationship between a pair of pixels $i,j$ is symmetric, i.e., that
$a_{i,j}=a_{j,i}$, which would preclude $\mathbf{D}_{\text{4DIR}}$ from being
used as a precision matrix. Symmetry can be enforced by applying the procedure

$\mathbf{A_{\text{4DIR}}}^{\prime}=\frac{\mathbf{A_{\text{4DIR}}}+\mathbf{A_{\text{4DIR}}}^{T}}{2},$

which effectively sets the value of $a_{i,j}$ and $a_{j,i}$ to the average of
the two. This correction has little effect on the orientations encoded in the
graph and allows the new matrix $\mathbf{L}_{\text{4DIR}}^{\prime}$ to be used
as a precision matrix $\mathbf{D}_{\text{4DIR}}$.

###  2.5 Spatial model with arbitrary orientations

The previously presented model is limited in angular resolution, as it can
only encode four orientations. In order to overcome this limitation we
consider a $3\times 3$ neighborhood around every pixel and determine the
weight for each of the $8$ possible neighbors by sampling from a continuous
function at discrete positions corresponding to the centers of the neighboring
pixels. We refer to this model as ANYDIR. The sampled function is

$a_{i,j}=\mathbf{f}(i,j)=\frac{\left|\sin(\phi_{\mathrm{pix}_{j}}-\phi_{\mathrm{tensor}_{i}})\right|^{\alpha}}{r_{\mathrm{pix}_{j}}^{\beta}},\quad\alpha,\beta>0,$

where $\phi_{\mathrm{pix}_{j}}$ is the angle of the line connecting the
central pixel $i$ with neighboring pixel $j$, $\phi_{\mathrm{tensor}_{i}}$ is
an angle representing the main orientation of the structure tensor at pixel
$i$, $r_{\mathrm{pix}_{j}}$ is the distance between pixels $i$ and $j$ (see
Figure 2). Additionally, $\alpha$ and $\beta$ are adjustable parameters taking
non-negative real values, where $\alpha$ controls the width of the
distribution around $\phi_{\mathrm{tensor}_{i}}$ and $\beta$ penalizes the
values in the diagonal neighbors with respect to the horizontal and vertical
ones. In this work we use $\alpha=12$ and $\beta=5$.

Fig. 2: Variables involved in weighting of the neighborhoods for the ANYDIR
model. The structure tensor is shown in the middle.

This neighborhood formulation generates neighborhoods that vary continuously
with respect to $\phi_{\mathrm{tensor}_{i}}$, which allows the representation
of arbitrary orientations in the structures around pixels, at the cost of
reduced sparsity in the precision matrix. As in the previous case the
adjacency matrix $\mathbf{A}_{\text{ANYDIR}}$ is not guaranteed to be
symmetric. The same kind of correction can be applied in order to generate a
graph Laplacian matrix $\mathbf{L}^{\prime}_{\text{ANYDIR}}$ suitable for use
as a precision matrix $\mathbf{D}_{\text{ANYDIR}}$.

##  3 Results

We compared the results produced by all three precision matrices in the
context of Bayesian fMRI analysis. We used preprocessed data from subject
$100307$ of the Human Connectome Project [20], specifically the
$T_{1}$-weighted volume and the motor task data. The analysis was performed
slice-by-slice using the SPM package for Matlab together with the extension
developed by Sidén et al. [13]. The only modification to the code was the
addition of the proposed precision matrices. Our modified code can be found at
(https://github.com/DavidAbramian/adaptiveBayesianPrior).

As a preprocessing step, the $T_{1}$-weighted volume was downsampled to match
the resolution of the fMRI data. This is necessary, since the precision matrix
has to be estimated from the $T_{1}$-weighted volume and later applied to the
fMRI data.

Due to the time-consuming nature of the MCMC analysis it was carried out on a
single representative axial slice. The Gibbs sampling algorithm was iterated
$10,000$ times, with $1,000$ warmup iterations, and with a thinning factor of
$5$.

###  3.1 Anatomical adaptiveness

Figure 3 illustrates the orientation of the pixel neighborhoods generated by
the two proposed approaches, which determine the conditional dependencies
encoded in prior precision matrix $\mathbf{D}$. In both cases the models have
successfully adapted to the anatomical structure given by the $T_{1}$-weighted
volume, as spatial prior dependence is placed along lines and edges and not
across them. The 4DIR method shows limited angular resolution, as it can only
represent neighborhoods in four orientations, while the ANYDIR produces
neighborhoods in arbitrary orientations.

Fig. 3: Comparison of neighborhood structures implied by both of the proposed
models. Lines represent the orientation of spatial prior dependencies at the
given point. Both methods adapt these dependencies at each point in accordance
with the anatomical structure. The 4DIR method only makes use of four possible
orientations, while with ANYDIR the orientations vary continuously.

###  3.2 Functional MRI results

Figure 4 shows, for all three models, the regression coefficients obtained for
a right hand motor task, as well as posterior probability maps (PPMs)
quantifying the probability of the effect of said task exceeding $0.2\%$ of
the global mean signal and thresholded at $0.8$. The regression coefficients
for both of the proposed models clearly reflect anatomical spatial patterns
absent from the UGL results. The patterns are highly angular for the 4DIR
model as a result of the limited angular resolution, while the ANYDIR model
results in more natural curved patterns.

The PPMs from all three methods are similar, showing large activation in the
motor cortex and close to the central sulcus, and smaller activations in the
somatosensory cortex. However, the activations detected by both of the
proposed methods are slightly narrower and extend further along anatomical
lines, indicating that the priors respect anatomical boundaries.

Fig. 4: Bayesian GLM regression results obtained using MCMC. Left: regression
coefficient corresponding to right hand motor task. Right: PPMs for the
probability of effect of the right hand motor task exceeding $0.2\%$ of the
global mean signal, thresholded at $0.8$. Top: UGL model. Middle: 4DIR model.
Bottom: ANYDIR model.

##  4 Discussion

We have proposed two new Bayesian spatial priors for fMRI analysis that allow
for a locally anisotropic spatial dependence over voxels. These priors were
used to encode anatomical structure, resulting in anatomically-adaptive
smoothing for fMRI data. The priors can be easily incorporated into the
existing framework for Bayesian fMRI analysis, requiring minimal
modifications.

While our presentation is centered on 2D priors, both of the proposed
approaches can be extended to 3D without requiring significant modifications.
The priors can also be improved by incorporating additional information from
the structure tensor to, for example, use the UGL prior in areas with uniform
intensity. These adaptive spatial models can also be applied in Bayesian
frameworks for diffusion MRI data [21].

## References

  * [1] Ola Friman, Magnus Borga, Peter Lundberg, and Hans Knutsson,  “Adaptive analysis of fMRI data,”  NeuroImage, vol. 19, no. 3, pp. 837–845, 2003.
  * [2] Rajesh R Nandy and Dietmar Cordes,  “Novel nonparametric approach to canonical correlation analysis with applications to low CNR functional MRI data,”  Magnetic Resonance in Medicine, vol. 50, no. 2, pp. 354–365, 2003\.
  * [3] Joakim Rydell, Hans Knutsson, and Magnus Borga,  “Bilateral filtering of fMRI data,”  IEEE Journal of Selected Topics in Signal Processing, vol. 2, no. 6, pp. 891–896, 2008.
  * [4] Anders Eklund, Mats Andersson, and Hans Knutsson,  “Fast random permutation tests enable objective evaluation of methods for single-subject fMRI analysis,”  International journal of biomedical imaging, vol. 2011, 2011.
  * [5] Hamid Behjat, Nora Leonardi, Leif Sörnmo, and Dimitri Van De Ville,  “Anatomically-adapted graph wavelets for improved group-level fMRI activation mapping,”  NeuroImage, vol. 123, pp. 185–199, 2015.
  * [6] Xiaowei Zhuang, Zhengshi Yang, Tim Curran, Richard Byrd, Rajesh Nandy, and Dietmar Cordes,  “A family of locally constrained CCA models for detecting activation patterns in fMRI,”  NeuroImage, vol. 149, pp. 63–84, 2017.
  * [7] Gabriele Lohmann, Johannes Stelzer, Eric Lacosse, Vinod J Kumar, Karsten Mueller, Esther Kuehn, Wolfgang Grodd, and Klaus Scheffler,  “LISA improves statistical analysis for fMRI,”  Nature communications, vol. 9, no. 1, pp. 4014, 2018.
  * [8] Will Penny, Stefan Kiebel, and Karl Friston,  “Variational Bayesian inference for fMRI time series,”  NeuroImage, vol. 19, no. 3, pp. 727–741, 2003.
  * [9] Will Penny and Guillaume Flandin,  “Bayesian analysis of fMRI data with spatial priors,”  in Proceedings of the Joint Statistical Meeting (JSM). American Statistical Association. Citeseer, 2005.
  * [10] Will Penny, Guillaume Flandin, and Nelson Trujillo-Barreto,  “Bayesian comparison of spatially regularised general linear models,”  Human brain mapping, vol. 28, no. 4, pp. 275–293, 2007.
  * [11] Lee M Harrison, W Penny, Jean Daunizeau, and Karl J Friston,  “Diffusion-based spatial priors for functional magnetic resonance images,”  Neuroimage, vol. 41, no. 2, pp. 408–423, 2008.
  * [12] Lee M Harrison, Will Penny, Guillaume Flandin, Christian C Ruff, Nikolaus Weiskopf, and Karl J Friston,  “Graph-partitioned spatial priors for functional magnetic resonance images,”  NeuroImage, vol. 43, no. 4, pp. 694–707, 2008.
  * [13] Per Sidén, Anders Eklund, David Bolin, and Mattias Villani,  “Fast Bayesian whole-brain fMRI analysis with spatial 3D priors,”  NeuroImage, vol. 146, pp. 211–225, 2017.
  * [14] Per Sidén, Finn Lindgren, David Bolin, Anders Eklund, and Mattias Villani,  “Spatial 3D Matern priors for fast whole-brain fMRI analysis,”  arXiv preprint arXiv:1906.10591, 2019.
  * [15] Fan RK Chung,  Spectral graph theory,  American Mathematical Society, 1997.
  * [16] Gösta H Granlund and Hans Knutsson,  Signal processing for computer vision,  Kluwer Academic Publishers, 1995.
  * [17] Hans Knutsson and Mats Andersson,  “What’s so good about quadrature filters?,”  in Proceedings 2003 International Conference on Image Processing (Cat. No. 03CH37429). IEEE, 2003, vol. 3, pp. III–61.
  * [18] Hans Knutsson,  “Representing local structure using tensors,”  in Scandinavian Conference on Image Analysis, 1989, pp. 244–251.
  * [19] Hans Knutsson, Carl-Fredrik Westin, and Mats Andersson,  “Representing local structure using tensors II,”  in Scandinavian Conference on Image Analysis. Springer, 2011, pp. 545–556.
  * [20] David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al.,  “The WU-Minn human connectome project: an overview,”  Neuroimage, vol. 80, pp. 62–79, 2013.
  * [21] Xuan Gu, Per Sidén, Bertil Wegmann, Anders Eklund, Mattias Villani, and Hans Knutsson,  “Bayesian diffusion tensor estimation with spatial priors,”  in International Conference on Computer Analysis of Images and Patterns. Springer, 2017, pp. 372–383.