# Unveiling the Quantum Toroidal Dipole in Nanosystems: Quantization,
Interaction Energy, and Measurement

Mircea Dolineanu, Alexandru-Lucian Nastasia , and Dragoş-Victor Anghel
Institutul National de Cercetare-Dezvoltare pentru Fizica si Inginerie
Nucleara Horia Hulubei; University of Bucharest, Doctoral School of Physics,
<EMAIL_ADDRESS>National de Cercetare-Dezvoltare
pentru Fizica si Inginerie Nucleara Horia Hulubei; University of Bucharest,
Faculty of Physics<EMAIL_ADDRESS>JINR, Dubna,
Moscow region, 141980, Russia; Institutul National de Cercetare-Dezvoltare
pentru Fizica si Inginerie Nucleara Horia Hulubei; corresponding author
<EMAIL_ADDRESS>
###### Abstract

We investigate the physical properties of a charged quantum particle confined
to a toroidal surface in the presence of a filiform current along the system’s
rotational axis. Our analysis reveals that the interaction between the
particle and the current induces a non-zero toroidal dipole in the particle’s
stationary states. We demonstrate that the differences between the toroidal
dipole projections for different energy levels can be quantized in units of
$\hbar R/(4m_{p})$ (where $R$ is the major radius of the torus and $m_{p}$ is
the particle mass), suggesting the existence of toroidal dipole quanta.
Furthermore, we find that both the toroidal dipole projection and the energy
eigenvalues exhibit periodic behavior with respect to the current intensity,
with a period that depends solely on the torus’s aspect ratio $R/r$, where r
is the minor radius. This periodicity opens up the possibility of using the
current intensity to manipulate and measure the toroidal dipole projection. We
also observe abrupt changes in the toroidal dipole projection and energy
eigenvalues around integer multiples of the current half-period. These changes
provide further evidence for the quantization of the toroidal dipole in such
systems. The interaction energy between the particle and the current follows
the classical electrodynamics form, suggesting a potential method for
measuring and manipulating the toroidal dipole projection along the current
axis. The quantization rules we have identified represent hallmarks of the
quantum toroidal dipole in nanosystems and could lead to the development of
novel devices based on this fundamental property.

Keywords: toroidal dipole operator; quantum observables; nanosystems;
metamaterials.

##  1 Introduction

In 1957, Zeldovich introduced the notion of ”anapole” as a multipole moment
different from the electric and magnetic ones, which could correspond to the
parity non-conserving term introduced by him for describing the $\beta$-decay
[1],

$\hat{{\cal H}}_{\beta}\propto\mathbf{S}\cdot\mathbf{J}^{\rm
ext}=\mathbf{S}\cdot(\gradient\crossproduct H^{\rm ext});$ (1)

in Eq. (1) ${\bf S}$ is the spin of the particle and ${\bf J}^{\rm ext}$ is
the external current corresponding to the external magnetic field $H^{ext}$.
The toroidal dipole was further studied by Dubovik and Cheskov in the
framework of classical electrodynamics, who showed that it is related to a
third class of multipoles, called toroidal multipoles [2, 3], which are
distinct from the usual electric and magnetic multipole families by their
symmetry properties under parity and time reversal transformations. A
comprehensive review of the multipole expansion in electrodynamics with a
clear focus on toroidal moments can be found in the book by Nanz [4]. As the
toroidal moments and the anapole have been studied in many contexts, from
particle physics to macroscopic physics [5, 6], these concepts were used to
describe systems at all scales, like nuclear physics [7], particle physics [8,
9, 10, 11], condensed matter physics [12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22], and metamaterials [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39].

The toroidal dipole, with the expression [13]

${\bf T}=\frac{1}{10}\int_{V}\Big{[}{\bf r}({\bf r}\cdot{\bf j})-2r^{2}{\bf
j}\Big{]}d^{3}{\bf r},$ (2)

where $\mathbf{j}$ is the current density and $\mathbf{r}$ is the position
vector, is the lowest order toroidal multipole and the easiest to observe.
Toroidal structures of currents, spins, charge distributions or magnetic field
distributions, which produce non-zero toroidal dipoles, appear in nuclei,
molecules (for example, in benzene, proteins, toroidal carbon cages),
biological cells, nanosystems, specially designed metamaterials, etc. (see [5,
14, 24] for examples). Also, toroidal polarizability and transitions could be
observed in hydrogen-like atoms [40, 14]. But, despite the obvious differences
between toroidal and electrical dipoles, their radiation patterns are
identical, which makes them indistinguishable at large distances. This
property made many authors question the physical significance of the toroidal
dipole and stimulated its study in the near field (see, for example, [36]).
The near field properties are important not only to justify the existence of
this third class of multipoles, but also because they are responsible for a
host of interesting effects, like non-radiating charge-current configurations,
nonreciprocal interactions, temperature and refractive index sensing, negative
refractive index, cloaking, superlensing, etc. (see [26, 39] and citations
therein).

In most of the previously mentioned papers, the toroidal moments are treated
in the context of classical electrodynamics. This may be justified for
condensed matter, solid state, and metamaterial applications, if the
microscopic toroidal structures (the elementary cells in metamaterials) are
large enough and quantum effects may be neglected. However, in order to
understand what would happen at smaller scales, a quantum approach to toroidal
moments is needed, especially, considering that there are possible
applications of toroidal moments, for example in dark matter particle models
and WISP detection, quantum computing and nanoscale imaging [5, 6]. In this
regard, the toroidal dipole operators

$\hat{T}_{i}\equiv\frac{1}{10m_{p}}\sum_{j=1}^{3}\left(x_{i}x_{j}-2r^{2}\delta_{ij}\right)\hat{p}_{j},$
(3)

with $x_{i}$ and $\hat{p}_{i}$ being the Cartesian position and conjugate
momentum operators ($i=1,2,3$), were introduced in [40] in order to study the
toroidal polarizability of the hydrogen atom. Later, the mathematical
properties of this operator were analyzed in [41], where a first attempt was
made at solving the eigenproblem. When the operator was defined on the whole
space, it was found to be hypermaximal, and the possibility of it being a
quantum observable is questionable. Recently, the operator was particularized
for the case of a particle on a torus-shaped nanostructure (see Fig. 1) in the
absence of electromagnetic field and a general formalism was established for
such systems [42, 43, 44]. The self-adjointness of the operator was proved,
therefore the toroidal dipole may be a quantum observable [43]. The eigenvalue
problem for this system was also solved, but the results revealed that this
operator’s eigenstates are outside the Hilbert space of (square integrable)
wavefunctions, raising another question regarding to the observability of the
toroidal dipole, since the wavefunction cannot collapse on such an eigenvector
after an eventual measurement [44].

Once established that the toroidal dipole operator is self-adjoint, we have to
understand how this quantum property can be observed or even measured. In the
previous studies, no method to measure the toroidal dipole was used or
proposed. Instead, some region of the system was studied by some (commercially
available, in general) simulation software and the magnetic or current field
lines were visually analyzed to identify specific toroidal configurations.
This would represent the proof of the presence of toroidal dipoles, which may
then be correlated with other specific phenomena, like the cancellation of the
electromagnetic radiation at large distances.

Here, we propose a method to identify and measure the toroidal dipole in
quantum systems. We start from the interaction Hamiltonian

$H_{I}\equiv-\mu_{0}q\left[{\bf J}+\dot{{\bf D}}\right]\cdot{\bf T},$ (4)

which represents the hallmark of the presence of a toroidal dipole in the
system, where ${\bf J}$ and $\dot{{\bf D}}$ are charge and displacement
currents, respectively, $\mu_{0}=1.256637\times 10^{-6}\ {\rm
kg\,m/(s^{2}A^{2})}$ is the vacuum permittivity, and $q$ is the charge of the
particle. The interaction (4) of the system with an external current leads to
the modification of the energy levels, proportional to the product between the
current intensity and the expectation value of the toroidal dipole.
Furthermore, if the toroidal dipole expectation values on the Hamiltonian
eigenstates are quantized (as we shall see it happens here in most cases), the
modification of the energy levels is also quantized, which constitutes an
additional verification of the model and of the quantum nature of the toroidal
dipole. This method is analog to the detection and measurement of the magnetic
moment of a particle, by its interaction with an external magnetic field. In
that case also, the energy levels change by amounts proportional to the
intensity of the field multiplied to the magnetic moment of the particle.
Furthermore, since the projection of the magnetic moment on the field
direction is quantized, the splitting of the energy levels is also quantized
for a fixed value of the intensity of the field.

Figure 1: The system of curvilinear coordinates $(\theta,\phi,\xi)$ on a torus
made of a thin layer of material. The inner and outer surfaces of the material
are two tori with the same major radius $R$, but minor radii $r$ and
$r+\xi_{max}$, respectively, where $r+\xi_{max}<R$ and $\xi_{max}\ll r$ – that
is, $\xi=0$ on the interior surface and $\xi=\xi_{max}$ on the exterior
surface; we denote $a\equiv R/r>1$.

The system we study is depicted in Fig. 1. A quantum particle of charge $q$ is
confined on a torus, in the presence of a current $I$ which flows along the
torus’ rotation axis $z$. We solve the eigenvalues and eigenvectors problem
for the Hamiltonian operator and we calculate the expectation values of the
projection $\hat{T}_{3}$ (3) of the toroidal dipole operator on the stationary
states. We show how the interaction between the current and toroidal dipole
induced in the system modify the energy levels, providing in this way a method
to measure the toroidal dipole and to show its quantum nature.

The paper is organized as follows. In Section 2 we introduce the formalism and
write down the Hamiltonian, together with all the necessary formulas that we
shall use in the following sections. In Section 3, we calculate
analytically–in the low current limit–the energy levels, the energy
eigenstates, and the expectation values of $\hat{T}_{3}$, proving the validity
of Eq. (4), together with the quantization of the projection of the toroidal
dipole expectation values, which further leads to the quantization of the
change of the energy levels. In Section 4 we solve the problem numerically for
arbitrary values of $I$, to find the range of validity of Eq. (4) and of the
quantization of the toroidal dipole expectation values. Section 5 is reserved
for conclusions and the appendices for some analytical calculations necessary
in the main body of the manuscript, but which are not important for the
understanding of the physical phenomena presented here.

##  2 The particle on a torus surface in the presence of an external
electromagnetic field

The particle is confined in the region between two tori, as in Fig. 1. Both
tori have cylindrical symmetry around the $z$ axis, common origin $O$ (the
origin of the Cartesian system of coordinates), the same major radius $R$, and
the minor radii $r$ and $r+\xi_{max}$, where $\xi_{max}\ll r$. We also
consider that $r+\xi_{max}<R$, so there is a hole in the middle of the
external torus. As in [43], we introduce the coordinates $(\theta,\phi,\xi)$
($\xi$ replaces here the notation for the coordinate $q$ of [43]), where
$\phi$ is the azimuthal angle, $\xi$ is the coordinate along the minor radius
between the inner torus (at $\xi=0$) and the external torus (at
$\xi=\xi_{max}$), whereas $\theta$ is the angle formed by the minor radius
with the $(x,y)$ plane, as shown in Fig. 1. Then, the Cartesian coordinates
$(x,y,z)$ may be written as $x=\left[R+(r+\xi)\cos{\theta}\right]\cos{\phi}$,
$y=\left[R+(r+\xi)\cos{\theta}\right]\sin{\phi}$, $z=(r+\xi)\sin{\theta}$, and
the corresponding scale factors are $h_{\theta}=r+\xi$,
$h_{\phi}=R+(r+\xi)\cos{\theta}$, and $h_{\xi}=1$. The determinant of the
Jacobian of this transformation is
$\operatorname{det}(J)=-(r+\xi)\left[R+(r+\xi)\cos{\theta}\right]$ [43].

###  2.1 The Hamiltonian

If a spinless particle, of charge $q$ and mass $m_{p}$, is in an external
electromagnetic field of vector and scalar potentials ${\bf A}$ and $\phi$,
respectively, then its Hamiltonian is

$\hat{H}=\frac{1}{2m_{p}}\left(\hat{\mathbf{p}}-q\mathbf{A}\right)^{2}+q\Phi.$
(5)

Writing $\hat{H}$ in the coordinates $(\theta,\phi,\xi)$ and taking into
consideration that the system is tightly confined in the $\xi$ direction, then
only the coordinates $\theta$ and $\phi$ remain relevant and the Hamiltonian
is reduced to the curved 2D surface of the torus [43]:

