# Interaction of Diatomic Molecules with Photon Angular Momentum

Thomas B. Bahder  Aviation and Missile Research, Development, and Engineering
Center,
US Army RDECOM, Redstone Arsenal, AL 35898, U.S.A.

###### Abstract

The interaction of a diatomic molecule with photons carrying well-defined
angular momentum and parity is investigated to determine whether photon
absorption can induce molecular rotational transitions between states having
angular momentum $\Delta J>1$. A transformation from laboratory coordinates to
coordinates with origin at the center-of-mass of the nuclei is used to obtain
the interaction between the photons and the molecule’s center-of-mass,
electronic, and rotational degrees of freedom. For molecules making
transitions between rotational levels, there is a small parameter, $ka\ll 1$,
where $k$ is the photon wave vector and $a$ is the size of the molecule, which
enters into the $Ej$ and $Mj$ photon absorption probabilities. For electric
photons having arbitrary angular momentum $j\hbar$, the probability of
absorbing an $E(j+1)$ photon divided by the probability of absorbing an $Ej$
photon, scales as $(ka)^{2}/(2j+1)^{2}$. The probability of absorbing an $Mj$
photon, divided by the probability of absorbing and $Ej$ photon scales
according to the same factor.

##  I Introduction

A molecule is significantly more complicated than an atom because, in addition
to center of mass motion, a molecule has vibrational and rotational degrees of
freedom. Work on molecular spectra has a long history. A significant milestone
was made in 1927 by Born and Oppenheimer Born and Oppenheimer (1927), who
assumed that the lighter electrons in a molecule adjust adiabatically to the
slower motion of the heavier nuclei, which remain at all times in their
instantaneous ground state. Further progress in understanding molecular
spectra was made by treating a rotating molecule as a symmetric top or
asymmetric top Wang (1929), sometimes called a rigid rotor Landau and Lifshitz
(1984); Davydov (1965), see Ref. Van Vleck (1929) and references contained
therein. In the rigid rotor treatment of molecules, vibrational and center-of-
mass motion are essentially neglected.

Solving the complete quantum mechanical many body problem consisting of
electrons and nuclei for a given type of molecule will in-principle provide
the energy spectrum of the molecule. In practice, this is a complicated task
which usually results in heavy numerical work. For this reason, much analytic
work has been done on simple diatomic molecules Landau and Lifshitz (1984);
Davydov (1965). For modern reviews on the theory of the molecular
Hamiltonians, see Meyer Meyer (2002), and for an overview about diatomic
molecules, see Brown and Carrington Brown and Carrington (2003).

While the spectrum of a molecule can in-principle be obtained by solving the
corresponding quantum mechanical many body problem, in experiments we usually
probe molecules with electromagnetic fields. In order to interpret
experiments, we must calculate selection rules and probabilities of transition
between molecular quantum states. This involves investigating the coupling of
a molecule to the electromagnetic field. The fundamental theory of a second
quantized electromagnetic radiation field coupled to charges is discussed in
Ref. Akheizer and Berestetskii (1965); Berestetskii et al. (1982). Detailed
treatments dealing with radiation coupling to molecules are given by Atkins
and Wooley Atkins and Woolley (1970), Woolley Woolley (1971, 1975), Craig and
Thirunamachandran Craig and Thirunamachandran (1998) and a recent detailed
treatment is given by Sindelka Sindelka and Moiseyev (2006).

In the past two decades, a completely different line of investigation has
dealt with the angular momentum carried by the electromagnetic field. Beth
Beth (1936) made the first experimental observation that light carries angular
momentum. Beth demonstrated that if a beam of right-hand circularly polarized
light passed through a birefringent plate, and the beam was converted into a
left-hand polarized beam, the plate experienced a torque of $2\hbar$ for each
photon in the beam. This experiment is cited as a demonstration that a light
beam carries spin angular momentum. Years later, in 1992, Allen et al. showed
that laser light with a Laguerre-Gaussian amplitude distribution has well-
defined orbital angular momentum Allen et al. (1992). Many theoretical and
experimental investigations followed, see Ref. Torres and Torner (2011). For
paraxial beams, i.e., in the paraxial approximation, it has been demonstrated
that the total angular momentum can be separated into a sum of spin angular
momentum and orbital angular momentum Allen et al. (1999). However, for a
general electromagnetic radiation field, spin and orbital angular momentum
cannot be separately defined, and only the total angular momentum is a well-
defined quantity Akheizer and Berestetskii (1965); Davydov (1965);
Berestetskii et al. (1982). The interaction of atoms with light carrying
angular momentum were explored theoretically in a number of works Huang
(1994); Jáuregui (2004); Grinter (2008).

Within the rigid rotor model of a molecule, one can easily see that an
electromagnetic radiation field can induce (dipole, quadrupole, etc…)
rotational transitions in a diatomic molecule through its multipole moments,
see for example p. 601 in Section 136 of DavydovDavydov (1965). However, the
rigid rotor treatment neglects the molecule center-of-mass degree of freedom.
In the absence of a radiation field, the molecular center-of-mass motion
decouples from the other degrees of freedom and there is no issue to consider
further, see Eqs. (30) and (38). However, in the presence of a radiation
field, the molecular center of mass plays a key role in the molecule-field
interaction, see Eq. (42) and Ref. Thomas (1970a); Thomas and Joy (1970);
Thomas (1970b, 1971); Sindelka and Moiseyev (2006). References Thomas (1970a);
Thomas and Joy (1970); Thomas (1970b, 1971) discuss the non-trivial role of
the center of molecular mass when a classical plane polarized electromagnetic
wave (corresponding to photons having $j=\hbar$) interacts with a molecule. In
order to investigate in detail molecular absorption of photons having
arbitrary angular momentum ($j>n\hbar)$, where $n>1$, one has to revisit the
interaction Hamiltonian between the molecule and radiation field. As already
remarked, this is complicated because a molecule has several degrees of
freedom: center-of-mass motion, electronic motion, vibrational motion, and
rotational motion.

In this manuscript, I consider a simple diatomic molecule interacting with
photons carrying arbitrary angular momentum, $j\hbar$. I investigate the
mechanism that leads to transitions between molecular rotational levels, with
emphasis on the role of the molecular center-of-mass. Of particular interest
are probabilities of transitions between rotational levels whose angular
momentum differs by more than one unit, $\Delta J>1$. I write down a
Hamiltonian for a diatomic molecule, for the radiation field, and for their
interaction. By making successive transformations from the laboratory
coordinates into the center-of-mass coordinates, and then to center-of-
molecule coordinates, I obtain interaction terms that show the coupling of a
diatomic molecule to the radiation field. For the molecule, there is a small
parameter, $ka\ll 1$, where $k$ is the wave vector and $a$ is the size of the
molecule. I find how the probability of the molecule absorbing an $Ej$ and
$Mj$ multipole photon scales with $ka$. In this treatment, I neglect all
relativistic interactions containing electron and nuclear spins.

##  II Molecular and Field Hamiltonian

I consider a closed system consisting of a diatomic molecule interacting with
a quantized radiation field. The Hamiltonian can be described as a sum of
three terms Davydov (1965); Berestetskii et al. (1982); Akheizer and
Berestetskii (1965); Atkins and Woolley (1970); Woolley (1971, 1975):

${{\hat{H}}}={\hat{H}}_{molecule}+{{\hat{H}}_{field}}+{\hat{H}}_{\mathop{\rm
int}}(t)$ (1)

where ${\hat{H}}_{molecule}$ is the Hamiltonian for the molecule alone (with
no field present), ${{\hat{H}}_{field}}$ is the Hamiltonian for the field
(with no molecule present) and ${\hat{H}}_{\mathop{\rm int}}$ is the
interaction Hamiltonian between the radiation field and the molecule. I assume
the diatomic molecule consists of $n$ electrons and two nuclei. I use a
laboratory system of Cartesian coordinates, ${\bf R}_{i}=(X_{i},Y_{i},Z_{i})$,
$i=1,2,\ldots,n$, for the electrons and ${\bf
R}_{\alpha}=(X_{\alpha},Y_{\alpha},Z_{\alpha})$, $\alpha=1,2$, for the nuclei.
I take the basis vectors of the Cartesian coordinate system to be $({\bf
e}_{x},{\bf e}_{y},{\bf e}_{z})\equiv({\bf e}_{1},{\bf e}_{2},{\bf e}_{3})$.
The position of the origin of the laboratory system of coordinates is
arbitrary. In what follows, I distinguish between coordinates, e.g.,
coordinates of the $i$th electron, ${\bf R}_{i}=(X_{i},Y_{i},Z_{i})$, and the
geometric object, $(X_{i}{\bf e}_{x},Y_{i}{\bf e}_{y},Z_{i}{\bf e}_{z},)$,
which represents the position of the $i$th electron and is independent of the
coordinate system. Note that these two quantities, ${\bf
R}_{i}=(X_{i},Y_{i},Z_{i})$ and $(X_{i}{\bf e}_{x},Y_{i}{\bf e}_{y},Z_{i}{\bf
e}_{z},)$, transform differently under change of coordinate systems. The first
one transforms as coordinate components and the second is a geometric quantity
that does not transform.

The Hamiltonian for the diatomic molecule can be written as Akheizer and
Berestetskii (1965); Meyer (2002)

${{\hat{H}}_{molecule}}=\hat{T}+{{\hat{V}}_{e-e}}+{{\hat{V}}_{n-n}}+{{\hat{V}}_{e-n}}$
(2)

The kinetic energy operator for electrons and nuclei is given by

$\hat{T}=\sum\limits_{i=1}^{n}{\frac{{{\bf{\hat{P}}}_{i}^{2}}}{{2m}}}+\sum\limits_{\alpha=1,2}{\frac{{{\bf{\hat{P}}}_{\alpha}^{2}}}{{2{M_{\alpha}}}}}$
(3)

where $m$ is the mass of an electron and $M_{1}$ and $M_{2}$ are the masses of
the nuclei, and $n$ is the number of electrons. Note that I have not assumed
that the molecule is electrically neutral. The electron and nuclear momentum
operators, ${\bf{\hat{P}}}_{i}=\frac{\partial}{{\partial{{\bf{R}}_{i}}}}$ and
${\bf{\hat{P}}}_{\alpha}=\frac{\partial}{{\partial{{\bf{R}}_{\alpha}}}}$, are
conjugate to the coordinates, ${\bf R}_{i}=(X_{i},Y_{i},Z_{i})$ and ${\bf
R}_{\alpha}=(X_{\alpha},Y_{\alpha},Z_{\alpha})$. The Coulomb interaction terms
for electrons, nuclei, and electrons and nuclei, respectively, are given by

$\begin{array}[]{l}{{\hat{V}}_{e-e}}=\frac{{\,{e^{2}}}}{{4\pi{\varepsilon_{0}}}}\sum\limits_{i<j}{\frac{1}{{\left|{{{\bf{R}}_{i}}-{{\bf{R}}_{j}}}\right|}}}\\\
{{\hat{V}}_{n-n}}=\frac{{{Z_{1}}\,{Z_{2}}\,{e^{2}}}}{{4\pi{\varepsilon_{0}}}}\frac{1}{{\left|{{{\bf{R}}_{2}}-{{\bf{R}}_{1}}}\right|}}\\\
{{\hat{V}}_{e-n}}=-\frac{{\,{e^{2}}}}{{4\pi{\varepsilon_{0}}}}\sum\limits_{\alpha=1,2}{\sum\limits_{i=1}^{n}{\frac{{{Z_{\alpha}}}}{{\left|{{{\bf{R}}_{i}}-{{\bf{R}}_{\alpha}}}\right|}}}}\\\
\end{array}$ (4)

where $e$ is the electronic charge, $\varepsilon_{0}$ is the permitivity of
vacuum, $Z_{\alpha}$ is the atomic number and $M_{\alpha}$ is the mass of
nucleus $\alpha$.

##  III Photons with Well-Defined Angular Momentum

For the the radiation field, I can choose Laguerre-Gaussian modes Allen et al.
(1999); Romero et al. (2002), which are solutions of the paraxial wave
equation, and in which a separation can be made between spin angular momentum
and orbital angular momentum. These types of photons carry angular momentum
$\hbar l$. Such modes are convenient when dealing with cylindrical symmetry,
however, they are only approximate solutions of Maxwell equations.
Alternatively, we can choose vector spherical harmonic modes Davydov (1965);
Berestetskii et al. (1982), which are exact solutions of Maxwell’s equations,
but then spin angular momentum and orbital angular momentum cannot be
separated. This is not a flaw of vector spherical harmonics, rather it is a
general feature of electromagnetic radiation fields. The vector spherical
harmonic modes are a complete set of vector eigenfunctions, so any field can
be expressed in terms of these modes. For these reasons, I take the radiation
field to be described by vector spherical harmonic modes. The Hamiltonian for
the radiation field, see Eq. (1), is given by

