main.tex

# Charge density model for the interaction of molecules with vortex beams

Mikhail Maslov<EMAIL_ADDRESS>Institute of Science and Technology
Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria  Georgios M.
Koutentakis<EMAIL_ADDRESS>Institute of Science and
Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria  Mateja
Hrast  Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400
Klosterneuburg, Austria  Oliver H. Heckl  Christian Doppler Laboratory for
Mid-IR Spectroscopy and Semiconductor Optics, Faculty Center for Nano
Structure Research, Faculty of Physics, University of Vienna, Boltzmanngasse
5, 1090 Vienna, Austria  Mikhail Lemeshko<EMAIL_ADDRESS>Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400
Klosterneuburg, Austria

###### Abstract

The interaction of molecules with the orbital angular momentum of light has
long been argued to benefit structural studies and quantum control of
molecular ensembles. We derive a general description of the light-matter
interaction in terms of the coupling between spherical gradients of the
electric field and an effective molecular charge density that exactly
reproduces molecular multipole moments. Our model can accommodate for an
arbitrary complexity of the molecular structure and is applicable to any
electric field, with the exception of tightly focused beams. Within this
framework, we derive the general mechanism of angular momentum exchange
between the spin and orbital angular momenta of light, molecular rotation and
its center-of-mass motion. We demonstrate that vortex beams strongly enhance
certain ro-vibrational transitions that are considered forbidden in the case
of a non-helical light. Finally, we discuss the experimental requirements for
the observation of novel transitions in state-of-the-art spatially-resolved
spectroscopy measurements.

##  I Introduction

Ro-vibrational spectroscopy is one of the most precise non-invasive methods
for studying the properties of molecules [1, 2, 3]. It is used in biochemistry
[4, 5, 6], medicine [7, 8], studies of the atmosphere [9, 10] and interstellar
medium [11]. In molecular physics, it has been applied, e.g., to identify
possible isomer structures of water [12, 13], study van der Waals interactions
[14, 15] and hydrogen bonds [16, 17]. In a typical infrared spectrum,
molecular rotational transitions are constrained by selection rules that stem
from the angular momentum (AM) conservation. Each photon of the probe light
carries at most one quantum of the spin AM [18, 19], which limits the range of
transitions that can be observed using a single-photon absorption
spectroscopy.

This range can be extended either by employing multi-photon schemes [20] or by
increasing the AM carried by the single photon. In particular, in addition to
the spin AM, light can also have orbital angular momentum (OAM) in the form of
the helical beam phase [18, 19]. Twisted photons with large values of OAM can
be generated by combining fundamental laser modes [21, 22]: using a
diffraction grating [23], a spiral phase plate [24] or a metasurface [25].
Outside the regime, where the spin and OAM are strongly coupled [26], they
represent two conceptually different characteristics of light. The former is
related to the polarization of the electric field, while the latter is induced
by its spatial gradient [18, 19]. Molecular spectroscopy with a non-twisted
(spin) light is a highly developed research field [1, 2, 3], but the
interaction of molecules with the OAM is a subject of only recent theoretical
[27, 28, 29, 30, 31, 32, 33] and experimental [34, 35] research.

Aside from non-trivially modifying selection rules and extending the
capabilities of ro-vibrational spectroscopy, vortex beams may also be a
promising tool for quantum control and metrology. For instance, they could
enable the creation of complex rotational molecular states using a single
laser instead of a convoluted multiple-beam setup. Relatively small gas-phase
molecules offer several advantages for such applications, compared to readily
available condensed matter systems, based, e.g., on liquid crystals [36],
micro and nanoparticles [35, 37], metamaterials [38, 39], molecular excitons
[40] and microwave plasmonic resonators [41]. First, as opposed to fabricated
systems, all molecules of the same kind are fundamentally identical, which
offers an advantage in using them as qubits [42]. On the other hand, gas-phase
molecules are almost entirely decoupled from their environment, which paves
the way to cool them to ultracold temperatures while mantaining full control
over their individual quantum states [43, 44], as well as to make use of very
narrow spectroscopic linewidths for precision measurements and fundamental
tests [45, 46].

Many of the previous theories on the interaction of molecules with vortex
beams [27, 28, 29, 30] consider a system of point charges that emulates the
molecule. There are two major drawbacks to this approach. First, the required
number of effective charges and their amplitudes are difficult to justify. In
particular, one needs to account for the redistribution of the electron
density resulting from chemical bonding [47, 48], as well as the shielding of
the nuclear charge by core electrons [49, 50, 51]. The second drawback is the
mathematical complexity and narrow applicability of the resulting approach.
Various implementations of point-charge models [27, 28, 29, 30] cannot be
straightforwardly extended to more complex molecules, e.g., tri-atomic
molecules, like water, which is argued to have rich physics associated with
its bending mode [52, 53].

To overcome these limitations, we model the molecule with a continuous charge
distribution. Similar to the Wigner-Eckart theorem in atomic spectroscopy
[54], our model absorbs most of the theoretical complexity of the molecular
structure into macroscopic molecular multipole moments (dipole, quadrupole,
etc). They can be measured experimentally [55, 56, 57, 58, 59] or extracted
from quantum chemistry calculations [60, 61, 62, 63, 64]. The charge density
is set to exactly reproduce the multipole moments, while embedding various
internal molecular degrees of freedom. For simplicity, we consider only
molecular rotations and a single vibrational mode, which corresponds to the
case of a diatomic molecule. Nevertheless, our framework can scale to more
complicated molecular geometries. It may also be extended to describe
molecular chirality, thus being relevant for the recent studies on resolving
enantiomers using a twisted light [32, 33, 34, 35, 65]. Finally, our model
clearly illustrates the transfer of the spin and OAM of light to the internal
and external (center-of-mass) molecular degrees of freedom.

The paper is organized as follows. In Sec. II, we introduce the molecular
charge density model. We explain how various degrees of freedom can be
embedded into the model. In Sec. III, we derive the general expression for the
light-matter interaction that, provided the validity of the multipole
expansion, is valid for any molecular structure and any electric field
profile, with the exception of tightly focused beams. Our analysis of this
expression reveals general selection rules for the angular momentum exchange.
In Sec. IV, we consider the Laguerre-Gaussian electric field profile, and
calculate contributions of various molecular rotational transitions to the
total Hamiltonian. We infer scenarios in which vortex beams offer an advantage
over the non-twisted light. In Sec. V, we calculate transition rates of ro-
vibrational transitions. In Sec. VI, we suggest an experimental scheme that
may be capable of measuring the enhancement of ro-vibrational transitions,
induced by the OAM of light. We conclude the paper with the discussion of
possible applications and extensions of our model in Sec. VII.

##  II Charge density model

Reasonably sized molecules behave like point sources when interacting with
infrared or optical radiation. This results from the separation of
characteristic length scales of the electric field and molecular charge
density. For instance, in the case of the carbon monosulfide (CS), the
wavelength corresponding to the lowest vibrational transition $(\nu=0\to 1)$
is $\lambda_{\text{vib}}\approx 8\,\text{µm}$, while the bond length is
$a\approx 1.5\,\text{\r{A}}$ [59]. This large difference enables the well-
defined multipole expansion [66, 67]. In particular, one can assume that the
electric field interacts with multipole moments of the whole molecule rather
than individual constituent charges. This allows us to pay less attention to
microscopic details of the molecule, which have a minor effect on its
interaction with radiation. Instead, we develop an effective description of
the molecule by defining the dependence of its multipole moments on the
internal molecular degrees of freedom.

We consider a reference frame, co-rotating with the molecule, with the origin
at its center of mass (herewith called the molecular frame). In this frame, we
suggest the effective model for the physical molecular charge density
$\rho_{\text{mol}}(\bm{r}^{\prime})$. For simplicity, our implementation
describes only the dependence of $\rho_{\text{mol}}(\bm{r}^{\prime})$ on
vibrations and rotations of the molecule and omits other field-induced
displacements of intramolecular charges, e.g., it is adiabatic with respect to
electronic transitions. Nevertheless, other degrees of freedom can be embedded
into the model following our example. For instance, in Sec. VII, we discuss a
possible extension of our model to describe the molecular chirality.

We consider a charge distribution $\rho(\bm{r}^{\prime})$ on a sphere of an
infinitesimally small radius $\chi\to 0$, centered at the molecular frame
origin (see Fig. 1(a)). One can write down the expansion of
$\rho(\bm{r}^{\prime})$ in terms of real-valued spherical harmonics

$\rho_{\text{mol}}(\bm{r}^{\prime})\sim\rho(\bm{r}^{\prime})=\sum\limits_{\lambda,\mu}\alpha_{\lambda,\mu}\mathcal{Y}_{\lambda,\mu}(\Omega^{\prime})\delta(r^{\prime}-\chi)/\chi^{\lambda+2}\,,$
(1)

where the coordinates
$\bm{r}^{\prime}=\\{r^{\prime},\Omega^{\prime}\\}=\\{r,\theta^{\prime},\phi^{\prime}\\}$
and real-valued harmonics $\mathcal{Y}_{\lambda,\mu}(\Omega^{\prime})$ are
defined as

$\mathcal{Y}_{\lambda,\mu}(\Omega^{\prime})=\begin{cases}i\big{(}Y_{\lambda,\mu}(\Omega^{\prime})-Y^{\ast}_{\lambda,\mu}(\Omega^{\prime})\big{)}/\sqrt{2}&\mu<0\\\
Y_{\lambda,0}(\Omega^{\prime})&\mu=0\\\
\big{(}Y_{\lambda,-\mu}(\Omega^{\prime})+Y^{\ast}_{\lambda,-\mu}(\Omega^{\prime})\big{)}/\sqrt{2}&\mu>0\end{cases}\,,$
(2)

where $Y_{\lambda,\mu}(\Omega^{\prime})$ are complex-valued harmonics with the
Condon-Shortley phase [68].