$\displaystyle\hat{{\cal H}}^{2D}_{elmg}$ $\displaystyle=$
$\displaystyle\hat{{\cal
H}^{2D}_{\text{free}}}+i\frac{q\hbar}{m_{p}}\frac{1}{r(R+r\cos\theta)}\left\\{(R+r\cos\theta)A_{\theta}\partialderivative{\theta}+rA_{\phi}\partialderivative{\phi}+\frac{1}{2}\partialderivative{[A_{\theta}(R+r\cos\theta)]}{\theta}+\frac{r}{2}\partialderivative{(A_{\phi})}{\phi}\right\\}$
(6a) $\displaystyle+\frac{q^{2}}{2m_{p}}(A_{\theta}^{2}+A_{\phi}^{2})+q\Phi,$
where $\displaystyle\hat{{\cal H}}^{2D}_{\text{free}}$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}a^{2}}{2m_{p}R^{2}}\left(\partialderivative[2]{\theta}-\frac{\sin{\theta}}{\left(a+\cos{\theta}\right)}\partialderivative{\theta}+\frac{1}{\left(a+\cos{\theta}\right)^{2}}\partialderivative[2]{\phi}+\frac{a^{2}}{4\left(a+\cos{\theta}\right)^{2}}\right)\equiv\frac{\hbar^{2}}{2m_{p}R^{2}}\tilde{{\cal
H}}^{2D}_{\text{free}}$ (6b)

is the free Hamiltonian, which was analyzed in [43]. In units
$\hbar^{2}/(2m_{p}R^{2})$, $\hat{{\cal H}}^{2D}_{elmg}$ becomes

$\displaystyle\hat{{\cal H}}^{2D}_{\text{elmg}}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{2m_{p}R^{2}}\left(\tilde{{\cal
H}}^{2D}_{\text{free}}+i\frac{a}{a+\cos\theta}\left\\{2(a+\cos\theta)\tilde{A}_{\theta}\partialderivative{\theta}+2\tilde{A}_{\phi}\partialderivative{\phi}+\partialderivative{[\tilde{A}_{\theta}(a+\cos\theta)]}{\theta}+\partialderivative{(\tilde{A}_{\phi})}{\phi}\right\\}\right.$
(7)
$\displaystyle\left.+(\tilde{A}_{\theta}^{2}+\tilde{A}_{\phi}^{2})+q\tilde{\Phi}\right)\equiv\frac{\hbar^{2}}{2m_{p}R^{2}}\tilde{{\cal
H}}^{2D}_{\text{elmg}},$

where $\tilde{{\bf A}}\equiv(qR/\hbar){\bf A}$ and
$\tilde{\Phi}\equiv(2m_{p}R^{2}/\hbar^{2})\Phi$.

Let’s assume that through the hole of the torus, along the $z$ axis, passes a
thin straight wire of length $L$, which carries a constant current $I$ (if $L$
is big enough, the circuit may be closed at large distances with additional
wires and a current source which would not influence the torus). The vector
potential corresponding to this setup is:

$\mathbf{A}=\left(0,0,\frac{\mu_{0}I}{2\pi}\ln{\left(\frac{2L}{\rho}\right)}\right)\,,$
(8)

where $\rho=R+r\cos\theta$ is the distance from the $z$-axis and $\mu_{0}$ is
the vacuum permeability introduced in Eq. (4). Considering that
$A_{\phi}=\mathbf{A}\cdot\hat{\mbox{\boldmath$\phi$}}$,
$A_{\theta}=\mathbf{A}\cdot\hat{\mbox{\boldmath$\theta$}}$, and
$A_{\xi}=\mathbf{A}\cdot\hat{\mbox{\boldmath$\xi$}}$, we can rewrite the
vector potential in the $(\theta,\phi,\xi)$ coordinate system:

$A_{\phi}=0,\quad
A_{\theta}=\frac{\mu_{0}I}{2\pi}\ln{\left(\frac{2L}{\rho}\right)}\cos\theta,\quad
A_{\xi}=\frac{\mu_{0}I}{2\pi}\ln{\left(\frac{2L}{\rho}\right)}\sin\theta.$ (9)

Plugging (9) into (7), we obtain

$\displaystyle\tilde{{\cal H}}^{(2D)}_{\text{elmg}}$ $\displaystyle=$
$\displaystyle\tilde{{\cal H}}^{(2D)}_{\text{free}}+\tilde{{\cal
H}}^{(2D)}_{I},\quad\tilde{{\cal
H}}^{(2D)}_{I}\equiv\frac{iI}{I_{0}}a\log{\left(\frac{2a}{a+\cos\theta}\frac{L}{R}\right)}\cos\theta\partialderivative{\theta}+\frac{iI}{2I_{0}}\frac{a\sin\theta}{a+\cos\theta}\Bigg{[}\cos\theta$
(10)
$\displaystyle-(a+2\cos\theta)\log{\left(\frac{2a}{a+\cos\theta}\frac{L}{R}\right)}\Bigg{]}+\frac{I^{2}\cos^{2}\theta}{4I_{0}^{2}}\log^{2}{\left(\frac{2a}{a+\cos\theta}\frac{L}{R}\right)},$

where $I_{0}\equiv\pi\hbar/(\mu_{0}qR)$ will be our unit of current.

###  2.2 The momentum operator and the vector basis

Following Refs. [43, 44], we introduce the momentum operators in the
coordinates $(\theta,\phi,\xi)$. The momentum along the $\phi$ direction is
related to the $z$ component of the angular momentum operator ($\hat{L}_{3}$),
whereas the momentum along the $\xi$ direction is not self-adjoint (because of
the boundary conditions), so it is not an observable. The momentum operator
along the $\theta$ direction is

$\hat{p}^{(\theta)}\equiv\frac{-i\hbar}{R}a\left\\{\frac{\partial}{\partial\theta}-\frac{\sin\theta}{2(a+\cos\theta)}\right\\}.$
(11)

The eigenvectors of $\hat{p}^{(\theta)}$, combined with the eigenvectors of
$\hat{L}_{3}$,

${\cal F}^{(\Lambda)}_{n,m}\equiv\frac{a}{2\pi
R}\frac{e^{i(n\theta+m\phi)}}{\sqrt{a+\cos\theta}}\quad({\rm where}\ n,m\ {\rm
are\ integers}),$ (12)

form a basis in the Hilbert space $H$ of wave functions of variables
$(\theta,\phi)$ defined on the surface of the torus.

In Ref. [43] we used the basis (12) to numerically diagonalize the operators
$\hat{{\cal H}}^{2D}_{\text{free}}$ and $\hat{T}^{(\theta)}_{3}$. In this
basis we have the matrix elements

$\langle{\cal F}^{(\Lambda)}_{n_{1},m_{1}}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\Lambda)}_{n_{2},m_{2}}\rangle=a^{2}\left[\left(n_{1}^{2}-\frac{1}{4}\right)\delta_{n_{1}n_{2}}+\frac{\left({|\delta
n|}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{|\delta
n|}}}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left(m_{1}^{2}-\frac{1}{4}\right)\right]\delta_{m_{1}m_{2}}$
(13a) where $\delta n\equiv n_{1}-n_{2}$. From Eq. (13a) we observe that the
diagonal elements $\langle{\cal F}^{(\Lambda)}_{n,m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\Lambda)}_{n,m}\rangle=a^{2}\left[\left(n^{2}-\frac{1}{4}\right)+\frac{a}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left(m^{2}-\frac{1}{4}\right)\right]=\langle{\cal
F}^{(\Lambda)}_{-n,m}|\tilde{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\Lambda)}_{-n,m}\rangle$ (13b) increase linearly with $n^{2}$, whereas
the non-diagonal elements ($n_{1}\neq n_{2}$) $\langle{\cal
F}^{(\Lambda)}_{n_{1},m}|\hat{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\Lambda)}_{n_{2},m}\rangle=\frac{a^{2}\left({|\delta
n|}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{|\delta
n|}}}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left(m^{2}-\frac{1}{4}\right)$ (13c)

decrease roughly in geometrical progression with $|\delta n|$ because of the
factor $\left(\sqrt{a^{2}-1}-a\right)^{{|\delta n|}}$, since
$|\sqrt{a^{2}-1}-a|<1$.

####  2.2.1 Symmetric and antisymmetric basis vectors

Exploiting the symmetry of the diagonal matrix elements to the change of $n$
into $-n$ (13b) and the fast decrease of non-diagonal elements with $|\delta
n|$, we introduce the symmetric ($+$) and antisymmetric ($-$) basis vectors

${\cal F}^{(+)}_{0,m}\equiv{\cal F}^{(\Lambda)}_{0,m}\quad{\rm and}\quad{\cal
F}^{(\pm)}_{n,m}\equiv\frac{1}{\sqrt{2}}\left[{\cal
F}^{(\Lambda)}_{n,m}\pm{\cal F}^{(\Lambda)}_{-n,m}\right],$ (14a) for any
$n\in\mathbb{N}^{*}$ (positive integer), which satisfy the conditions
$\langle{\cal F}^{(+)}_{n_{1},m_{1}}|{\cal
F}^{(+)}_{n_{2},m_{2}}\rangle=\langle{\cal F}^{(-)}_{n_{1},m_{1}}|{\cal
F}^{(-)}_{n_{2},m_{2}}\rangle=\delta_{n_{1},n_{2}}\delta_{m_{1},m_{2}},\quad\langle{\cal
F}^{(+)}_{n_{1},m_{1}}|{\cal F}^{(-)}_{n_{2},m_{2}}\rangle=0.$ (14b)

In the basis (14) we have

$\displaystyle\langle{\cal F}^{(\pm)}_{n,m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(\pm)}_{n,m}\rangle$ $\displaystyle=$
$\displaystyle
a^{2}\left\\{\left[\left(n^{2}-\frac{1}{4}\right)+\frac{a}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left(m^{2}-\frac{1}{4}\right)\right]\pm\frac{\left({2n}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{2n}}}{\left(a^{2}-1\right)^{\frac{3}{2}}}\right.$
(15a) $\displaystyle\left.\times\left(m^{2}-\frac{1}{4}\right)\right\\},$
$\displaystyle\langle{\cal F}^{(\pm)}_{n_{1},m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(\pm)}_{n_{2},m}\rangle$ $\displaystyle=$
$\displaystyle
a^{2}\left\\{\frac{\left({|n_{1}-n_{2}|}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{|n_{1}-n_{2}|}}}{\left(a^{2}-1\right)^{\frac{3}{2}}}\right.$
$\displaystyle\left.\pm\frac{\left({|n_{1}+n_{2}|}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{|n_{1}+n_{2}|}}}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left(m^{2}-\frac{1}{4}\right)\right\\},\quad{\rm
for}\quad n_{1}\neq n_{2},\quad{\rm and}$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle\langle{\cal F}^{(+)}_{n_{1},m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(-)}_{n_{2},m}\rangle=\langle{\cal
F}^{(-)}_{n_{1},m}|\tilde{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(+)}_{n_{2},m}\rangle.$ (15c)

We observe again the fast increase of the diagonal elements with $n$, whereas
the non-diagonal elements decrease fast with $|n_{1}-n_{2}|$ and
$|n_{1}+n_{2}|$. For this reason, the basis (14a) represents a good
approximation for the eigenfunctions of $\hat{{\cal H}}^{(2D)}_{\text{free}}$
[43] (the higher the energy level $n$, the better the approximation is).

Equations (15) imply that the Hilbert space $H$ may be split in two: $H\equiv
H^{(+)}\oplus H^{(-)}$, where $H^{(+)}$ and $H^{(-)}$ are spanned by the bases
$\left\\{{\cal F}^{(+)}_{n,m}\right\\}$ and $\left\\{{\cal
F}^{(-)}_{n,m}\right\\}$, respectively. Then, the eigenfunctions of
$\hat{{\cal H}}^{(2D)}_{\text{free}}$ belong either to $H^{(+)}$ or to
$H^{(-)}$, but they are never a sum of vectors from both spaces.

####  2.2.2 The matrix elements of the toroidal dipole in the absence of the
electromagnetic field

In the absence of the electromagnetic fields, the toroidal dipole operator has
the expression (3), which, in the variables $(\theta,\xi)$ has the hermitic
components [43]