${{\hat{H}}_{field}}=\sum\limits_{\lambda,k,j,m}{\hbar{\omega_{k}}\left({{{\hat{a}}^{\dagger}}_{\lambda
kjm}\,{{\hat{a}}_{\lambda kjm}}+\frac{1}{2}}\right)\,}$ (5)

where ${{{\hat{a}}^{\dagger}}_{\lambda kjm}}$ and ${{{\hat{a}}}_{\lambda
kjm}}$ are the creation and destruction operators for electric or magnetic
multipole photons, labeled by $\lambda=(E,M)$, respectively, having wave
vector magnitude $k=\omega_{k}/c$, and carrying angular momentum $j\hbar$ and
angular momentum projection $m\hbar$ on z-axis . The creation and annihilation
operators satisfy the commutation relations Davydov (1965)

$\left[{{{\hat{a}}_{\lambda
kjm}},{{\hat{a}}^{\dagger}}_{\lambda^{\prime}{\kern 1.0pt}k^{\prime}{\kern
1.0pt}j^{\prime}{\kern 1.0pt}m^{\prime}}}\right]={\delta_{\lambda{\kern
1.0pt}\lambda^{\prime}}}{\delta_{k{\kern 1.0pt}k^{\prime}}}\,{\delta_{j{\kern
1.0pt}j^{\prime}}}\,{\delta_{m{\kern 1.0pt}m^{\prime}}}$ (6)

Each field state with well-defined photon number has well-defined energy. Each
photon in a state labeled by quantum numbers ($\lambda,k,j,m$) has wave vector
$k$, square of angular momentum $j(j+1)\hbar^{2}$, z-component of angular
momentum $m\hbar$, and parity $(-1)^{j}$ for $\lambda=E$ and parity
$(-1)^{j+1}$ for $\lambda=M$. Photons are of two types: an $Ej$ photon is an
electric photon and a $Mj$ photon is a magnetic photon, each carrying $j\hbar$
units of angular momentum Davydov (1965); Akheizer and Berestetskii (1965);
Berestetskii et al. (1982). The vector potential operator is given by

${\bf{\hat{A}}}({\bf{r}},t)=\sum\limits_{\lambda,k,j,m}{\left[{{{\hat{a}}_{\lambda
kjm}}\,{\bf{A}}_{kjm}^{\lambda}({\bf{r}},t)+{{\hat{a}}^{\dagger}}_{\lambda
kjm}\,{\bf{A}}_{kjm}^{\lambda\;*}({\bf{r}},t)}\right]}$ (7)

where ${{\bf{A}}_{kjm}^{E}({\bf{r}},t)}$ is the photon wave function for
electric 2j pole radiation and ${{\bf{A}}_{kjm}^{M}({\bf{r}},t)}$ is the
photon wave function for magnetic 2j pole radiation Davydov (1965); Akheizer
and Berestetskii (1965); Berestetskii et al. (1982). In a large sphere of
radius $\cal{R}$, the magnetic and electric photon wave functions are given by
Davydov (1965)

${\bf{A}}_{kjm}^{M}({\bf{r}},t)=\frac{1}{{\sqrt{V}}}{f_{j}}(kr)\,{\bf{Y}}_{jjm}(\theta,\phi)\,{e^{-i{\omega_{k}}t}}$
(8)

and

$\displaystyle{\bf{A}}_{kjm}^{E}({\bf{r}},t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{V}\,\sqrt{2j+1}}\left[\sqrt{j}\,f_{j+1}(kr){\bf{Y}}_{j,j+1,m}(\theta,\phi)\right.$
(9)
$\displaystyle\left.-\sqrt{j+1}\,{f_{j-1}}(kr){\bf{Y}}_{j,j-1,m}(\theta,\phi)\right]\,{e^{-i{\omega_{k}}t}}$

where photon frequency is $\omega_{k}=ck$ and the spherical quantization
volume is $V$. The radial functions, $f_{j}(kr)$, are given in terms of Bessel
functions of the first kind Davydov (1965), $J_{\nu+\frac{1}{2}}(x)$

${f_{l}}(k\,r)=\frac{{{{\left({2V}\right)}^{1/2}}}}{R}\frac{1}{{{J_{l+\frac{3}{2}}}\left({{x_{l+\frac{1}{2},n}}}\right)}}\frac{{{J_{l+\frac{1}{2}}}\left({k\,r}\right)}}{{{r^{1/2}}}}$
(10)

and the allowed set of discrete wave vectors, $k_{ln}$, where $l$ and $n$ are
positive integers, are given by

${J_{l+\frac{1}{2}}}\left({{k_{l{\kern 1.0pt}n}}\,\cal{R}}\right)=0$ (11)

These discrete wave vectors determine the photon energies,
$\hbar\omega_{k}=\hbar ck_{l\,n}$, in a sphere of radius $R$. For a given
value of $l$, the spectrum of allowed wave vectors, $\\{k_{l\,n}\\}$ is the
same for electric and magnetic photons.

The three vector wave functions (vector spherical harmonics), ${\bf
Y}_{j,l,m}(\theta,\phi)$, $j=l+1,l,|l-1|$, are defined as direct products of
two different irreducible representations of the three-dimensional rotation
group, ${\bf e}_{\mu}$, and $Y_{lm}(\theta,\phi)$, coupled by Clebsch-Gordon
coeficients:

${{\bf{Y}}_{j,l,m}}(\theta,\phi)=\sum\limits_{\mu}{\left\langle{{1,l,\mu,m-\mu}}\mathrel{\left|{\vphantom{{1,l,\mu,m-\mu}{jm}}}\right.\kern-1.2pt}{{jm}}\right\rangle}\,{{\bf{e}}_{\mu}}\,{Y_{l,m-\mu}}$
(12)

Here ${\bf e}_{\mu}$ are three spin $S=1$ eigenfunctions,
$\chi_{\mu}(\sigma)$, whose components are ${\bf
e}_{\mu}=\\{\chi_{\mu}(+1),\chi_{\mu}(0),\chi_{\mu}(-1)\\}$, and
$\sigma={+1,0,-1}$ is the spin variable Davydov (1965); Landau and Lifshitz
(1984). The ${\bf e}_{\mu}$, $\mu=\\{-1,0,+1\\}$, are the spherical basis
vectors, which are simultaneous eigenfunctions of the spin angular momentum
operators $\hat{{\bf S}}^{2}$ and $\hat{S}_{z}$:

$\displaystyle\hat{{\bf S}}^{2}\,{\bf e}_{\mu}$ $\displaystyle=$
$\displaystyle 2{\kern 1.0pt}{\bf e}_{\mu}$ $\displaystyle\hat{S}_{z}\,{\bf
e}_{\mu}$ $\displaystyle=$ $\displaystyle\mu\,{\bf e}_{\mu}$ (13)

The spherical basis vectors are related to the Cartesian basis vectors,
$\\{{\bf e}_{x},{\bf e}_{y},{\bf e}_{z}\\}$, by

$\displaystyle{{\bf{e}}_{+}}$ $\displaystyle=$
$\displaystyle-\frac{1}{{\sqrt{2}}}\left({{{\bf{e}}_{x}}+i{{\bf{e}}_{y}}}\right)$
$\displaystyle{{\bf{e}}_{-}}$ $\displaystyle=$
$\displaystyle\frac{1}{{\sqrt{2}}}\left({{{\bf{e}}_{x}}-i{{\bf{e}}_{y}}}\right)$
$\displaystyle{{\bf{e}}_{0}}$ $\displaystyle=$ $\displaystyle{{\bf{e}}_{z}}$
(14)

In terms of the Cartesian basis vectors, $\\{{\bf e}_{x},{\bf e}_{y},{\bf
e}_{z}\\}$, the action of the spin operator is given by

$\hat{S}_{i}{\bf e}_{k}=i\sum_{l=1,2,3}\,e_{ikl}\,{\bf e}_{l}$ (15)

where $e_{ikl}$ is the totally antisymmetric Levi-Civita tensor, where
$e_{123}=+1$.

The $Y_{l,m}$ in Eq. (12) are the usual scalar spherical Harmonics that are
eigenfunctions of the orbital angular momentum operators $\hat{{\bf L}}^{2}$
and $\hat{L}_{z}$:

$\displaystyle\hat{{\bf L}}^{2}\,Y_{lm}$ $\displaystyle=$ $\displaystyle
l(l+1)Y_{lm}$ $\displaystyle\hat{L}_{z}\,Y_{lm}$ $\displaystyle=$
$\displaystyle mY_{lm}$ (16)

According to the definitions in Eq.(12)–(16), the product functions,
${{\bf{e}}_{\mu}}\,{Y_{l,m}}$, are simultaneous eigenfunctions of $\hat{{\bf
L}}^{2}$, $\hat{L}_{z}$, $\hat{{\bf S}}^{2}$ and $\hat{S}_{z}$, with
eigenvalues $l(l+1)$, $m$, $2$, and $\mu$, respectively.

The spherical vector harmonics, ${{\bf{Y}}_{j,l,m}}(\theta,\phi)$, form an
irreducible representaion of the three-dimensional rotation group on the unit
sphere, and are simultaneous eigenfunctions of $\hat{J}^{2}$, $\hat{J}_{z}$,
$\hat{{\bf L}}^{2}$, and $\hat{{\bf S}}^{2}$:

$\displaystyle{{\hat{J}}_{z}}{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$
$\displaystyle=$ $\displaystyle m{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$
(17)
$\displaystyle{{\hat{J}}^{2}}{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$
$\displaystyle=$ $\displaystyle
j(j+1){{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$ (18)
$\displaystyle{{\hat{L}}^{2}}{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$
$\displaystyle=$ $\displaystyle
l(l+1){{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$ (19)
$\displaystyle{{\hat{S}}^{2}}{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$
$\displaystyle=$ $\displaystyle 2{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right)$
(20)

The vector spherical harmonics, ${{\bf{Y}}_{j,l,m}}(\theta,\phi)$, satisfy the
orthogonality relation

$\int{{d^{2}}\Omega}\,{\bf{Y}}_{jlm}^{*}\left({\theta,\phi}\right)\cdot{\bf{Y}}_{j^{\prime}l^{\prime}m^{\prime}}\left({\theta,\phi}\right)={\delta_{j\,j^{\prime}}}\,{\delta_{l\,l^{\prime}}}\,{\delta_{m\,m^{\prime}}}$
(21)

where the integration is over the unit sphere, and completeness relation

$\sum\limits_{j\,l\,m}{{{\bf{Y}}_{j,l,m}}\left({\theta,\phi}\right){{\bf{Y}}_{j,l,m}}\left({\theta^{\prime},\phi^{\prime}}\right)}={\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\leftrightarrow$}}\over{I}}}}\,{\delta^{2}}\left({\Omega-\Omega^{\prime}}\right)$
(22)

Under the parity operator (spatial inversion), the basis vectors ${\bf
e}_{\mu}$ change sign, and we take the scalar spherical harmonics, $Y_{lm}$,
to be defined so that they are multiplied by $(-1)^{l}$, and therefore the
vector spherical harmonics ${\bf Y}_{j,l,m}(\theta,\phi)$ are eigenfunctions
of parity with eigenvalues $(-1)^{l+1}$.

From Eq. (8) we see that magnetic photons have well-defined spin and orbital
angular momentum, since ${\bf{A}}_{kjm}^{M}({\bf{r}},t)$ is an eigenfunction
of spin and orbital angular momentum. However, from Eq. (9) we see that
electric photons with quantum numbers $(j,m)$ are a superposition of orbital
angular momentum $l=j-1$ and $l=j+1$, and do not have well-defined orbital
angular momentum.

The number operator for photons of type $\lambda$ is given by

$\hat{{n}}_{\lambda}=\sum\limits_{k,j,m}{{{\hat{a}}^{\dagger}}_{\lambda
kjm}\,{{\hat{a}}_{\lambda kjm}}\,}$ (23)

##  IV Molecule-Field Interaction

The electromagnetic field couples to the charges in the molecule, which
consist of electrons and nuclei. The interaction Hamiltonian describing the
coupling of the radiation field to the molecule is given by Atkins and Woolley
(1970); Woolley (1971, 1975); Berestetskii et al. (1982); Meyer (2002);
Sindelka and Moiseyev (2006)

${{\hat{H}}_{{\mathop{\rm
int}}}}(t)=\frac{e}{m}\sum\limits_{i=1}^{n}{{\bf{\hat{A}}}({{\bf{R}}_{i}},t)\cdot{{{\bf{\hat{P}}}}_{i}}}-e\sum\limits_{\alpha=1}^{2}{\frac{{{Z_{\alpha}}}}{{{M_{\alpha}}}}{\bf{\hat{A}}}({{\bf{R}}_{\alpha}},t)\cdot{{{\bf{\hat{P}}}}_{\alpha}}}$
(24)

where the vector potential operator is in laboratory coordinates, and is given
in Eq. (7). In Eq. (24), the first term is the coupling of the field to
electrons and the second term is the coupling of the field to the nuclei. The
radiation field, ${\bf{\hat{A}}}({\bf r},t)$, is quantized in the laboratory
coordinates, ${\bf R}$. The electron and nuclear momentum operators in Eq.
(24), ${\bf P}_{i}$ and ${\bf P}_{\alpha}$, are also in laboratory
coordinates. I have dropped terms that are quadratic in the vector potential
because they are responsible for photon scattering and are small for photon
absorption and emission. Also, I have dropped small relativistic terms that
couple the electron and nuclear spins to the field Landau and Lifshitz (1984);
Davydov (1965); Berestetskii et al. (1982). Furthermore, in Eq. (24), I have
assumed the Coulomb gauge for the vector potential, so that ${\rm
div}\,\,{\bf{\hat{A}}}({\bf r},t)=0$ and the scalar potential is zero.

In many experiments, the wavelength of the electromagnetic radiation field is
long compared to the dimensions of the physical system, such as the molecule.
We can make a Taylor series expansion of the vector potential terms
${\bf{\hat{A}}}({\bf R}_{i},t)$ and ${\bf{\hat{A}}}({\bf R}_{\alpha},t)$ about
the center of mass of the molecule ${\bf R}_{o}$. Keeping only the first term,
I make the approximations

${\bf{\hat{A}}}({{\bf{R}}_{i}},t)\approx{\bf{\hat{A}}}({{\bf{R}}_{\alpha}},t)\approx{\bf{\hat{A}}}({{\bf{R}}_{o}},t)$
(25)

##  V Transformation to Center of Mass of the Molecule

In order to explore the transfer of angular momentum from photons to
molecules, we must look at degrees of freedom that can ”absorb” angular
momentum. To this end, I transform the molecular Hamiltonian in Eq. (2) from
the laboratory coordinates to center-of-mass-of-molecule coordinates. I take
the Cartesian basis vectors in laboratory coordinates to be collinear with the
basis vectors in center-of-mass-of-molecule coordinates. I follow the notation
used by Brown and Carrington Brown and Carrington (2003). The position of the
center of mass of the molecule, in laboratory coordinates, is given by

${{\bf{R}}_{o}}=\frac{m}{M}\sum\limits_{i=1}^{n}{{{\bf{R}}_{i}}+\frac{1}{M}\sum\limits_{\alpha=1}^{2}{{M_{\alpha}}{{\bf{R}}_{\alpha}}}}$
(26)

where the total molecular mass is given by

$M=\sum\limits_{i=1}^{n}{m_{i}}+\sum\limits_{\alpha=1}^{2}{{M_{\alpha}}}=nm+M_{1}+M_{2}$
(27)

where $n$ is the number of electrons in the molecule. The transformation of
electron and nuclear coordinates, from laboratory coordinates,
${{\bf{R}}_{i}}$ and ${{\bf{R}}_{\alpha}}$, to the center of mass of molecule
coordinates, ${{{\bf{R^{\prime}}}}_{i}}$ and ${{{\bf{R^{\prime}}}}_{\alpha}}$,
is given by

$\displaystyle{{{\bf{R^{\prime}}}}_{i}}$ $\displaystyle=$
$\displaystyle{{\bf{R}}_{i}}-{{\bf{R}}_{o}}$
$\displaystyle{{{\bf{R^{\prime}}}}_{\alpha}}$ $\displaystyle=$
$\displaystyle{{\bf{R}}_{\alpha}}-{{\bf{R}}_{o}}$ (28)

Note that the transformations in Eq. (28) are transformations of vector
components, see comment above Eq. (2). I introduce the inter-nuclear vector,
${\bf R}$, in terms of the coordinates of the two nuclei (in laboratory
coordinates)

${\bf R}={\bf R}_{2}-{\bf R}_{1}$ (29)

Using the above definitions of the transformation to the center-of-mass-of-
molecule coordinates, the transformation of the kinetic energy of the
molecule, given in Eq. (3), to the center-of-mass-of-molecule coordinates,
gives Brown and Carrington (2003)

$T=\frac{1}{{2M}}{\bf{P}}_{o}^{2}+\frac{1}{{2\mu}}{\bf{P}}_{\bf{R}}^{2}+\frac{1}{{2m}}\sum\limits_{i=1}^{n}{{\bf{P^{\prime}}}_{i}^{2}-}\frac{1}{{2M}}\sum\limits_{j,k=1}^{n}{{{{\bf{P^{\prime}}}}_{i}}\cdot}{{{\bf{P^{\prime}}}}_{k}}$
(30)

where the momentum operators

${{\bf{P}}_{o}}=\frac{\hbar}{i}\frac{\partial}{{\partial{{\bf{R}}_{o}}}}$ (31)
${{\bf{P}}_{\bf{R}}}=\frac{\hbar}{i}\frac{\partial}{{\partial{\bf{R}}}}$ (32)
${{\bf{P}}_{j}^{\prime}}=\frac{\hbar}{i}\frac{\partial}{{\partial{{\bf{R}}_{j}^{\prime}}}}$
(33)

are conjugate to the coordinates, ${\bf R}_{o}$, ${\bf R}$ and ${\bf
R}_{j}^{\prime}$, and where $\mu$ is the reduced nuclear mass

$\frac{1}{\mu}=\frac{1}{{{M_{1}}}}+\frac{1}{{{M_{2}}}}$ (34)

The first term in Eq. (30) is the center of mass kinetic energy of the
molecule, expressed in laboratory coordinates. The second term in Eq. (30) is
the relative kinetic energy of the nuclei, expressed in terms of the inter-
nuclear vector ${\bf R}$. The third term in Eq. (30) is the kinetic energy of
the electrons. The last term in Eq. (30) is the so-called mass polarization
terms and results from fluctuations of electron position about the massive
nuclei.

The transformation of the Coulomb interaction terms, given in Eq. (4), from
the laboratory coordinates to center-of-mass-of-molecule coordinates, is done
similarly.

##  VI Transformation to Center of Mass of the Nuclei

In order to discuss rotation of the molecule about the line joining the
nuclei, I define the position of the center of mass of the nuclei (in center-
of-mass-of-molecule coordinates):

${{{\bf{R^{\prime}}}}_{N}}=\frac{1}{{{M_{1}}+{M_{2}}}}\left({{M_{1}}{{{\bf{R^{\prime}}}}_{1}}+{M_{2}}{{{\bf{R^{\prime}}}}_{2}}}\right)$
(35)

The positions of the electrons in the center-of-mass-of-nuceli coordinates are
given by

$\displaystyle{{{\bf{R^{\prime\prime}}}}_{i}}$ $\displaystyle=$
$\displaystyle{{{\bf{R^{\prime}}}}_{i}}-{{{\bf{R^{\prime}}}}_{N}}$ (36)
$\displaystyle=$
$\displaystyle{{{\bf{R^{\prime}}}}_{i}}+\frac{1}{{{M_{1}}+{M_{2}}}}\sum\limits_{j=1}^{n}{{m_{j}}}{{{\bf{R^{\prime}}}}_{j}}$

In the center-of-mass-of-nuceli coordinates, the momentum operator conjugate
to the electron position, ${{{\bf{R^{\prime\prime}}}}_{i}}$, is given by

${{{\bf{P^{\prime\prime}}}}_{j}}=\frac{\hbar}{i}\frac{\partial}{{\partial{{{\bf{R^{\prime\prime}}}}_{j}}}}$
(37)

The transformation of nuclear coordinates, ${\bf R}_{\alpha}$, to center-of-
mass-of-nuceli coordinates is similar.

The total kinetic energy of the molecule, given in equation Eq. (30),
transformed to center-of-nuclei coordinates, is given by Judd (1975); Brown
and Carrington (2003)

$\displaystyle T$ $\displaystyle=$
$\displaystyle\frac{1}{2M}{\bf{P}}_{o}^{2}+\frac{1}{{2\mu}}{\bf{P}}_{\bf{R}}^{2}+\frac{1}{{2m}}\sum\limits_{i=1}^{n}{\bf{P^{\prime\prime}}}_{i}^{2}$
(38) $\displaystyle+$
$\displaystyle\frac{1}{{2\left({{M_{1}}+{M_{2}}}\right)}}\sum\limits_{j,k=1}^{n}{{{{\bf{P^{\prime\prime}}}}_{j}}\cdot}{{{\bf{P^{\prime\prime}}}}_{k}}$

Transforming the Coulomb terms in Eq. (2) and (4) to the center-of-mass-of-
the-nuclei coordinates (which are designated by two primed superscripts),
these terms become

$\begin{array}[]{l}{{\hat{V}}_{e-e}}=\frac{{\,{e^{2}}}}{{4\pi{\varepsilon_{0}}}}\sum\limits_{i<j}{\frac{1}{{\left|{{{{\bf{R^{\prime\prime}}}}_{i}}_{j}}\right|}}}\\\
{{\hat{V}}_{n-n}}=\frac{{{Z_{1}}\,{Z_{2}}\,{e^{2}}}}{{4\pi{\varepsilon_{0}}}}\frac{1}{R}\\\
{{\hat{V}}_{e-n}}=-\frac{{\,{e^{2}}}}{{4\pi{\varepsilon_{0}}}}\sum\limits_{\alpha=1,2}{\sum\limits_{i=1}^{n}{\frac{{{Z_{\alpha}}}}{{\left|{{{{\bf{R^{\prime\prime}}}}_{i\alpha}}}\right|}}}}\\\
\end{array}$ (39)

where $|{{{\bf R^{\prime\prime}}}}_{ij}|=\left|{{{\bf
R^{\prime\prime}}}}_{i}-{{{\bf R^{\prime\prime}}}}_{j}\right|$ is the distance
between electrons, $|{{{\bf R^{\prime\prime}}}}_{i\alpha}|=\left|{{{{\bf
R^{\prime\prime}}}}_{i}}-{{{\bf R^{\prime\prime}}}}_{\alpha}\right|$ is the
distance between the $i$th electron and $\alpha$th nucleus, and ${\bf R}$ is
the inter-nuclear vector defined in Eq. (29), with magnitude $R={\bf R}$. Note
that the ${{{\bf R^{\prime\prime}}}}_{i\alpha}$ are expressible in terms of
the electron coordinates, ${{{\bf R^{\prime\prime}}}}_{i}$, and the inter-
nuclear vector, ${\bf R}$, so the molecular wavefunction depends on the center
of mass coordinate, ${\bf R}_{o}$, the inter-nuclear coordinate, ${\bf R}$,
and the electron coordinates ${{{\bf R^{\prime\prime}}}}_{i}$, for
$i=1,2,\cdots,n$.

In summary, when the molecule Hamiltonian in laboratory coordinates, given in
Eq. (2), is transformed into the center-of-mass-of-nuceli coordinates, it is
given by the sum of Eq. (38) and (39). Note that in the absence of the
electromagnetic field, the center-of-mass motion of the molecule decouples
from the internal motion of the molecule, i.e., the molecular Hamiltonian does
not contain mixed terms with center-of-mass and internal degrees of freedom.
So, in the absence of the radiation field, the total molecular wavefunction is
the product:

${\psi_{QLM}}\left({{{\bf{R}}_{o}}}\right)\,\phi\left({{\bf{R}},{{{\bf{R^{\prime\prime}}}}_{1}},{{{\bf{R^{\prime\prime}}}}_{2}},\cdots,{{{\bf{R^{\prime\prime}}}}_{n}},}\right)$
(40)

where the center-of-mass wavefunction is of the form

${\psi_{QLM}}\left({{{\bf{R}}_{o}}}\right)={R_{QL}}\left({Q{R_{o}}}\right)\,{Y_{LM}}\left({{\theta_{o}},{\phi_{o}}}\right)$
(41)

and
$\phi\left({{\bf{R}},{{{\bf{R^{\prime\prime}}}}_{1}},{{{\bf{R^{\prime\prime}}}}_{2}},\cdots,{{{\bf{R^{\prime\prime}}}}_{n}},}\right)$
is the wave function for internal degrees of freedom.

##  VII Interaction Hamiltonian in CM of Nuclei Coordinates

As previously stated, in the absence of the radiation field, the molecule is
described by a Hamiltonian that is the sum of Eq. (38) and (39). In this case,
the center of mass motion, described by coordinate ${\bf R}_{o}$, decouples
from the internal motion, which is described by coordinate ${\bf R}$. However,
when an external radiation field is present, see the interaction Hamiltonian
given in Eq. (42), the radiation field introduces a coupling between the
center of mass motion of the molecule and its internal motions Sindelka and
Moiseyev (2006).

Next, I assume that the molecule interacts with radiation, according to the
interaction Hamiltonian in Eq. (24). It is possible to re-write the molecule-
field interaction in a form that allows more insight than that given by the
interaction Hamiltonian in Eq. (24). In order to do this, I make the long
wavelength approximation, given in Eq. (25), in the interaction Hamiltonian in
Eq. (24). Then, I transform the electron and nuclear momentum operators to
center-of-mass-of-nuclei coordinates. The molecule-field interaction
Hamiltonian ${\hat{H}_{{\mathop{\rm int}}}}(t)$ in Eq. (24) then becomes

$\displaystyle{\hat{H}_{{\mathop{\rm int}}}}(t)$ $\displaystyle=$
$\displaystyle\frac{e}{M}(n-{Z_{1}}-{Z_{2}})\,{\bf{\hat{A}}}({{\bf{R}}_{o}},t)\cdot\,{{{\bf{\hat{P}}}}_{o}}+$
(42) $\displaystyle+$ $\displaystyle{\kern
1.0pt}e\left({\frac{1}{m}+\frac{{{Z_{1}}+{Z_{2}}}}{{{M_{1}}+{M_{2}}}}}\right){\kern
1.0pt}{\bf{\hat{A}}}({{\bf{R}}_{o}},t)\cdot\,\sum\limits_{i=1}^{n}{{{{\bf{P^{\prime\prime}}}}_{i}}}$
$\displaystyle+$ $\displaystyle{\kern
1.0pt}e\left({\frac{{{Z_{1}}}}{{{M_{1}}}}-\frac{{{Z_{2}}}}{{{M_{2}}}}}\right){\kern
1.0pt}{\bf{\hat{A}}}({{\bf{R}}_{o}},t)\cdot\,{{{\bf{\hat{P}}}}_{R}}$

where ${{{\bf{\hat{P}}}}_{o}}$ is the relative momentum operator for the
center of mass of the molecule (in laboratory coordinates) that is defined in
Eq. (31), ${{{{\bf{P^{\prime\prime}}}}_{i}}}$ is the momentum operator of the
$i$th electron (in center-of-mass-of-nuclei coordinates) given in Eq. (37),
and ${{{\bf{\hat{P}}}}_{R}}$ is the momentum operator that is conjugate to the
inter-nuclear vector ${\bf R}$ and is given by Eq. (29).

The total Hamiltonian for the molecule, the field, and their interaction is
given by the sum of Eq. (38), (39) and (42). The coordinates in the total
Hamiltonian are the center-of-mass position (in laboratory coordinates), ${\bf
R}_{o}$, the inter-nuclear position (nuclear relative coordinate) ${\bf R}$,
and the electron positions in the center-of-mass-of-nuclei coordinates,
${{\bf{R^{\prime\prime}}}_{i}}$. The basis states for the combined system of
molecule and field are of the form

$\Psi\left({{{\bf{R}}_{o}},{\bf{R}},{{{\bf{R^{\prime\prime}}}}_{1}},\ldots,{{{\bf{R^{\prime\prime}}}}_{n}}}\right)\,\left|{{{({n_{1}})}_{{\alpha_{1}}}}\,,{{({n_{2}})}_{{\alpha_{2}}}}\,,\cdots}\right\rangle$
(43)

where
$\Psi\left({{{\bf{R}}_{o}},{\bf{R}},{{{\bf{R^{\prime\prime}}}}_{1}},\ldots,{{{\bf{R^{\prime\prime}}}}_{n}}}\right)$
are the molecular states that are eigenstates of the sum of Eq. (38) and (39),
and
$\left|{{{({n_{1}})}_{{\alpha_{1}}}}\,,{{({n_{2}})}_{{\alpha_{2}}}}\,,\cdots}\right\rangle$
are the photon states that are eigenstates of Eq. (5), where there are $n_{i}$
photons in mode $\alpha_{i}$ and I use the short-hand notation
$\alpha=\\{\lambda,k,j,m\\}$.

The first term in the interaction Hamiltonian in Eq. (42) is an ionic term,
and is proportional to the pre-factor $n-{Z_{1}}-{Z_{2}}$. This pre-factor is
zero for neutral molecules, whose electron number is equal to the sum of the
atomic numbers of the two nuclei, $n={Z_{1}}+{Z_{2}}$. For the case of a
molecule that is not neutral (an ion) this first term couples the momentum of
the center of mass of the molecule to the radiation field. This term leads to
center of mass motion of the (charged) molecule as a whole.

The second term in the interaction Hamiltonian in Eq. (42) couples the sum of
electron momenta (in center-of-mass-of-nuclei coordinates) to the radiation
field. This term allows a transfer of angular momentum from the radiation
field to the electrons. The pre-factor for this second term is never zero.

The third term in the interaction Hamiltonian in Eq. (42) couples the
momentum,

${{\bf{P}}_{\bf{R}}}=-i\hbar\left[{\frac{\partial}{{\partial{R_{X}}}}{{\bf{e}}_{x}}+\frac{\partial}{{\partial{R_{Y}}}}{{\bf{e}}_{y}}+\frac{\partial}{{\partial{R_{Z}}}}{{\bf{e}}_{z}}}\right]$
(44)

which is conjugate to the inter-nuclear vector,
${\bf{R}}=\left({{R_{X}},{R_{Y}},{R_{Z}}}\right)$, to the electromagnetic
radiation field, ${\bf{\hat{A}}}({{\bf{R}}_{o}},t)$. Here, $({\bf e}_{x},{\bf
e}_{y},{\bf e}_{z})$ are the basis vectors for the space-fixed coordinates of
${\bf{R}}=\left({{R_{X}},{R_{Y}},{R_{Z}}}\right)$. The pre-factor for this
term is zero for homo-nuclear diatomic molecules Sindelka and Moiseyev (2006),
which have $M_{1}=M_{2}$ and $Z_{1}=Z_{2}$. Homo-nuclear diatomic molecules
have an inversion symmetry about the nuclear center of mass. Molecules which
have $Z_{1}=Z_{2}$ but $M_{1}\neq M_{2}$, (e.g., different isotopes, such as
HD), lack this inversion symmetry, and then this third term in Eq. (42) is
generally non-zero. This third term contains a coupling of the radiation field
to the rotation of the molecule about the line joining the nuclei. This can be
seen by considering the unitary transformation, $\cal{M}$, from space-fixed
coordinates with origin at the center-of-nuclear mass,
${\bf{X}}=\left({{X},{Y},{Z}}\right)$, to molecule-fixed coordinates,
${\bf{x}}=\left({{x},{y},{z}}\right)$, where the $z$-axis is pointing along
the inter-molecular vector ${\bf R}_{2}-{\bf R}_{1}$:

$\left({\begin{array}[]{*{20}{c}}X\\\ Y\\\ Z\\\
\end{array}}\right)=\cal{M}\left({\phi,\theta,\chi}\right)\cdot\left[{\begin{array}[]{*{20}{c}}x\\\
y\\\ z\\\ \end{array}}\right]$ (45)

where $\left({\phi,\theta,\chi}\right)$ are the Euler angles between
coordinates ${\bf{X}}=\left({{X},{Y},{Z}}\right)$ and
${\bf{x}}=\left({{x},{y},{z}}\right)$. I take the sequence of rotations
through angles $\\{\phi,\theta,\chi\\}$ to be along the $z$, $y$, $z$ axes, so
the rotation matrix is given by Judd (1975); Zare (1988); Brown and Carrington
(2003)

$\cal{M}\left({\phi,\theta,\chi}\right)=\left({\begin{array}[]{*{20}{c}}{\cos\phi\cos\theta\cos\chi-\sin\phi\sin\chi}&{-\sin\phi\cos\chi-\cos\phi\cos\theta\sin\chi}&{\cos\phi\sin\theta}\\\
{\sin\phi\cos\theta\cos\chi+\cos\phi\sin\chi}&{\cos\phi\cos\chi-\sin\phi\cos\theta\sin\chi}&{\sin\phi\sin\theta}\\\
{-\sin\theta\cos\chi}&{\sin\theta\sin\chi}&{\cos\theta}\\\
\end{array}}\right)$ (46)

For internuclear vector ${\bf R}=\left({{R_{X}},{R_{Y}},{R_{Z}}}\right)$, the
transformation from space-fixed coordinates to molecule-fixed coordinates,
$\left({{R_{X}},{R_{Y}},{R_{Z}}}\right)\rightarrow\left({R,\theta,\phi,\chi}\right)$,
is given by

$\left({\begin{array}[]{*{20}{c}}{{R_{X}}}\\\ {{R_{Y}}}\\\ {{R_{Z}}}\\\
\end{array}}\right)=\cal{M}\left({\phi,\theta,\chi}\right)\cdot\left[{\begin{array}[]{*{20}{c}}0\\\
0\\\ R\\\ \end{array}}\right]$ (47)

For a diatomic molecule, since the vector
${\bf{R}}=\left({{R_{X}},{R_{Y}},{R_{Z}}}\right)$ is parallel to the molecule-
fixed $z$-axis, the components $\left({{R_{X}},{R_{Y}},{R_{Z}}}\right)$ do not
depend on $\chi$, and we are free to choose its value for convenience Judd
(1975); Zare (1988); Brown and Carrington (2003). Using the transformation in
Eq. (47), the relative momemtum in Eq. (44) becomes Brown and Carrington
(2003):

${{\bf{P}}_{\bf{R}}}=\hbar\left[{-\frac{1}{R}{\bf{n}}\times{\bf{N}}-i{\bf{n}}\frac{\partial}{{\partial
R}}}\right]$ (48)

where ${\bf n}$ is the inter-nuclear unit vector, given in space-fixed
coordinates by

${\bf{n}}=\frac{{\bf R}_{2}-{\bf R}_{1}}{|{\bf R}_{2}-{\bf
R}_{1}|}=\sin\theta\,\cos\phi\,{{\bf{e}}_{x}}+\sin\theta\,\sin\phi\,{{\bf{e}}_{y}}+\cos\theta\,\,{{\bf{e}}_{z}}$
(49)

and ${\bf N}$ is the angular momentum operator in space-fixed coordinates with
origin at center-of-mass of the nuclei

$\displaystyle{\bf{N}}$ $\displaystyle=$ $\displaystyle
i\left[{\cot\theta\cos\phi\frac{\partial}{{\partial\phi}}+\sin\phi\frac{\partial}{{\partial\theta}}-{\mathop{\rm
c}\nolimits}\sec\theta\cos\phi\frac{\partial}{{\partial\chi}}}\right]{{\bf{e}}_{x}}$
(50) $\displaystyle+$ $\displaystyle
i\left[{\cot\theta\sin\phi\frac{\partial}{{\partial\phi}}-\cos\phi\frac{\partial}{{\partial\theta}}-{\mathop{\rm
c}\nolimits}\sec\theta\sin\phi\frac{\partial}{{\partial\chi}}}\right]{{\bf{e}}_{y}}$
$\displaystyle-$ $\displaystyle
i\frac{\partial}{{\partial\phi}}{{\bf{e}}_{z}}$

Using the form of ${\bf P}_{\bf R}$ in Eq. (48) in the third term in Eq. (42),
gives

${\bf{\hat{A}}}({{\bf{R}}_{o}})\cdot\,{{{\bf{\hat{P}}}}_{R},t}=-{\bf{\hat{A}}}({{\bf{R}}_{o}},t)\cdot\left[{\frac{1}{R}{\bf{n}}\times{\bf{N}}+i{\bf{n}}\frac{\partial}{{\partial
R}}}\right]$ (51)

which shows that the radiation field couples to both molecular rotations and
molecular vibrations. The first term in Eq. (51) is a coupling of the
radiation field to molecular rotations, through the molecular angular momentum
${\bf N}$. The second term in Eq. (51) is a coupling of the radiation field to
molecular vibrations through the length of the inter-nuclear distance $R$.

Babiker et al. investigated the interaction of a hydrogenic bound system
consisting of a nucleus and a bound electron, interacting with a linearly
polarized radiation field carrying orbital angular momentum Babiker et al.
(2002). They found that an exchange of orbital angular momentum in an electric
dipole transition occurs only between light and the center of mass motion of
the atom—they concluded that the internal electronic-type motion does not
participate in any exchange of orbital angular momentum in a dipole
transition. The presence of the second term in Eq. (42) suggests that this is
not the case for a diatomic molecule.

Alexandrescu et al. studied a charged diatomic molecule (an ion, with two
nuclei) interacting with a Laguerre-Gaussian beam carrying orbital angular
momentum Alexandrescu et al. (2006). They seem to disagree with the
conclusions of Babiker et al., stating that the electronic motion and the
electromagnetic field can exchange one unit of OAM within the electronic
dipole interaction.

I note that the interaction Hamiltonian in Eq. (42) is not of the form $-{\bf
d}\cdot{\bf E}$, which is a dipole coupling to the radiation field ${\bf E}$.
The interaction Hamiltonian in Eq. (42) contains all multipole couplings with
the radiation field. The interaction Hamiltonian in Eq. (42) allows
investigation of transfer of angular momentum from the radiation field to the
molecule, including transfer of angular momentum to the rotation of the
nuclei. Regarding the dipole moment, consider the dipole moment of the
molecule, expressed in terms of variables ${\bf R^{\prime\prime}}_{i}$ and
${\bf R^{\prime\prime}}_{\alpha}$ in the center-of-mass-of-nuclei coordinates:

${\bf{d}}=-{\bf{e}}\sum\limits_{i=1}^{n}{{{{\bf{R^{\prime\prime}}}}_{i}}}+e\sum\limits_{\alpha=1}^{2}{{Z_{\alpha}}{{{\bf{R^{\prime\prime}}}}_{\alpha}}}$
(52)

When the dipole moment is transformed to the set of coordinates $\\{{\bf
R}_{o},{\bf R},{\bf R^{\prime\prime}}_{1},\ldots{\bf
R^{\prime\prime}}_{n}\\}$, it becomes

$\displaystyle{\bf{d}}$ $\displaystyle=$
$\displaystyle-{\bf{e}}\left[{1+\frac{m}{M}\left({{Z_{1}}+{Z_{2}}}\right)}\right]\sum\limits_{i=1}^{n}{{{{\bf{R^{\prime\prime}}}}_{i}}}$
(53) $\displaystyle+$ $\displaystyle
e\left({{Z_{1}}+{Z_{2}}}\right){{\bf{R}}_{o}}-e\left({\frac{{{Z_{1}}}}{{{M_{1}}}}-\frac{{{Z_{2}}}}{{{M_{2}}}}}\right)\mu{\kern
1.0pt}{\bf{R}}$

In order to isolate the rotational part of the Hamiltonian, contained in ${\bf
P}_{\bf R}$, it was necessary to transform to coordinates $\\{{\bf R}_{o},{\bf
R},{\bf R^{\prime\prime}}_{1},\ldots{\bf R^{\prime\prime}}_{n}\\}$. However,
after performing his transformation, the Hamiltonian in Eq. (42) is not of the
form of a dipole coupling to the radiation field, $-{\bf d}\cdot{\bf E}$. This
is not a result of the transformation itself, rather it is due to the added
rotational degrees of freedom of a molecule as compared to an atom.

##  VIII Transition Probabilities

The form of the molecule-radiation interaction in Eq. (42) shows that there is
a clear coupling of the radiation field to molecular rotations. The center of
mass coordinate, ${\bf R}_{o}$, plays a central role in this coupling. As
stated earlier, in order to determine the effect of this coupling on
observable quantities, transition matrix elements and probabilities for photon
absorption must be calculated. In order to estimate these matrix elements, we
must have knowledge of the form of the wavefunctions for the diatomic
molecule, which are schematically shown in Eq. (43). We must also assume some
form for the modes of the radiation field, for example, vector spherical
harmonics or Laguerre-Gaussian modes. Group theoretic methods based on tensor
operators can then be used to determine ratios of transition matrix elements
Judd (1975); Zare (1988); Brown and Carrington (2003). This work is in
progress. However, some statements can be made about parity selection rules
and how the transition matrix elements scale.

Exact selection rules based on conservation of parity can be written down.
Assume that under inversion of coordinates (under the parity operator), the
parity of the initial and final molecular wavefunctions are $P_{i}$ and
$P_{f}$, respectively, and the parity of the emitted photon is $P_{\gamma}$,
where parity eigenvalues takes values $\pm 1$. If the initial state of the
whole system has one photon, and the final state has no photons, then parity
conservation requires $P_{i}\,P_{\gamma}=P_{f}$, or, since each parity
eigenvalue is $\pm 1$, parity conservation requires Berestetskii et al. (1982)

$P_{i}\,P_{f}=P_{\gamma}$ (54)

Electric $Ej$ photons have parity $P^{Ej}_{\gamma}=(-1)^{j}$ and magnetic $Mj$
photons have parity $P^{Mj}_{\gamma}=(-1)^{j+1}$. For absorbing one photon of
a given type, parity conservation requires

$P_{i}\,P_{f}=(-1)^{j}\qquad{\rm electric\,\,\,{\it Ej}\,\,\,photon}$ (55)

or

$P_{i}\,P_{f}=(-1)^{j+1}\qquad{\rm magnetic\,\,\,{\it Mj}\,\,\,photon}$ (56)

An estimate of the scaling of transition probabilities can be made by
considering the molecule as a system of charges interacting with a radiation
field. Consider a molecule in an initial state with energy $E_{i}$, absorbing
a single photon of energy $\hbar\omega_{k}=\hbar ck$, and making a transition
to a final state with energy $E_{f}$. The transition matrix element contains
the quantity $kr$, which enters in the vector potential in Eqs. (8) and (9).
The quantity $r$ has a characteristic scale of a diatomic molecule, on the
order of $a\sim 10^{-10}$ m. For transitions between rotational states, I take
the frequency $f\sim 300$ GHz. Therefore, there is a small dimensionless
quantity, $kr\sim ka\sim 2\pi a/\lambda\sim 6\times 10^{-7}$ in the transition
matrix elements. The vector potentials in Eqs. (8) and (9) depend on the
function $f_{l}(kr)$, which in turn depend on the Bessel function
$J_{l+\frac{1}{2}}(kr)$. Using the small argument, $kr\ll 1$, expansion of the
Bessel function, in Eq. (10), the radial function, $f_{l}(kr)$, scales as (up
to a constant):

${f_{l}}\left({k\,r}\right)\sim{\left({k{\kern
1.0pt}r}\right)^{l+\frac{1}{2}}}\frac{1}{{1\cdot 3\cdot 5\cdots(2l+1)}}$ (57)

The vector potential for $Ej$ photons scales as
${\bf{A}}_{kjm}^{E}\left({\bf{r}}\right)\sim{f_{j-1}}\left({k\,r}\right)$, and
the vector potential for $Mj$ photons scales as
${\bf{A}}_{kjm}^{M}\left({\bf{r}}\right)\sim{f_{j}}\left({k\,r}\right)$. The
probability of absorbing a photon scales as the square of the vector potential
(square of the matrix element). Therefore, the ratio of the probability for
absorbing an $Mj$ photon, $P(Mj)$, divided by the probability of absorbing an
$Ej$ photon, $P(Ej)$, scales as

$\frac{{P\left({Mj}\right)}}{{P\left({Ej}\right)}}\sim\frac{{{{\left({k{\kern
1.0pt}a}\right)}^{2}}}}{{{{\left({2j+1}\right)}^{2}}}}\sim\frac{{{{4\times
10}^{-13}}}}{{{{\left({2j+1}\right)}^{2}}}}$ (58)

For a given angular momentum $j$, the probability of absorption of a magnetic
photon is much less than absorption of an electric photon.

The ratio of the probability for absorbing an $E(j+1)$ photon, $P(E(j+1)$,
divided by the probability of absorbing an $Ej$ photon, $P(Ej)$, also scales
as

$\frac{{P\left({E(j+1)}\right)}}{{P\left({Ej}\right)}}\sim\frac{{{{\left({k{\kern
1.0pt}a}\right)}^{2}}}}{{{{\left({2j+1}\right)}^{2}}}}\sim\frac{{{{4\times
10}^{-13}}}}{{{{\left({2j+1}\right)}^{2}}}}$ (59)

Form Eq. (59), we see that the probability of absorbing higher angular
momentum photons rapidly decreases with increasing $j$. As the simplest
example, the ratio of probability of absorbing an $E2$ photon to the
probability of absorbing an $E1$ photon is:

$\frac{{P\left({E2}\right)}}{{P\left({E1}\right)}}\sim{4\times 10^{-14}}$ (60)

The scaling of the absorption probabilities for molecules, given in Eqs.
(58)–(60), makes it a challenging proposition to observe absorption of photons
having angular momentum $j>1$.

##  IX Summary

I have investigated the interaction between a diatomic molecule and photons
with well-defined angular momentum. I have exploited the transformation from
laboratory coordinates to coordinates with origin at the center of mass of the
nuclei to show the nuclear rotational degrees of freedom explicitly (through
the vector ${\bf R}$). For the photon modes, I have assumed spherical vector
harmonics in a large sphere of radius $R$. These modes are exact solutions of
Maxwell equations, and consequently, any other modes, such as the Laguerre-
Gaussian modes, can be expanded in terms of them. However, in general, spin
and orbital angular momentum are not separable in vector spherical harmonic
modes. I used the standard minimal coupling between elementary charges and
radiation field to couple the radiation field to the molecule, neglecting
quadratic terms in the vector potential and neglecting relativistic
interaction with spins. For molecules, which absorb at millimeter wavelengths,
there is a small parameter, $k\,a$ where $k$ is a wave vector and $a$ is the
size of the molecule. This small parameter was used to determine the scaling
of absorption probabilities of $Ej$ and $Mj$ type photons. For transitions
between rotational levels of the molecule having $\Delta J>1$, the probability
of absorption of photons with angular momentum $j>1$ is small, because the
factor $ka\ll 1$, where wave vector $k=2\pi/\lambda$ corresponds to photons in
the millimeter wavelength region, see Section VIII. Note that the small
parameter, $k\,a$, also exists for atoms Huang (1994); Jáuregui (2004);
Grinter (2008), however, it is larger for atoms because $k$ corresponds to
visible wavelengths, whereas for rotational transitions in molecules $k$
corresponds to millimeter wavelengths.

For electric photons having arbitrary angular momentum $j\hbar$, the
probability of rotational transition by absorbing an electric $E(j+1)$ photon
divided by the probability of absorbing an electric $Ej$ photon, scales as
$(ka)^{2}/(2j+1)^{2}$. The probability of absorbing a magnetic $Mj$ photon is
much smaller than the probability of absorbing an electric $Ej$ photon. The
probability of absorbing a magnetic $Mj$ photon, divided by the probability of
absorbing an electric $Ej$ photon, also scales as $(ka)^{2}/(2j+1)^{2}$.

###### Acknowledgements.

The author acknowledges interesting discussions with Henry Everitt and Matt
Goodman.

## References

  * Born and Oppenheimer (1927) M. Born and J. R. Oppenheimer, Ann. Phys. 84, 457 484 (1927).
  * Wang (1929) S. C. Wang, Phys. Rev. 34, 243 (1929).
  * Landau and Lifshitz (1984) L. D. Landau and E. M. Lifshitz, _QuantumMechanics–Non-Relativistic Theory_ (Pergamon Press, New York, N.Y. USA, 1984), 3rd ed.
  * Davydov (1965) A. S. Davydov, _Quantum Mechanics_ (Pergamon Press, New York, N.Y., USA, 1965), 2nd ed.
  * Van Vleck (1929) J. H. Van Vleck, Phys. Rev. 33, 467 (1929), URL http://link.aps.org/doi/10.1103/PhysRev.33.467.
  * Meyer (2002) H. Meyer, Annu. Rev. Phys. Chem. 53, 141 (2002).
  * Brown and Carrington (2003) J. Brown and A. Carrington, _Rotational Spectroscopy of Diatomic Molecules_ (Cambridge University Press, New York, N. Y., USA, 2003).
  * Akheizer and Berestetskii (1965) A. I. Akheizer and V. B. Berestetskii, _Quantum Electrodynamics_ (John Wiley & Sons, New York, N.Y., USA, 1965).
  * Berestetskii et al. (1982) V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitaevskii, _Quantum Electrodynamics_ (Pergamon Press, Elmsford, New York, USA, 1982), 2nd ed.
  * Atkins and Woolley (1970) P. W. Atkins and R. G. Woolley, Proc. Roy. Soc. Lond. A 319, 549 (1970).
  * Woolley (1971) R. G. Woolley, Proc. Roy. Soc. Lond. A 321, 557 (1971).
  * Woolley (1975) R. G. Woolley, in _Volume 33_ , edited by I. Progogine and S. A. Rice (University of Chicago, Chicago, Illinois, 1975), Advances in Chemical Physics, pp. 153–233.
  * Craig and Thirunamachandran (1998) D. P. Craig and T. Thirunamachandran, _Molecular Quantum Electrodynamics_ (Dover Publications, Inc., Mineola, N.Y., USA, 1998), dover edition ed.
  * Sindelka and Moiseyev (2006) M. Sindelka and N. Moiseyev, J. Phys. Chem. A 110, 5561 5571 (2006), URL http://arxiv.org/abs/quant-ph/0601195.
  * Beth (1936) R. A. Beth, Phys. Rev. 50, 115 (1936).
  * Allen et al. (1992) L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Phys. Rev. A 45, 8185 (1992), URL http://link.aps.org/doi/10.1103/PhysRevA.45.8185.
  * Torres and Torner (2011) J. P. Torres and L. Torner, eds., _Twisted Photons: Applications of Light With Orbital Angular Momentum_ (WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2011).
  * Allen et al. (1999) L. Allen, M. J. Padgett, and M. Babiker, in _Progress on Optics, Volume 39_ , edited by E. Wolf (Elsevier, New York, N.Y. USA, 1999), Progress on Optics.
  * Huang (1994) K.-N. Huang, Am. J. Phys. 62, 73 (1994).
  * Jáuregui (2004) R. Jáuregui, Phys. Rev. A 70, 033415 (2004), URL http://link.aps.org/doi/10.1103/PhysRevA.70.033415.
  * Grinter (2008) R. Grinter, J. Phys. B: At. Mol. Opt. Phys. 41, 1 (2008).
  * Thomas (1970a) I. L. Thomas, Phys. Rev. A 2, 72 (1970a), URL http://link.aps.org/doi/10.1103/PhysRevA.2.72.
  * Thomas and Joy (1970) I. L. Thomas and H. W. Joy, Phys. Rev. A 2, 1200 (1970), URL http://link.aps.org/doi/10.1103/PhysRevA.2.1200.
  * Thomas (1970b) I. L. Thomas, Phys. Rev. A 2, 1675 (1970b), URL http://link.aps.org/doi/10.1103/PhysRevA.2.1675.
  * Thomas (1971) I. L. Thomas, Phys. Rev. A 3, 1022 (1971), URL http://link.aps.org/doi/10.1103/PhysRevA.3.1022.
  * Romero et al. (2002) L. C. D. Romero, D. L. Andrews, and M. Babiker, J. Opt. B: Quantum Semiclass. Opt. 4, S66 (2002).
  * Judd (1975) B. R. Judd, _Angular Momentum Theory for Diatomic Molecules_ (Academic Press, New York, USA, 1975).
  * Zare (1988) R. N. Zare, _Angular Momentum: Understanding Spatial Aspecta in Chemistry and Physics_ (Wiley-Interscience Publication, New York, USA, 1988).
  * Babiker et al. (2002) M. Babiker, C. R. Bennett, D. L. Andrews, and L. C. Dávila Romero, Phys. Rev. Lett. 89, 143601 (2002), URL http://link.aps.org/doi/10.1103/PhysRevLett.89.143601.
  * Alexandrescu et al. (2006) A. Alexandrescu, D. Cojoc, and E. D. Fabrizio, Phys. Rev. Lett. 96, 243001 (2006).


# Geometrical optics of large deviations of Brownian motion in inhomogeneous
media

Tal Bar<EMAIL_ADDRESS>Baruch Meerson<EMAIL_ADDRESS>Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904,
Israel

###### Abstract

Geometrical optics provides an instructive insight into Brownian motion,
“pushed” into a large-deviations regime by imposed constraints. Here we extend
geometrical optics of Brownian motion by accounting for diffusion
inhomogeneity in space. We consider three simple model problems of Brownian
motion on the line or in the plane in situations where the diffusivity of the
Brownian particle depends on one spatial coordinate. One of our results
describes “Brownian refraction”: an analog of refraction of light passing
through a boundary between two media with different refraction indices.

##  I Introduction

Large deviations of Brownian motion have recently attracted much attention
from statistical physicists Grosberg ; Ikeda ; Schuss ; Nechaev ; SM2019 ;
M2019 ; MS2019a ; M2019b ; MM2020 ; Agranovetal ; Vladimirov ; M2020 ; Valov .
One setting where large deviations arise naturally is stochastic searching of
a distant target by multiple random searchers, such as receptors searching for
a reaction site on the surface of a living cell Coombs , or sperm cells
searching for an oocyte Eisenbach . In such systems the multiple searchers
essentially compete among themselves for reaching the target first and
performing a biological function. When the number of searchers is very large
(for example, there are about $3\times 10^{8}$ sperm cells which attempt to
reach the oocyte after copulation in humans), the searcher which succeeds in
arriving first does it in an unusually short time. As a result, this
searcher’s trajectory looks quite differently from _typical_ trajectories of a
single searcher conditioned on reaching the target at _any_ time Grosberg ;
Ikeda ; Schuss ; SM2019 ; M2019 ; MS2019a ; M2019b ; MM2020 ; Agranovetal ;
M2020 .

For a single Brownian searcher, the probability of an extremely fast arrival
at the target is very small; it corresponds to a short-time tail of the
complete probability distribution of the arrival times. Short-time tails of
this type can be efficiently described theoretically by the optimal
fluctuation method which, in the context of Brownian motion, is called
“geometrical optics of Brownian motion” Grosberg ; Ikeda ; Schuss ; SM2019 ;
M2019 ; MS2019a ; M2019b ; MM2020 ; Agranovetal ; M2020 . An important
advantage of geometrical optics is that it allows to account, in a simple
manner, for additional constraints. Furthermore, apart from providing the
distribution tail of interest, geometrical optics predicts the optimal path of
the system, that is the most likely trajectory of the Brownian particle which
dominates the probability distribution tail. The optimal path gives a valuable
insight into the physics of the large deviation in question.

Previous works on geometrical optics of Brownian motion assumed that the
Brownian motion, “pushed” into a large-deviation regime by imposed
constraints, occurs in a homogeneous medium (possibly with some impenetrable
regions) Grosberg ; Ikeda ; Schuss ; SM2019 ; M2019 ; MS2019a ; M2019b ;
MM2020 ; Agranovetal ; M2020 . In some applications, however, the homogeneity
assumption is unrealistic. For example, the diffusivity of a reactant
macromolecule, exploring a living cell can vary quite significantly due to a
varying physical and chemical composition of different cell regions Bressloff
. Here we consider large deviations of a Brownian particle with a space-
dependent diffusivity, as described by the Langevin equation

$\dot{\mathbf{x}}=\sqrt{2D[\mathbf{x}(t)]}\,\boldsymbol{\xi}(t)\,.$ (1)

Here $\mathbf{x}(t)$ is the displacement of the particle, $D=D(\mathbf{x})$ is
the space-dependent diffusivity, and $\boldsymbol{\xi}(t)$ is the Gaussian
white noise with unit norm and zero mean. Using the Hänggi-Klimontovich
interpretation of the multiplicative noise in Eq. (1) Sokolov , one can write
down the diffusion equation

$\frac{\partial P(\mathbf{x},t)}{\partial
t}=\nabla\cdot\left[D(\mathbf{x})\nabla P(\mathbf{x},t)\right]$ (2)

for the probability $P(\mathbf{x},t)$ of observing the particle at position
$\mathbf{x}$ at time $t$. In one spatial dimension Eqs. (1) and (2) were
studied in Ref. key-20  in connection to phenomena accompanying anomalous
diffusion, and in Ref. Farago2014  in the context of fluctuation-dissipation
relation. Here we are interested in large deviations in this system. The
starting point of our analysis will be the Wiener’s action

$-\ln P[\mathbf{x}(t)]\simeq
S=\int_{0}^{T}\frac{\dot{\mathbf{x}}^{2}(t)}{4D\left[\mathbf{x}(t)\right]}\,dt\,,$
(3)

which gives the probability $\mathcal{P}[\mathbf{x}(t)]$ of a Brownian path
$\mathbf{x}\left(t\right)$ up to pre-exponential factors. A large-deviation
regime corresponds to a very large action, $S\gg 1$. It is this regime where
the (tail of the) probability distribution of interest is dominated by a
single optimal path $\mathbf{x}(t)$: a deterministic trajectory which
minimizes the Wiener’s action (3) over all possible trajectories subject to
additional, problem-specific constraints.

This is the essence of the method of geometrical optics of large deviations of
Brownian motion Grosberg ; Ikeda ; Schuss ; SM2019 ; M2019 ; MS2019a ; M2019b
; MM2020 ; Agranovetal ; M2020 , previously developed for homogeneous media.
For $D=\text{const}$ (and in the absence of any integral constraints) the
action (3) is minimized by a particle motion with constant speed,
$|\dot{\mathbf{x}}|=\text{const}$, along a geodesic: the shortest path between
the initial and final position of the particle, subject to additional
constraints, if any. The action along such a path is

$S=\frac{1}{4D}\int_{0}^{T}\left(\frac{\mathcal{L}}{T}\right)^{2}\,dt=\frac{\mathcal{L}^{2}}{4DT},$
(4)

where $\mathcal{L}$ is the path’s length. For $D=\text{const}$ the problem
therefore reduces to minimizing $\mathcal{L}$ subject to the additional
constraints MS2019a .

Here we are interested in inhomogeneous diffusion as described by Eqs.
(1)-(3). To this end we consider three simple model problems where $D$ depends
only on one spatial coordinate $x$. As we will see, the diffusion
inhomogeneity can lead to interesting effects. One of them is the “Brownian
refraction”, which emphasizes the geometrical nature of large deviations of
Brownian motion and provides an additional justification to the term
“geometrical optics of diffusion”.

Here is a layout of the remainder of the paper. In Sec. II we use the
geometrical optics to determine the short-time asymptotic of the propagator of
the Brownian motion on the line with a piecewise constant diffusivity $D(x)$.
In this case the propagator can be also found exactly by solving Eq. (2) (see
the Appendix). This allows us to verify geometrical optics’ predictions.
Importantly, the latter also include the optimal path of the system, which is
not directly available from the exact solution. In Sec. III we find the short-
time asymptotic of the propagator of the Brownian motion on the plane with a
piecewise constant $D(x)$. It is here where we encounter the Brownian
refraction. In Sec. IV the geometrical optics is used to determine the short-
time asymptotic of the propagator of the Brownian motion on the line with a
continuously varying diffusivity $D(x)$. A brief discussion of our results and
of their possible extensions is included in Sec. V.

##  II Brownian motion on the line with a piecewise constant $D(x)$

Our first example is Brownian motion on a line with a piecewise constant
diffusivity:

$D\left(x\right)=\begin{cases}D_{1},&x<\ell\,,\\\
D_{2}\equiv\kappa^{2}D_{1},&x>\ell,\end{cases}$ (5)

where $\ell>0$. The particle starts at the origin at $t=0$, and we are
interested in the short-time, $T\to 0$, asymptotic of the probability density
$P(L,T)$ of observing the particle at a point $x=L>\ell$ at time $T>0$.

On each of the intervals $0<x<\ell$ and $\ell<x<L$ the diffusivity is
constant. Therefore, in view of the above-mentioned geometrical optics results
for $D=\text{const}$, the optimal path $x(t)$ must be a piecewise linear
function of time. Let us denote by $\tau$ the (a priori unknown) time of
crossing the boundary $x=\ell$ between the two constant-$D$ regions (we assume
that this time is unique.) Demanding that $x(t)$ obey the conditions $x(0)=0$
and $x(T)=L$, we obtain

$x\left(t\right)=\begin{cases}\frac{\ell t}{\tau},&0<t<\tau\,,\\\
\frac{\ell(T-t)+L(t-\tau)}{T-\tau},&\tau<t<T\,.\end{cases}$ (6)

The action on this composite trajectory is the sum of contributions, described
by Eq. (4), from each of the two intervals:

$S\simeq\frac{\ell^{2}}{4D_{1}\tau}+\frac{\left(L-\ell\right)^{2}}{4D_{2}\left(T-\tau\right)}=\frac{1}{4}\left(\frac{t_{1}}{\tau}+\frac{t_{2}}{T-\tau}\right)\,,$
(7)

where

$t_{1}=\frac{\ell^{2}}{D_{1}}\quad\mbox{and}\quad
t_{2}=\frac{\left(L-\ell\right)^{2}}{D_{2}}=\frac{\left(L-\ell\right)^{2}}{\kappa\sqrt{D_{1}}}$
(8)

are the characteristic diffusion times through each of the two media. The
minimum action is achieved for

$\tau=\frac{\sqrt{t_{1}}}{\sqrt{t_{1}}+\sqrt{t_{2}}}\,T=\frac{T}{\frac{1}{\kappa}\left(\frac{L}{\ell}-1\right)+1}\,,$
(9)

and we obtain the short-time asymptotic

$-\ln P(L,T\to 0)\simeq
S_{\text{min}}=\frac{\left(\sqrt{t_{1}}+\sqrt{t_{2}}\right)^{2}}{4T}=\frac{\left[L-\left(1-\kappa\right)\ell\right]^{2}}{4D_{2}T}\,.$
(10)

The same $T\to 0$ asymptotic (10) follows from exact solution of this problem
(see the Appendix for details):

$\displaystyle P_{\text{exact}}\left(L,T\right)$
$\displaystyle=\frac{1}{\sqrt{\pi
T}\left(\sqrt{D_{1}}+\sqrt{D_{2}}\right)}e^{-\frac{\left[L-\left(1-\kappa\right)\ell\right]^{2}}{4D_{2}T}}\,.$
(11)

Indeed, when $T\to 0$, the pre-exponential factor in Eq. (11) (which is missed
by the leading-order geometrical optics calculations) becomes subleading, and
Eqs.(10) and (11) coincide in the leading order.

Figure 1 a and b gives two examples of simulated Brownian paths in this
setting. The particle starts at the origin and arrives at $x=L$ at two
different times $T$. The figure also shows the piecewise linear optimal path
$x(t)$ at $T\to 0$, predicted by Eqs. (6) and (9) for the same $D_{1}$ and
$D_{2}$. As $T$ decreases, the simulated paths start approaching the predicted
short-time optimal path.

Figure 1: Solid lines: two simulated Brownian paths for the setting of Sec.
II. The parameters are $L=6$, $\ell=3$, $D_{1}=10$, $D_{2}=3$, and $T=6$ (a)
and $T=0.45$ (b). The dashed lines show the piecewise linear optimal path
predicted by geometrical optics, Eqs. (6) and (9).

Notice that the geometrical optics result (10) gives a correct leading-order
asymptotic $P(L,T\to 0)$ for a whole class of settings dealing with a Brownian
particle starting at the origin and arriving at $x=L$ at $t=$. One such
setting – when the Brownian motion is defined on the whole line – is solved
exactly in the Appendix. But one could actually consider a Brownian motion on
any interval $(a,b)$ which includes in itself the interval $0<x<L$, with
either reflecting or absorbing boundary conditions at the boundaries $x=a$ and
$x=b$. In this case the exact solution would of course be different from Eq.
(11), but its leading-order $T\to 0$ asymptotic will still be described by Eq.
(10). This is explained by the character of the optimal path: when conditioned
on reaching $x=L$ in a very short time, the particle cannot spend any time on
exploring regions beyond the interval $(0,L)$, rendering them irrelevant for
the leading-order calculations.

##  III Brownian motion on the plane with a piecewise constant $D(x)$:
Brownian refraction

Now let us consider a Brownian motion in the $(x,y)$ plane with a piecewise
constant diffusivity

$D\left(x\right)=\begin{cases}D_{1}\,,&x<\ell\,,\\\
D_{2}\equiv\kappa^{2}D_{1}\,,&x>\ell\,,\end{cases}$ (12)

for any $y$. The Brownian particle starts at the origin at $t=0$, and we are
interested in the short-time asymptotic, $T\to 0$, of the probability
$P(L_{x},L_{y},T)$ for the particle to reach the point $(L_{x},L_{y})$, so
that $L_{x}>\ell$, at time $t=T$, see Fig. 2.

Figure 2: A simulated Brownian path for the setting of Sec. III with
$\kappa>1$. The parameters are $L_{x}=L_{y}=10$, $\ell=5$, $D_{1}=1$,
$D_{2}=3$, and $T=4.17$ (solid line). The dashed lines show the piecewise
linear optimal path predicted by the geometrical optics in the limit of $T\to
0$. This path obeys the Snell’s law (16).

We assume that the particle crosses the boundary $x=\ell$ between the two
regions with different diffusivities at a unique time $0<\tau<T$, and the
$y$-coordinate of the particle at this moment is $y_{0}$. The action on this
trajectory is the sum of two terms, described by Eq. (4), which come from each
of the two media:

$S(\tau,y_{0})=\frac{1}{4}\left[\frac{\ell^{2}+y_{0}^{2}}{D_{1}\tau}+\frac{\left(L_{x}-\ell\right)^{2}+\left(L_{y}-y_{0}\right)^{2}}{D_{2}\left(T-\tau\right)}\right]=\frac{1}{4}\left(\frac{t_{1}}{\tau}+\frac{t_{2}}{T-\tau}\right)\,,$
(13)

with obvious definitions of the characteristic diffusion times $t_{1}$ and
$t_{2}$ for each of the two media. Now we should minimize the action (13) with
respect to the two auxiliary parameters $\tau$ and $y_{0}$. It is convenient
to do it sequentially. Notice that Eqs. (13) and (7) have the same form in
terms of $\tau$, $t_{1}$ and $t_{2}$. Therefore, at fixed $y_{0}$, the minimum
of the action with respect to $\tau$ is achieved at
$\tau=\frac{\sqrt{t_{1}}}{\sqrt{t_{1}}+\sqrt{t_{2}}}T$ [see Eq. (9)], and we
obtain

$S(y_{0})=\frac{1}{4T}\left(\sqrt{t_{1}}\\!+\\!\sqrt{t_{2}}\right)^{2}\\!=\\!\frac{1}{4T}\\!\left[\sqrt{\frac{\ell^{2}\\!+\\!y_{0}^{2}}{D_{1}}}+\sqrt{\frac{\left(L_{x}\\!-\\!\ell\right)^{2}\\!+\\!\left(L_{y}\\!-\\!y_{0}\right)^{2}}{D_{2}}}\right]^{2}.$
(14)

The minimization with respect to the remaining parameter $y_{0}$ demands
$dS(y_{0})/dy=0$, which leads to the relation

$\frac{L_{y}-y_{0}}{\sqrt{D_{2}}\sqrt{\left(L_{x}-l\right)^{2}+\left(L_{y}-y_{0}\right)^{2}}}=\frac{y_{0}}{\sqrt{D_{1}}\sqrt{\ell^{2}+y_{0}^{2}}}$
(15)

To better understand this result we notice that $y_{0}$ and $L_{y}-y_{0}$ can
be written as $y_{0}=\ell\tan\theta_{1}$, and
$L_{y}-y_{0}=\left(L_{x}-\ell\right)\tan\theta_{2}$, respectively, where
$\theta_{1}$ ($\theta_{2}$) is the angle between the particle trajectory in
the left (right) region and the straight line $y=y_{0}$. In other words,
$\theta_{1}$ and $\theta_{2}$ are the angles of incidence and refraction,
respectively. In the new notation Eq. (15) becomes, remarkably,

$\frac{\sin\theta_{1}}{\sin\theta_{2}}=\frac{\sqrt{D_{1}}}{\sqrt{D_{2}}}=\frac{1}{\kappa}\,,$
(16)

providing an unexpected Brownian-motion analog of Snell’s law of optics (the
latter describes the refraction of light or other waves passing through a
boundary between two media with different refraction indices Crawford ;
Feynman ). Notice that the role of the refraction index $n$ in optics is
played here by $1/\sqrt{D}$. The “Brownian refraction”, described by Eq. (16),
provides a further justification to the term “geometrical optics of (large
deviations of) Brownian motion” Grosberg ; Ikeda ; Schuss ; SM2019 ; M2019 ;
MS2019a ; M2019b ; MM2020 ; Agranovetal ; M2020 111To remind the reader, the
Snell’s law of refraction can be derived by minimizing the arrival time of the
ray to the target point Feynman . For the Brownian motion the arrival time is
fixed, and one minimizes the Wiener’s action (3). In the light of these
differences the Brownian refraction phenomenon is somewhat surprising..

Let us complete the solution. The relation (15) yields a quartic equation for
$y_{0}$, the analytical solutions of which is too bulky to work with. It is
more convenient to present the solution in a parametric form. Using Eq. (15),
we can express $L_{y}$ via $y_{0}$:

$L_{y}=y_{0}\left[1+\frac{\left(L_{x}-l\right)\kappa}{\sqrt{l^{2}+\left(1-\kappa^{2}\right)y_{0}^{2}}}\right].$
(17)

Plugging this expression into Eq. (14), we obtain

$S=\frac{\ell^{2}+y_{0}^{2}}{4D_{2}T}\left[\kappa+\frac{\left(L_{x}-\ell\right)}{\sqrt{\ell^{2}+\left(1-\kappa^{2}\right)y_{0}^{2}}}\right]^{2}\,.$
(18)

Equations (17) and (18) yield $-\ln P\simeq S$ (which we are after) in a
parametric form, where $y_{0}$ serves as the parameter.

It is instructive to look at the limits of $S$ at $\kappa\to 0$ and
$\kappa\to\infty$ at fixed $D_{2}$:

$S\left(\kappa\right)\simeq\begin{cases}\frac{\left(L_{x}-\ell\right)^{2}}{4D_{2}T}\,,&\kappa\to
0\,,\\\ \frac{\ell^{2}}{4D_{1}T}\,,&\kappa\to\infty\,.\end{cases}$ (19)

The corresponding optimal paths in these two limits are the following. At
$\kappa\to 0$ we have $y_{0}\to L_{y}$, and $\tau\to 0$. The form of the
trajectory is

$y(x,\kappa\to 0)\simeq\begin{cases}\frac{L_{y}}{\ell}\,x\,,&0<x<\ell\,,\\\
L_{y}\,,&\ell<x<L_{x}\,.\end{cases}$ (20)

Here the diffusion in media 1 is fast, and the particle reaches the point
$(\ell,L_{y})$ almost immediately so as to minimize the distance it has to
travel in media 2, where the diffusivity is finite. The particle spends almost
all of the allocated time $T$ in media 2.

As $\kappa\to\infty$, $y_{0}$ goes to zero, and $\tau$ approaches $T$. Now the
form of the trajectory is

$y(x,\kappa\to\infty)\simeq\begin{cases}0\,,&0<x<\ell\,,\\\
\frac{L_{y}(x-\ell)}{L_{x}-\ell}\,,&\ell<x<L_{x}\,.\end{cases}$ (21)

Here the diffusion in media 1 is slow. Therefore, the particle moves along the
$x$-axis so as to minimize the distance it has to travel here, and its spends
almost all of the allocated time $T$ to arrive at $x=\ell$. Then it almost
immediately “jumps” to the point $(L_{x},L_{y})$ along the straight line.

##  IV Brownian motion on the line with a continuously varying $D(x)$

Here we return to one spatial dimension and suppose that a Brownian particle
is released at $t=0$ at $x=0$ in a medium with a continuously varying
diffusivity $D\left(x\right)$.

We are interested in the short-time asymptotic of the probability $P(L,T)$ of
observing the particle at $x=L$ at time $t=T$. Using geometrical optics, we
identify the Lagrangian

$\mathfrak{L}\left(x,\dot{x}\right)=\frac{\dot{x}^{2}}{4D(x)}$ (22)

of the Wiener’s action (3). Since this Lagrangian does not depend explicitly
on time, there is an integral of motion which, in this case, coincides with
the Lagrangian itself, $\mathfrak{L}\left(x,\dot{x}\right)=\text{const}$.
Therefore, we can write

$\dot{x}=C\sqrt{\ D\left(x\right)}\,,$ (23)

where $C$ is a constant to be determined. A comparison of Eq. (23) with the
Langevin equation (1) shows that the optimal realization of the rescaled
Gaussian white noise $\xi(t)$, conditioned on a very fast arrival of the
Brownian particle at $x=L$, is simply a constant, equal to $C/\sqrt{2}$.

Equation (23) can be solved in quadratures for any $D(x)$, and we obtain

$Ct=\int_{0}^{x}\frac{dy}{\sqrt{D(y)}}\,,$ (24)

where we used the boundary condition $x(0)=0$. Using the second boundary
condition $x(T)=L$, we determine the constant $C$:

$C=\frac{1}{T}\int_{0}^{L}\frac{dx}{\sqrt{D(x)}}\,.$ (25)

Equations (24) and (25) yield the optimal path $x(t)$ in an implicit form for
any $D(x)$. Further, plugging Eq. (23) into Eq. (3) and using Eq. (25), we
obtain a general expression for the action:

$S=\int_{0}^{T}\frac{\dot{x}^{2}}{4D\left(x\right)}\,dt=\frac{1}{4T}\left[\int_{0}^{L}\frac{dx}{\sqrt{D(x)}}\right]^{2}\,.$
(26)

Figure 3: The function $f(\beta)$ from Eq. (28).

As a simple example, consider a linear diffusivity profile,
$D(x)=D_{0}(1+\beta x/L)$, where we must assume $\beta>-1$ in order to keep
$D(x)$ positive in the relevant region $0<x<L$. Using Eqs. (24), we obtain

$x(t)=L\frac{t}{T}\left[\frac{t}{T}+\frac{2\left(\sqrt{\beta+1}-1\right)}{\beta}\left(1-\frac{t}{T}\right)\right]\,.$
(27)

The action (26) takes the form

$S=\frac{L^{2}}{D_{0}T}\,f(\beta)\,,\;\mbox{where}\;\>f(\beta)=\frac{\left(\sqrt{\beta+1}-1\right)^{2}}{\beta^{2}}\,.$
(28)

It gives the short-time asymptotic of $P(L,T)$ via the relation $-\ln
P(L,T)\simeq S$. A graph of the function $f(\beta)$ is shown in Fig. 3. As one
can see, the action decreases monotonically with an increase of $\beta$. It
reaches its maximum value $L^{2}/(DT)$ at $\beta=-1$. In this case the
diffusivity vanishes at $x=L$. As $\beta$ goes to infinity, $f(\beta)$ goes
down as $1/\beta$. Finally, at $\beta\to 0$ we reproduce the constant-
diffusivity action $S=L^{2}/(4D_{0}T)$.

The optimal path $x(t)$, as described by Eq. (27), is shown in Fig. 4 for
three different values of $\beta$. In general, $x(t)$ is a parabola: concave
or convex depending on whether $\beta$ is positive or negative, respectively.
For $\beta=0$ (a constant diffusivity) $x(t)$ is a straight line as to be
expected.

Figure 4: The optimal path $x(t)$, as described by Eq. (27), for the linear
diffusivity profile $D(x)=D_{0}(1+\beta x/L)$ with $\beta=10$ (blue) and $-1$
(magenta). The black line shows $x(t)$ for the constant diffusivity,
$\beta=0$.

##  V Discussion

We considered three simple model problems, illustrating the role of spatial
variations of the diffusivity in large deviations of Brownian motion. One of
our predictions is what can be called “Brownian refraction”: a close analog of
refraction of light passing through a boundary between two media with
different refraction indices. This phenomenon appears when one conditions the
Brownian particle on arrival at a distant point in a very short time.

It would be interesting to extend the geometrical-optics approach by
accommodating additional local or integral constraints, or constraints in the
form of inequalities, like it was done for the homogeneous diffusion Grosberg
; Ikeda ; Schuss ; SM2019 ; M2019 ; MS2019a ; M2019b ; MM2020 ; Agranovetal ;
M2020

When dealing, in Sec. III, with a Brownian motion on the plane, we assumed
that the diffusivity $D$ depends only on one spatial coordinate $x$, and this
dependence is piecewise constant. This theory can be extended to a
continuously varying $D(x)$. Indeed, the Lagrangian of Eq. (3),

$\mathfrak{L}(x,\dot{x},\dot{y})=\frac{\dot{x}^{2}}{4D(x)}+\frac{\dot{y}^{2}}{4D(x)}$
(29)

describes an effective classical mechanics with two degrees of freedom. This
system possesses two integrals of motion. The first of them is the Lagrangian
itself,

$\frac{\dot{x}^{2}}{4D(x)}+\frac{\dot{y}^{2}}{4D(x)}=E=\text{const}\,,$ (30)

An additional integral of motion, $\dot{y}/D[x(t)]=\text{const}$, follows from
the fact that the Lagrangian (29) does not depend on $y$. The presence of two
integrals of motion allows one to solve the problem in quadratures for any
$D(x)$.

Extensions of the theory to a more general case, where $D(x,y)$ depends both
on $x$, and on $y$, encounter an immediate difficulty coordinates , because in
this case the geometrical-optics description of the system involves only one
integral of motion: the Lagrangian $\mathfrak{L}(x,y,\dot{x},\dot{y})$ itself.
The lack of additional independent integral of motion makes the geometrical-
optics formulation non-integrable analytically Tabor . Numerical solutions,
however, are certainly possible. In addition, perturbative analytical
treatments can be possible if there is a small parameter in the coordinate
dependence $D(x,y)$. It would be interesting to explore these possibilities in
a future work.

Finally, geometrical optics of large deviations of Brownian motions present
only one, although important, example of a class of systems which can be
studied with the optimal fluctuation method (OFM). In recent years the OFM has
been employed for analysis of large deviations in additional Markovian
Touchette  and non-Markovian M2019 ; M2022 ; MO2022 ; M2023  stochastic
processes. Another well-known large-deviation technique, based on the Feynman-
Kac formula DV ; EG , applies to a certain class of large deviations of _time-
averaged_ quantities. There is also a group of large-deviation methods, in
different areas of physics, which rely on the “single-big-jump” principle, see
Ref. bigjump  and references therein. A skillful use of all these methods (and
sometimes of their combinations Smith ) should help develop a better
understanding of many important, fascinating and, quite often, non-intuitive
large deviations and rare events.

## Acknowledgment

We thank Oded Farago for useful comments. BM was supported by the Israel
Science Foundation (Grant No. 1499/20).

##  Appendix A Exact solution for Sec. II

Here we present exact solution of the equation

$\frac{\partial P(x,t)}{\partial t}=\frac{\partial}{\partial
x}\left[D\left(x\right)\frac{\partial P}{\partial x}\right]$ (31)

for the probability density $P(x,t)$ on an infinite line $|x|<\infty$, subject
to the initial condition

$P\left(x,t=0\right)=\delta\left(x\right)\,.$ (32)

The piecewise constant diffusivity $D(x)$ is given by

$D\left(x\right)=\begin{cases}D_{1}\,,&x<\ell\,,\\\
D_{2}\,,&x>\ell\,,\end{cases}$ (33)

so we can solve the problem separately in the two regions,

$\frac{\partial P}{\partial t}=\begin{cases}D_{1}\frac{\partial^{2}P}{\partial
x^{2}}\,,&-\infty<x<\ell\,,\\\ D_{2}\frac{\partial^{2}P}{\partial
x^{2}}\,,&\ell<x<\infty\,,\end{cases}$ (34)

and properly match the solutions at $x=\ell$. Let us introduce the functions
$f(x,t)$ and $g(x,t)$ as follows:

$P\left(x,t\right)=\begin{cases}f\left(x,t\right)\,,&x<\ell\,,\\\
g\left(x,t\right)\,,&x>\ell\,.\end{cases}$ (35)

The continuity of the probability density $P(x,t)$ and of the probability flux
$-D(x)\frac{\partial P(x,t)}{\partial x}$ at $x=\ell$ yield two conditions:

$f\left(\ell,t\right)=g\left(\ell,t\right)\,,\quad D_{1}\frac{\partial
f}{\partial x}\Big{|}_{x=l}=D_{2}\frac{\partial g}{\partial
x}\Big{|}_{x=l}\,.$ (36)

The problem can be solved by Laplace transform. Alternatively, one can look
for the solutions in the two regions as sums of two rescaled and shifted
Gaussians. As one can check explicitly, the resulting solution (which is not
new, see e.g. Ref. Farago2020 ) is the following:

$P\left(x,t\right)\\!=\\!\\!\begin{cases}\frac{1}{\sqrt{4\pi
D_{1}t}}\left[\\!e^{-\frac{x^{2}}{4D_{1}t}}\\!+\\!\frac{\sqrt{D_{1}}-\sqrt{D_{2}}}{\sqrt{D_{1}}+\sqrt{D_{2}}}e^{-\frac{\left(x-2l\right)^{2}}{4D_{1}t}}\right]\\!,&\\!x<\ell,\\\
\frac{1}{\sqrt{\pi
t}\left(\sqrt{D_{1}}+\sqrt{D_{2}}\right)}e^{-\frac{\left[x-\ell\left(1-\sqrt{\frac{D_{2}}{D_{1}}}\right)\right]^{2}}{4D_{2}t}}\\!,&\\!x>\ell.\end{cases}$
(37)

In Sec. II we compare the exact result (37) for $x>\ell$ with our prediction
from geometrical optics, see Eq. (11).

## References

  * (1) A. Grosberg and H. Frisch, J. Phys. A: Math. Gen. 36, 8955 (2003).
  * (2) N. Ikeda and H. Matsumoto, in Memoriam Marc Yor – Séminaire de Probabilités XLVII (Lecture Notes in Mathematics, vol 2137), edited by C. Donati-Martin et al (Cham: Springer, Berlin, 2015), p. 497.
  * (3) K. Basnayake, A. Hubl, Z. Schuss, and D. Holcman, Phys. Lett. A 382, 3449 (2018).
  * (4) S. Nechaev, K. Polovnikov, S. Shlosman, A. Valov, and A. Vladimirov, Phys. Rev. E 99 012110 (2019).
  * (5) N. R. Smith and B. Meerson, J. Stat. Mech. (2019) 023205.
  * (6) B. Meerson, J. Stat. Mech. (2019) 013210.
  * (7) B. Meerson and N. R. Smith, Phys. A: Math. Theor. 52, 415001 (2019).
  * (8) B. Meerson, Int. J. Mod. Phys. B 33, 1950172 (2019).
  * (9) T. Agranov, P. Zilber, N. R. Smith, T. Admon, Y. Roichman, and B. Meerson, Phys. Rev. Res. 2, 013174 (2020).
  * (10) A. Vladimirov, S. Shlosman, and S. Nechaev, Phys. Rev. E 102, 012124 (2020).
  * (11) S. N. Majumdar and B. Meerson, J. Stat. Mech. (2020) 023202.
  * (12) B. Meerson, J. Stat. Mech. (2020) 103208.
  * (13) S. Nechaev and A. Valov, J. Phys. A: Math. Theor. 54, 465001 (2021).
  * (14) D. Coombs, R. Straube, and M.Ward, SIAM J. Appl. Math. 70, 302 (2009).
  * (15) M. Eisenbach and L. C. Giojalas, Nat. Rev. Mol. Cell Biol. 7, 276 (2006).
  * (16) P. C. Bressloff and J. M. Newby, Rev. Mod. Phys. 85, 135 (2013).
  * (17) I. M. Sokolov, Chem. Phys. 375, 359 (2010).
  * (18) A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys. 15, 20220 (2013).
  * (19) O. Farago and N. Grønbech-Jensen, J. Stat. Phys. 156, 1093 (2014).
  * (20) F. S. Crawford Jr., Waves (Berkeley Physics Course, Vol. 3) (Mcgraw-Hill, New York, 1968).
  * (21) R. P. Feynman, R. B. Leighton, and M. Sands, The Feynman Lectures on Physics, vol. 1 (Addison-Wesley, Reading, MA, 1963), Chapter 27.
  * (22) Integrability is preserved if one can perform a coordinate transformation so that $D$ depends only on one of the new coordinates.
  * (23) M. Tabor, Chaos and Integrability in Nonlinear Dynamics (Wiley, New York, 1989).
  * (24) D. Nickelsen and H. Touchette, Phys. Rev. Lett. 121, 090602 (2018).
  * (25) B. Meerson, Phys. Rev. E 100, 042135 (2019).
  * (26) B. Meerson, Phys. Rev. E 105, 034106 (2022).
  * (27) B. Meerson and G. Oshanin, Phys. Rev. E 105, 064137 (2022).
  * (28) B. Meerson, Phys. Rev. E 107, 064122 (2023).
  * (29) M. D. Donsker and S. R. S. Varadhan, Commun. Pure Appl. Math. 28, 1 (1975); 28, 279 (1975); 29, 389 (1976); 36, 183 (1983).
  * (30) J. Gärtner, Theory Probab. Appl. 22, 24 (1977); R. S. Ellis, Ann. Probab. 12, 1 (1984).
  * (31) A. Vezzani, E. Barkai, and R. Burioni, Sci. Rep. 10, 2732 (2020).
  * (32) N. R. Smith, Phys. Rev. E 105, 014120 (2022).
  * (33) O. Farago, J. Comp. Phys. 423, 109802 (2020).