The choice of charge distribution (1) is not accidental. The harmonics
$\mathcal{Y}_{\lambda,\mu}(\Omega^{\prime})$ of a degree $\lambda$ provide a
complete basis for the irreducible representation of the SO$(3)$ group of the
dimension $(2\lambda+1)$ [68, 69, 70]. Therefore, the expansion coefficients
$\alpha_{\lambda,\mu}$, called spherical multipole moments of the order
$\lambda$ [66], reflect point group symmetry of the molecule [71]. In an
experiment, it is more practical to define axial multipole moments, like
$d_{x}$, $d_{y}$ and $d_{z}$. They can be mapped to $\alpha_{\lambda,\mu}$ by
integrating the charge density $\rho(\bm{r}^{\prime})$ with the Cartesian
tensor of a rank $\lambda$. For instance, for the given dipole moment vector
$\bm{d}=\\{d_{x},d_{y},d_{z}\\}$, the condition

$\begin{pmatrix}d_{x}\\\ d_{y}\\\ d_{z}\end{pmatrix}=\lim_{\chi\to
0}\int\rho(\bm{r}^{\prime})\begin{pmatrix}x^{\prime}\equiv
r^{\prime}\sin{\theta^{\prime}}\cos{\phi^{\prime}}\\\ y^{\prime}\equiv
r^{\prime}\sin{\theta^{\prime}}\sin{\phi^{\prime}}\\\ z^{\prime}\equiv
r^{\prime}\cos{\theta^{\prime}}\end{pmatrix}\mathrm{d}^{3}\bm{r}^{\prime}\,,$
(3)

yields:
$\\{\alpha_{1,-1},\alpha_{1,0},\alpha_{1,1}\\}=\sqrt{\frac{3}{4\pi}}\\{d_{y},d_{z},d_{x}\\}$.
Similar mappings can be derived for all higher-order multipole moment tensors,
e.g., quadrupole, as shown in Fig. 1(b). When calculating the corresponding
integrals, which are similar to Eq. (3), the denominator $1/\chi^{\lambda+2}$
in Eq. (1) balances out the Cartesian tensor and Jacobian and assures that the
integrals are finite in the limit $\chi\to 0$.

Besides describing the molecular point group symmetry [71], multipole moments
also reflect the longitudinal geometry of the molecule [60], including its
vibrations. In particular, numerous studies have shown that the dipole [55,
56] and quadrupole [64] moments are different in the ground and excited
vibrational states. Unlike the dependence of multipoles on molecule’s
rotational symmetry, their vibrational dependence is difficult to formalize in
the case of an arbitrary polyatomic molecule. One possible approach is to
model vibrations with the modes of the quantum harmonic oscillator and define
transition multipole moments between vibrational states. If one considers a
molecule with a single vibrational degree of freedom, e.g., diatomic molecule,
multipole moments can be expanded in Taylor series with respect to the
vibrational coordinate $\hat{q}$. For instance, for the quadrupole moment
matrix $Q_{i,j}$, the expansion reads

$\displaystyle Q_{i,j}(\hat{q})$
$\displaystyle=Q_{i,j}(\hat{q}=0)+\bigg{[}\frac{\partial Q_{i,j}}{\partial
q}\bigg{]}_{\hat{q}=0}\hat{q}+...$
$\displaystyle=\sum\limits_{\nu,\nu^{\prime}}Q_{i,j}^{\nu,\nu^{\prime}}(\hat{a}^{\dagger})^{\nu^{\prime}}(\hat{a})^{\nu}\,,$
(4)

where $i,j\in\\{x,y,z\\}$, $\nu^{(\prime)}\in\mathbb{Z}^{+}$ and
$\hat{q}=q_{0}(\hat{a}^{\dagger}+\hat{a})$, where $\hat{a}^{(\dagger)}$ is the
(creation) annihilation operator of the vibrational mode and $q_{0}$ is the
characteristic distance. The coefficient $Q_{i,j}^{\nu,\nu^{\prime}}$ is the
quadrupole moment matrix associated with the transition from the vibrational
state $\ket{\nu}$ to the state $\ket{\nu^{\prime}}$. This approach can be
straightforwardly generalized to the case of multiple contributing vibrational
modes $\hat{q}_{n}$.

Figure 1: Summary of the effective charge density model. (a) Physical
molecular charge distribution $\rho_{\text{mol}}(\bm{r}^{\prime})$ can be
emulated by the charge density $\rho(\bm{r}^{\prime})$ on a sphere with radius
$\chi\to 0$ in the molecular frame (left). The corresponding laboratory-frame
charge density $\rho(\bm{r},\hat{\Omega}_{\text{mol}})$ can be obtained using
the rotational operator $\hat{\mathcal{D}}(\hat{\Omega}_{\text{mol}})$, where
$\hat{\Omega}_{\text{mol}}$ is the molecular orientation operator (right). (b)
Spherical multipole moments $\alpha_{\lambda,\mu}$ – coefficients in the
expansion of $\rho(\bm{r}^{\prime})$ over real-valued spherical harmonics
$\mathcal{Y}_{\lambda,\mu}(\Omega^{\prime})$, can be directly mapped to axial
multipole moments, as exemplified here by the quadrupole moment matrix
$Q_{i,j}$, where ${i,j\in\\{x,y,z\\}}$. (c) Transition quadrupole moments
$Q_{i,j}^{\nu,\nu^{\prime}}$ are calculated numerically using density
functional theory. For the CS molecule, we obtain the adiabatic potential
energy curve $E^{\text{vib}}(q)$ (markers) and the associated eigenenergies
$E^{\text{vib}}_{\nu}$ (dashed lines) and eigenstates $\psi_{\nu}(q)$
(shades), where $q$ is the vibrational coordinate (main panel). Both the
spectrum and eigenstates are asymmetric with respect to the equilibrium
internuclear distance $q_{\text{eq}}$. We calculate the coordinate dependence
of the quadrupole moment $Q_{zz}(q)$ (inset). To obtain the transition
quadrupole moment, we calculate the expectation value of $Q_{zz}(q)$ with
respect to the eigenstates $\psi_{\nu(\nu^{\prime})}(q)$ (see Eq. (5)). If the
potential energy curve $E^{\text{vib}}(q)$ is only slightly anharmonic,
$Q_{zz}^{\nu,\nu^{\prime}}$ can be approximated with the
$(\nu+\nu^{\prime})$’th derivative of $Q_{zz}(q)$ at $q=0$.

Another approach to obtaining the transition multipole moments is to use
various numerical ab-initio methods [61, 62]. In Fig. 1(c), we present our
calculation for the CS molecule based on the StoBe-deMon implementation [72]
of density functional theory [73]. By using the local Perdew-Wang exchange-
correlation potential of Ref. [74] and constraining the Kohn-Sham orbitals
[75, 73] to be rotationally invariant around the interatomic axis, we obtain
the adiabatic potential energy curve $E^{\text{vib}}(q)$ (markers). Here the
vibrational coordinate $q$ is the shift of the internuclear distance from its
equilibrium value $q_{\text{eq}}\approx 1.552\,\text{\AA}$. The resulting
curve can be reasonably approximated by the anharmonic Morse potential (solid
line). We also calculate the vibrational levels $E^{\text{vib}}_{\nu}$ (dashed
lines) and the corresponding wave functions $\psi_{\nu}(q)$ (shades), which
are asymmetric with respect to the equilibrium coordinate $q_{\text{eq}}$. In
the inset of Fig. 1(c), we plot the dependence of the quadrupole moment
$Q_{zz}$ on the vibrational coordinate $q$ (markers), which has non-linear
profile (cf. solid line). The transition quadrupole moment matrix can be
straightforwardly obtained by calculating the integral

$Q_{i,j}^{\nu,\nu^{\prime}}=\int\psi^{\ast}_{\nu^{\prime}}(q)Q_{i,j}(q)\psi_{\nu}(q)\mathrm{d}q\,.$
(5)

In the case of a harmonic energy spectrum $E^{\text{vib}}(q)$ and symmetric
wave functions $\psi_{\nu(\nu^{\prime})}(q)$, the transition quadrupole moment
can be approximated, in the first order, by the derivative:
$Q_{i,j}^{\nu,\nu^{\prime}}\approx\bigg{[}\dfrac{\partial^{\nu+\nu^{\prime}}Q_{i,j}}{\partial
q^{\nu+\nu^{\prime}}}\bigg{]}_{q=0}$.

The dependence of the charge density $\rho(\bm{r}^{\prime})$ on the molecular
rotation becomes apparent after casting the molecular-frame spherical
harmonics (with angle $\Omega^{\prime}$) in terms of the laboratory-frame
angles $\Omega$

$\displaystyle Y_{\lambda,\mu}(\Omega^{\prime})=$
$\displaystyle\hat{\mathcal{D}}\big{(}\hat{\Omega}_{\text{mol}}\big{)}Y_{\lambda,\mu}(\Omega)$
$\displaystyle=\sum\limits_{\zeta}D^{\lambda}_{\zeta,\mu}\big{(}\hat{\Omega}_{\text{mol}}\big{)}Y_{\lambda,\zeta}(\Omega)\,,$
(6)

where $\hat{\mathcal{D}}\big{(}\hat{\Omega}_{\text{mol}}\big{)}$ is the
rotation operator (see Fig. 1(a)),
$D^{\lambda}_{\zeta,\mu}\big{(}\hat{\Omega}_{\text{mol}}\big{)}$ are Wigner
$D$-matrices [68] and $\hat{\Omega}_{\text{mol}}$ is the molecular orientation
operator.

Finally, to illustrate our model better, we consider a general diatomic
heteronuclear molecule with the dipole moment vector
$\bm{d}=\\{0,0,\hat{d}^{\text{vib}}_{z}\\}$ and quadrupole moment matrix
$Q_{i,j}=-\hat{Q}^{\text{vib}}_{z,z}\,\text{diag}\\{1/2,1/2,-1\\}$, where
$\hat{d}^{\text{vib}}_{z}$ and $\hat{Q}^{\text{vib}}_{z,z}$ are operators with
respect to molecular vibrational states. Omitting all the higher multipole
moments, the effective charge density reads