$\displaystyle\hat{T}^{(\theta)}_{3}$ $\displaystyle=$
$\displaystyle\frac{\hbar
R}{10m_{p}}\left\\{i\frac{3a\cos(2\theta)+4\left(a^{2}+1\right)\cos\theta+5a}{2a}\frac{\partial}{\partial\theta}-\tilde{T}_{\theta}(\theta,\xi)\right\\},$
(16a) $\displaystyle\hat{T}^{(\xi)}_{3}$ $\displaystyle=$
$\displaystyle\frac{\hbar
R}{10m_{p}}\left\\{i\frac{2a^{2}+1+3a\cos\theta}{a^{2}}\sin\theta
R\frac{\partial}{\partial\xi}-\tilde{T}_{\xi}(\theta,\xi)\right\\},$ (16b)
such that $\tilde{T}_{\theta}(\theta,\xi)=-\tilde{T}_{\xi}(\theta,\xi)$, so
$\displaystyle\hat{T}_{3}=\hat{T}^{(\theta)}_{3}+\hat{T}^{(\xi)}_{3}.$ (16c)
In the limit $\xi_{max}\ll r$, $\tilde{T}_{\theta}(\theta,\xi)$ and
$\tilde{T}_{\xi}(\theta,\xi)$ are simplified to $\tilde{T}_{\theta}\equiv
i\frac{\left[9a\cos^{2}\theta+2\left(5a^{2}+2\right)\cos\theta+a\left(2a^{2}+3\right)\right]\sin\theta}{2a\left(a+\cos\theta\right)}=-\tilde{T}_{\xi}.$
(16d)

The matrix elements of $\hat{T}^{(\theta)}_{3}$ are [43]

$\langle{\cal F}^{(\Lambda)}_{n_{1},m_{1}}|\hat{T}_{3}^{(\theta)}|{\cal
F}^{(\Lambda)}_{n_{2},m_{2}}\rangle=-\frac{\hbar
R}{10m_{p}}\frac{n_{1}+n_{2}}{2a}\left\\{\frac{5a}{2}\delta_{n_{1}n_{2}}+\left(a^{2}+1\right)\left(\delta_{n_{1},n_{2}+1}+\delta_{n_{1},n_{2}-1}\right)+\frac{3a}{4}\left(\delta_{n_{1},n_{2}+2}+\delta_{n_{1},n_{2}-2}\right)\right\\}\delta_{m_{1}m_{2}},$
(17)

which, in the basis (14a), give

$\langle{\cal F}^{(+)}_{n_{1},m_{1}}|\hat{T}_{3}^{(\theta)}|{\cal
F}^{(+)}_{n_{2},m_{2}}\rangle=\langle{\cal
F}^{(-)}_{n_{1},m_{1}}|\hat{T}_{3}^{(\theta)}|{\cal
F}^{(-)}_{n_{2},m_{2}}\rangle=0$ (18)

for any $n_{1}>0$, $n_{2}>0$, $m_{1}$, and $m_{2}$.

To calculate the expectation values of the whole $\hat{T}_{3}$ operator (16c),
we first observe that since the eigenstates of the Hamiltonian are stationary
states, the probability current in the $\xi$ direction should be zero (no net
current flowing in the $\xi$ direction). Taking this into account, plus the
conditions (16c) and (16d), the expression for the expectation value of
$\hat{T}_{3}$ simplifies to

$\left\langle\Psi^{(0)}\left|\hat{T}_{3}\right|\Psi^{(0)}\right\rangle=\frac{\hbar
R}{10m_{p}}\left\langle\Psi^{(0)}\left|\tilde{T}_{3}^{({\bf
j},0)}\right|\Psi^{(0)}\right\rangle,\quad{\rm where}\quad\tilde{T}_{3}^{({\bf
j},0)}\equiv
i\frac{3a\cos(2\theta)+4\left(a^{2}+1\right)\cos\theta+5a}{2a}\frac{\partial}{\partial\theta}$
(19)

where $\Psi^{(0)}$ is an arbitrary eigenfunction of $\hat{{\cal
H}}^{(2D)}_{\text{free}}$. In Appendix B we calculate the matrix elements of
$\hat{T}_{3}^{({\bf j},0)}$ in the basis (12), using the definitions given in
Appendix A. Using these matrix elements we obtain a condition similar to (18),
namely

$\langle{\cal F}^{(+)}_{n_{1},m_{1}}|\tilde{T}_{3}^{({\bf j},0)}|{\cal
F}^{(+)}_{n_{2},m_{2}}\rangle=\langle{\cal
F}^{(-)}_{n_{1},m_{1}}|\tilde{T}_{3}^{({\bf j},0)}|{\cal
F}^{(-)}_{n_{2},m_{2}}\rangle=0$ (20)

for any $n_{1}>0$, $n_{2}>0$, $m_{1}$, and $m_{2}$.

Equations (18), (19), and (20) imply that the expectation values of both,
$\hat{T}_{3}^{(\theta)}$ and $\hat{T}_{3}$ on the eigenstates of $\hat{{\cal
H}}^{(2D)}_{\text{free}}$ are zero.

####  2.2.3 The matrix elements of the interaction Hamiltonian

In the basis (12), the matrix elements of the interaction Hamiltonian (10) are

$\displaystyle\left\langle{\cal F}^{(\Lambda)}_{n_{2}-n,m}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal
F}^{(\Lambda)}_{n_{2},m}\right\rangle=\frac{I}{I_{0}}a\left\\{-\log\left(\frac{2aL}{R}\right)\left[\frac{n_{2}}{2}(\delta_{n+1}+\delta_{n-1})+\frac{I_{n+2}-I_{n-2}}{8}\right]-n_{2}\left[\frac{\log\frac{2}{a+\sqrt{a^{2}-1}}}{2}\right.\right.$
$\displaystyle\left.\times\left(\delta_{n+1}+\delta_{n-1}\right)+\left.\frac{I_{n+2}-I_{n}}{4(n+1)}\right|_{n+1\neq
0}+\left.\frac{I_{n}-I_{n-2}}{4(n-1)}\right|_{n-1\neq
0}\right]+\frac{I_{n+2}-I_{n-2}}{8}-\frac{(n+1)I^{(ln^{2})}_{n+1}+(n-1)I^{(ln^{2})}_{n-1}}{8}$
$\displaystyle\left.-\frac{a}{4}\left[\log\left(\frac{2aL}{R}\right)\left(I_{n+1}-I_{n-1}\right)+nI^{(ln^{2})}_{n}\right]\right\\}+\frac{I^{2}}{16I_{0}^{2}}\left\\{\log^{2}\left(\frac{2aL}{R}\right)(\delta_{n+2}+2\delta_{n}+\delta_{n-2})\right.$
$\displaystyle+\log\left(\frac{2aL}{R}\right)\left[+2\log\frac{2}{a+\sqrt{a^{2}-1}}\left(\delta_{n+2}+2\delta_{n}+\delta_{n-2}\right)+\left.\frac{I_{n+3}-I_{n+1}}{n+2}\right|_{n+2\neq
0}+2\left.\frac{I_{n+1}-I_{n-1}}{n}\right|_{n\neq 0}\right.$
$\displaystyle\left.\left.+\left.\frac{I_{n-1}-I_{n-3}}{n-2}\right|_{n-2\neq
0}\right]+\left(I^{(ln^{2})}_{n+2}+2I^{(ln^{2})}_{n}+I^{(ln^{2})}_{n-2}\right)\right\\}$
(21)

where (see Appendix A)

$\displaystyle I_{n}(a)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}\frac{e^{in\theta}}{a+\cos\theta}d\theta=\frac{\left(-a+\sqrt{a^{2}-1}\right)^{|n|}}{\sqrt{a^{2}-1}},$
(22a) $\displaystyle I_{n}^{(ln^{2})}(a)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}e^{in\theta}\ln^{2}(a+\cos\theta)\,d\theta=\frac{1}{2\pi}\int_{0}^{2\pi}\cos(n\theta)\ln^{2}(a+\cos\theta)\,d\theta.$
(22b)

Using the expression (15) and (21) we write the Hamiltonian matrix in the
basis $\left\\{{\cal F}^{(\pm)}_{n,m}\right\\}$ (14), which we then
diagonalize numerically to obtain the eigenvalues and eigenvectors, as
explained in Section 4. Using this basis is especially important in the case
$I=0$, since it preserves the (anti)symmetry of the eigenvectors.

###  2.3 The toroidal dipole in the presence of electromagnetic fields

In the presence of electromagnetic fields, the definition (3) of the toroidal
dipole operator should be changed by taking into consideration the
contribution of the fields to the current density operator:

${\bf
j}=\frac{1}{2m_{p}}\left(\Psi^{*}\hat{\bm{p}}\Psi-\Psi\hat{\bm{p}}\Psi^{*}-2q\bm{A}|\Psi|^{2}\right).$
(23)

With this expression for ${\bf j}$, the components of the toroidal dipole
operator change into

$\displaystyle\hat{T}_{i}=\frac{1}{10m_{p}}\sum_{j=1}^{3}(x_{i}x_{j}-2r^{2}\delta_{ij})(\hat{p}_{j}-qA_{j}).$
(24)

The $z$ component written in toroidal coordinates $(\theta,\xi)$ is

$\displaystyle\hat{T}_{3}=-\frac{i\hbar}{10m_{p}}[z\rho\hat{\bm{\rho}}-(2\rho^{2}+z^{2})\hat{\bm{z}}]\cdot\left[\hat{\bm{\theta}}\left(\frac{1}{r}\frac{\partial}{\partial\theta}-\frac{iq}{\hbar}A_{\theta}\right)+\hat{\bm{\xi}}\left(\frac{\partial}{\partial\xi}-\frac{iq}{\hbar}A_{\xi}\right)\right],$
(25)

where $A_{\theta}$ and $A_{\phi}$ are the components of ${\bf A}$ along the
$\hat{\mbox{\boldmath$\theta$}}$ and $\hat{\mbox{\boldmath$\xi$}}$ directions.
We split $\hat{T}_{3}$ into the hermitic components $\hat{T}_{3}^{(\theta)}$
and $\hat{T}_{3}^{(\xi)}$ according to the recipe given in [43] (as we did in
Section 2.2.2 for the case without electromagnetic fields), to obtain

$\displaystyle\hat{T}_{3}^{(\theta)}$ $\displaystyle=$
$\displaystyle\frac{\hbar
R}{10m_{p}}\left\\{\frac{\left[z\rho\hat{\bm{\rho}}-(2\rho^{2}+z^{2})\hat{\bm{z}}\right]\cdot\hat{\bm{\theta}}}{R^{2}}\left(-ia\frac{\partial}{\partial\theta}-\frac{qR}{\hbar}A_{\theta}\right)-\tilde{T}_{\theta}\right\\},$
(26) $\displaystyle\hat{T}_{3}^{(\xi)}$ $\displaystyle=$
$\displaystyle\frac{\hbar
R}{10m_{p}}\left\\{\frac{\left[z\rho\hat{\bm{\rho}}-(2\rho^{2}+z^{2})\hat{\bm{z}}\right]\cdot\hat{\bm{\xi}}}{R^{2}}\left(-iR\frac{\partial}{\partial\xi}-\frac{qR}{\hbar}A_{\xi}\right)-\tilde{T}_{\xi}\right\\},$
(27)

where $\tilde{T}_{\theta}$ and $\tilde{T}_{\xi}$ are given in Eqs. (16).
Taking into consideration the confinement in the $\xi$ direction, the
component of $\hat{T}_{3}$ which corresponds to the probability current along
$\xi$–that is, $\partial/\partial\xi-iqA_{q}/\hbar$–should be zero on any
eigenstate of the Hamiltonian. Taking this into consideration together with
(16d), by calculating
$\hat{T}_{3}\equiv\hat{T}_{3}^{(\theta)}+\hat{T}_{3}^{(\xi)}$ we remain only
with the term proportional to the probability current along the $\theta$
direction:

$\left\langle\Psi^{(I)}\left|\hat{T}_{3}\right|\Psi^{(I)}\right\rangle\equiv\left\langle\Psi^{(I)}\left|\hat{T}_{3}^{({\bf
j})}\right|\Psi^{(I)}\right\rangle,$ (28)

where $|\Psi^{(I)}\rangle$ is an eigenstate of the Hamiltonian (10) and