$\displaystyle\rho_{\text{diatomic}}\big{(}\bm{r},\hat{d}^{\text{vib}}_{z},\hat{Q}^{\text{vib}}_{zz},\hat{\Omega}_{\text{mol}}\big{)}$
$\displaystyle=\sum\limits_{\mu}\hat{d}^{\text{vib}}_{z}Y_{1,\mu}\big{(}\hat{\Omega}_{\text{mol}}\big{)}Y_{1,-\mu}(\Omega)\frac{\delta(r-\chi)}{\chi^{3}}$
$\displaystyle+\sum\limits_{\mu}\frac{3}{2}\hat{Q}^{\text{vib}}_{z,z}Y_{2,\mu}\big{(}\hat{\Omega}_{\text{mol}}\big{)}Y_{2,-\mu}(\Omega)\frac{\delta(r-\chi)}{\chi^{4}}\,.$
(7)

Multipole moment tensors $\bm{d}$ and $Q_{i,j}$ have a more complicated
structure for polyatomic (non-linear) molecules. Moreover, they depend on more
than one vibrational degree of freedom. Modelling such molecules is an active
challenge in quantum chemistry and is outside the scope of this manuscript.
Nevertheless, provided that the analytical or numerical dependence of
multipole moment tensors on internal degrees of freedom is known, and provided
that the multipole expansion is valid, our charge density model is valid for
any molecular geometry.

##  III Light-matter interaction and angular momentum exchange

After introducing the effective description of molecular properties, we
proceed to define the Hamiltonian of the light-matter interaction.

Our starting point is the non-relativistic multipolar QED Hamiltonian, which
can be obtained from the Coulomb gauge Hamiltonian by performing the unitary
Power-Zienau-Wooley (PZW) transformation [76, 77, 78]. Both the aforementioned
Hamiltonians are fully identical for the on-energy-shell processes, while PZW
Hamiltonian gives, arguably, a better representation for the off-energy-shell
processes [76, 79]. We consider a nonmagnetic molecule, thus we can safely
omit its interaction with the magnetic field further on. The resulting
Hamiltonian reads: $\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{\text{int}}$,
where the light-matter interaction can be expressed as

$\mathcal{H}_{\text{int}}=-\int\big{(}\bm{\mathcal{P}}(\bm{r}_{\text{EF}})\cdot\bm{E}(\bm{r}_{\text{EF}})\big{)}\mathrm{d}^{3}\bm{r}_{\text{EF}}+\text{H.c.}\,,$
(8)

where $\bm{r}_{\text{EF}}$ is a coordinate in the laboratory frame. The
Hamiltonian $\mathcal{H}_{0}$ describes the “unperturbed” system, i.e., the
electric field and molecule in the non-interacting regime. In particular, it
includes the kinetic energy of the center of mass, rotational and vibrational
molecular energies, which are discussed in Sec. VI, as well as the energy of
the electric field

$\mathcal{E}_{\text{EF}}=\frac{\varepsilon_{0}}{2}\int|\bm{E}(\bm{r}_{\text{EF}})|^{2}\mathrm{d}^{3}\bm{r}_{\text{EF}}\,,$
(9)

where $\varepsilon_{0}$ is the vacuum permittivity. Following the common
practice, we also put the divergent integral
$\int|\bm{\mathcal{P}}(\bm{r}_{\text{EF}})|^{2}\mathrm{d}^{3}\bm{r}_{\text{EF}}$
into $\mathcal{H}_{0}$. This particular divergence, its origin and remedy are
thoroughly discussed in Ref. [79].

To define the polarization field $\bm{\mathcal{P}}(\bm{r}_{\text{EF}})$, we
begin with the polarization
$\bm{\mathcal{P}}_{\text{PZW}}(\bm{r}_{\text{EF}})$ of the PWZ Hamiltonian. It
provides the exact description of the generic system of $N$ point-charges

$\bm{\mathcal{P}}_{\text{PZW}}(\bm{r}_{\text{EF}})=\sum\limits_{i=1}^{N}q_{i}\bm{r}_{i}\int\limits_{0}^{1}\mathrm{d}\eta\,\delta\big{(}\bm{r}_{\text{EF}}-\bm{R}-\eta\bm{r}_{i}\big{)}\,,$
(10)

where $q_{i}$ and $\bm{r}_{i}$ are the charges and their coordinates in the
center-of-mass reference frame, $\bm{R}$ is the center-of-mass position in the
laboratory frame. This relation can be straightforwardly generalized to
describe the polarization field $\bm{\mathcal{P}}(\bm{r}_{\text{EF}})$,
generated by the continuous charge distribution $\rho(\bm{r})$, introduced in
Sec. II

$\bm{\mathcal{P}}(\bm{r}_{\text{EF}})=\lim_{\chi\to
0}\int\limits\mathrm{d}^{3}\bm{r}\,\rho(\bm{r})\bm{r}\int\limits_{0}^{1}\mathrm{d}\eta\,\delta\big{(}\bm{r}_{\text{EF}}-\bm{R}-\eta\bm{r}\big{)}\,.$
(11)

Substituting the polarization field into Eq. (8), we obtain

$\mathcal{H}_{\text{int}}=-\lim_{\chi\to
0}\int\mathrm{d}^{3}\bm{r}\,\rho(\bm{r})\int\limits_{0}^{1}\mathrm{d}\eta\,\big{(}\bm{r}\cdot\bm{E}(\bm{R}+\eta\bm{r})\big{)}+\text{H.c.}\,$
(12)

We consider the electric field of a non-tightly focused beam. It allows us to
use the small-angle approximation and split the electric field into a
polarization vector and a spatial field profile:
${\bm{E}(\bm{R}+\eta\bm{r})=\bm{E}_{0}\,E(\bm{R}+\eta\bm{r})}$. In the case of
a tightly focused beam, the spin and OAM of light, associated respectively
with its polarization and spatial profile, are strongly coupled [80, 26],
which complicates the analysis of the OAM transfer to the molecule.
Nevertheless, tight focusing was shown to aid the OAM transfer, as discussed
in Sec. VII, rendering the description of such light sources an important
future extension to our theory. Without the loss of generality, we consider a
beam propagating along the $z$-axis. For such a beam, $E_{0,z}=0$, which
yields: $(\bm{r}\cdot\bm{E}_{0})=\sqrt{4\pi/3}\,r\sum_{\sigma=\pm
1}\epsilon_{\sigma}Y_{1\sigma}(\Omega_{\bm{r}})$, where
$\epsilon_{\sigma}=(E_{0,x}+i\sigma E_{0,y})/\sqrt{2}$ expresses $\bm{E}_{0}$,
in terms of the circular polarization $\sigma=\pm 1$.

To proceed, we expand the spatial electric field profile
$E(\bm{R}+\eta\bm{r})$ around the center of mass position $\bm{R}$, where
$r\ll R$. Unlike the previous studies [27, 28, 29, 30], that considered the
Cartesian Taylor expansion, we employ the spherical expansion [81, 82, 83]

$\displaystyle
E(\bm{R}+\eta\bm{r})=\exp(\eta\bm{r}\cdot\bm{\nabla}_{\bm{R}})E(\bm{R})$
$\displaystyle=\sum\limits_{nlm}c_{n,l}(\eta
r)^{2n+l}Y^{\ast}_{l,m}(\Omega_{\bm{r}})\big{[}\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}\,,$
(13)

where $c_{n,l}=\frac{\pi 2^{l+2}\kappa^{2n}(l+n)!}{n!(2l+2n+1)!}$, $n,l\geq
0$, $m\in[-l,l]$ and $\kappa$ is the wavenumber of the electric field.
$\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})$ are the solid harmonics of the
gradient operator $\bm{\nabla}_{\bm{R}}$, also known as spherical tensor
gradient operators. Their detailed overview is the subject of Refs. [82, 81,
83]. For the sake of simplicity, we would refer to the term
$\big{[}\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}$ as the
spherical gradient. For differentiable electric field profiles $E(\bm{R})$ and
specific values of $l$ and $m$, it can be calculated analytically, after
expanding the solid harmonics in Cartesian coordinates

$\displaystyle\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})=\sqrt{\frac{(2l+1)(l+m)!(l-m)!}{4\pi}}$
$\displaystyle\times\sum\limits_{k}\frac{(-\partial_{X}-i\partial_{Y})^{m+k}(\partial_{X}-i\partial_{Y})^{k}\partial_{Z}^{l-m-2k}}{2^{m+2k}(m+k)!k!(l-m-2k)!}\,,$
(14)

where $\partial_{X}$, $\partial_{Y}$ and $\partial_{Z}$ are the Cartesian
components of the gradient vector $\bm{\nabla}_{\bm{R}}$, $\max(-m,0)\leq
k\leq\lfloor\frac{l-m}{2}\rfloor$, with $\lfloor x\rfloor$ being the floor
function yielding the largest integer less than $x$.

We further substitute the electric field expansion (13) into the Hamiltonian
(12). We consider the molecular-frame charge density, given by Eq. (1), and
apply the rotation rule (6) to obtain the charge density $\rho(\bm{r})$ in the
laboratory coordinates. After integrating over $\eta$ and $\bm{r}$ and taking
the limit $\chi\to 0$, we obtain the final relation for the interaction
Hamiltonian

$\mathcal{H}_{\text{int}}=\sum\limits_{\begin{subarray}{c}nlm\\\
\sigma\mu\end{subarray}}\gamma_{n,l,m,\sigma}\hat{\alpha}_{\lambda,\mu}\epsilon_{\sigma}\big{[}\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}\Big{|}_{\bm{R}=\hat{\bm{R}}}\begin{cases}D^{\lambda}_{m-\sigma,\mu}\big{(}\hat{\Omega}\big{)}-(-1)^{\mu}D^{\lambda}_{m-\sigma,-\mu}\big{(}\hat{\Omega}\big{)}&\mu<0\\\
\sqrt{2}D^{\lambda}_{m-\sigma,0}(\hat{\Omega})&\mu=0\\\
D^{\lambda}_{m-\sigma,-\mu}\big{(}\hat{\Omega}\big{)}+(-1)^{-\mu}D^{\lambda}_{m-\sigma,\mu}\big{(}\hat{\Omega}\big{)}&\mu>0\end{cases}+\text{H.c.}\,,$
(15)

where
$\gamma_{n,l,m,\sigma}=c_{n,l}C^{\lambda,0}_{l,0;1,0}C^{l,m}_{\lambda,m-\sigma;1,\sigma}/(\lambda\sqrt{2})$,
where the multipole order $\lambda=2n+l+1$ and
$C^{L,M}_{l,m;l^{\prime},m^{\prime}}$ are the Clebsch-Gordan coefficients
[68]. Operators $\hat{\bm{R}}$, $\hat{\alpha}_{\lambda,\mu}$ and
${\hat{\Omega}\equiv\hat{\Omega}_{\text{mol}}}$ act respectively on the
center-of-mass, molecular vibrational and rotational states.

The Hamiltonian (15) is the major result of our model. Provided the validity
of the multipole expansion, it holds for any molecular charge density and any
profile of the electric field, except for tightly focused beams.
$\mathcal{H}_{\text{int}}$ already reveals the angular momentum exchange and
the corresponding selection rules. To illustrate this fact, we consider
cylindrical coordinates ${\bm{R}=\\{R,\Phi,Z\\}}$ and expand the electric
field profile into Fourier series with respect to the polar angle $\Phi$:
$E(\bm{R})=\sum_{M}E_{M}(R,Z)e^{iM\Phi}$. Each Fourier component $E_{M}(R,Z)$
is an eigenvalue of the $\hat{L}_{z,\Phi}=-i\hbar\,\partial/\partial\Phi$
operator and can be associated with the magnetic quantum number M. The
spherical gradient can be straightforwardly calculated, using Eq. (14)

$\big{[}\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}=\sum\limits_{M}\tilde{E}_{M,l}(R,Z)e^{i(M+m)\Phi}\,.$
(16)

As a result of applying the spherical gradient operator, Fourier amplitudes
change: ${E_{M}(R,Z)\to\tilde{E}_{M,l}(R,Z)}$, and magnetic quantum numbers
shift: $M\to M+m$. The non-zero spherical gradient for the certain $l$ and $m$
represents the transfer of $M+m$ quanta of angular momentum from the beam to
the center-of-mass motion $(\hat{\bm{R}})$. The ability to directly associate
a derivative from the Taylor series to the angular momentum transfer is the
main benefit of using the spherical expansion in Eq. (13) instead of the
Cartesian one. Similarly, Wigner D-matrices
$D^{\lambda}_{m-\sigma,\pm\mu}\big{(}\hat{\Omega}_{\text{mol}}\big{)}$
describe the transfer of $\sigma-m$ quanta of AM from the electric field to
the molecular rotation $(\hat{\Omega}_{\text{mol}})$. Together, these
transfers satisfy the angular momentum conservation: $M+\sigma$ quanta of AM
of the electric field get redistributed into $M+m$ quanta of center-of-mass AM
and $\sigma-m$ quanta of molecular AM. Apart from the angular momentum
exchange, Eq. (15) infers two major corollaries.

If one considers only the electric-dipole component of the light-matter
interaction, the only remaining terms in Eq. (15) are those with $\lambda=1$,
thus $n=l=m=0$. The spherical gradient
$\big{[}\mathcal{R}_{0,0}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}\propto
E(\bm{R})$ is independent of summation indices and can be extracted outside
the sum. This means, that the spatial profile of the electric field
$E(\bm{R})$ is simply an amplitude modifier to the interaction Hamiltonian. It
is unable to affect molecular rotational transitions. In this case, molecular
transitions are controlled only by the circular polarization $\sigma$, which
was shown to be associated with the spin angular momentum of light [18, 19].
This result is in agreement with the previous studies [27, 29]. In particular,
for a diatomic molecule, described by the charge density (7), the contribution
to the light-matter interaction from the electric dipole reads

$\mathcal{H}_{\text{int}}^{\text{dip}}=-\sqrt{\frac{4\pi}{3}}\hat{d}^{\text{vib}}_{z}E(\hat{\bm{R}})\sum\limits_{\sigma}\epsilon_{\sigma}Y_{1,\sigma}(\hat{\Omega})\,,$
(17)