$\displaystyle\hat{T}_{3}^{({\bf j})}$ $\displaystyle=$
$\displaystyle-\frac{i\hbar}{10m_{p}}\left[z\rho\hat{\bm{\rho}}-(2\rho^{2}+z^{2})\hat{\bm{z}}\right]\cdot\hat{\bm{\theta}}\left(\frac{1}{r}\frac{\partial}{\partial\theta}-\frac{iq}{\hbar}A_{\theta}\right)\equiv\frac{\hbar
R}{10m_{p}}\left(\tilde{T}_{3}^{({\bf j},0)}+\tilde{T}_{3}^{({\bf
j},I)}\right),$ $\displaystyle{\rm with}$ $\displaystyle\tilde{T}_{3}^{({\bf
j},I)}\equiv\frac{3a\cos(2\theta)+4\left(a^{2}+1\right)\cos\theta+5a}{2a^{2}}\tilde{A}_{\theta},$
(29)

whereas $A_{\theta}$ and $\tilde{T}_{3}^{({\bf j},0)}$ are defined in Eqs. (9)
and (19), respectively, and
$\tilde{A}_{\theta}\equiv\left(qR/\hbar\right)A_{\theta}$. The matrix elements
of $\tilde{T}_{3}^{({\bf j},0)}$ and $\tilde{T}_{3}^{({\bf j},I)}$ are
calculated in Appendix B and the final result is

$\displaystyle\langle{\cal F}^{(\Lambda)}_{n_{2}-n,m_{1}}|\tilde{T}_{3}^{({\bf
j})}|{\cal
F}^{(\Lambda)}_{n_{2},m_{2}}\rangle=-\delta_{m_{1},m_{2}}\left(\frac{3}{4}\left(n_{2}-\frac{1}{2}\right)\delta_{n+2}+\frac{3}{4}\left(n_{2}+\frac{1}{2}\right)\delta_{n-2}+\frac{1}{4a}\left[4n_{2}\left(a^{2}+1\right)+a^{2}-2\right]\delta_{n+1}\right.$
$\displaystyle+\frac{1}{4a}\left[4n_{2}\left(a^{2}+1\right)-a^{2}+2\right]\delta_{n-1}+\frac{5n_{2}}{2}\delta_{n}-{\rm
sgn}(n)\frac{(a^{2}-1)}{2}\left(\sqrt{a^{2}-1}-a\right)^{|n|}-\frac{I}{4a^{2}I_{0}}\left\\{\log(\frac{2aL}{R})\right.$
$\displaystyle\times\left[\frac{3a}{4}\left(\delta_{n+3}+\delta_{n-3}\right)+\left(a^{2}+1\right)\left(\delta_{n+2}+\delta_{n-2}\right)+\frac{13a}{4}\left(\delta_{n+1}+\delta_{n-1}\right)+2\left(a^{2}+1\right)\delta_{n}\right]$
(30)
$\displaystyle\left.\left.-\left[\frac{3a}{4}\left(I^{(\ln)}_{n+3}(a)+I^{(\ln)}_{n-3}(a)\right)+\left(a^{2}+1\right)\left(I^{(\ln)}_{n+2}(a)+I^{(\ln)}_{n-2}(a)\right)+\frac{13a}{4}\left(I^{(\ln)}_{n+1}(a)+I^{(\ln)}_{n-1}(a)\right)+2\left(a^{2}+1\right)I^{(\ln)}_{n}(a)\right]\right\\}\right).$

After we diagonalize the Hamiltonian, we use the expression (30) to calculate
the expectation values of $\hat{T}_{3}$ on the Hamiltonian eigenfunctions, by
using the property (28). The results are reported in the next section,
together with the diagonalization of the Hamiltonian. In the next sections we
shall see that the external electromagnetic fields mix the vector sub-spaces
$H^{(+)}$ and $H^{(-)}$ in the Hamiltonian eigenfunctions, so the expectation
values of $\hat{T}_{3}^{(\theta)}$ on these eigenfunctions is not zero for
$I\neq 0$.

##  3 Analytical approximations and the measurement of the toroidal dipole

At large $n$ or $a$, the diagonal elements of the free particle Hamiltonian in
the basis $\left\\{{\cal F}^{(+)}_{n,m},{\cal F}^{(-)}_{n,m}\right\\}$ are
much bigger than the non-diagonal ones (see Eqs. 15). In such cases, the
Hamiltonian eigenfunctions are well approximated by the basis functions [45]
and the energy levels are almost double degenerate (15):

$\displaystyle\left|\langle{\cal F}^{(+)}_{n,m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(+)}_{n,m}\rangle-\langle{\cal
F}^{(-)}_{n,m}|\tilde{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(-)}_{n,m}\rangle\right|=\frac{2\left({2n}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{2n}}}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left|m^{2}-\frac{1}{4}\right|$
$\displaystyle\ll\langle{\cal F}^{(\pm)}_{n,m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(\pm)}_{n,m}\rangle-\langle{\cal
F}^{(\pm)}_{n-1,m}|\tilde{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\pm)}_{n-1,m}\rangle=(2n-1)\left[1\mp\frac{2a^{2}\left(a-\sqrt{a^{2}-1}\right)^{2n-1}}{\sqrt{a^{2}-1}}\right],$
(31)

since $0<a-\sqrt{a^{2}-1}<1$. In such a case, at small values of $I$, the
contribution of $\tilde{{\cal H}}^{2D}_{I}$ to the Hamiltonian may be treated
as a perturbation in the 2D space of basis $\left({\cal F}^{(+)}_{n,m},{\cal
F}^{(-)}_{n,m}\right)$. Since a non-zero current $I$ breaks the symmetry of
the system, it becomes more convenient to work in the basis $\left({\cal
F}^{(\Lambda)}_{n,m},{\cal F}^{(\Lambda)}_{-n,m}\right)$, which spans the same
2D space. In the new basis, the free Hamiltonian matrix is given by Eqs. (13)
and has the symmetric form (69) discussed in Appendix C,

$A\left(\begin{array}[]{cc}1&\epsilon\\\ \epsilon&1\end{array}\right){\bf
v}=A\alpha{\bf v},$ (32a) where
$A=a^{2}\left[\left(n^{2}-\frac{1}{4}\right)+\frac{a}{\left(a^{2}-1\right)^{\frac{3}{2}}}\left(m^{2}-\frac{1}{4}\right)\right]\quad\text{and}\quad\epsilon=\frac{\left({2n}\sqrt{a^{2}-1}+a\right)\left(\sqrt{a^{2}-1}-a\right)^{{2n}}}{\left(n^{2}-\frac{1}{4}\right)\left(a^{2}-1\right)^{\frac{3}{2}}+a\left(m^{2}-\frac{1}{4}\right)},$
(32b)

with the obvious eigenvalues $\langle{\cal F}^{(+)}_{n,m}|\tilde{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(+)}_{n,m}\rangle$ and $\langle{\cal
F}^{(-)}_{n,m}|\tilde{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(-)}_{n,m}\rangle$ and eigenvectors ${\cal F}^{(+)}_{n,m}$ and ${\cal
F}^{(-)}_{n,m}$, respectively.

The matrix elements of the interaction Hamiltonian are

$\displaystyle\left\langle{\cal F}^{(\Lambda)}_{nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal
F}^{(\Lambda)}_{nm}\right\rangle=\frac{an(a-\sqrt{a^{2}-1})I}{I_{0}}+\frac{I^{2}}{8I_{0}^{2}}\log\left(\frac{2aL}{R}\right)\left\\{\log\left(\frac{2aL}{R}\right)+2\log\frac{2}{a+\sqrt{a^{2}-1}}\right.$
$\displaystyle\left.+\left(a-\sqrt{a^{2}-1}\right)^{2}+I^{(ln^{2})}_{2}+I^{(ln^{2})}_{0}\right\\},$
(33a) $\displaystyle\left\langle{\cal F}^{(\Lambda)}_{-nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal
F}^{(\Lambda)}_{-nm}\right\rangle=-\frac{an(a-\sqrt{a^{2}-1})I}{I_{0}}+\frac{I^{2}}{8I_{0}^{2}}\log\left(\frac{2aL}{R}\right)\left\\{\log\left(\frac{2aL}{R}\right)+2\log\frac{2}{a+\sqrt{a^{2}-1}}\right.$
$\displaystyle\left.+\left(a-\sqrt{a^{2}-1}\right)^{2}+I^{(ln^{2})}_{2}+I^{(ln^{2})}_{0}\right\\},$
(33b) $\displaystyle\left\langle{\cal F}^{(\Lambda)}_{-nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal
F}^{(\Lambda)}_{nm}\right\rangle=\frac{I}{I_{0}}a\left(\sqrt{a^{2}-1}-a\right)^{2n}\left\\{\log\left(\frac{2aL}{R}\right)\frac{a}{2}+\frac{n}{2}\left(\frac{a-\sqrt{a^{2}-1}}{2n+1}+\frac{a+\sqrt{a^{2}-1}}{2n-1}\right)\right.$
$\displaystyle\left.-\frac{a}{2}-\frac{a}{2}\log\left(\frac{2aL}{R}\right)\right\\}+\frac{I^{2}}{16I_{0}^{2}}\log\left(\frac{2aL}{R}\right)\left(\sqrt{a^{2}-1}-a\right)^{2n}\left[\frac{\left(a-\sqrt{a^{2}-1}\right)^{2}}{n+1}+\frac{2}{n}+\frac{\left(a+\sqrt{a^{2}-1}\right)^{2}}{n-1}\right]$
$\displaystyle-\frac{I}{I_{0}}a\frac{(2n+1)I^{(ln^{2})}_{2n+1}+(2n-1)I^{(ln^{2})}_{2n-1}+4naI^{(ln^{2})}_{2n}}{8}+\frac{I^{2}}{16I_{0}^{2}}\left(I^{(ln^{2})}_{2n+2}+2I^{(ln^{2})}_{2n}+I^{(ln^{2})}_{2n-2}\right),$
(33c)

where we used the expression (22a) for $I_{n}$ and in the calculation of the
non-diagonal elements we assumed $n\geq 2$, so that $2n-3\geq 0$ (since we are
interested, in general, in large $n$ and $a$). The matrix elements (33) are
added to the matrix elements of the free Hamiltonian and we obtain a matrix of
the form (71),

$A_{I}\left(\begin{array}[]{cc}1+\delta&\epsilon_{I}\\\
\epsilon_{I}&1-\delta\end{array}\right){\bf v}_{I}=A_{I}\alpha_{I}{\bf
v}_{I}.$ (34a) where $\displaystyle A_{I}$ $\displaystyle=$ $\displaystyle
A+\frac{I^{2}}{8I_{0}^{2}}\log\left(\frac{2aL}{R}\right)\left\\{\log\left(\frac{2aL}{R}\right)-2\log(a)+I^{(ln^{2})}_{2}+I^{(ln^{2})}_{0}\right\\},$
(34b) $\displaystyle\delta_{I}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{I}}\frac{nI}{2I_{0}},$ (34c)
$\displaystyle\epsilon_{I}$ $\displaystyle=$ $\displaystyle\frac{\epsilon
A+\left\langle{\cal F}^{(\Lambda)}_{-nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal F}^{(\Lambda)}_{nm}\right\rangle}{A_{I}}$ (34d)

Equations (34) have the solutions

$\alpha_{I}^{(\pm)}=1\pm\sqrt{\epsilon_{I}^{2}+\delta_{I}^{2}},\ {\bf
v}_{I}^{(\pm)}=N_{I}^{(\pm)}\left(\begin{array}[]{c}1\\\
-\frac{\delta_{I}}{\epsilon_{I}}\pm\sqrt{1+\frac{\delta_{I}^{2}}{\epsilon_{I}^{2}}}\end{array}\right),\
N_{I}^{(\pm)}=\frac{1}{\sqrt{2\sqrt{1+\frac{\delta_{I}^{2}}{\epsilon_{I}^{2}}}\left(\sqrt{1+\frac{\delta_{I}^{2}}{\epsilon_{I}^{2}}}\mp\frac{\delta_{I}}{\epsilon_{I}}\right)}}.$
(35)

In Appendix D we calculated analytic approximations for $I_{n}$ and
$I_{n}^{(\ln^{2})}$ in the case $a\gg 1$, where we have also

$a-\sqrt{a^{2}-1}=a\left(1-\sqrt{1-1/a^{2}}\right)\approx 1/(2a)\ll 1.$ (36)