which is identical to the result of Ref. [29]. In contrast to this, if one
considers the quadrupole interaction term, with $\lambda=2$, thus $n=0$ and
$l=1$, the Hamiltonian (15) couples the spherical gradient
$\big{[}\mathcal{R}_{1,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}$ to the
D-matrix of the molecular rotation angle $\hat{\Omega}$. For a diatomic
molecule, the electric-quadrupole interaction term is given by

$\displaystyle\mathcal{H}_{\text{int}}^{\text{quad}}=\pi\sqrt{2/3}\,\hat{Q}^{\text{vib}}_{z,z}$
$\displaystyle\times\sum\limits_{m\sigma}\epsilon_{\sigma}C^{1,m}_{2,m-\sigma;1,\sigma}\big{[}\mathcal{R}_{1,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}\Big{|}_{\bm{R}=\hat{\bm{R}}}Y^{\ast}_{2,m-\sigma}(\hat{\Omega})\,.$
(18)

It clearly illustrates the coupling between the center-of-mass coordinate
$\hat{\bm{R}}$ and the molecular rotation angle $\hat{\Omega}$, mediated by
the electric field (via summation index $m$).

##  IV Laguerre-Gaussian beams and enhancement of rotational transitions

To analyze the angular momentum exchange revealed by the Hamiltonian (15) in
more detail, one needs to specify the spatial profile of the electric field
$E(\bm{R})$. We consider an electric field that can be expressed as a
combination of Laguerre-Gaussian (LG) modes. These modes are solutions to the
paraxial Helmholtz equation in cylindrical coordinates
$\bm{R}=\\{R,\Phi,Z\\}$. The magnitude of the electric field in a LG mode is
given by [22]

$\displaystyle
E_{PM}(\bm{R})=\gamma_{PM}\bigg{(}\frac{R}{\omega_{Z}}\bigg{)}^{|M|}\exp[-\frac{R^{2}}{\omega^{2}_{Z}}]e^{-iM\Phi}e^{-i\kappa
Z}$
$\displaystyle\times\Bigg{\\{}\frac{\omega_{0}}{\omega_{Z}}\exp[-i\bigg{(}\frac{\kappa
R^{2}}{2\mathbb{R}_{Z}}-\psi_{Z}\bigg{)}]\mathcal{L}^{|M|}_{P}\bigg{(}\frac{2R^{2}}{\omega^{2}_{Z}}\bigg{)}\Bigg{\\}}\,,$
(19)

where $\kappa=2\pi n/\lambda_{\text{beam}}$ is the wavenumber, $n$ is the
refraction index, $\lambda_{\text{beam}}$ is the wavelength of light, and
$\mathcal{L}^{j}_{i}(x)$ with $i\in\mathbb{Z}^{+},\,j\in\mathbb{R}$ are the
generalized Laguerre polynomials. The beam is focused at $Z=0$, with the waist
function $\omega_{Z}=\omega_{0}\sqrt{1+(Z/Z_{R})^{2}}$ along the $z$-axis,
where $\omega_{0}$ is the waist at the focus and
${Z_{R}=\pi\omega_{0}^{2}n/\lambda_{\text{beam}}}$ is the Rayleigh length. The
radius of the wavefront curvature $\mathbb{R}_{Z}=Z(1+(Z_{R}/Z)^{2})$ and the
Gouy phase $\psi_{Z}=(2P+|M|+1)\atan(Z/Z_{R})$. The electric field (19) is
normalized in the sense of the Dirac delta-function with respect to the axial
coordinate $Z$, and to unity with respect to coordinates $R$ and $\Phi$, with
the normalization constant $\gamma_{PM}=\sqrt{2^{|M|+1}P!/(\pi(P+|M|)!)}$. The
electric field $E_{PM}(\bm{R})$ is an eigenfunction of the $\hat{L}_{z,\Phi}$
operator, defined in Sec. III, with the magnetic quantum number $M$. Eq. (19)
with $M=0$ describes a Gaussian beam without OAM, and serves within this
manuscript as a benchmark representation of the non-helical light.

The action of the spherical gradient on the electric field $E_{PM}(\bm{R})$ is
readily given by Eq. (16), where the transformed Fourier amplitude can be
calculated using Eq. (14). However, the resulting expression is cumbersome and
is therefore omitted here. Instead, to illustrate our findings in a simple
way, we consider an approximation to the electric field (19)

$E^{\text{foc}}_{M}(\bm{R})=\gamma_{0M}\left(R/\omega_{0}\right)^{|M|}e^{-R^{2}/\omega^{2}_{0}}e^{-iM\Phi}e^{-i\kappa
Z}\,.$ (20)

This corresponds to a LG beam with $P=0$ within the in-focus approximation,
which is valid, when $Z\ll Z_{R}$, thus $\omega_{Z}\to\omega_{0}$, $1/R_{Z}\to
0$ and $\psi_{Z}\to 0$. Previous studies [28, 29] have already considered a
similar spatial profile

$E^{\text{foc,cent}}_{M}(\bm{R})=\frac{1}{\sqrt{|M|!}}\left(\frac{R}{\omega_{0}}\right)^{|M|}e^{-iM\Phi}e^{-ikZ}\,,$
(21)

which corresponds to Eq. (20) with the additional in-center approximation,
which is valid when $R\ll\omega_{0}$. The major drawback of the approximation
(21) is the divergence of the electric field at large distances $R$. While
this problem can be solved by introducing a cutoff
$\omega_{\text{cut}}\sim\omega_{0}$, the value of this cutoff was argued to
affect the result of numerical calculations [84]. Moreover, the light-matter
interaction was shown to depend on the average position of molecule’s center
of mass with respect to the cutoff distance $\omega_{\text{cut}}$ [29].
Therefore, to ensure the validity of our findings, we choose the cutoff-free
approximation to the electric field (20).

Using the approximation $E^{\text{foc}}_{M}(\bm{R})$, we examine the OAM
transfer to the molecule. In the Hamiltonian (15), we choose the lowest-
multipole term that couples the gradient of the electric field to the
molecular rotation. Namely, we consider the electric-quadrupole interaction,
i.e., the general case of Eq. (18). We calculate the spherical gradients,
which describe the OAM transfer, i.e., the gradients with $m=\pm 1$, as the
$m=0$ gradient describes purely the spin transfer. The gradients read

$\displaystyle\big{[}\mathcal{R}_{1,\pm
1}(\bm{\nabla}_{\bm{R}})E^{\text{foc}}_{M}(\bm{R})\big{]}=$
$\displaystyle=-\sqrt{\frac{3}{8\pi}}e^{\pm
i\Phi}\bigg{[}\frac{M}{R}\mp\bigg{(}\frac{|M|}{R}-\frac{2R}{\omega_{0}^{2}}\bigg{)}\bigg{]}E^{\text{foc}}_{M}(\bm{R})\,.$
(22)

In what follows, the absolute value of the spherical gradient (22) serves as a
quantitative measure for the OAM transfer. The three terms in the square
brackets of Eq. (22) correspond respectively to derivatives of the spiral beam
phase $(\sim e^{iM\Phi})$, radial distribution for $R\ll\omega_{0}$
$(\sim(R/\omega_{0})^{|M|})$ and exponential decay for $R\gg\omega_{0}$ $(\sim
e^{-R^{2}/\omega_{0}^{2}})$.

One can also employ the approximation $E^{\text{foc,cent}}_{M}(\bm{R})$ to
calculate the spherical gradient. In this case, the term in Eq. (22)
associated with the exponential decay is missing. Hence, for the beam with
$M=0$, the amplitude of the OAM transfer is zero, as also shown in Ref. [29].
The possible conclusion is that a non-helical beam cannot transfer any
additional angular momentum, apart from the spin, to the molecule. This
conclusion, however, does not hold when a more accurate profile, like
$E^{\text{foc}}_{M}(\bm{R})$, is taken into account. As revealed by Eq. (22)
with $M=0$, even a non-vortex Gaussian beam has a small impact on the
molecular rotations, in agreement with Ref. [27]. This effect has a simple
explanation. For instance, in the case of a diatomic molecule, the two nuclei,
in general, experience a slightly different electric field, due to its spatial
gradient, which creates an effective disbalance among different rotational
states.

After showing that the OAM transfer occurs both in the case of a non-twisted
and twisted light, we examine how vortex beams can enhance the transfer,
compared to non-vortex beams. As revealed by Eq. (22), there are two
contributions to the OAM transfer. Namely, the derivative of the spiral beam
phase (first term in square brackets) and the derivative of the radial profile
(round brackets). While the latter contribution is present for a Laguerre-
Gaussian beam with any value of $M$, the former contribution describes
exclusively the OAM transfer stimulated by the non-zero vorticity. This
transfer channel becomes the dominant one in the $R\to 0$ regime. In
particular, the amplitude of the OAM transfer for a beam with $M\neq 0$, which
is $\sim M/R$, can be significantly larger than the transfer amplitude for the
light with $M=0$, which is $\sim R/\omega_{0}^{2}$. This indicates that the
experimental observation of such an enhancement would require the molecules to
be placed around the center of the beam, as further discussed in Sec. VI.

Our numerical analysis also reveals that, in the out-of-focus regime, i.e.,
when $Z\gg Z_{R}$, besides the angular momentum exchange arising due to radial
and phase variation of the electric field, there is an AM exchange induced by
the non-zero curvature of the wavefront. This can be intuitively understood as
follows, the curved wavefront effectively deforms the radial distribution of
the electric field, which, in turn, affects the magnitude of AM transfer.

##  V Ro-vibrational transition rates

To verify the major argument of Sec. IV, namely, that the OAM of light can
significantly enhance certain rotational transitions, we calculate the
transition rates of ro-vibrational transitions.

We assume that the molecular state is the direct product:
${\ket{\Psi}=\ket{\psi^{\text{vib}}_{\nu}}\ket{\psi^{\text{rot}}_{J,N_{1},N_{2}}}\ket{\psi^{\text{CM}}_{\bm{R}_{\text{CM}}}}}$,
i.e., we neglect weak correlation effects, e.g., correlations stemming from
the ro-vibrational coupling. We leave the vibrational state
$\ket{\psi^{\text{vib}}_{\nu}}$ implicit, since we directly calculate the
transition multipole moments using numerical methods, as discussed in Sec. II.
The rotational state $\ket{\psi^{\text{rot}}_{J,N_{1},N_{2}}}$ is an
eigenstate of the kinetic energy of a rigid rotor:
$\mathcal{H}_{0,\text{rot}}=B_{x}\hat{J}_{x}^{2}+B_{y}\hat{J}_{y}^{2}+B_{z}\hat{J}_{z}^{2}$,
where $\hat{\bm{J}}$ is the angular momentum of the molecular rotation. The
corresponding wave function reads

$\bra{\Omega_{\text{mol}}}\ket{\psi^{\text{rot}}_{J,N_{1},N_{2}}}=D^{J}_{N_{1},N_{2}}(\Omega_{\text{mol}})\,,$
(23)

where $\ket{\Omega_{\text{mol}}}$ are eigenstates of the molecular angle
operator $\hat{\Omega}_{\text{mol}}$.

For most molecules at relevant experimental temperatures, the thermal de-
Broglie wavelength is significantly smaller than the wavelength of the probe
light. For instance, for a CS molecule at $T=20\,\text{K}$: the de-Broglie
wavelength is $\lambda_{\text{dB}}=h\sqrt{\beta/(2\pi m)}\approx
0.6\,\text{\AA}$, where $m$ is the molecular mass, $\beta=1/(k_{\text{B}}T)$
is the Boltzmann factor with the Boltzmann constant $k_{\text{B}}$. At the
same time, the wavelength corresponding to the lowest vibrational transition
is $\lambda_{\text{beam}}\sim 8\,\text{µm}$. Provided that, for a non-tightly
focused beam, the beam waist is larger than the wavelength, we obtain:
$\omega_{0}>\lambda_{\text{beam}}\gg\lambda_{\text{dB}}$. Therefore, the wave
function that describes the center-of-mass position inside the beam can be
approximated by the three-dimensional $\delta$-function

$\bra{\bm{R}}\ket{\psi^{\text{CM}}_{\bm{R}_{\text{CM}}}}=\delta^{(3)}(\bm{R}-\bm{R}_{\text{CM}})\,,$
(24)

where $\ket{\bm{R}}$ is an eigenstate of the center-of-mass position operator
$\hat{\bm{R}}$.

The transition amplitude between the initial $\ket{\Psi_{i}}$ and final
$\ket{\Psi_{f}}$ molecular states with respect to the full light-matter
interaction Hamiltonian (15) reads

$\displaystyle\mathcal{M}_{i\to
f}=\bra{\Psi_{f}}\mathcal{H}_{\text{int}}\ket{\Psi_{i}}$
$\displaystyle=\sum\limits_{nlm\mu}\mathcal{I}^{\text{vib};\,n,l,\mu}_{\nu,\nu^{\prime}}\,\mathcal{I}^{\text{CM};\,l,m}_{\bm{R}_{\text{CM}},\bm{R}^{\prime}_{\text{CM}}}\,\mathcal{I}^{\text{rot};\,n,l,m,\mu}_{J,N_{1},N_{2},J^{\prime},N^{\prime}_{1},N^{\prime}_{2}}\,,$
(25)

where states $\ket{\Psi_{i(f)}}$ are characterised by sets of quantum numbers:
$\\{\nu^{(\prime)},\bm{R}^{(\prime)}_{\text{CM}},J^{(\prime)},N_{1}^{(\prime)},N_{2}^{(\prime)}\\}$
and the vibrational integral,
$\mathcal{I}^{\text{vib};\,n,l,\mu}_{\nu,\nu^{\prime}}=\bra{\psi^{\text{vib}}_{\nu^{\prime}}}\alpha_{\lambda,\mu}\ket{\psi^{\text{vib}}_{\nu}}$,
with $\lambda=2n+l+1$, can be expressed via the transition multipole moments,
as discussed in Sec. II. The rotational and center-of-mass integrals are

$\mathcal{I}^{\text{rot};\,n,l,m,\mu}_{J,N_{1},N_{2},J^{\prime},N^{\prime}_{1},N^{\prime}_{2}}=\bigg{[}\sum\limits_{\sigma}\frac{8\pi^{2}}{2\lambda+1}\gamma_{n,l,m,\sigma}\epsilon_{\sigma}C^{J^{\prime},N_{1}^{\prime}}_{\lambda,m-\sigma;J,N_{1}}\bigg{]}\begin{cases}C^{J^{\prime},N_{2}^{\prime}}_{\lambda,\mu;J,N_{2}}-(-1)^{\mu}C^{J^{\prime},N_{2}^{\prime}}_{\lambda,-\mu;J,N_{2}}&\mu<0\\\
\sqrt{2}C^{J^{\prime},N_{2}^{\prime}}_{\lambda,0;J,N_{2}}&\mu=0\\\
C^{J^{\prime},N_{2}^{\prime}}_{\lambda,-\mu;J,N_{2}}+(-1)^{-\mu}C^{J^{\prime},N_{2}^{\prime}}_{\lambda,\mu;J,N_{2}}&\mu>0\end{cases}\,,$
(26)
$\mathcal{I}^{\text{CM};\,l,m}_{\bm{R}_{\text{CM}},\bm{R}^{\prime}_{\text{CM}}}=\bra{\psi^{\text{CM}}_{\bm{R}^{\prime}_{\text{CM}}}}\big{[}\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}})E(\bm{R})\big{]}\Big{|}_{\bm{R}=\hat{\bm{R}}}\ket{\psi^{\text{CM}}_{\bm{R}_{\text{CM}}}}=\delta^{(3)}(\bm{R}_{\text{CM}}-\bm{R}_{\text{CM}}^{\prime})\big{[}\mathcal{R}_{l,m}(\bm{\nabla}_{\bm{R}_{\text{CM}}})E(\bm{R}_{\text{CM}})\big{]}\,.$
(27)

Provided the transition amplitude $\mathcal{M}_{i\to f}$, we obtain the
transition rate using Fermi’s golden rule

$\Gamma_{i\to f}(\omega)=\frac{2\pi}{\hbar}\big{|}\mathcal{M}_{i\to
f}\big{|}^{2}\delta_{E}(\hbar\omega-\Delta E_{fi})\,,$ (28)

where $\omega$ is the angular frequency of the photon that drives the
transition and $\Delta E_{fi}=E_{f}-E_{i}$ is the energy difference between
the initial and final states. The $\delta$-function in Eq. (28) denotes the
density of energy states and has therefore the units of inverse energy.

##  VI Possible experimental scheme

After defining transition rates of molecular ro-vibrational transitions, we
proceed with a suggestion for the generic proof-of-principle experimental
scheme. The proposed setup may be capable of revealing the enhancement of
rotational transition amplitudes, induced by the non-zero vorticity of the
probe light.

We consider a vacuum chamber of the characteristic length $L$ (along the
$z$-axis) that contains molecules in the gas phase under pressure $P$ and at
temperature $T$. As discussed in Sec. IV, the OAM-induced rotational
enhancement depends on the position of the molecular center of mass with
respect to the beam axis. For this reason, we suggest measuring the absorbance
of light in a spatially-resolved manner. In particular, we propose measuring
the ratio between the absorbed and incident power for the molecules, which
reside within the optical path of light that is collected by an adjustable
circular aperture $\tilde{\bm{R}}_{0}$, placed outside the chamber and
centered on the beam axis. For small chambers $(L<Z_{R})$, such molecules are
approximately confined to the cylinder, described by the effective aperture
$\bm{R}_{0}$ (see Fig. 2(a)). Note that due to the small average thermal
center-of-mass velocity $\expectationvalue{v_{T}}$ of molecules, their motion
within the chamber is negligible on the time-scale of the photon propagation
$(\approx L/c)$. For instance, $\expectationvalue{v_{T}}L/c\approx 20~{}{\rm
nm}$ for $L=12~{}{\rm cm}$ and $T=20~{}{\rm K}$.

The absorbance $\mathcal{A}(\omega)$ of photons with the angular frequency
$\omega$ can be calculated using the Beer-Lambert’s law:
$\mathcal{A}(\omega)=1-\exp(-\chi_{\text{tot}}(\omega))$, where the total
attenuation $\chi_{\text{tot}}(\omega)=\sum_{i,f}\chi_{i\to f}(\omega)$ and
$\chi_{i\to f}(\omega)$ is the attenuation, associated with the single ro-
vibrational transition $i\to f$. Note that here, unlike Sec. V, indices $i$
and $f$ denote only the vibrational and rotational quantum numbers, i.e.,
$i(f)\equiv\\{\nu^{(\prime)},J^{(\prime)},N_{1}^{(\prime)},N_{2}^{(\prime)}\\}$.
The single-channel attenuation: $\chi_{i\to
f}(\omega)=\int_{-L/2}^{L/2}\tilde{\chi}_{i\to f}(\omega,Z)\mathrm{d}Z$. The
local attenuation at the position $Z$ is the ratio:
$\tilde{\chi}(\omega,Z)_{i\to f}=\mathbb{P}^{i\to
f}_{\text{abs}}(Z)/\mathbb{I}_{\text{inc}}(Z)$, where
$\mathbb{I}_{\text{inc}}(Z)$ is the energy transfer rate of the incident
electric field, given by

$\mathbb{I}_{\text{inc}}(Z)=\frac{c\varepsilon_{0}}{2}\int\limits_{0}^{\bm{R}_{0}}|E(\bm{R})|^{2}\,\mathrm{d}^{2}\bm{R}\,,$
(29)

and $\mathbb{P}^{i\to f}_{\text{abs}}(Z)$ is the energy flux, i.e., energy
transfer rate per unit length, absorbed by molecules, undergoing the
transition $i\to f$ and residing in a thin section of the chamber
$(\mathrm{d}Z)$. The absorbed energy flux is given by

$\displaystyle\mathbb{P}^{i\to f}_{\text{abs}}(Z)$
$\displaystyle=\rho_{0}\frac{\rho_{\text{B}}(E_{i},T)}{\mathcal{Z}}$
$\displaystyle\times\int\limits_{0}^{R_{0}}\mathrm{d}^{2}\bm{R}^{i}_{\text{CM}}\int\limits_{0}^{\infty}\mathrm{d}^{3}\bm{R}^{f}_{\text{CM}}\Big{(}\hbar\omega\Gamma_{i\to
f}(\omega)\Big{)}\,,$ (30)

where $\rho_{0}=\beta P$ is the equilibrium density of molecules,
$\rho_{\text{B}}(E_{i},T)=\exp(-\beta E_{i})$ is the Boltzmann probability to
occupy the initial state, ${\mathcal{Z}=\sum_{i}\rho_{\text{B}}(E_{i},T)}$ is
the canonical partition function. The integral over the initial center-of-mass
state is two-dimensional and covers the effective aperture $\bm{R}_{0}$ (see
Fig. 2(a)). The integral over the final state is three-dimensional and covers
the whole space. At a finite temperature $T$, spectral lines are broadened, so
instead of a $\delta$-function, used in the definition of the transition rate
(28), we employ the Doppler-broadened line profile

$\displaystyle\delta_{E}(\hbar\omega-\Delta
E_{fi}))\to\rho_{\text{DB}}(\omega,\Delta E_{fi},T)$
$\displaystyle=\sqrt{\frac{\beta mc^{2}}{2\pi(\Delta
E_{fi})^{2}}}\exp(-\frac{\beta mc^{2}(\hbar\omega-\Delta E_{fi})^{2}}{2(\Delta
E_{fi})^{2}})\,.$ (31)

In our numerical calculation, we consider CS molecules with the mass
$\mathbb{M}=44.076\,\text{amu}$ at the pressure $P=1\,\text{mbar}$ and
temperature $T=20\,\text{K}$. The rotational state of a linear molecule, like
CS, can be described by two quantum numbers $\\{l,m\\}$ instead of three
$\\{J,N_{1},N_{2}\\}$, as required in the general case (see Sec. V). For the
molecular state, defined in Sec. V, the difference in energy between the
initial and final states $\Delta E_{fi}$ is the difference in rotational and
vibrational energies. The combined rotational-vibrational energy for a linear
molecule can be calculated using the Dunham expansion [85] with coefficients
for CS provided in Ref. [86]. We consider a Laguerre-Gaussian beam (19) with a
width $\omega_{0}=300\,\text{µm}$ that is tuned to the second vibrational
overtone of CS, i.e., it excites the vibrational transition
$\nu=0\to\nu^{\prime}=2$ with the characteristic angular frequency
$\omega_{\text{vib}}=2488\,\text{cm}^{-1}$. Note that, compared to Eq. (19),
the Rayleigh range in a real experimental setup is given by:
$Z_{R}=\pi\omega_{0}^{2}/(\lambda\mathfrak{M}^{2})$, where $\mathfrak{M}^{2}$
is the beam quality factor [87]. For a Laguerre-Gaussian beam:
$\mathfrak{M}^{2}\sim 2|M|+1$. We calculate the corresponding transition
dipole $d^{0,2}_{z}=9.14\cdot 10^{-3}\,\text{D}$ and quadrupole
$Q^{0,2}_{zz}=-7.83\cdot 10^{-3}\,\text{D}\cdot\text{\AA}$ moments using
density functional theory, as described in Sec. II. In a real experiment,
transition multipole moments can be also measured directly. Our numerical
analysis has shown that the chamber length $L=12\,\text{cm}$ and the effective
aperture $R_{0}=140\,\text{µm}$ maximize the overall transition amplitudes.

We plot single-channel attenuation matrices $\chi_{i\to f}(\omega)$ for LG
beams with $M=0$ and $M=3$ in Fig. 2(b). To visually distinguish “forbidden”
and “allowed” transitions, we set the cutoff for the attenuation: $10^{-6}$.
Despite the fact, that matrices for non-helical and helical light have very
similar structure, our calculation reveals additional “forbidden” rotational
transitions, enabled by the OAM of light. To examine the effect in more
detail, we calculate the attenuation for different $\Delta m$-subchannels
within the $\Delta l=2$ rotational transition (see Fig. 2(c)). In this case,
we do not impose any numerical cutoff and reveal that the transitions with
$|\Delta m|=2$ are substantially enhanced by the OAM of light, confirming the
prediction of Sec. IV.

From Fig. 2(b), one can notice that the selection rules on the azimuthal
quantum number $l$ remain the same for $M=0$ and $M=3$. In other words,
$\Delta l$ transitions are not modified by the OAM of light. This can be
explained from the semiclassical perspective. The electric field in a vortex
beam is an eigenstate of the $\hat{L}_{z}$ operator, thus, it has the
cylindrical symmetry, whereas quantum number $l$ is the eigenvalue of the
$\hat{L}^{2}$ operator and describes the spherical symmetry of the molecular
state. Since $[\hat{L}_{z},\hat{L}^{2}]=0$, the beam is incapable of modifying
the spherical symmetry of a molecular state.

Figure 2: (a) Scheme of the generic proof-of-principle experiment. A probe
beam with the width profile $\omega_{z}$ propagates through the molecular
chamber of the characteristic length $L$. The transmitted light is collected
on the aperture $\tilde{R}_{0}$ outside the cavity. For small cavities
$(L<Z_{R})$, the domain of molecules that interact with the light, collected
on the aperture, can be approximated by the cylinder, defined by the effective
aperture $R_{0}$ (dotted line). (b) Single-channel attenuation $\chi_{i\to
f}(\omega)$ for rotational transitions
$\ket{l,m}\to\ket{l^{\prime},m^{\prime}}$ within the vibrational overtone
$(\nu=0\to\nu^{\prime}=2)$ of the CS molecule. Compared to the non-helical
light $(M=0)$, twisted light $(M=3)$ enables more transitions with $\chi_{i\to
f}(\omega)>10^{-6}$. (c) Single-channel attenuation $\chi_{i\to f}(\omega)$
for different $\Delta m$-subchannels within the $\Delta l=2$ rotational
transition, that corresponds to electric-quadrupole interaction. OAM of light
substantially enhances $|\Delta m|=2$ transitions that can be resolved in the
presence of a Stark shift in the system.