Using the expressions (36) and (77), for $n\geq 2$, Eqs. (33) simplify to

$\displaystyle\left\langle{\cal F}^{(\Lambda)}_{nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal F}^{(\Lambda)}_{nm}\right\rangle$ $\displaystyle=$
$\displaystyle\frac{nI}{2I_{0}}\left\\{1+\frac{I}{8nI_{0}}\log\left(\frac{2aL}{R}\right)\log\left(\frac{2L}{R}\right)\right\\},$
(37a) $\displaystyle\left\langle{\cal F}^{(\Lambda)}_{-nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal F}^{(\Lambda)}_{-nm}\right\rangle$ $\displaystyle=$
$\displaystyle-\frac{nI}{2I_{0}}\left\\{1-\frac{I}{4nI_{0}}\log\left(\frac{2aL}{R}\right)\log\left(\frac{2L}{R}\right)\right\\},$
(37b) $\displaystyle\left\langle{\cal F}^{(\Lambda)}_{-nm}\left|\tilde{{\cal
H}}^{2D}_{I}\right|{\cal F}^{(\Lambda)}_{nm}\right\rangle$ $\displaystyle=$
$\displaystyle\frac{I}{I_{0}}\frac{1}{8a^{2n-2}}\left(\frac{1}{2^{2n-3}(2n-1)}-\left(\sum_{k=1}^{2n-2}\frac{1}{k}\right)-\frac{2(3n-1)}{n(2n-1)}\right.$
$\displaystyle\left.+\frac{I}{4(n-1)I_{0}}\left\\{\log\left(\frac{2L}{R}\right)\frac{1}{2^{2n-3}}+\ln
a\left[2+\left(\sum_{k=1}^{2n-3}\frac{1}{k}\right)\right]\right\\}\right).$ In
the same approximation, the matrix elements of the free Hamiltonian (13)
become $\langle{\cal F}^{(\Lambda)}_{n,m}|\hat{{\cal
H}}^{(2D)}_{\text{free}}|{\cal F}^{(\Lambda)}_{n,m}\rangle=\langle{\cal
F}^{(\Lambda)}_{-n,m}|\hat{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\Lambda)}_{-n,m}\rangle\approx
a^{2}\left(n^{2}-\frac{1}{4}\right)\quad\text{and}\quad\langle{\cal
F}^{(\Lambda)}_{-n,m}|\hat{{\cal H}}^{(2D)}_{\text{free}}|{\cal
F}^{(\Lambda)}_{n,m}\rangle=\frac{\left({2n}+1\right)}{(2a)^{2n}}\left(m^{2}-\frac{1}{4}\right).$
(37d)

Plugging Eqs. (37) into (34) we obtain

$\displaystyle A_{I}$ $\displaystyle=$ $\displaystyle
a^{2}\left(n^{2}-\frac{1}{4}\right)+\frac{I^{2}}{8I_{0}^{2}}\log\left(\frac{2aL}{R}\right)\log\left(\frac{2L}{R}\right),$
(38a) $\displaystyle\delta_{I}$ $\displaystyle=$
$\displaystyle\frac{1}{A_{I}}\frac{nI}{2I_{0}}=\frac{nI}{2\left(n^{2}-\frac{1}{4}\right)I_{0}a^{2}}\frac{1}{1+\frac{I^{2}}{8\left(n^{2}-\frac{1}{4}\right)I_{0}^{2}a^{2}}\log\left(\frac{2aL}{R}\right)\log\left(\frac{2L}{R}\right)},$
(38b) $\displaystyle\epsilon_{I}$ $\displaystyle=$
$\displaystyle\frac{\left({2n}+1\right)\left(m^{2}-\frac{1}{4}\right)(2a)^{-2n}+\left\langle{\cal
F}^{(\Lambda)}_{-nm}\left|\tilde{{\cal H}}^{2D}_{I}\right|{\cal
F}^{(\Lambda)}_{nm}\right\rangle}{A_{I}},$ (38c)
$\displaystyle\frac{\delta_{I}}{\epsilon_{I}}$ $\displaystyle=$
$\displaystyle\frac{k_{1}I/I_{0}}{1+k_{2}(I/I_{0})+k_{3}(I/I_{0})^{2}},\quad\text{where}$
(38d) $\displaystyle k_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{n(2a)^{2n}}{2\left({2n}+1\right)\left(m^{2}-\frac{1}{4}\right)},$
(38e) $\displaystyle k_{2}$ $\displaystyle\equiv$
$\displaystyle-\frac{2^{2n-3}a^{2}}{\left(m^{2}-\frac{1}{4}\right)}\frac{\left(\sum_{k=1}^{2n-2}\frac{1}{k}\right)+\frac{2(3n-1)}{n(2n-1)}-\frac{1}{2^{2n-3}(2n-1)}}{{2n}+1},$
(38f) $\displaystyle k_{3}$ $\displaystyle\equiv$
$\displaystyle\frac{2^{2n-5}a^{2}}{\left(m^{2}-\frac{1}{4}\right)}\frac{\left[2+\left(\sum_{k=1}^{2n-3}\frac{1}{k}\right)\right]\ln
a+\frac{1}{2^{2n-3}}\log\left(\frac{2L}{R}\right)}{({2n}+1)(n-1)},\quad\text{and}\quad
n\geq 2.$ (38g)

Figure 2: The solution $x_{\text{d}}\equiv I_{\text{d}}/I_{0}$ of the equation
$D(x)=0$, in the case $m=0$.

Figure 3: The ratio $|\delta_{I}/\epsilon_{I}|$ in the case $m=0$ and
$n=2,3,5$, as indicated in each plot (for $m=0$, $\delta_{I}/\epsilon_{I}<0$
for small values of $I/I_{0}$).

In the limit $I\to 0$, $A_{I}\to
A,\delta_{I}\to\delta,\epsilon_{I}\to\epsilon$, and the system (34) transforms
into (32), discussed above. Let’s analyze now what happens as $I$ increases
from zero (because of the symmetry of the system, we shall consider only the
case $I\geq 0$). The denominator in the expression (38d) is a second order
polynomial that we shall denote by $D(x)\equiv 1+k_{2}x+k_{3}x^{2}$, with
$x\equiv I/I_{0}$. If $m=0$, then $k_{2}>0$ whereas $k_{1},k_{3}<0$.
Calculating the derivative $D^{\prime}(x)=k_{2}+2k_{3}x$, we observe that
both, $D(0)$ and $D^{\prime}(0)$ are positive at $x=0$. Since
$\lim_{x\to\infty}D(x)=-\infty$, this implies that for $m=0$, $D(x)$ has one
and only one solution at $x>0$, namely

$\displaystyle x_{\text{d}}$ $\displaystyle\equiv$
$\displaystyle\frac{I_{\text{d}}}{I_{0}}=\frac{-k_{2}-\sqrt{k_{2}^{2}-4k_{3}}}{2k_{3}}.$
(39)

(see Fig. 2). Calculating

$\frac{d(\delta_{I}/\epsilon_{I})}{d(I/I_{0})}=\frac{k_{1}\left[1-k_{3}(I/I_{0})^{2}\right]}{\left[k_{3}(I/I_{0})^{2}+k_{2}I/I_{0}+1\right]^{2}},$
(40)

we observe that for $m=0$, $d(\delta_{I}/\epsilon_{I})/d(I/I_{0})<0$ for any
$I$. Since $\delta_{I}/\epsilon_{I}|_{I=0}<0$, Eq. (40) implies that
$\delta_{I}/\epsilon_{I}|_{I<I_{\text{d}}}<0$ and
$\delta_{I}/\epsilon_{I}|_{I>I_{\text{d}}}>0$, whereas
$\delta_{I}/\epsilon_{I}|_{I=i_{\text{d}}}$ is divergent. In Fig. 3 we plot
the $\delta_{I}/\epsilon_{I}$ for $n=2,3,5$. We see that in the vicinity of
$I=0$, $\delta_{I}/\epsilon_{I}$ increases very fast with $I$, especially when
$a$ or $n$ are big.

For $m\neq 0$, $k_{1},k_{3}>0$ and $k_{2}<0$. The function
$\delta_{I}/\epsilon_{I}$ is positive at $I=0$, whereas its derivative (40),
also positive at $I=0$, is equal to zero at

$\frac{I_{\text{max}}}{I_{0}}=\frac{1}{\sqrt{k_{3}}},$ (41)

where $\delta_{I}/\epsilon_{I}$ reaches its local maximum

$\left(\frac{\delta_{I}}{\epsilon_{I}}\right)_{\text{max}}=\frac{k_{1}}{k_{2}+2\sqrt{k_{3}}}.$
(42)

The ratio $I_{\text{max}}/I_{0}$ is plotted in Fig. 4 and in Fig. 5 we plot
$\delta_{I}/\epsilon_{I}$ for $m=1$ and $n=2,3,5$.

Figure 4: The ratio $I/I_{0}$ at which $\delta_{I}/\epsilon_{I}$ reaches its
maximum for $m\neq 0$, $n\geq 2$, and $a\geq 2$.

Figure 5: The ratio $\delta_{I}/\epsilon_{I}$ for $m=1$ and $n=2,3,5$, as
indicated in each plot. Figure 6: The absolute value of the slope
$d(\delta_{I}/\epsilon_{I})/d(I/I_{0})$ at $I=0$ and $m=0$.

From Eqs. (38) and Figs. (3) and (5) we observe that if $n\geq 2$ and $a\gg 1$
(in practice, $a>2$ is enough in most cases) the ratio
$\delta_{I}/\epsilon_{I}$ increases very fast with $I$ in the vicinity of
$I=0$. From Eq. (40) we find the slope

$\lim_{I\to
0}\left[\frac{d(\delta_{I}/\epsilon_{I})}{d(I/I_{0})}\right]=k_{1},$ (43)

plotted in Fig. 6. Therefore, we may find a range for $I$, such that
$I/I_{0}\ll 1\ll k_{1}I/I_{0}$ and $\delta_{I}/\epsilon_{I}\gg 1$. In such a
case, as we showed in the end of Appendix C, the eigenvectors of the
Hamiltonian converge to $\left({\cal F}^{(\Lambda)}_{n,m},{\cal
F}^{(\Lambda)}_{-n,m}\right)$, with the corresponding eigenvalues

$\tilde{E}^{(\pm)}_{n}=\left\langle{\cal F}^{(\Lambda)}_{\pm
n,m}\left|\tilde{{\cal H}}^{2D}_{\text{elmg}}\right|{\cal F}^{(\Lambda)}_{\pm
n,m}\right\rangle=A_{I}\alpha_{I}^{(\pm)}=A_{I}\left(1\pm\delta_{I}\right)=A\pm\frac{nI}{2I_{0}}.$
(44a) We observe, therefore, a splitting of the quasi-degenerate energy levels
by $\Delta\tilde{E}_{n}^{(I)}=n\frac{I}{I_{0}}.$ (44b)

Further, the expectation values of $\tilde{T}_{3}$ on the Hamiltonian
eigenstates (see Eq. 30) have very simple values,

$\langle\tilde{T}_{3}\rangle_{n}=\left\langle{\cal
F}^{(\Lambda)}_{n,m}\left|\tilde{T}_{3}\right|{\cal
F}^{(\Lambda)}_{n,m}\right\rangle=-\frac{5n}{2}=-\left\langle{\cal
F}^{(\Lambda)}_{-n,m}\left|\tilde{T}_{3}\right|{\cal
F}^{(\Lambda)}_{-n,m}\right\rangle\equiv-nt_{3}^{(0)}.$ (45)

Notice that $t_{3}^{(0)}=5/2$ is different from the quanta $t_{3}$ of the
eigenvalues of the operator $\hat{T}_{3}^{(\theta)}$ obtain in Ref. [44].
Combining Eqs. (44) and (45) we obtain

$\Delta\tilde{E}_{n}^{(I)}=-\frac{2}{5}\frac{I}{I_{0}}\langle\tilde{T}_{3}\rangle_{n}.$
(46)

Including the dimensional scales for energy (7) and toroidal dipole (16), we
obtain the linear relation between the energy deviation and toroidal dipole

$\frac{\Delta
E_{n}^{(I)}}{2}\equiv\frac{\hbar^{2}}{2m_{p}R^{2}}\frac{\Delta\tilde{E}_{n}^{(I)}}{2}=-\mu_{0}q\frac{I}{\pi
R^{2}}\langle\tilde{T}_{3}\rangle_{n},$ (47)