The inability of the OAM to affect the selection rules on the azimuthal
quantum number $l$ necessitates an additional requirement for the experiment.
In particular, in the aforementioned setup molecular ro-vibrational energies
are degenerate with respect to the magnetic quantum number $m$. This means
that in an absorption spectroscopy, one can only measure a weighted sum of
$\chi_{i\to f}(\omega)$ for all possible $\Delta m$. The transition with
$|\Delta m|=1$, which corresponds to the electric-dipole interaction with
light, is much stronger than lower-order quadrupole transitions. Thus, the
degeneracy renders it impossible to distinguish $M=0$ and $M=3$ beams. The
straightforward solution is to lift the degeneracy of molecular energies with
respect to $m$. For instance, this can be done by applying a strong static
electric field along the beam axis. Our numerical calculation indicates that
the Stark field $E_{\text{Stark}}\sim 8\,\text{kV}/\text{cm}$ is necessary to
be able to resolve the Doppler-broadened spectral lines in the absorption
spectroscopy of the CS.

##  VII Conclusion and outlook

Motivated by the potential benefits of using vortex beams for structural
studies and quantum control of molecules, we developed a general analytical
framework to describe the interaction of molecules with twisted light. Unlike
the previous studies, which involved various point-charge models [27, 28, 29,
30], we modeled the molecule using a continuous charge distribution. We
considered molecules interacting with long-wavelength, infrared or optical
light, which allowed us to describe light-matter interaction in terms of
multipole moments and avoid describing the details of the molecular structure.
Using the prototypical case of a single vibrational degree of freedom coupled
to the molecular rotations, we showed how various internal molecular
characteristics can be embedded into the analytical definition of the
multipole moments.

We used the multipolar QED Hamiltonian, also known as the PZW Hamiltonian [76,
77, 78], which was proven practical in molecular spectroscopy applications
[79]. We applied the spherical Taylor expansion [82, 81, 83] and derived the
general light-matter interaction Hamiltonian that describes the coupling
between spherical gradients of the electric field and molecular multipole
moments. This Hamiltonian is applicable to any molecular structure, provided
the validity of the multipole expansion, and to any electric field, with the
notable exception of tightly focused beams.

From the ground up, our analytical framework is built around the notion of the
angular momentum: we make use of spherical coordinates and rotational point
symmetry. Therefore, our theory unambiguously describes the angular momentum
exchange between the spin and OAM of light, molecular rotation and its center-
of-mass motion. In particular, we confirm a well-known result [27, 29] that,
within the electric-dipole approximation, the OAM of light cannot couple to
the molecular rotation irrespective of the beam profile. At the quadrupole
order, however, we demonstrate that vortex beams strongly enhance certain ro-
vibrational transitions, which are considered “forbidden” in the case of non-
helical light.

The enhancement strongly depends on the mean position of the molecular center
of mass with respect to the beam axis. Recent years have seen a rapid
development in the field of optical trapping and laser cooling of molecules
[88, 89, 90, 91, 92]. Therefore, it is realistic to suggest measuring the ro-
vibrational enhancement in these new setups. Another pathway is to study the
gas-phase absorption spectroscopy in a spatially and frequency-resolved
manner. We provide proof-of-principle calculations by suggesting an experiment
that may be capable of probing the enhancement. Based on our findings, we
discuss the experimental requirements for the observation of the enhanced ro-
vibrational transitions.

Our analytical framework is rather general, and should be applicable to a
broad range of problems that concern angular momentum exchange induced by the
light-matter interaction. Nevertheless, it relies on a few approximations.
Lifting some of these approximations will be an important future step. In
particular, one could extend the theory to beams in the regime of tight
focusing, which was shown to alter and enhance the OAM interaction with chiral
molecules and nanostructures [65]. In such electric fields, the spin and OAM
of light are strongly coupled [26], which makes the analysis of the OAM
transfer to the molecule more complicated.

Another important extension concerns adapting our model to account for more
complex molecular structures. In particular, polyatomic molecules of different
point groups can be addressed by generalizations of the charge distribution,
allowing for the coupling of multiple degrees of freedom to twisted light. A
very interesting perspective along this direction is the extension of our
approach to chiral molecules. In this context our model can address the
minimal requirements for the generation of helical dichroism [31, 32, 33] that
is a promissing experimental technique for the detection of chirality. Recent
studies have shown that an important aspect of the chirality-related physics
is the interplay of magnetic multipole transitions with electric multipole
ones. Such transitions can also be expressed in terms of a model Hamiltonian,
similar to the electric ones, by incorporating an effective magnetization
density for the molecule.

## Acknowledgements

We are grateful to Emilio Pisanty and Philipp Lunt for valuable discussions.
G.M.K. gratefully acknowledges funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie Grant
Agreement No. 101034413. M.L. acknowledges support by the European Research
Council (ERC) Starting Grant No. 801770 (ANGULON). O.H.H. acknowledges support
by the Austrian Science Fund (FWF) [P 36040, M 2561]. Furthermore, the
financial support by the Austrian Federal Ministry for Digital and Economic
Affairs, the National Foundation for Research, Technology and Development and
the Christian Doppler Research Association is gratefully acknowledged.

## References

  * Townes and Schawlow [2013] C. Townes and A. Schawlow,  _Microwave Spectroscopy_, Dover books on physics (Dover Publications, 2013).
  * Gordy  _et al._ [1984] W. Gordy, R. L. Cook, and A. Weissberger,  _Microwave Molecular Spectra_ , Vol. 18 (Wiley New York, 1984).
  * Hollas [2004] J. M. Hollas,  _Modern spectroscopy_ (John Wiley & Sons, 2004).
  * Barth [2007] A. Barth, Biochimica et Biophysica Acta (BBA) - Bioenergetics 1767, 1073 (2007).
  * Matthäus  _et al._ [2008] C. Matthäus, B. Bird, M. Miljković, T. Chernenko, M. Romeo, and M. Diem, Methods in cell biology 89, 275 (2008).
  * Movasaghi  _et al._ [2008] Z. Movasaghi, S. Rehman, and D. I. ur Rehman, Applied Spectroscopy Reviews 43, 134 (2008).
  * Parker [2012] F. Parker,  _Applications of infrared spectroscopy in biochemistry, biology, and medicine_ (Springer Science & Business Media, 2012).
  * Waynant  _et al._ [2001] R. W. Waynant, I. K. Ilev, and I. Gannot, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 359, 635 (2001).
  * Chance  _et al._ [1989] K. V. Chance, D. G. Johnson, and W. A. Traub, Journal of Geophysical Research: Atmospheres 94, 11059 (1989).
  * Rinsland  _et al._ [1989] C. P. Rinsland, D. W. Johnson, A. Goldman, and J. S. Levine, Nature 337, 535 (1989).
  * Markus  _et al._ [2016] C. R. Markus, J. N. Hodges, A. J. Perry, G. S. Kocheril, H. S. P. Müller, and B. J. McCall, The Astrophysical Journal 817, 138 (2016).
  * Pérez  _et al._ [2012] C. Pérez, M. T. Muckle, D. P. Zaleski, N. A. Seifert, B. Temelso, G. C. Shields, Z. Kisiel, and B. H. Pate, Science 336, 897 (2012).
  * Perakis  _et al._ [2016] F. Perakis, L. De Marco, A. Shalit, F. Tang, Z. R. Kann, T. D. Kühne, R. Torre, M. Bonn, and Y. Nagata, Chemical Reviews 116, 7590 (2016).
  * Hutson [1990] J. M. Hutson, Annual Review of Physical Chemistry 41, 123 (1990).
  * Cohen and Saykally [1992] R. C. Cohen and R. J. Saykally, The Journal of Physical Chemistry 96, 1024 (1992).
  * Legon and Millen [1982] A. C. Legon and D. J. Millen, Faraday Discuss. Chem. Soc. 73, 71 (1982).
  * Legon [1990] A. C. Legon, Chem. Soc. Rev. 19, 197 (1990).
  * Bliokh and Nori [2015] K. Y. Bliokh and F. Nori, Physics Reports 592, 1 (2015), transverse and longitudinal angular momenta of light.
  * Andrews and Babiker [2012] D. L. Andrews and M. Babiker,  _The angular momentum of light_ (Cambridge University Press, 2012).
  * Ashfold and Howe [1994] M. N. R. Ashfold and J. D. Howe, Annual Review of Physical Chemistry 45, 57 (1994).
  * Tamm and Weiss [1990] C. Tamm and C. O. Weiss, J. Opt. Soc. Am. B 7, 1034 (1990).
  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Phys. Rev. A 45, 8185 (1992).
  * Bazhenov  _et al._ [1992] V. Bazhenov, M. Soskin, and M. Vasnetsov, Journal of Modern Optics 39, 985 (1992).
  * Schemmel  _et al._ [2014] P. Schemmel, G. Pisano, and B. Maffei, Opt. Express 22, 14712 (2014).
  * Karimi  _et al._ [2014] E. Karimi, S. A. Schulz, I. De Leon, H. Qassim, J. Upham, and R. W. Boyd, Light: Science & Applications 3, e167 (2014).
  * Bliokh  _et al._ [2015] K. Y. Bliokh, F. J. Rodríguez-Fortuño, F. Nori, and A. V. Zayats, Nature Photonics 9, 796 (2015).
  * Babiker  _et al._ [2002] M. Babiker, C. R. Bennett, D. L. Andrews, and L. C. Dávila Romero, Phys. Rev. Lett. 89, 143601 (2002).
  * Alexandrescu  _et al._ [2006] A. Alexandrescu, D. Cojoc, and E. D. Fabrizio, Phys. Rev. Lett. 96, 243001 (2006).
  * Mondal  _et al._ [2014] P. K. Mondal, B. Deb, and S. Majumder, Phys. Rev. A 89, 063418 (2014).
  * Mukherjee  _et al._ [2017] K. Mukherjee, S. Majumder, P. K. Mondal, and B. Deb, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 015004 (2017).
  * Forbes and Andrews [2018] K. A. Forbes and D. L. Andrews, Opt. Lett. 43, 435 (2018).
  * Forbes and Andrews [2019] K. A. Forbes and D. L. Andrews, Phys. Rev. Res. 1, 033080 (2019).
  * Forbes and Jones [2021] K. A. Forbes and G. A. Jones, Phys. Rev. A 103, 053515 (2021).
  * Araoka  _et al._ [2005] F. Araoka, T. Verbiest, K. Clays, and A. Persoons, Phys. Rev. A 71, 055401 (2005).
  * Brullot  _et al._ [2016] W. Brullot, M. K. Vanbel, T. Swusten, and T. Verbiest, Science Advances 2, e1501349 (2016).
  * Müllner  _et al._ [2022] S. Müllner, F. Büscher, A. Möller, and P. Lemmens, Phys. Rev. Lett. 129, 207801 (2022).
  * Porfirev  _et al._ [2023] A. Porfirev, S. Khonina, and A. Kuchmizhak, Progress in Quantum Electronics 88, 100459 (2023).
  * Ren and Gu [2018] H. Ren and M. Gu, Laser Photonics Rev. 12, 1700255 (2018).
  * Arikawa  _et al._ [2020] T. Arikawa, T. Hiraoka, S. Morimoto, F. Blanchard, S. Tani, T. Tanaka, K. Sakai, H. Kitajima, K. Sasaki, and K. Tanaka, Sci. Adv. 6, eaay1977 (2020).
  * Zang and Lusk [2017] X. Zang and M. T. Lusk, Phys. Rev. A 96, 013819 (2017).
  * Zhang and Cui [2020] X. Zhang and T. J. Cui, ACS Photonics 7, 3291 (2020).
  * Albert  _et al._ [2020] V. V. Albert, J. P. Covey, and J. Preskill, Phys. Rev. X 10, 031050 (2020).
  * Langen  _et al._ [2023] T. Langen, G. Valtolina, D. Wang, and J. Ye, Quantum state manipulation and science of ultracold molecules (2023), arXiv:2305.13445 [cond-mat.quant-gas] .
  * Koch  _et al._ [2019] C. P. Koch, M. Lemeshko, and D. Sugny, Rev. Mod. Phys. 91, 035005 (2019).
  * Kobayashi  _et al._ [2019] J. Kobayashi, A. Ogino, and S. Inouye, Nature Communications 10, 3771 (2019).
  * Verma  _et al._ [2020] M. Verma, A. M. Jayich, and A. C. Vutha, Phys. Rev. Lett. 125, 153201 (2020).
  * Bader  _et al._ [2004] R. F. W. Bader, W. H. Henneker, and P. E. Cade, The Journal of Chemical Physics 46, 3341 (2004).
  * Bader and Bandrauk [2003] R. F. W. Bader and A. D. Bandrauk, The Journal of Chemical Physics 49, 1653 (2003).
  * Slater [1930] J. C. Slater, Phys. Rev. 36, 57 (1930).
  * Clementi and Raimondi [1963] E. Clementi and D. L. Raimondi, The Journal of Chemical Physics 38, 2686 (1963).
  * Waldron  _et al._ [2001] K. A. Waldron, E. M. Fehringer, A. E. Streeb, J. E. Trosky, and J. J. Pearson, Journal of Chemical Education 78, 635 (2001).
  * Seki  _et al._ [2020] T. Seki, K.-Y. Chiang, C.-C. Yu, X. Yu, M. Okuno, J. Hunger, Y. Nagata, and M. Bonn, The Journal of Physical Chemistry Letters 11, 8459 (2020).
  * Devlin  _et al._ [2001] J. P. Devlin, J. Sadlej, and V. Buch, The Journal of Physical Chemistry A 105, 974 (2001).
  * Corney [1978] A. Corney,  _Atomic and laser spectroscopy_ (Clarendon Press Oxford, 1978).
  * Winnewisser and Cook [1968] G. Winnewisser and R. L. Cook, Journal of Molecular Spectroscopy 28, 266 (1968).
  * Patel  _et al._ [1979] D. Patel, D. Margolese, and T. R. Dykea, The Journal of Chemical Physics 70, 2740 (1979).
  * Buckingham and Disch [1963] A. Buckingham and R.-L. Disch, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 273, 275 (1963).
  * Buckingham  _et al._ [1968] A. D. Buckingham, R. L. Disch, and D. A. Dunmur, Journal of the American Chemical Society 90, 3104 (1968).
  * Johnson III [2022] R. D. Johnson III, NIST computational chemistry comparison and benchmark database, NIST standard reference database number 101 (2022).
  * Green and Jones [2006] N. G. Green and T. B. Jones, Journal of Physics D: Applied Physics 40, 78 (2006).
  * Pulay  _et al._ [1979] P. Pulay, G. Fogarasi, F. Pang, and J. E. Boggs, Journal of the American Chemical Society 101, 2550 (1979).
  * Pulay  _et al._ [1981] P. Pulay, G. Fogarasi, and J. E. Boggs, The Journal of Chemical Physics 74, 3999 (1981).
  * Buckingham [1959] A. Buckingham, Quarterly Reviews, Chemical Society 13, 183 (1959).
  * Truhlar [1972] D. G. Truhlar, International Journal of Quantum Chemistry 6, 975 (1972).
  * Woźniak  _et al._ [2019] P. Woźniak, I. D. Leon, K. Höflich, G. Leuchs, and P. Banzer, Optica 6, 961 (2019).
  * Jackson [1999] J. D. Jackson, Classical electrodynamics (1999).
  * Scully and Zubairy [1999] M. O. Scully and M. S. Zubairy, Quantum optics (1999).
  * Varshalovich  _et al._ [1988] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii,  _Quantum theory of angular momentum_ (World Scientific, 1988).
  * Hall [2013a] B. C. Hall, Radial potentials and the hydrogen atom, in  _Quantum Theory for Mathematicians_ (Springer New York, New York, NY, 2013) pp. 393–418.
  * Hall [2013b] B. C. Hall, Lie groups, lie algebras, and representations, in  _Quantum Theory for Mathematicians_ (Springer New York, New York, NY, 2013) pp. 333–366.
  * Gelessus  _et al._ [1995] A. Gelessus, W. Thiel, and W. Weber, Journal of Chemical Education 72, 505 (1995).
  * Hermann  _et al._ [2014] K. Hermann, L. Pettersson, M. Casida, C. Daul, A. Goursot, A. Koester, E. Proynov, A. St-Amant, D. Salahub, _et al._ , STOBE software (2014).
  * Jones [2015] R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
  * Perdew and Wang [1992] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Power and Zienau [1959] E. A. Power and S. Zienau, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 251, 427 (1959).
  * Atkins and Woolley [1970] P. Atkins and R. Woolley, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 319, 549 (1970).
  * Woolley [1971] R. G. Woolley, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 321, 557 (1971).
  * Andrews  _et al._ [2018] D. L. Andrews, G. A. Jones, A. Salam, and R. G. Woolley, The Journal of Chemical Physics 148, 040901 (2018).
  * Zhao  _et al._ [2007] Y. Zhao, J. S. Edgar, G. D. M. Jeffries, D. McGloin, and D. T. Chiu, Phys. Rev. Lett. 99, 073901 (2007).
  * Weniger [2002] E. J. Weniger, International Journal of Quantum Chemistry 90, 92 (2002).
  * Weniger [2000] E. J. Weniger, International Journal of Quantum Chemistry 76, 280 (2000).
  * Weniger [2005] E. J. Weniger, Collection of Czechoslovak chemical communications 70, 1225 (2005).
  * Wu  _et al._ [2022] Y. Wu, S. Gargiulo, F. Carbone, C. H. Keitel, and A. Pálffy, Phys. Rev. Lett. 128, 162501 (2022).
  * Dunham [1932] J. L. Dunham, Phys. Rev. 41, 721 (1932).
  * Ram  _et al._ [1995] R. Ram, P. Bernath, and S. Davis, Journal of Molecular Spectroscopy 173, 146 (1995).
  * Siegman [1998] A. E. Siegman, in  _DPSS (Diode Pumped Solid State) Lasers: Applications and Issues_ (Optica Publishing Group, 1998) p. MQ1.
  * Anderegg  _et al._ [2018] L. Anderegg, B. L. Augenbraun, Y. Bao, S. Burchesky, L. W. Cheuk, W. Ketterle, and J. M. Doyle, Nature Physics 14, 890 (2018).
  * Mitra  _et al._ [2020] D. Mitra, N. B. Vilas, C. Hallas, L. Anderegg, B. L. Augenbraun, L. Baum, C. Miller, S. Raval, and J. M. Doyle, Science 369, 1366 (2020).
  * Zhang  _et al._ [2020] J. T. Zhang, Y. Yu, W. B. Cairncross, K. Wang, L. R. B. Picard, J. D. Hood, Y.-W. Lin, J. M. Hutson, and K.-K. Ni, Phys. Rev. Lett. 124, 253401 (2020).
  * Yu  _et al._ [2021] Y. Yu, K. Wang, J. D. Hood, L. R. B. Picard, J. T. Zhang, W. B. Cairncross, J. M. Hutson, R. Gonzalez-Ferez, T. Rosenband, and K.-K. Ni, Phys. Rev. X 11, 031061 (2021).
  * Cairncross  _et al._ [2021] W. B. Cairncross, J. T. Zhang, L. R. B. Picard, Y. Yu, K. Wang, and K.-K. Ni, Phys. Rev. Lett. 126, 123402 (2021).