which is identical to (4), if we make the natural identification for the
current density, $j\equiv I/(\pi R^{2})$.

Equation (47) gives us a direct method to observe the toroidal dipole and to
measure its expectation values. In this way one may observe also the
quantization (45), as a proof of the quantum nature of $T_{3}$.

##  4 Numerical study of stationary states and toroidal moment values

Figure 7: The energies (in dimensionless units) of the first eleven
Hamiltonian eigenstates a functions of $I$ (in $I_{0}$ units) for $a=1.2,2,3$
and $m=0,1,2$, as indicated in each plot. The legend, shown in (a), is the
same in every plot: $0=$ ground state, $1=$ first excited state, and so on.

Figure 8: The expectation vales of $\tilde{T}_{3}$ on the Hamiltonian
eigenstates vs. $I$ (in $I_{0}$ units) for $a=1.2,2,3$ and $m=0,1,2$, as
indicated in each plot. The legend, shown only once, is the same in every
plot: $0=$ ground state, $1=$ first excited state, $2=$ second excited state,
$\ldots$.

Figure 9: The expectation vales of $\hat{T}_{3}$ on the Hamiltonian
eigenstates vs. $I$ (in $I_{0}$ units) for $a=1.2,2,3$ and $m=0,1,2$, as
indicated in each plot. The legend, shown only once, is the same in every
plot: $0=$ ground state, $1=$ first excited state, $2=$ second excited state,
$\ldots$.

Figure 10: The relation between the energy level splitting and the product
$(I/I_{0})\langle\tilde{T}_{3}\rangle_{n}$, at $a=3$ and $m=0,1,2$. We see
that
$\Delta\tilde{E}_{n}^{(I)}\approx(2/5)(I/I_{0})\langle\tilde{T}_{3}\rangle_{n}$
and the approximation gets better and better as $\eta$ increases.

Figure 11: The coefficients of the Hamiltonian eigenstates in the $\\{{\cal
F}^{(\Lambda)}_{n,m}\\}$ basis.

For numerical computations we use $L=2000$ nm and $R=250$ nm, whereas the
filiform current $I$ (flowing along the $z$ axis) and the torus’ minor radius
$r$, characterized by the parameter $a\equiv R/r$, are varying (see Fig. 1).
The Hamiltonian of the system is given by Eqs (10) and (6b). Similarly to the
procedure detailed in [43], we numerically computed the matrix elements
$\langle{\cal F}^{(\pm)}_{n_{1},m_{1}}|\tilde{{\cal H}}^{(2D)}|{\cal
F}^{(\pm)}_{n_{2},m_{2}}\rangle$ using the expressions (15) and (21), and we
performed numerical diagonalization. The energy levels as functions of $I$ are
shown in Fig. 7 for different values of $a$ and $m$. In connection to these,
we calculated the expectation values of $\tilde{T}_{3}$ on the Hamiltonian
eigenstates (using Eq. 30), which we plot in Figs. 8 and 9. We index the
energy levels by $\eta=0,1,2,\ldots$, starting from the ground state and going
up in energy. Therefore, the Hamiltonian eigenvalues and eigenstates will be
denoted by $E_{\eta}(I)$ and $|\Psi_{\eta}(I)\rangle$, such that
$E_{\eta}(I)\leq E_{\eta^{\prime}}(I)$ if $\eta<\eta^{\prime}$. For the
expectation values of $\tilde{T}_{3}$ we shall use the simpler notation
$\langle\tilde{T}_{3}\rangle_{\eta}(I)\equiv\langle\Psi_{\eta}(I)|\tilde{T}_{3}|\Psi_{\eta}(I)\rangle$.
To keep the notations simple, the dependence of $|\Psi_{\eta}(I)\rangle$,
$E_{\eta}(I)$, and $\langle\tilde{T}_{3}\rangle_{\eta}(I)$ on $a$ and $m$ will
be considered implicit. Even more, sometimes we may consider implicit also
dependence of these quantities on $I$. Comparing this notation with the one
used in Eqs. (44), we observe that, if the approximations of Section 3 apply,
then

$E_{\eta=2n-1}\approx E_{n}^{(-)}\quad\text{and}\quad E_{\eta=2n}\approx
E_{n}^{(+)}.$ (48a) In Eqs. (44) and (48b) we have $n>0$, since $E_{n}^{(+)}$
and $E_{n}^{(-)}$ refer to $n$ and $-n$, respectively. To avoid confusion,
from now on we shall use the correspondence
$\eta\equiv\begin{cases}2n,&\text{if}\quad n>0,\\\ -2n-1,&\text{if}\quad
n<0.\end{cases}$ (48b)

At $I=0$, the energy levels (see Fig. 7) are the same as the ones calculated
in [43] and the expectation values of $\tilde{T}_{3}$ are zero, as proved in
Section 2.2.2 (Fig. 9). In the presence of the external current, the energy
levels change and the energy eigenstates acquire non-zero toroidal dipole. The
cause of the energy levels change is the interaction between the toroidal
dipoles and the current–analogously to the modification of the energy levels
of a particle with nonzero magnetic moment, in the presence of an external
magnetic field. Both, $E_{\eta}(I)$ and
$\langle\tilde{T}_{3}\rangle_{\eta}(I)$ show a periodic patterns which depend
only on $a$ (and are independent of $\eta$ and $m$).

To understand this better, let’s look at Fig. 9. If $a$ or $\eta$ are
sufficiently large, $\langle\tilde{T}_{3}\rangle_{\eta}(I)$ varies rapidly
with $I/I_{0}$, approaching the value $-5n/2$ (48b), in accordance with the
analytical calculations of Section 3. It is not necessary to give here a
precise definition for when “$a$ or $\eta$ are sufficiently large” and we
shall give it only a qualitative significance. For example, we may roughly say
that this condition is satisfied for the levels 7 to 10 of panel (a), levels 5
to 10 of panels (d) and (e), levels 7 to 10 of panel (f), levels 3 to 10 of
panel (g), and levels 5 to 10 of panels (h) and (i).

On a larger interval of currents (see Fig. 8), if $a$ or $\eta$ are
sufficiently large, we observe that $\langle\tilde{T}_{3}\rangle_{\eta}(I)$
decreases linearly with $I$, having abrupt changes of sign at regular
intervals $lI_{\text{s}}$, where $l$ is an integer and $I_{\text{s}}$ depends
on $a$, but not on $\eta$ or $m$. In each interval of linear variation
included between the limits $lI_{\text{s}}$ and $(l+1)I_{\text{s}}$, the slope
$d\langle\tilde{T}_{3}\rangle_{\eta}(I)/d(I/I_{0})=-(5/4)(I_{0}/I_{\text{s}})$,
is the same for any $m$, $l$, $a$, and $\eta$. Therefore, the decrease of the
linear approximation of $\langle\tilde{T}_{3}\rangle_{\eta}(I)$ on any
interval $(lI_{\text{s}},(l+1)I_{\text{s}})$ takes the universal value of
$5/4$, independent of current, angular momentum, shape of the torus or energy
level, as long as $a$ or $\eta$ are sufficiently large. Moreover, differences
between expectation values take also universal forms

$\langle\tilde{T}_{3}\rangle_{\eta}(I)-\langle\tilde{T}_{3}\rangle_{\eta-1}(I)=\frac{(-1)^{\eta-1}5}{2}\,\eta\quad\text{and}\quad\langle\tilde{T}_{3}\rangle_{\eta}(I)-\langle\tilde{T}_{3}\rangle_{\eta-2}(I)=\frac{(-1)^{\eta-1}5}{2},$
(49)

which are valid for any current $I$, torus shape $a$, angular momentum $m$, or
energy level $\eta$, if $a$ and $\eta$ are large enough and the current is in
the region where $\langle\tilde{T}_{3}\rangle_{\eta}(I)$ varies linearly, as
seen in Fig. 8. For all values of $a$ and $m$ (not only when $a$ and $\eta$
are sufficiently large), the functions $\langle\tilde{T}_{3}\rangle_{\eta}(I)$
are periodic, of period $2I_{\text{s}}$.

In Fig. 7 we plot the energy levels $E_{\eta}(I)$ as functions of $I/I_{0}$
for a few values of $a$ and $m$. We clearly see same periodicity, of period
$2I_{\text{s}}$, as in Fig. 7. For the levels where $a$ and $\eta$ are
sufficiently large, the variation of $E_{\eta}(I)$ is proportional to
$I\langle\tilde{T}_{3}\rangle_{\eta}(I)$, as we proved in Section 3 (Eq. 47)
and we can see in Fig. 10. The numerical calculations show that the validity
of the proportionality relation (47) extends over a larger interval of
currents than that on which the approximations of Section 3 are valid. At
$I=lI_{\text{s}}$ ($l>0$), the functions $E_{\eta}(I)$ are continuous, but
their derivatives change abruptly.

In Fig. 11 we plot the first coefficients of the expansion of the energy
eigenstates in the basis formed by the momentum eigenstates $\\{{\cal
F}^{(\Lambda)}_{n,m}\\}$. In all the plots, at $I=0$ we obtain again the
coefficients from Ref. [43], but, as $I$ increases, the wavefunctions
significantly change.

##  5 Conclusions

In this work, we studied the quantum properties of a particle on a torus in
the presence of a filiform current flowing along the rotational axis of the
system. We observed the quantization of the energy levels and their dependence
on the current intensity and we calculated the expectation values of the
projection of the toroidal dipole operator on the rotation axis in the
stationary states. We showed that the expectation values are zero if the
current is zero, but they become non-zero in the presence of non-zero current,
providing a method to observe and measure the quantum toroidal dipole.

Our analytical and numerical calculations revealed several important features
of the system. First, we found that the energy levels and the expectation
values of the toroidal dipoles are periodic functions of the current
intensity, with period $2I_{\text{s}}$, where $I_{\text{s}}$ depends only on
the aspect ratio of the torus. Second, we observed that the derivative of the
energy with respect to the current intensity, as well as the expectation
values of the toroidal dipole exhibit sharp changes of sign at integer
multiples of $I_{\text{s}}$. Third, we found that the toroidal dipole moment
decreases linearly with current outside of these narrow regions, with a
universal change in its scaled, dimensionless value equal to $-5/4$ over half-
period. Fourth, we showed that the differences between the toroidal dipole
moments of consecutive and second consecutive energy levels exhibit universal
forms. Fifth (and maybe the most significant), in general, the changes, due to
transitions between energy levels, of the expectation values of the toroidal
dipole are quantized, the scaled, dimensionless value of the quanta being
$5/2$.

Our results provide a theoretical framework for understanding the interaction
between a particle with a toroidal dipole moment and an external current. They
also pave the way for the experimental observation and measurement of toroidal
dipoles in quantum systems.

##  6 Acknowledgments

The work was financially supported by the projects PN 23210101/2023 and ELI-RO
81-44/2020.

##  Appendix A Analytical expression for $I_{n}(a)$

We want find analytical expressions for the integrals

$\displaystyle I_{n}(a)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}\frac{e^{in\theta}}{a+\cos\theta}d\theta.$
(50)

For this, we write

$\displaystyle\frac{1}{2}(I_{n+1}+I_{n-1})+aI_{n}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{e^{in\theta}\frac{1}{2}(e^{i\theta}+e^{-i\theta})+e^{in\theta}a}{a+\cos(\theta)}\differential\theta=\frac{1}{2\pi}\int_{0}^{2\pi}e^{in\theta}\differential\theta=\delta_{n}$
(51)

We obtain two cases: Case i) $n=0$

$\displaystyle\frac{1}{2}(I_{1}(a)+I_{-1}(a))+aI_{0}(a)=1$ (52)

Case ii) $n\neq 0$

$\displaystyle\frac{1}{2}(I_{n+1}(a)+I_{n-1}(a))+aI_{n}(a)=0$ (53)

From cases i) and ii) we get a recurrence relation

$I_{n+1}(a)+2aI_{n}(a)+I_{n-1}(a)=0$ (54)

which has the characteristic equation for $n>0$

$\displaystyle\lambda^{2}+2a\lambda+1=0.$ (55)

From Eq. (55) we obtain

$I_{n}(a)=C_{1}(a)\left(-a-\sqrt{a^{2}-1}\right)^{n}+C_{2}(a)\left(-a+\sqrt{a^{2}-1}\right)^{n}.$

We need $C_{1}(a)$ and $C_{2}(a)$ in order to find a general relation for
$I_{n}(a)$. For this, we make the following observations:

$-a-\sqrt{a^{2}-1}=\frac{1}{-a+\sqrt{a^{2}-1}}\quad{\rm and}\quad
I_{n}=I_{n}^{*}=I_{-n}.$ (56)

Then, we have two possibilities:

$\displaystyle I_{n}(a)=C(a)\left(-a-\sqrt{a^{2}-1}\right)^{n},\quad{\rm or}$
(57a) $\displaystyle I_{n}(a)=C(a)\left(-a+\sqrt{a^{2}-1}\right)^{n}.$ (57b)

Plugging the forms (57) into the condition (52) we get $C(a)=\pm
1/\sqrt{a^{2}-1}$, where the “$-$” and “$+$” signs correspond to the
possibilities (57a) and (57b), respectively. Since $I_{0}(a)>0$, it follows
that

$C(a)=\frac{1}{\sqrt{a^{2}-1}}\quad{\rm and}\quad
I_{n}(a)=\frac{\left(-a+\sqrt{a^{2}-1}\right)^{|n|}}{\sqrt{a^{2}-1}}.$ (58)

##  Appendix B The matrix elements of $\hat{T}_{3}$ in the momentum basis

We compute the matrix elements of the $z$ component of the toroidal dipole
operator $\hat{T}_{3}$ (29). Using the definitions (19) and (29), we obtain

$\displaystyle\hat{T}_{3}^{({\bf j},0)}$ $\displaystyle\equiv$
$\displaystyle\frac{i\hbar
R}{10m_{p}}\left[\frac{3a\cos(2\theta)+4\left(a^{2}+1\right)\cos\theta+5a}{2a}\frac{\partial}{\partial\theta}\right],$
$\displaystyle\hat{T}_{3}^{({\bf j},I)}$ $\displaystyle\equiv$
$\displaystyle-\frac{i\hbar
R}{10m_{p}}\left[\frac{3a\cos(2\theta)+4\left(a^{2}+1\right)\cos\theta+5a}{2a}\frac{iqR}{\hbar
a}A_{\theta}\right],$ $\displaystyle\left\langle{\cal
F}^{(\Lambda)}_{n_{1},m_{1}}\left|\tilde{T}_{3}^{({\bf j},0)}\right|{\cal
F}^{(\Lambda)}_{n_{2},m_{2}}\right\rangle\equiv
i\delta_{m_{1},m_{1}}\left[K_{1}(n_{2},n_{2}-n_{1})+K_{2}(n_{2}-n_{1})\right],\quad{\rm
where}$ (59) $\displaystyle K_{1}(n_{2},n)=\frac{1}{4\pi
a}\int_{0}^{2\pi}d\theta
e^{in\theta}\left\\{in_{2}\left[3a\cos(2\theta)+4\left(a^{2}+1\right)\cos\theta+5a\right]+(3a\cos\theta-a^{2}+2)\sin\theta\right\\}$
$\displaystyle K_{2}(n)=\frac{a^{2}-1}{4\pi}\int_{0}^{2\pi}d\theta
e^{in\theta}\frac{\sin\theta}{(a+\cos\theta)}d\theta$

(in the second and third line we used $n\equiv n_{2}-n_{1}$). For $K_{1}$ we
obtain immediately

$\displaystyle K_{1}(n_{2},n)$ $\displaystyle=$
$\displaystyle\frac{3i}{4}\left(n_{2}-\frac{1}{2}\right)\delta_{n+2,0}+\frac{3i}{4}\left(n_{2}+\frac{1}{2}\right)\delta_{n-2,0}+\frac{i}{4a}\left[4n_{2}\left(a^{2}+1\right)+a^{2}-2\right]\delta_{n+1,0}$
(60)
$\displaystyle+\frac{i}{4a}\left[4n_{2}\left(a^{2}+1\right)-a^{2}+2\right]\delta_{n-1,0}+\frac{5in_{2}}{2}\delta_{n,0}.$

For $K_{2}(n)$, we first observe that it is imaginary and

$K_{2}(-n)=K_{2}^{*}(n)=-K_{2}(n).$ (61)

Then, we may write

$\displaystyle aK_{2}(n)+\frac{1}{2}[K_{2}(n-1)+K_{2}(n+1)]$ $\displaystyle=$
$\displaystyle i\frac{a^{2}-1}{4}(\delta_{n-1,0}-\delta_{n+1,0}).$

Directly computing

$K_{2}(0)=0,\quad
K_{2}(1)=-i\frac{(a^{2}-1)}{2}\left(\sqrt{a^{2}-1}-a\right),\quad
K_{2}(2)=-i\frac{(a^{2}-1)}{2}\left(\sqrt{a^{2}-1}-a\right)^{2},$ (62)

we can check that Eq. (B) is satisfied. Then, for any $n>1$ we have the
recurrence relation

$a\cdot K_{2}(n)+\frac{1}{2}[K_{2}(n-1)+K_{2}(n+1)]=0,$ (63)

with the solutions

$K_{2}(n)=\begin{cases}C_{1}\left(-a-\sqrt{a^{2}-1}\right)^{n}\quad{\rm or}\\\
C_{2}\left(-a+\sqrt{a^{2}-1}\right)^{n}.\end{cases}$ (64)

Comparing (64) with (62) we finally obtain

$K_{2}(n)=-i\,{\rm
sgn}(n)\frac{(a^{2}-1)}{2}\left(\sqrt{a^{2}-1}-a\right)^{|n|},$ (65)

where we consider that ${\rm sgn}(0)=0$. Plugging Eqs. (60) and (65) into (59)
we obtain

$\displaystyle\langle{\cal F}^{(\Lambda)}_{n_{2}-n,m_{1}}|\tilde{T}_{3}^{({\bf
j},0)}|{\cal
F}^{(\Lambda)}_{n_{2},m_{2}}\rangle=-\delta_{m_{1},m_{2}}\left\\{\frac{3}{4}\left[n_{2}(\delta_{n+2,0}+\delta_{n-2,0})-\frac{1}{2}(\delta_{n+2,0}-\delta_{n-2,0})\right]+n_{2}\frac{a^{2}+1}{a}(\delta_{n+1,0}+\delta_{n-1,0})\right.$
$\displaystyle\left.+\frac{a^{2}-2}{4a}(\delta_{n+1,0}-\delta_{n-1,0})+\frac{5n_{2}}{2}\delta_{n,0}-{\rm
sgn}(n)\frac{(a^{2}-1)}{2}\left(\sqrt{a^{2}-1}-a\right)^{|n|}\right\\}$ (66)

For the electromagnetic contribution we have

$\displaystyle\langle{\cal F}^{(\Lambda)}_{n_{2}-n,m_{1}}|\tilde{T}_{3}^{({\bf
j},I)}|{\cal
F}^{(\Lambda)}_{n_{2},m_{2}}\rangle=\frac{I}{4a^{2}I_{0}}\delta_{m_{1},m_{2}}\left\\{\log(\frac{2aL}{R})\left[\frac{3a}{4}\left(\delta_{n+3,0}+\delta_{n-3,0}\right)+\left(a^{2}+1\right)\left(\delta_{n+2,0}+\delta_{n-2,0}\right)\right.\right.$
$\displaystyle\left.+\frac{13a}{4}\left(\delta_{n+1,0}+\delta_{n-1,0}\right)+2\left(a^{2}+1\right)\delta_{n,0}\right]-\left[\frac{3a}{4}\left(I^{(\ln)}_{n+3}(a)+I^{(\ln)}_{n-3}(a)\right)+\left(a^{2}+1\right)\left(I^{(\ln)}_{n+2}(a)+I^{(\ln)}_{n-2}(a)\right)\right.$
$\displaystyle\left.\left.+\frac{13a}{4}\left(I^{(\ln)}_{n+1}(a)+I^{(\ln)}_{n-1}(a)\right)+2\left(a^{2}+1\right)I^{(\ln)}_{n}(a)\right]\right\\},$
(67)

with the notation

$\displaystyle
I^{(ln)}_{n}(a)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{in\theta}\log(a+\cos\theta)d\theta=\left.\frac{I_{n-1}(a)-I_{n+1}(a)}{n}\right|_{n\neq
0}-\log(\frac{2}{a+\sqrt{a^{2}-1}})\delta_{n}.$ (68)

##  Appendix C A simple 2D model

The 2D matrix equation

$\left(\begin{array}[]{cc}1&\epsilon\\\ \epsilon&1\end{array}\right){\bf
v}=\alpha{\bf v},$ (69)

has two solutions, one symmetric ($+$) and one antisymmetric ($-$),

${\bf v}^{(\pm)}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\ \pm
1\end{array}\right)\quad{\rm with}\quad\alpha^{(\pm)}=1\pm\epsilon.$ (70)

We assume that $\epsilon$ is a small and positive variable and we add a
perturbation $\pm\delta$ to the diagonal elements of the matrix (including the
variables $\epsilon$ and $\delta$ in the notations of the eigenvalues and
eigenvectors),

$\left(\begin{array}[]{cc}1+\delta&\epsilon\\\
\epsilon&1-\delta\end{array}\right){\bf
v}(\epsilon,\delta)=\alpha(\epsilon,\delta){\bf v}(\epsilon,\delta).$ (71)

Equation (71) gives the solutions

$\alpha^{(\pm)}(\epsilon,\delta)=1\pm\sqrt{\epsilon^{2}+\delta^{2}},\ {\bf
v}^{(\pm)}(\epsilon,\delta)=N^{(\pm)}(\epsilon,\delta)\left(\begin{array}[]{c}1\\\
-\frac{\delta}{\epsilon}\pm\sqrt{1+\frac{\delta^{2}}{\epsilon^{2}}}\end{array}\right),\
N^{(\pm)}(\epsilon,\delta)=\frac{1}{\sqrt{2\sqrt{1+\frac{\delta^{2}}{\epsilon^{2}}}\left(\sqrt{1+\frac{\delta^{2}}{\epsilon^{2}}}\mp\frac{\delta}{\epsilon}\right)}}.$
(72)

We represent the vectors ${\bf v}^{(\pm)}(\epsilon,\delta)$ as complex numbers

$\zeta^{(\pm)}(\epsilon,\delta)\equiv
N^{(\pm)}(\epsilon,\delta)\left[1+i\left(-\frac{\delta}{\epsilon}\pm\sqrt{1+\frac{\delta^{2}}{\epsilon^{2}}}\right)\right]\equiv
e^{i\theta^{(\pm)}(\epsilon,\delta)},\quad{\rm
where}\quad\tan\left(\theta^{(\pm)}(\epsilon,\delta)\right)\equiv-\frac{\delta}{\epsilon}\pm\sqrt{1+\frac{\delta^{2}}{\epsilon^{2}}},$
(73)

and in Fig. 12 we plot the angles $\theta^{(\pm)}(\epsilon,\delta)$. At
$\delta=0$, we recover the symmetric and antisymmetric solutions,
$\theta^{(+)}(\epsilon,0)=\pi/4$ and $\theta^{(-)}(\epsilon,0)=-\pi/4$,
respectively, whereas at $\delta/\epsilon\to\pm\infty$ we have
$\lim_{(\delta/\epsilon)\to\infty}\theta^{(+)}(\delta,\epsilon)=0$,
$\lim_{(\delta/\epsilon)\to\infty}\theta^{(-)}(\delta,\epsilon)=-\pi/2$,
$\lim_{(\delta/\epsilon)\to-\infty}\theta^{(+)}(\delta,\epsilon)=\pi/2$, and
$\lim_{(\delta/\epsilon)\to-\infty}\theta^{(-)}(\delta,\epsilon)=0$.

Figure 12: The angles $\theta^{(+)}$ (red line) and $\theta^{(-)}$ (blue
line), as functions of $\delta/\epsilon$. At $\delta=0$, $\theta^{(+)}$
corresponds to the symmetric solution, whereas $\theta^{(-)}$ corresponds to
the antisymmetric one.

##  Appendix D Taylor expansions for $I_{n}(a)$ and $I_{n}^{(\ln^{2})}(a)$

It is useful to calculate the Taylor expansions in the lowest order of
$I_{n}(a)$ and $I_{n}^{(ln^{2})}(a)$, when $a\gg 1$. For $I_{n}(a)$ we simply
obtain

$I_{n}(a)\stackrel{{\scriptstyle a\gg
1}}{{\approx}}\frac{1}{a}\left(\frac{-1}{2a}\right)^{|n|}$ (74)

The calculation for $I_{n}^{(\ln^{2})}(a)$ is more complicated. First, we
write

$\displaystyle I_{n}^{(ln^{2})}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}e^{in\theta}\left[\ln^{2}a+\frac{\cos^{2}\theta}{a^{2}}+\frac{\cos^{4}\theta}{(2a^{2})^{2}}+\frac{\cos^{6}\theta}{(3a^{3})^{2}}+\ldots+2\ln
a\left(\frac{\cos\theta}{a}-\frac{\cos^{2}\theta}{2a^{2}}+\frac{\cos^{3}\theta}{3a^{3}}-\frac{\cos^{4}\theta}{4a^{4}}+\ldots\right)\right.$
(75)
$\displaystyle\left.+\sum_{m=3}^{\infty}(-1)^{m-2}\frac{\cos^{m}\theta}{a^{m}}\sum_{k=1}^{m-1}\frac{1}{(m-k)}\frac{1}{k}\right]d\theta$

Then, in the lowest order of approximation for each term $m$, we retain

$\cos^{m}\theta=\left(\frac{e^{i\theta}+e^{-i\theta}}{2}\right)^{m}\approx\frac{e^{im\theta}+e^{-im\theta}}{2}$
(76)

since the other terms give higher orders of approximation. Plugging (76) into
(75) and using the identity

$\sum_{k=1}^{m-1}\frac{m}{k(m-k)}=\sum_{k=1}^{m-1}\left(\frac{1}{k}+\frac{1}{m-k}\right)=2\sum_{k=1}^{m-1}\frac{1}{k}$

we get

$\displaystyle I_{n}^{(ln^{2})}(a)$ $\displaystyle\stackrel{{\scriptstyle a\gg
1}}{{\approx}}$
$\displaystyle\ln^{2}a\,\delta_{n}+\sum_{m=1}^{\infty}\frac{\delta_{n+2m}+\delta_{n-2m}}{2m^{2}a^{2m}}+\ln
a\left(\frac{\delta_{n+1}+\delta_{n-1}}{a}-\frac{\delta_{n+2}+\delta_{n-2}}{2a^{2}}\right)$
(77)
$\displaystyle+\sum_{m=3}^{\infty}(-1)^{m}\left[\left(\sum_{k=1}^{m-1}\frac{1}{k}\right)-\ln
a\right]\frac{\delta_{n+m}+\delta_{n-m}}{ma^{m}}.$

Therefore,

$I_{0}^{(\ln^{2})}\stackrel{{\scriptstyle a\gg 1}}{{\approx}}\ln a,\quad
I_{1}^{(\ln^{2})}\stackrel{{\scriptstyle a\gg 1}}{{\approx}}\frac{\ln
a}{a},\quad I_{2}^{(\ln^{2})}\stackrel{{\scriptstyle a\gg
1}}{{\approx}}-\frac{\ln a}{2a^{2}}+\frac{1}{2a^{2}}$ (78)

## References

  * [1] Ya. B. Zeldovich.  Electromagnetic interaction with parity violation.  Sov. Phys.-JETP, 6:1184, 1958.
  * [2] V M Dubovik and A A Cheshkov.  Form-factors and multipoles in electromagnetic interactions.  Sov. Phys. JETP, 24, 1965.
  * [3] V M Dubovik and A A Cheshkov.  Multipole expansion in classical and quantum field theory and radiation.  Sov. J. Particles Nucl., 5:318, 1974.
  * [4] S. Nanz.  Toroidal Multipole Moments in Classical Electrodynamics: An Analysis of their Emergence and Physical Significance.  BestMasters. Springer Fachmedien Wiesbaden, 2016.
  * [5] Nikitas Papasimakis, V.A. Fedotov, Vassili Savinov, T.A. Raybould, and N.I. Zheludev.  Electromagnetic toroidal excitations in matter and free space.  Nature Mat., 15:263, 2016.
  * [6] V. V. Flambaum.  Nuclear anapole moment and tests of the standard model.  AIP Conference Proceedings, 477(1):14, 1999.
  * [7] C. S. Wood, S. C. Bennett, D. Cho, B. P. Masterson, J. L. Roberts, C. E. Tanner, and C. E. Wieman.  Measurement of parity nonconservation and an anapole moment in cesium.  Science, 275(5307):1759, 1997.
  * [8] C. M. Ho and R. J. Scherrer.  Anapole dark matter.  Phys. Lett. B, 722:341, 2013.
  * [9] Luis G. Cabral-Rosetti, Myriam Mondragón, and Esteban Reyes-Pérez.  Anapole moment of the lightest neutralino in the cMSSM.  Nucl. Phys. B, 907:1, 2016.
  * [10] E. E. Radescu.  On the electromagnetic properties of majorana fermions.  Phys. Rev. D, 32:1266, 1985.
  * [11] Dubovik Vladimir M. and Valentin E. Kuznetsov.  The toroid dipole moment of the neutrino.  International Journal of Modern Physics A, 13:5257, 1998.
  * [12] Charles Kittel.  Theory of the structure of ferromagnetic domains in films and small particles.  Phys. Rev., 70:965–971, 1946.
  * [13] V.M. Dubovik and V.V. Tugushev.  Toroid moments in electrodynamics and solid-state physics.  Phys. Rep., 187:145, 1990.
  * [14] A. Ceulemans, L. F. Chibotaru, and P. W. Fowler.  Molecular anapole moments.  Phys. Rev. Lett., 80:1861, 1998.
  * [15] Daniel Khomskii.  Classifying multiferroics: Mechanisms and effects.  Physics, 2:20, 2009.
  * [16] Aleksandr Pyatakov and Anatoly Zvezdin.  Magnetoelectric and multiferroic media.  Physics-Uspekhi, 55:557, 2012.
  * [17] N. Tolstoi and A. Spartakov.  Aromagnetism: A new type of magnetism.  J. Exp. Theor. Phys. Lett., 52:161, 01 1990.
  * [18] V.A. Fedotov, K. Marinov, Allan Boardman, and N.I. Zheludev.  On the aromagnetism and anapole moment of anthracene nanocrystals.  New J. Phys., 9:95, 04 2007.
  * [19] Pierre Tolédano, Dmitry Khalyavin, and Laurent Chapon.  Spontaneous toroidal moment and field-induced magnetotoroidic effects in $Ba_{2}CoGe_{2}O_{7}$.  Phys. Rev. B, 84:094421, 2011.
  * [20] Takahiro Shimada, Yuuki Ichiki, Gen Fujimoto, Le Lich, Tao Xu, Jie Wang, and Hiroyuki Hirakata.  Ferrotoroidic polarons in antiferrodistortive $SrTiO_{3}$.  Phys. Rev. B, 101, 2020.
  * [21] Yu Popov, A. Kadomtseva, G. Vorob’ev, VA Timofeeva, D. Ustinin, Anatoly Zvezdin, and M. Tegeranchi.  Magnetoelectric effect and toroidal ordering in $Ga_{2-x}Fe_{x}O_{3}$.  Journal of Experimental and Theoretical Physics, 87:146, 07 1998\.
  * [22] Jannis Lehmann, Claire Donnelly, Peter Derlet, Laura Heyderman, and Manfred Fiebig.  Poling of an artificial magneto-toroidal crystal.  Nature Nanotechnology, 14:141, 2019.
  * [23] A. Ahmadivand, B. Gerislioglu, and Z. Ramezani.  Toroidal Metamaterials.  Springer, Cham, Springer Nature Switzerland AG.
  * [24] I. I. Naumov, L. Bellaiche, and H. Fu.  Unusual phase transitions in ferroelectric nanodisks and nanorods.  432:737, 2004.
  * [25] K. Marinov, A. D. Boardman, V. A. Fedotov, and N. Zheludev.  Toroidal metamaterial.  New J. Phys., 9:324, 2007.
  * [26] T. Kaelberer, V. A. Fedotov, N. Papasimakis, D. P. Tsai, and Zheludev N. I.  Toroidal dipolar response in a metamaterial.  Science, 330:1510, 2010.
  * [27] Alexey A. Basharin, Maria Kafesaki, Eleftherios N. Economou, Costas M. Soukoulis, Vassili A. Fedotov, Vassili Savinov, and Nikolay I. Zheludev.  Dielectric metamaterials with toroidal dipolar response.  Phys. Rev. X, 5:011036, Mar 2015.
  * [28] Y. Yang and S. I. Bozhevolnyi.  Nonradiating anapole states in nanophotonics: from fundamentals to applications.  Nanotechnology, 30:204001, 2019.
  * [29] V Savinov, N Papasimakis, D P Tsai, and N I Zheludev.  Optical anapoles.  Communications Physics, 2:10, 2019.
  * [30] Nahid Talebi, Surong Guo, and Peter Aken.  Theory and applications of toroidal moments in electrodynamics: Their emergence, characteristics, and technological relevance.  Nanophotonics, 7, 1 2017.
  * [31] E. Gurvitz, K. Ladutenko, P. Dergachev, A. Evlyukhin, A. Miroshnichenko, and A. Shalin.  The high-order toroidal moments and anapole states in all‐dielectric photonics.  Laser and Photonics Reviews, 13:1800266, 2019.
  * [32] Alexandre Zagoskin, A. Chipouline, Evgeni Il’ichev, Robert Johansson, and Franco Nori.  Toroidal qubits: Naturally-decoupled quiet artificial atoms.  Scientific Reports, 5, 2014.
  * [33] Y. He, G. Guo, T. Feng, Y. Xu, and A. E. Miroshnichenko.  Toroidal dipole bound states in the continuum.  Phys. Rev. B, 98:161112, 2018.
  * [34] W. Huang, S. Yin, K. Wang, Y. Zhang, and J. Han.  Robust and broadband integrated terahertz coupler conducted with adiabatic following.  New J. Phys., 21:113004, 2019.
  * [35] Yu Liang, Fengchun Zhang, Xu-Guang Huang, and Huicheng Zheng.  Polarization-controlled triple-band absorption in all-metal nanostructures with magnetic dipoles and anapole responses.  Appl. Phys. Express, 12:062014, 2019.
  * [36] N. Pavlov, I. Stenishchev, A. Ospanova, P. Belov, P. Kapitanova, and A. Basharin.  Toroidal dipole mode observation in situ.  Physica Status Solidi (b), 257:1900406, 2020.
  * [37] M. Li, Q. Ma, A. Luo, and W. Hong.  Multiple toroidal dipole symmetry-protected bound states in the continuum in all-dielectric metasurfaces.  Optics & Laser Technology, 154:108252, 2022.
  * [38] T. Yezekyan, V. A. Zenin, J. Beermann, and S. I. Bozhevolnyi.  Anapole states in gap-surface plasmon resonators.  Nano Lett., 22:6098–6104, 2022.
  * [39] S. B. Saadatmand, V. Ahmadi, and S. M. Hamidi.  Quasi-bic based all-dielectric metasurfaces for ultra-sensitive refractive index and temperature sensing.  Sci. Rep., 13:20625, 2023.
  * [40] A Costescu and E.E Radescu.  Dynamic toroid polarizability of atomic hydrogen.  Ann. Phys., 209:13, 1991.
  * [41] D. V. Anghel.  Mathematical considerations regarding the toroidal momentum operator.  J. Phys. A: Math. Gen., 30:3515–3525, 1997.
  * [42] D. V. Anghel and A. T. Preda.  Quantized toroidal dipole eigenvalues in nano-systems.  J. Phys.: Conf. Ser., 2090:012151, 2021.
  * [43] M. Dolineanu, A. T. Preda, and D. V. Anghel.  The toroidal dipole operator in nanostructures.  Physica A, 598:127377, 2022.
  * [44] D. V. Anghel and M. Dolineanu.  The eigenvalues and eigenfunctions of the toroidal dipole operator.  Physica Scripta, 98:015223, 2022.
  * [45] M. Dolineanu, A. T. Preda, and D. V. Anghel.  Corrigendum to “the toroidal dipole operator in nanostructures” [physica a 598 (2022) 127377].  Physica A, 617:128684, 2023.