∎

11institutetext: Vienna Institute of Technology
11email<EMAIL_ADDRESS>
_Vienna University of Technology_

# On the Motion of Billiards in Ellipses

Hellmuth Stachel

###### Abstract

For billiards in an ellipse with an ellipse as caustic, there exist canonical
coordinates such that the billiard transformation from vertex to vertex is
equivalent to a shift of coordinates. A kinematic analysis of billiard motions
paves the way to an explicit canonical parametrization of the billiard and
even of the associated Poncelet grid. This parametrization uses Jacobian
elliptic functions to the numerical eccentricity of the caustic as modulus.

###### Keywords:

billiard billiard motion confocal conics, elliptic functions

###### MSC:

51N20 53A17 33E05, 22E30

††journal: European Journal of Mathematics

##  1 Introduction

Already for two centuries, billiards in ellipses have attracted the attention
of mathematicians, beginning with J.-V. Poncelet and A. Cayley. The assertion
that one periodic billiard inscribed in an ellipse $e$ and tangent to a
confocal ellipse $c$ implies a one-parameter family of such polygons, is known
as the standard example of a Poncelet porism. It was Cayley who derived
analytical conditions for such a pair $(e,c)$ of ellipses.

In 2005 S. Tabachnikov published a book on billiards from the viewpoint of
integrable systems Tabach . The book DR_Buch  and various papers by V.
Dragović and M. Radnović addressed billiards in conics and quadrics within the
framework of dynamical systems.

Computer animations carried out by D. Reznik, stimulated a new vivid interest
on this well studied topic, where algebraic and analytic methods are meeting.
Originally, Reznik’s experiments focused on billiard motions in ellipses,
i.e., on the variation of billiards with a fixed circumscribed ellipse $e$ and
a fixed caustic $c\,$. He published a list of more than 80 numerically
detected invariants in 80  and contributed, together with his coauthors R.
Garcia, J. Koiller and M. Helman, several proofs. Other authors like A.
Akopyan, M. Bialy, A. Chavez-Caliz, R. Schwartz, and S. Tabachnikov published
several proofs and found more invariants (e.g., in Ako-Tab ; Bialy-Tab ;
Chavez ).

For a long time, at least since Jacobi’s proof of the Poncelet theorem on
periodic billiards (see further references in (DR_russ, , p. 320)), it has
been wellknown that there is a tight connection between billiards and elliptic
functions (note also (Duistermaat, , Sect. 11.2) and Bobenko ). On the other
hand, S. Tabachnikov proved in his book Tabach  the existence of canonical
parameters on ellipses with the property that the billiard transformation
between consecutive vertices $P_{i}\mapsto P_{i+1}$ of a billiard acts like a
shift on the parameters.

The goal of this paper is to prove that Jacobian elliptic functions with the
numerical eccentricity of the caustic $c$ as modulus pave the way to canonical
coordinates on the ellipse $e\,$. This is a consequence of a kinematic
analysis of the billiard motion. It yields an infinitesimal transformation in
the plane which preserves a family of confocal ellipses while it permutes the
confocal hyperbolas as well as the tangents of the caustic. Integration
results in a group of transformations with a canonical parameter. In terms of
elliptic functions, we also obtain a mapping that sends a square grid together
with the diagonals to a Poncelet grid.111 Recently, parametrizations of
confocal conics in termins of elliptic functions were also presented in
Bobenko , but not from the viewpoint of billiards. The paper concludes with
applying the results of the velocity analysis to a few invariants of periodic
billiards. These invariants deal mainly with the distances which occur on each
side between the contact point with the caustic and the endpoints.

##  2 Confocal conics and billiards

At the begin, we recall a few properties of confocal conics. A family of
confocal central conics is given by

$\frac{x^{2}}{a^{2}+k}+\frac{y^{2}}{b^{2}+k}=1,\ \mbox{where}\
k\in{\mathbb{R}}\setminus\\{-a^{2},-b^{2}\\}$ (2.1)

serves as a parameter in the family. All these conics share the focal points
$F_{1,2}=(\pm d,0)$, where $d^{2}:=a^{2}-b^{2}$.

The confocal family sends through each point $P$ outside the common axes of
symmetry two orthogonally intersecting conics, one ellipse and one hyperbola
(Conics, , p. 38). The parameters $(k_{e},k_{h})$ of these two conics define
the elliptic coordinates of $P$ with

$-a^{2}<k_{h}<-b^{2}<k_{e}\,.$

If $(x,y)$ are the cartesian coordinates of $P$, then $(k_{e},k_{h})$ are the
roots of the quadratic equation

$k^{2}+(a^{2}+b^{2}-x^{2}-y^{2})k+(a^{2}b^{2}-b^{2}x^{2}-a^{2}y^{2})=0,$ (2.2)

while conversely

$x^{2}=\frac{(a^{2}+k_{e})(a^{2}+k_{h})}{d^{2}}\,,\quad
y^{2}=-\frac{(b^{2}+k_{e})(b^{2}+k_{h})}{d^{2}}\,.$ (2.3)

Suppose that $(a,b)$ in (2.1) are the semiaxes $(a_{c},b_{c})$ of the ellipse
$c$ with $k=0$. Then, for points $P$ on a confocal ellipse $e$ with semiaxes
$(a_{e},b_{e})$ and $k=k_{e}>0$, i.e., exterior to $c$, the standard
parametrization yields

$\begin{array}[]{c}P=(x,y)=(a_{e}\cos t,\,b_{e}\sin t),\ 0\leq
t<2\pi,\\\\[2.84526pt] \mbox{with}\ a_{e}^{2}=a_{c}^{2}+k_{e},\
b_{e}^{2}=b_{c}^{2}+k_{e}\,.\end{array}$ (2.4)

For the elliptic coordinates $(k_{e},k_{h})$ of $P$ follows from (2.2) that

$k_{e}+k_{h}=a_{e}^{2}\cos^{2}t+b_{e}^{2}\sin^{2}t-a_{c}^{2}-b_{c}^{2}.$

After introducing the respective tangent vectors of $e$ and $c$, namely

$\begin{array}[]{rcl}{\mathbf{t}}_{e}(t)&:=&(-a_{e}\sin t,\,b_{e}\cos
t),\\\\[2.27621pt] {\mathbf{t}}_{c}(t)&:=&(-a_{c}\sin t,\,b_{c}\cos
t),\end{array}\ \mbox{where}\
\|{\mathbf{t}}_{e}\|^{2}=\|{\mathbf{t}}_{c}\|^{2}+k_{e}\,,$ (2.5)

we obtain

$k_{h}=k_{h}(t)=-(a_{c}^{2}\sin^{2}t+b_{c}^{2}\cos^{2}t)=-\|{\mathbf{t}}_{c}(t)\|^{2}=-\|{\mathbf{t}}_{e}(t)\|^{2}+k_{e}$
(2.6)

and $\|{\mathbf{t}}_{e}(t)\|^{2}=k_{e}-k_{h}(t)\,$. Note that points on the
confocal ellipses $e$ and $c$ with the same parameter $t$ have the same
coordinate $k_{h}$. Consequently, they belong to the same confocal hyperbola
(Figure 1). Conversely, points of $e$ or $c$ on this hyperbola have a
parameter out of $\\{t,-t,\pi+t,\pi-t\\}$.

Let $\theta_{i}/2$ denote the angle between the tangents drawn from any point
$P_{i}\in e$ to $c$ and the tangent to $e$ at $P_{i}$ (Figures 2 or 3). Then
we obtain for $P_{i}=(a_{e}\cos t_{i},\,b_{e}\sin t_{i})$ with elliptic
coordinates $\left(k_{e},k_{h}(t_{i})\right)$

$\begin{array}[]{c}\sin^{2}{\displaystyle\frac{\theta_{i}}{2}}={\displaystyle\frac{k_{e}}{\|{\mathbf{t}}_{e}(t_{i})\|^{2}}}={\displaystyle\frac{k_{e}}{k_{e}-k_{h}(t_{i})}},\quad\tan{\displaystyle\frac{\theta_{i}}{2}}=\pm\sqrt{-{\displaystyle\frac{k_{e}}{k_{h}(t_{i})}}}\\\\[8.53581pt]
\mbox{and}\mskip
9.0mu\sin\theta_{i}=\pm{\displaystyle\frac{2\sqrt{-k_{e}k_{h}(t_{i})}}{k_{e}-k_{h}(t_{i})}}=\pm{\displaystyle\frac{2\|{\mathbf{t}}_{c}(t_{i})\|\sqrt{k_{e}}(\|{\mathbf{t}}_{e}(t_{i})\|^{2}}{\,}}.\end{array}$
(2.7)

For a proof see Sta_I . We can assume a counter-clockwise order of the
billiard. Hence, all exterior angles are positive.

Figure 1: Periodic billiard $P_{1}P_{2}\dots P_{9}$ in $e$ with the turning
number $\tau\\!=\\!1$ and the caustic $c$ as well as a part of the associated
Poncelet grid. The extended sides form a billiard with $\tau=2$ in $e^{(1)}$
and three triangles as billiards in $e^{(2)}$.

From (2.6) follows

$k_{h}=-\frac{a_{c}^{2}\tan^{2}t+b_{c}^{2}}{1+\tan^{2}t},\quad\mbox{hence}\quad\tan^{2}t(a_{c}^{2}+k_{h})=-b_{c}^{2}-k_{h}$

and furthermore

$\sin t\cos t=\frac{\tan
t}{1+\tan^{2}t}=\frac{\sqrt{-(b_{c}^{2}+k_{h})(a_{c}^{2}+k_{h})}}{a_{c}^{2}-b_{c}^{2}}=\frac{a_{h}\,b_{h}}{d^{2}}$
(2.8)

with $a_{h}$ and $b_{h}$ as semiaxes of the hyperbola corresponding to the
parameter $t$, i.e., $a_{h}^{2}=a_{c}^{2}+k_{h}$ and
$b_{h}^{2}=-(b_{c}^{2}+k_{h})$.

Let $\dots P_{1}P_{2}P_{3}\dots$ be a billiard in the ellipse $e$. Then the
extended sides intersect at points

$S_{i}^{(j)}:=\left\\{\begin{array}[]{rl}[P_{i-k-1},P_{i-k}]\cap[P_{i+k},P_{i+k+1}]&\mskip
9.0mu\mbox{for}\ j=2k,\\\\[1.70717pt]
[P_{i-k},P_{i-k+1}]\cap[P_{i+k},P_{i+k+1}]&\mskip 9.0mu\mbox{for}\
j=2k-1,\end{array}\right.$ (2.9)

where $i=\dots,1,2,3,\dots$ and $j=1,2,\dots$. These points are distributed on
different confocal conics: For fixed $j$, there are ellipses $e^{(j)}$ passing
through the points $S_{i}^{(j)}$. On the other hand, the points $S_{i}^{(2)}$,
$S_{i}^{(4)},\dots$ are located on the confocal hyperbola through $P_{i}$,
while $S_{i}^{(1)}$, $S_{i}^{(3)},\dots$ belong to the confocal hyperbola
through the contact point $Q_{i}$ between the side $P_{i}P_{i+1}$ and the
caustic $c\,$. This configuration is called the associated Poncelet grid
(Figure 1). For periodic billiards the sets of points $S_{i}^{(j)}$ and
associated conics are finite. The turning number $\tau$ of a periodic billiard
in $e$ with an ellipse as caustic counts how often one period of the billiard
surrounds the center $O$ of $e$ (note Figure 1).

For each billiard $P_{1}P_{2}\dots$ in $e$ with caustic $c$ there exists a
conjugate billiard $P_{1}^{\prime}P_{2}^{\prime}\dots$ in $e$ with the same
caustic. The axial scaling $c\to e$ maps the contact point $Q_{i}\in c$ of
$P_{i}P_{i+1}$ to $P_{i}^{\prime}$ while the inverse scaling brings $P_{i}$ to
the contact point $Q_{i-1}^{\prime}$ of $P_{i-1}^{\prime}P_{i}^{\prime}$ with
the caustic. The relation between these billiards is symmetric. For further
details see (Sta_I, , Sect. 3.2).

##  3 Velocity analysis

Let the first vertex of a billiard $P_{1}P_{2}\dots$ move smoothly along the
circumscribed ellipse $e$. Then this induces a continuous variation of all
other vertices along $e$ and also of the intersection points $S_{i}^{(j)}$
along $e^{(j)}$. We call this a billiard motion, though it neither preserves
angles or distances nor is an affine or projective motion.

According to Graves’s construction (Conics, , p. 47), we can conceive the
periodic billiard $P_{1}P_{2}\dots P_{N}$ as a flexible chain of fixed total
length $L_{e}$ and the caustic $c$ as a fixed non-circular chain wheel. The
vertices $P_{1},P_{2},\dots$ move along $e$ and relative to the chain such
that they keep the chain strengthened, while the chain contacts $c$ only at
the single points $Q_{1},Q_{2},\dots,Q_{N}$.

Let us pick out a single vertex $P_{2}$ (see Figure 2). In the language of
kinematics, the line spanned by the straight segment $Q_{1}P_{2}$ rolls at
$Q_{1}$ on $c$ (= fixed polode) while point $P_{2}$ moves along the line (=
moving polode) with the velocity vector ${\mathbf{v}}_{t_{1}}$ . The
instantaneous rotation about $Q_{1}$ with the angular velocity $\omega_{1}$
assigns to $P_{2}$ a velocity vector ${\mathbf{v}}_{n_{1}}$ orthogonal to
$Q_{1}P_{2}$ in order to keep the vector of absolute velocity of $P_{2}\,$,
namely ${\mathbf{v}}_{2}={\mathbf{v}}_{t_{1}}+{\mathbf{v}}_{n_{1}}$, tangent
to the ellipse $e$.

Figure 2: Graves’s string construction of an ellipse $e$ confocal to $c$.

Similarly, we have a second decomposition
${\mathbf{v}}_{2}={\mathbf{v}}_{t_{2}}+{\mathbf{v}}_{n_{2}}$, since at the
same time the line $[Q_{2},P_{2}]$ rotates about $Q_{2}$ with the angular
velocity $\omega_{2}$, while $P_{2}$ moves relative to this line. Due to the
constant length of the chain, the tangential components in these two
decompositions must be of equal lengths
$\|{\mathbf{v}}_{t_{2}}\|=\|{\mathbf{v}}_{t_{1}}\|$. Since the tangent $t_{P}$
to $e$ at $P_{2}$ bisects the exterior angle of $Q_{1}P_{2}Q_{2}$, the second
decomposition is symmetric w.r.t. $t_{P}$ to the first one. From
$\|{\mathbf{v}}_{n_{2}}\|=\|{\mathbf{v}}_{n_{1}}\|$ follows for the distances
$r_{2}:=\overline{P_{2}Q_{1}}$ and $l_{2}:=\overline{P_{2}Q_{2}}$

$l_{2}\,\omega_{2}=r_{2}\,\omega_{1},\quad\mbox{or}\quad\frac{\omega_{1}}{\omega_{2}}=\frac{l_{2}}{r_{2}},$
(3.1)

and similarly for all other vertices. If the billiard is $N$-periodic, then
the product of all ratios $l_{i}/r_{i}$ for $i=1,\dots,N$ yields

$\frac{l_{1}}{r_{1}}\cdot\frac{l_{2}}{r_{2}}\cdots\frac{l_{N}}{r_{N}}=\frac{\omega_{N}}{\omega_{1}}\cdot\frac{\omega_{1}}{\omega_{2}}\cdots\frac{\omega_{N-1}}{\omega_{N}}=1,$

which results in the equation

$l_{1}l_{2}\dots l_{N}=r_{1}r_{2}\dots r_{N}.$ (3.2)

listed as k116 in (80, , Table 2).

Figure 3 shows a graphical velocity analysis for the billiard motion of a
5-sided periodic billiard in $e$. We can begin this analysis by choosing an
arbitrary length for the arrow representing the velocity vector
${\mathbf{v}}_{2}$ of $P_{2}$. This defines the two components
${\mathbf{v}}_{t_{2}}$ and ${\mathbf{v}}_{n_{2}}$, where the latter determines
the angular velocity $\omega_{2}$ of the side $P_{2}P_{3}$ and furtheron the
absolute velocity ${\mathbf{v}}_{3}$ of $P_{3}$. This can be continued. From
now on, we denote the norms
$\|{\mathbf{v}}_{t_{1}}\|=\|{\mathbf{v}}_{t_{2}}\|$ and
$\|{\mathbf{v}}_{n_{1}}\|=\|{\mathbf{v}}_{n_{2}}\|$ of the respective
components of the velocity vector ${\mathbf{v}}_{i}$ of $P_{i}$ with $v_{t|i}$
and $v_{n|i}$.

Figure 3: Velocities of the vertices $P_{1},P_{2},\dots,P_{5}$ of a periodic
billiard in the ellipse $e$ with the caustic $c\,$.

In terms of the exterior angles $\theta_{1},\dots,\theta_{N}$ of the billiard,
we obtain from (3.1)

$\sin{\displaystyle\frac{\theta_{2}}{2}}={\displaystyle\frac{l_{2}\,\omega_{2}}{v_{2}}}={\displaystyle\frac{r_{2}\,\omega_{1}}{v_{2}}}\mskip
9.0mu\mbox{and}\mskip
9.0mu\cos{\displaystyle\frac{\theta_{2}}{2}}={\displaystyle\frac{v_{t|2}}{v_{2}}}\,,\mskip
9.0mu\mbox{where}\mskip 9.0muv_{2}:=\|{\mathbf{v}}_{2}\|.$ (3.3)

Let $R_{i}$ denote the pole of the line $[P_{i},P_{i+1}]$ with respect to
(w.r.t. in brief) $e\,$. Since the poles of a line $\ell$ w.r.t. confocal
conics lie on a line orthogonal to $\ell$, the side $P_{1}P_{2}$ is orthogonal
to $[Q_{1},R_{1}]$ (Figure 3), which means

$\overline{R_{1}Q_{1}}=l_{1}\tan\frac{\theta_{1}}{2}=r_{2}\tan\frac{\theta_{2}}{2}\,.$
(3.4)

From (3.4) and (3.3) follows

$l_{1}\tan\frac{\theta_{1}}{2}=l_{1}\frac{l_{1}\omega_{1}}{v_{t|1}}=r_{2}\tan\frac{\theta_{1}}{2}=r_{2}\frac{r_{2}\omega_{1}}{v_{t|2}}\mskip
9.0mu\mbox{and}\mskip
9.0mu\frac{v_{t|2}}{v_{t|1}}=\frac{r_{2}^{2}}{l_{1}^{2}}=\frac{\tan^{2}\frac{\theta_{1}}{2}}{\tan^{2}\frac{\theta_{2}}{2}}\,.$

This shows that, by virtue of (2.7), the products

$v_{t|1}\tan^{2}{\displaystyle\frac{\theta_{1}}{2}}=v_{t|2}\tan^{2}{\displaystyle\frac{\theta_{2}}{2}}=\dots=v_{t|N}\tan^{2}{\displaystyle\frac{\theta_{N}}{2}}=v_{t|i}{\displaystyle\frac{k_{e}}{\|{\mathbf{t}}_{c|i}\|^{2}}}$
(3.5)

for $i=1,2,\dots,N$ are equal along the billiard. We denote this quantity with
$C$. Instead of a free choice of $v_{2}\,$, it means no restriction of
generality to set $C=k_{e}\,$. Then we obtain for the point $P_{i}=(a_{e}\cos
t_{i},\,b_{e}\sin t_{i})$ of the ellipse $e$, by virtue of (2.7),

$\begin{array}[]{c}v_{t|i}=\|{\mathbf{t}}_{c}\|^{2}=-k_{h},\quad
v_{n|i}=v_{i}\sin\mbox{\small$\displaystyle\frac{\theta}{2}$}=\|{\mathbf{t}}_{c}\|\,\sqrt{k_{e}}=\sqrt{-k_{e}k_{h}}\\\\[5.69054pt]
v_{i}={\displaystyle\frac{\|{\mathbf{t}}_{c}\|^{2}}{\cos\frac{\theta}{2}}}=\|{\mathbf{t}}_{c}\|\,\|{\mathbf{t}}_{e}\|=\sqrt{k_{h}(k_{h}-k_{e})}\mskip
9.0mu\mbox{for}\mskip 9.0mut=t_{i}\mskip 9.0mu\mbox{and}\mskip
9.0muk_{h}=k_{h}(t_{i})\,.\end{array}$ (3.6)

##  4 Billiard motion and the underlying Lie group

Our specification of the quantity $C$ assigns to the vertex $P_{i}\in e$ with
parameter $t_{i}$ a non-vanishing velocity vector
${\mathbf{v}}_{i}=\|{\mathbf{t}}_{c}(t_{i})\|\,{\mathbf{t}}_{e}(t_{i})\,$.
This assignment can immediately be extended to all points of $e\,$. There
exists a parameter $u$ on $e$ such that the differentiation by $u$ gives the
said velocity vector. Let a dot indicate this differentiation. Then

${\mathbf{v}}(t)={\dot{\mathbf{p}}}(t)={\displaystyle\frac{\mathrm{d}\,{\mathbf{p}}(t)}{\mathrm{d}u}}=\|{\mathbf{t}}_{c}(t)\|\,{\mathbf{t}}_{e}(t)=\sqrt{a_{c}^{2}\sin^{2}t+b_{c}^{2}\cos^{2}t}\
{\mathbf{t}}_{e}(t).$ (4.1)

We can even extend this to all confocal ellipses of $c\,$. The assignment of a
velocity vector

${\mathbf{v}}(t)=\|{\mathbf{t}}_{c}(t)\|\,{\mathbf{t}}_{e}(t)=\dot{t}\;{\mathbf{t}}_{e}(t)$

to each point $P=(a_{e}\cos t,\,b_{e}\sin t)$ with
$a_{e}^{2}-b_{e}^{2}=a_{c}^{2}-b_{c}^{2}$ defines an ‘instant motion’ of the
plane, where

$\dot{t}=\frac{\mathrm{d}t}{\mathrm{d}u}=\|{\mathbf{t}}_{c}(t)\|=\sqrt{-k_{h}(t)}=\sqrt{a_{c}^{2}\sin^{2}t+b_{c}^{2}\cos^{2}t}\,.$
(4.2)

Figure 4: The infinitesimal motion assigns to each point of the Poncelet grid
a velocity vector such that the points on each tangent $t_{Q}$ to the caustic
$c$ remain aligned.

We prove below, that this ‘instant motion’ is compatible with the billiard and
the associated Poncelet grid. This means in particular, that the velocity
(3.6) is also valid for all points $S_{i}^{(j)}\in e^{(j)}$.

Figure 4 shows a portion of the Poncelet grid and the velocity vectors of a
couple of points, each represented by a scaled arrow. As indicated, for all
points on the tangent $t_{Q}$ to $c$ at any point $Q\in c$, the normal
components $\|{\mathbf{v}}_{n}\|$ of the respective velocity vectors
${\mathbf{v}}$ are proportional to the distance to $Q$. On the other hand, all
points on any confocal hyperbola share the tangential component
$\|{\mathbf{v}}_{t}\|$.

######  Theorem 1.

Let the billiard $P_{1}P_{2}\dots$ with the ellipse $c$ as caustic be moving
along the circumscribed ellipse $e\,$. Then the motion is the action of a one-
parameter Lie group $\Gamma$. Each transformation $\gamma(u)\in\Gamma$
preserves the confocal ellipses and permutes the confocal hyperbolas as well
as the tangents to $c\,$.

  1. 1.

If $(a_{c},b_{c})$ are the semiaxes of the caustic $c$ with the tangent
vectors ${\mathbf{t}}_{c}(t)=(-a_{c}\sin t,\,b_{c}\cos t)$, then for all
confocal ellipses $e$ with semiaxes $(a_{e},b_{e})$ the $\Gamma$ generating
instant motion is defined, up to a scalar, by the vector field

$(x,y)=(a_{e}\cos t,\,b_{e}\sin t)\,\
\mapsto\,\|{\mathbf{t}}_{c}\|\,{\mathbf{t}}_{e}=\sqrt{-k_{h}(t)}\,\left(-{\displaystyle\frac{a_{e}y}{b_{e}}},\,{\displaystyle\frac{b_{e}x}{a_{e}}}\right)$
(4.3)

with $a_{e}^{2}-b_{e}^{2}=a_{c}^{2}-b_{c}^{2}=d^{2}$ and
$a_{e}^{2}-a_{c}^{2}\geq 0$.

  2. 2.

If we parametrize the quadrant $x,y>0$ by elliptic coordinates as
${\mathbf{X}}(k_{e},k_{h})$, then the vector field can be expressed as

${\mathbf{X}}(k_{e},\,k_{h})\ \mapsto\
-2\sqrt{k_{h}(a_{c}^{2}+k_{h})(b_{c}^{2}+k_{h})}\;\frac{\partial{\mathbf{X}}}{\partial\,k_{h}}\,.$
(4.4)

The vector field defines a canonical parameter $u$ for the one-parameter Lie
group $\Gamma$, i.e., for transformations $\gamma(u)\in\Gamma$ holds
$\gamma(u_{2})\circ\gamma(u_{1})=\gamma(u_{1}+u_{2})$. At the same time, $u$
provides canonical coordinates222 Of course, we still obtain canonical
coordinates on the ellipses, when all velocity vectors are multiplied with any
constant $\lambda\in\mathbb{R}\setminus{0}$. on each confocal ellipse with the
property

$\gamma(2\varDelta u)\\!:\ P_{i}\mapsto P_{i+1},\ S_{i}^{(j)}\mapsto
S_{i+1}^{(j)}.$

###### Proof.

1\. The first derivative $\dot{t}$ in (4.2) is independent of the choice of
the ellipse $e$. Therefore $\gamma(u)$ permutes the confocal hyperboloids. On
the other hand, the representation
${\mathbf{v}}=\|{\mathbf{t}}_{c}\|\,{\mathbf{t}}_{e}$ reveals that all
confocal ellipses remain fixed. Furthermore, we verify that the position of
any point $P$ on the tangent $t_{Q}$ to $c$ at $Q$ (see Figure 4) is preserved
under the infinitesimal motion:

Given $P=\left(a_{e}\cos t,\,b_{e}\sin t\right)\in e$ and $Q=\left(a_{c}\cos
t^{\prime},\,b_{c}\sin t^{\prime}\right)$, the point $P$ lies on $t_{Q}$ if
and only if

$b_{c}a_{e}\cos t^{\prime}\cos t+a_{c}b_{e}\sin t^{\prime}\sin
t=a_{c}b_{c}\,.$ (4.5)

This is preserved under the infinitesimal motion if differentiation by $u$
based on (4.2) yields an identity, namely

$\begin{array}[]{l}\phantom{-}\|{\mathbf{t}}_{c}(t^{\prime})\|\left(-b_{c}a_{e}\sin
t^{\prime}\cos t+a_{c}b_{e}\cos t^{\prime}\sin
t\right)\phantom{.}\\\\[4.2679pt]
=-\|{\mathbf{t}}_{c}(t)\|\left(-b_{c}a_{e}\cos t^{\prime}\sin t+a_{c}b_{e}\sin
t^{\prime}\cos t\right).\end{array}$ (4.6)

In order to verify this, we square both sides and substitute from the squared
equation (4.5) the mixed term $2a_{c}b_{c}a_{e}b_{e}\sin t^{\prime}\cos
t^{\prime}\sin t\cos t\,$. After some computations, this yields for both sides

$d^{2}\left(\sin^{2}t-\sin^{2}t^{\prime}\right)\left(a_{c}^{2}b_{e}^{2}\sin^{2}t^{\prime}\sin^{2}t+b_{c}^{2}a_{e}^{2}\cos^{2}t^{\prime}\cos^{2}t-a_{c}^{2}b_{c}^{2}\right).$

The velocity analysis in (3.6) for the particular ellipse $e$ confirms, that
also the signs of both sides in (4.6) are equal.

2\. From (2.3) follows for $(x,y)={\mathbf{X}}(k_{e},\,k_{h})$

$2x\,\frac{\partial x}{\partial k_{h}}=\frac{a_{c}^{2}+k_{e}}{d^{2}}\,,\quad
2y\,\frac{\partial y}{\partial k_{h}}=-\frac{b_{c}^{2}+k_{e}}{d^{2}}$

and therefore

${\mathbf{X}}_{k_{h}}=\frac{\partial{\mathbf{X}}}{\partial\,k_{h}}=\frac{1}{2d^{2}}\left(\frac{a_{c}^{2}+k_{e}}{x},\,-\frac{b_{c}^{2}+k_{e}}{y}\right)=\frac{-1}{2d^{2}\sin
t\cos t}\,{\mathbf{t}}_{e}\,.$

This implies by (2.6) and (2.8)

$\|{\mathbf{t}}_{c}\|\,{\mathbf{t}}_{e}=\lambda{\mathbf{X}}_{k_{h}}\mskip
9.0mu\mbox{with}\mskip
9.0mu\lambda=-2a_{h}b_{h}\sqrt{-k_{h}}=-2\,\sqrt{k_{h}(a_{c}^{2}+k_{h})(b_{c}^{2}+k_{h})}\,,$

which confirms the claim in (4.4). ∎

In order to express the action of the transformation $\gamma(u)\in\Gamma$ on
an initial point $(a_{e}\cos t,\,b_{e}\sin t)$, we integrate the differential
equation (4.2)

$\dot{t}={\displaystyle\frac{\mathrm{d}t}{\mathrm{d}u}}=\sqrt{a_{c}^{2}\sin^{2}t+(a_{c}^{2}-d^{2})\cos^{2}t}=a_{c}\sqrt{1-m^{2}\cos^{2}t}$

with $m:=d/a_{c}<1$ as numerical eccentricity of the caustic $c$. The
substitution

$\varphi:=t-\mbox{\small$\displaystyle\frac{\pi}{2}$}$

results in

$\frac{\mathrm{d}\varphi}{\sqrt{1-m^{2}\sin^{2}\varphi}}=a_{c}\,\mathrm{d}u\,.$

The initial condition $\varphi=0$ for $u=0$ yields the unique solution

$a_{c}\,u(\varphi)=F(\varphi,m)=\int_{0}^{\varphi}\frac{\mathrm{d}\varphi}{\sqrt{1-m^{2}\sin^{2}\varphi}}$
(4.7)

with $F(\varphi,m)$ as the elliptic integral of the first kind with the
modulus $m$. The equation (4.7) shows the canonical coordinate $u$ in terms of
$\varphi$ with the quarter period

$K:=a_{c}\,u\left(\mbox{\small$\displaystyle\frac{\pi}{2}$}\right)=\int_{0}^{\pi/2}\frac{\mathrm{d}\varphi}{\sqrt{1-m^{2}\sin^{2}\varphi}}$

For the sake of simplicity, we introduce

$\tilde{u}(\varphi):=a_{c}\,u(\varphi)$ (4.8)

as a new canonical coordinate.

The inverse function of $\tilde{u}=F(\varphi,m)$, namely the Jacobian
amplitude $\varphi=\mathrm{am}\mskip 1.0mu(\tilde{u})$ leads to the Jacobian
elliptic functions, the elliptic sine

$\mskip 2.0mu\mathrm{sn}\mskip
2.0mu\tilde{u}=\sin(\mathrm{am}(\tilde{u}))=\sin\varphi=-\cos t$

with $\mskip 2.0mu\mathrm{sn}(-\tilde{u})=-\mskip 2.0mu\mathrm{sn}\mskip
2.0mu\tilde{u}\,$, the elliptic cosine

$\mskip 2.0mu\mathrm{cn}\mskip
2.0mu\tilde{u}=\cos(\mathrm{am}(\tilde{u}))=\cos\varphi=\sin t$

with $\mskip 2.0mu\mathrm{cn}(-\tilde{u})=\mskip 2.0mu\mathrm{cn}\mskip
2.0mu\tilde{u}\,$, and the delta amplitude

$\mskip 2.0mu\mathrm{dn}\mskip
2.0mu\tilde{u}=\sqrt{1-m^{2}\,\mathrm{sn}^{2}\tilde{u}}$

with $\mskip 2.0mu\mathrm{dn}(-\tilde{u})=\mskip 2.0mu\mathrm{dn}\mskip
2.0mu\tilde{u}\,$ as the third elliptic base function Hoppe . Moreover, for
$k\in\mathbb{Z}$ holds

$\begin{array}[]{c}\mskip 2.0mu\mathrm{sn}(\tilde{u}+2kK)=(-1)^{k}\mskip
2.0mu\mathrm{sn}\mskip 2.0mu\tilde{u}\,,\quad\mskip
2.0mu\mathrm{cn}(\tilde{u}+2kK)=(-1)^{k}\mskip 2.0mu\mathrm{cn}\mskip
2.0mu\tilde{u}\,,\\\\[2.84526pt] \mskip
2.0mu\mathrm{dn}(\tilde{u}+2kK)=(-1)^{k}\mskip 2.0mu\mathrm{dn}\mskip
2.0mu\tilde{u}\,.\end{array}$

This gives rise to the canonical parametrization of the ellipse $e$ with
semiaxes $(a_{e},b_{e})$ as

$\left(-a_{e}\mskip 2.0mu\mathrm{sn}\mskip 2.0mu\tilde{u},\,b_{e}\mskip
2.0mu\mathrm{cn}\mskip 2.0mu\tilde{u}\right)\quad\mbox{for}\mskip
9.0mu0\leq\tilde{u}<4K=4\tilde{u}\left(\mbox{\footnotesize$\displaystyle\frac{\pi}{2}$}\right).$

As an alternative, we can proceed with elliptic coordinates. From (4.4) and

$\frac{\mathrm{d}{\mathbf{X}}}{\mathrm{d}u}=\dot{k}_{e}\,\frac{\partial{\mathbf{X}}}{\partial\,k_{e}}+\dot{k}_{h}\,\frac{\partial{\mathbf{X}}}{\partial\,k_{h}}=-2\,\sqrt{k_{h}(a_{c}^{2}+k_{h})(b_{c}^{2}+k_{h})}\;\frac{\partial{\mathbf{X}}}{\partial\,k_{h}}$

follows for the orbits of the Lie group $\dot{k}_{e}=0$ and

$\dot{k}_{h}=-2\,\sqrt{k_{h}(a_{c}^{2}+k_{h})(b_{c}^{2}+k_{h})}\,.$

As expected, the orbits are confocal ellipses. In order to express the action
of $\gamma(u)\in\Gamma$ on an initial point ${\mathbf{X}}(k_{e|0},\,k_{h|0})$,
we have to integrate the differential equation

$\frac{\mathrm{d}k_{h}}{\sqrt{k_{h}(a_{c}^{2}+k_{h})(b_{c}^{2}+k_{h})}}=-2\,\mathrm{d}u$
(4.9)

with any initial condition. Again, we face an elliptic integral, this time in
the socalled Riemannian form.

######  Theorem 2.

  1. 1.

Let $c$ be the ellipse $c$ with semiaxes $(a_{c},b_{c})$ and linear
eccentricity $d=\sqrt{a_{c}^{2}-b_{c}^{2}}\,$. Then for all confocal ellipses
$e$ with semiaxes $(a_{e},b_{e})$, the inscribed billiards with the caustic
$c$ can be canonically parametrized using the Jacobian elliptic functions to
the modulus $m=d/a_{c}$ (= numerical eccentricity of $c$) as

$\left(-a_{e}\mskip 2.0mu\mathrm{sn}\mskip 2.0mu\tilde{u},\ b_{e}\mskip
2.0mu\mathrm{cn}\mskip 2.0mu\tilde{u}\right).$

This means that, if $b_{c}=b_{e}\mskip 2.0mu\mathrm{cn}(\varDelta\tilde{u})$,
then the vertices of the billiards in $e$ have the canonical parameters
$\tilde{u}=(\tilde{u}_{0}+2k\varDelta\tilde{u})$ for $k\in\mathbb{Z}$ and any
given initial $\tilde{u}_{0}\,$.

  2. 2.

Conversely, we obtain an ellipse $e$ for which the billiards with caustic $c$
are $N$-periodic with turning number $\tau$, where $\mathrm{gcd}(N,\tau)=1$,
by the choice

$\varDelta\tilde{u}=\frac{2\tau K}{N}$

with $K$ as the complete elliptic integral of the first kind to the modulus
$m$, provided that

$a_{e}=\frac{a_{c}\mskip 2.0mu\mathrm{dn}(\varDelta\tilde{u})}{\mskip
2.0mu\mathrm{cn}(\varDelta\tilde{u})}\mskip 9.0mu\mbox{and}\quad
b_{e}=\frac{b_{c}}{\mskip 2.0mu\mathrm{cn}(\varDelta\tilde{u})}\,.$ (4.10)

Figure 5: Dependence between the minor semiaxes $b_{e}$, $b_{c}$ and the
intervall $\varDelta\tilde{u}$.

###### Proof.

If the first vertex $P_{1}\in e$ of the billiard is chosen on the positive
$y$-axis, i.e., with canonical parameter $\tilde{u}=0$ (see Figure 5), then
the first contact point $Q_{1}$ has the parameter $\varDelta\tilde{u}$, and
the tangent to $c$ at $Q_{1}$ passes through $P_{1}=(0,b_{e})$. Hence, the
points $P_{1}$ and $Q_{1}$ are conjugate w.r.t. $c$, which means by (4.5) that
the product of the respective $y$-coordinates $b_{e}$ and $b_{c}\mskip
2.0mu\mathrm{cn}(\varDelta\tilde{u})$ equals $b_{c}^{2}\,$.

In view of the major semiaxis $a_{e}$ follows from
$\mathrm{dn}^{2}(\varDelta\tilde{u})=1-m^{2}\,\mathrm{sn}^{2}(\varDelta\tilde{u})$
and
$\mathrm{sn}^{2}(\varDelta\tilde{u})+\mathrm{cn}^{2}(\varDelta\tilde{u})=1$
that $\mskip 2.0mu\mathrm{dn}(\varDelta\tilde{u})=a_{e}\mskip
2.0mu\mathrm{cn}(\varDelta\tilde{u})/a_{c}\,$. ∎

######  Corollary 3.

If in the ellipse $e$ with semiaxes $(a_{e},b_{e})$ the billiard with caustic
$c$ is $N$-periodic with turning number $\tau=1$ and
$\varDelta\tilde{u}=\frac{2K}{N}\,$, then the associated Poncelet grid
contains the ellipses $e^{(1)},\,e^{(2)},\dots,\,e^{(k)}$,
$k=\left[\frac{N-3}{2}\right]$, with respective semiaxes

$\begin{array}[]{c}a_{e|1}={\displaystyle\frac{a_{c}\mskip
2.0mu\mathrm{dn}(2\mskip 1.0mu\varDelta\tilde{u})}{\mskip
2.0mu\mathrm{cn}(2\mskip 1.0mu\varDelta\tilde{u})}},\
b_{e|1}={\displaystyle\frac{b_{c}}{\mskip 2.0mu\mathrm{cn}(2\mskip
1.0mu\varDelta\tilde{u})}},\mskip 9.0mua_{e|2}={\displaystyle\frac{a_{c}\mskip
2.0mu\mathrm{dn}(3\mskip 1.0mu\varDelta\tilde{u})}{\mskip
2.0mu\mathrm{cn}(3\mskip 1.0mu\varDelta\tilde{u})}},\
b_{e|2}={\displaystyle\frac{b_{c}}{\mskip 2.0mu\mathrm{cn}(3\mskip
1.0mu\varDelta\tilde{u})}},\\\\[8.53581pt] \dots\ ,\
a_{e|k}={\displaystyle\frac{a_{c}\mskip 2.0mu\mathrm{dn}(k\mskip
1.0mu\varDelta\tilde{u})}{\mskip 2.0mu\mathrm{cn}(k\mskip
1.0mu\varDelta\tilde{u})}},\mskip
9.0mub_{e|k}={\displaystyle\frac{b_{c}}{\mskip 2.0mu\mathrm{cn}(k\mskip
1.0mu\varDelta\tilde{u})}}\,.\end{array}$

Figure 6: Canonical coordinates $u$ (blue) for the confocal hyperbolas and $v$
(red) for the confocal ellipses exterior to the caustic $c$ such that $u\pm
v=\mathrm{const.}$ represent the tangents of the caustic.

Corollary 3 reveals that $\varDelta\tilde{u}$ serves as a canonical coordinate
for confocal ellipses in the exterior of $c$. If $\varDelta\tilde{u}$
corresponds by (4.10) to the ellipse $e$ with semiaxes $(a_{e},b_{e})$, then
$2\mskip 1.0mu\varDelta\tilde{u}$ is the shift for the billiards in $e$ with
caustic $c\,$. If these billiards have the turning number 1, then increasing
the shift by $\varDelta\tilde{u}$ means to increase the turning number of the
billiard in a confocal ellipse by 1, while the caustic $c$ remains fixed
(Figure 6). The billiard $P_{1}P_{2}\dots$ and its conjugate
$P_{1}^{\prime}P_{2}^{\prime}\dots$ in $e$ (cf. (Sta_I, , Sect. 3.2))
intersect each other along the ellipse with the canonical coordinate
$\varDelta\tilde{u}/2\,$. Note that the ellipses $e^{(j)}$ and $e^{(N-2-j)}$
coincide while the corresponding $\varDelta\tilde{u}$’s differ in their signs.
For even $N$, the points $S_{i}^{(N/2-1)}$ are at infinity, and the line at
infinity as a limit of a confocal ellipse corresponds to
$\varDelta\tilde{u}=K$.

The following formulas express the elliptic coordinates $(k_{e},k_{h})$ of the
point $P=(-a_{e}\mskip 2.0mu\mathrm{sn}\mskip 2.0mu\tilde{u},\,b_{e}\mskip
2.0mu\mathrm{cn}\mskip 2.0mu\tilde{u})$ of $e$ in terms of the canonical
coordinate $\tilde{u}$ on $e$ and the shift $\varDelta\tilde{u}$ corresponding
to $e\,$.

$k_{e}=k_{e}(\varDelta\tilde{u})=\frac{a_{c}^{2}\,\mathrm{sn}^{2}\varDelta\tilde{u}}{\mathrm{cn}^{2}\varDelta\tilde{u}}(1-m^{2}),\quad
k_{h}=k_{h}(\tilde{u})=-a_{c}^{2}\,\mathrm{dn}^{2}\tilde{u}\,.$ (4.11)

This follows from

$k_{e}=a_{e}^{2}-a_{c}^{2}=a_{c}^{2}\,\frac{\mathrm{dn}^{2}\varDelta\tilde{u}-\mathrm{cn}^{2}\varDelta\tilde{u}}{\mathrm{cn}^{2}\varDelta\tilde{u}}=a_{c}^{2}(1-m^{2})\frac{\mathrm{sn}^{2}\varDelta\tilde{u}}{\mathrm{cn}^{2}\varDelta\tilde{u}}$

and

$k_{h}=-a_{c}^{2}\,\mathrm{cn}^{2}\mskip
2.0mu\tilde{u}-b_{c}^{2}\,\mathrm{sn}^{2}\mskip
2.0mu\tilde{u}=-a_{c}^{2}+d^{2}\,\mathrm{sn}^{2}\mskip
2.0mu\tilde{u}=-a_{c}^{2}+m^{2}a_{c}^{2}\,\mathrm{sn}^{2}\mskip
2.0mu\tilde{u}\,.$

Note that $k_{h}=k_{h}(\tilde{u})$ is a solution of (4.9).

In Izmestiev , an unordered pair of coordinates $(r,s)$ is proposed for each
point $P$ in the exterior of $c$, namely with $r$ and $s$ as canonical
coordinates of the tangency points for the tangent lines from $P$ to $c$ (see
also (MonGeom, , p. 358)). This means for $P=(-a_{e}\mskip
2.0mu\mathrm{sn}\mskip 2.0mu\tilde{u},\,b_{e}\mskip 2.0mu\mathrm{cn}\mskip
2.0mu\tilde{u})$ that

$r=\tilde{u}-\varDelta\tilde{u}\,,\quad s=\tilde{u}+\varDelta\tilde{u}\,,$

where $\varDelta\tilde{u}$ corresponds to $e$ according to (4.10).

If we keep the sum $\tilde{u}+\varDelta\tilde{u}$ or the difference
$\tilde{u}-\varDelta\tilde{u}$ constant, the corresponding point $P$ runs
along a tangent of the caustic $c$ (compare with (Bobenko, , Prop. 8.3)). For
the sake of simplicity, we replace in the summary below $\varDelta\tilde{u}$
by $\tilde{v}\,$.

Figure 7: The injective mapping ${\mathbf{Y}}$ sends the square grid of points
$Q_{i},P_{i}$ and $S_{i}^{(j)}$, $i=1,\dots,9\,$, $j=1,\dots,3\,$, to the
vertices and the diagonals to the confocal conics of the Poncelet grid
depicted in Figure 1.

######  Theorem 4.

Referring to the notation in Theorem 2, the injective mapping

$\begin{array}[]{rl}{\mathbf{Y}}\\!:&U\times
V\to{\mathbb{R}}^{2},\quad(\tilde{u},\,\tilde{v})\,\mapsto\,\left(-a_{c}{\displaystyle\frac{\mskip
2.0mu\mathrm{sn}\mskip 2.0mu\tilde{u}\,\mskip 2.0mu\mathrm{dn}\mskip
2.0mu\tilde{v}}{\mskip 2.0mu\mathrm{cn}\mskip 2.0mu\tilde{v}}},\
b_{c}{\displaystyle\frac{\mskip 2.0mu\mathrm{cn}\mskip 2.0mu\tilde{u}}{\mskip
2.0mu\mathrm{cn}\mskip 2.0mu\tilde{v}}}\right)\\\\[8.53581pt]
&\mathrm{for}\mskip 9.0muU:=\\{\tilde{u}\,|\,0\leq\tilde{u}<4K\\},\
V:=\\{\tilde{v}\,|\,0\leq\tilde{u}<K\\}\end{array}$

parametrizes the exterior of the caustic $c$ with semiaxes $(a_{c},b_{c})$ in
such a way, that the lines $\tilde{u}=\mathrm{const.}$ are branches of
confocal hyperbolas; $\tilde{v}=\mathrm{const.}$ are confocal ellipses and
$\tilde{u}\pm\tilde{v}=\mathrm{const.}$ tangents of $c\,$.

The domain of the mapping ${\mathbf{Y}}$ (Figure 7) can be extended to
${\mathbb{R}}^{2}$ and satisfies

${\mathbf{Y}}((\tilde{u}+4K),\,\tilde{v})={\mathbf{Y}}(\tilde{u},\,(\tilde{v}+2K))={\mathbf{Y}}(\tilde{u},\,-\tilde{v})={\mathbf{Y}}(\tilde{u},\,\tilde{v})$

and therefore
${\mathbf{Y}}(\tilde{u},\,(K+\tilde{v}))={\mathbf{Y}}(\tilde{u},\,(K-\tilde{v}))$
(Figure 7). The Liegroup $\Gamma$ mentioned in Theorem 1 is the
${\mathbf{Y}}$-transform of the group of translations along the
$\tilde{u}$-axis.

##  5 More about invariants of periodic billiards

In this section we study how the infinitesimal motion induced by the vector
field in (4.3) affects distances and angles of periodic billiards. As before,
the dot means differentiation by the canonical parameter $u$, and we call the
billiard motion canonical when it is parametrized by a canonical parameter
like $u\,$.

######  Lemma 5.

Let $P_{1}P_{2}\dots$ be a billiard in the ellipse $e$ with
$Q_{1},Q_{2},\dots$ as contact points with its caustic, the ellipse $c$. If
$t_{i}$ is the parameter of $P_{i}$ and $t_{i}^{\prime}$ that of $Q_{i}$, then

$r_{i}=\overline{Q_{i-1}P_{i}}={\displaystyle\frac{\|{\mathbf{t}}_{c}(t_{i-1}^{\prime})\|\,\|{\mathbf{t}}_{c}(t_{i})\|\sqrt{k_{e}}}{a_{c}b_{c}}},\quad
l_{i}=\overline{P_{i}Q_{i}}={\displaystyle\frac{\|{\mathbf{t}}_{c}(t_{i})\|\,\|{\mathbf{t}}_{c}(t_{i}^{\prime})\|\sqrt{k_{e}}}{a_{c}b_{c}}}.$

The canonical motion of the billiard induces for the side $P_{i}P_{i+1}$ the
instant angular velocity

$\omega_{i}=\frac{a_{c}b_{c}}{\|{\mathbf{t}}_{c}(t_{i}^{\prime})\|}=\,\frac{a_{c}b_{c}}{\sqrt{-k_{h}(t_{i}^{\prime})}}\,.$

Figure 8: Velocities ${\mathbf{v}}_{2}$ of the vertex $P_{2}$,
${\mathbf{v}}_{1}^{\prime}$ of the affine image $Q_{1}^{\prime}$ with
$\|{\mathbf{v}}_{1}^{\prime}\|=\|{\mathbf{v}}_{n_{1}}\|=\|{\mathbf{v}}_{n_{2}}\|$,
of the contact points $Q_{i}$ (with $Q_{i}^{*}$ as respective centers of
curvature), and of vertices $P_{i}^{\prime}$ of the conjugate billiard,
$i=1,2\,$.

###### Proof.

Referring to Figure 8, if the tangent $[Q_{1},P_{2}]$ rolls on $c$, then the
vertex $P_{2}$ receives the velocity vector ${\mathbf{v}}_{n_{2}}$ satisfying
(3.6), while the point of contact $Q_{1}$ moves with the velocity
$v_{c}(t_{1}^{\prime})$ along $c\,$. We can express this velocity in terms of
the radius of curvature $\rho_{c}(t_{1}^{\prime})$ of $c$ as

$v_{c}(t_{1}^{\prime})=\omega_{1}\,\rho_{c}(t_{1}^{\prime}),$

where $\rho_{c}(t)=\|{\mathbf{t}}_{e}(t)\|^{3}/a_{c}b_{c}$ by (Conics, , p.
79). On the other hand, from ${\mathbf{v}}_{c}=\dot{t}\,{\mathbf{t}}_{c}$ and
(4.2) follows
$v_{c}(t_{1}^{\prime})=\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|^{2}$. This yields

$\omega_{1}=\frac{v_{c}(t_{1}^{\prime})}{\rho_{c}(t_{1}^{\prime})}=\frac{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|^{2}\,a_{c}b_{c}}{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|^{3}}=\frac{a_{c}b_{c}}{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|}=\frac{a_{c}b_{c}}{\sqrt{-k_{h}(t_{1}^{\prime})}}$

by (2.6). Thus, we obtain for the velocity $v_{n|2}$ of $P_{2}$ by (3.6)

$r_{2}\,\omega_{1}=v_{n|2}=\|{\mathbf{t}}_{c}(t_{2})\|\sqrt{k_{e}},\mskip
9.0mu\mbox{hence}\mskip
9.0mur_{2}=\frac{\|{\mathbf{t}}_{c}(t_{2})\|\sqrt{k_{e}}}{a_{c}b_{c}}\,\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|.$

Similarly follows from $l_{2}\,\omega_{2}=v_{n|2}$ the stated expression for
the distance $l_{2}\,$. Note that
$t_{1},t_{1}^{\prime},t_{2},t_{2}^{\prime},t_{3},\dots$ is the sequence of
consecutive parameters of the points $P_{1},Q_{1},P_{2},Q_{2},P_{3},\dots$ The
formulas for $r_{i}$ and $l_{i}$ reveal again that the same distances appear
at the conjugate billiard. ∎

The angular velocity of the tangent to $e$ at $P_{2}$ equals the arithmetic
mean $(\omega_{1}+\omega_{2})/2$ (Figure 8). On the other hand, it is defined
by the radius of curvature $\rho_{e}$ of $e$ at $P_{2}$ and the velocity
$v_{2}$ by (3.6), since

$v_{2}=\rho_{e}(t_{2})\frac{\omega_{1}+\omega_{2}}{2}\,.$

This means by Lemma 5

$\|{\mathbf{t}}_{c}(t_{2})\|\,\|{\mathbf{t}}_{e}(t_{2})\|={\displaystyle\frac{\|{\mathbf{t}}_{e}(t_{2})\|^{3}}{a_{e}b_{e}}}\,{\displaystyle\frac{a_{c}b_{c}}{2}}\left({\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|}}+{\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{2}^{\prime})\|}}\right)$

and results by (2.7) in

$\frac{1}{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|}+\frac{1}{\|{\mathbf{t}}_{c}(t_{2}^{\prime})\|}=\frac{2a_{e}b_{e}}{a_{c}b_{c}}\,\frac{\|{\mathbf{t}}_{c}(t_{2})\|}{\|{\mathbf{t}}_{e}(t_{2})\|^{2}}=\frac{a_{e}b_{e}}{a_{c}b_{c}\sqrt{k_{e}}}\,\sin\theta_{2}\,.$
(5.1)

######  Theorem 6.

The exterior angles $\theta_{i}$ of an $N$-periodic billiard in an ellipse and
with an ellipse as caustic satisfy for $N\equiv 0\pmod{2}$

$\sum_{i=1}^{N}(-1)^{i}\sin\theta_{i}=0\mskip 9.0mu\mbox{and for}\mskip
9.0muN\equiv 0\mskip-15.0mu\pmod{4}\mskip
9.0mu\sum_{i=1}^{N/2}(-1)^{i}\sin\theta_{i}=0\,.$

###### Proof.

By virtue of (Sta_I, , Corollary 4.2), periodic billiards with even $N=2n$ are
centrally symmetric, which implies $\theta_{i}=\theta_{i+n}\,$. For $N\equiv
2\pmod{4}$ the sum from $1$ to $N$ must vanish since $(-1)^{i}=-(-1)^{(i+n)}$.
In the remaining case $N\equiv 0\pmod{4}$ follows from (5.1)

$\sin\theta_{i}=\frac{a_{c}b_{c}\sqrt{k_{e}}}{a_{e}b_{e}}\left(\frac{1}{\|{\mathbf{t}}_{c}(t_{i-1}^{\prime})\|}+\frac{1}{\|{\mathbf{t}}_{c}(t_{i}^{\prime})\|}\right)$

and further

$\begin{array}[]{l}\displaystyle\sum_{i=1}^{N/2}\,\sin\theta_{i}={\displaystyle\frac{a_{c}b_{c}\sqrt{k_{e}}}{a_{e}b_{e}}}\,\cdot\\\
\phantom{m}\left({\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{N}^{\prime})\|}}+{\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|}}-{\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{1}^{\prime})\|}}-{\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{2}^{\prime})\|}}+-\dots-{\displaystyle\frac{1}{\|{\mathbf{t}}_{c}(t_{n}^{\prime})\|}}\right).\end{array}$

This sum vanishes, since
${\mathbf{t}}_{c}(t_{n}^{\prime})=-{\mathbf{t}}_{c}(t_{N}^{\prime})\,$, due to
the odd turning number $\tau$ because of $\mbox{gcd}(N,\tau)=1\,$. ∎

At the same token, from

$\dot{\theta}_{i}=\omega_{i}-\omega_{i-1}=a_{c}b_{c}\left(\frac{1}{\|{\mathbf{t}}_{c}(t_{i}^{\prime})\|}-\frac{1}{\|{\mathbf{t}}_{c}(t_{i-1}^{\prime})\|}\right)$
(5.2)

and (5.1) follows

$\frac{\mathrm{d}}{\mathrm{d}u}\cos\theta_{i}=-\dot{\theta}_{i}\sin\theta_{i}=\frac{a_{c}^{2}b_{c}^{2}\sqrt{k_{e}}}{a_{e}b_{e}}\left(\frac{1}{\|{\mathbf{t}}_{c}(t_{i-1}^{\prime})\|^{2}}-\frac{1}{\|{\mathbf{t}}_{c}(t_{i}^{\prime})\|^{2}}\right).$

This shows that
$\frac{\mathrm{d}}{\mathrm{d}u}\left(\sum_{1}^{N}\cos\theta_{i}\right)$
vanishes and, therefore, $\sum_{1}^{N}\cos\theta_{i}$ is invariant against
billiard motions, which was first proved in Ako-Tab .

For the sake of completeness, we also focus on the variation of the side
lengths under the canonical billiard motion. We obtain

$\begin{array}[]{rcl}{\displaystyle\frac{\mathrm{d}}{\mathrm{d}u}}\,\overline{P_{i}P_{i+1}}&=&v_{t|i+1}-v_{t|i}\\\\[4.2679pt]
&=&\|{\mathbf{t}}_{c}(t_{i+1})\|^{2}-\|{\mathbf{t}}_{c}(t_{i})\|^{2}=d^{2}(\sin^{2}t_{i+1}-\sin^{2}t_{i}).\end{array}$
(5.3)

The vanishing sum over all $i$ confirms again the constant perimeter. We like
to recall that already in Bialy-Tab  some proofs for invariants were based on
differentiation.

Finally we concentrate on the effects showing up when the vertex $P_{i}$
traverses a quarter of the full period along the billiard.

######  Lemma 7.

As before, let
$t_{1},t_{1}^{\prime},t_{2},t_{2}^{\prime},\dots,t_{N}^{\prime}$ be the
sequence of parameters of an $N$-periodic billiard in an ellipse $e$ with an
ellipse $c$ as caustic. Then holds for $N\equiv 0\pmod{2}$:

$\left.\begin{array}[]{rl}\mbox{if}\mskip
9.0muN=4n:&\|{\mathbf{t}}_{c}(t_{i})\|\,\|{\mathbf{t}}_{c}(t_{i+n})\|=\sqrt{k_{h}(t_{i})\,k_{h}(t_{i+n})}\\\\[2.84526pt]
\mbox{if}\mskip
9.0muN=4n\\!+\\!2:&\|{\mathbf{t}}_{c}(t_{i})\|\,\|{\mathbf{t}}_{c}(t_{i+n}^{\prime})\|=\sqrt{k_{h}(t_{i})\,k_{h}(t_{i+n}^{\prime})}\end{array}\right\\}=a_{c}\,b_{c}\,,$

and the same after the parameter shift $t_{i}\mapsto t_{i}^{\prime}$ and
$t_{i}^{\prime}\mapsto t_{i+1}\,$.

###### Proof.

Based on the canonical parametrization by $\tilde{u}$, a quarter of the period
$4K$ corresponds to a shift by $K$. In the case $N=4n$ this shift effects
$t_{i}\mapsto t_{i+n}$ and $t_{i}^{\prime}\mapsto t_{i+n}^{\prime}$. If
$N=4n+2$, then $t_{i}\mapsto t_{i+n}^{\prime}$ and $t_{i}^{\prime}\mapsto
t_{i+n+1}$.
According to (4.11) holds $k_{h}=-a_{c}^{2}\,\mathrm{dn}^{2}\tilde{u}$ and by
(2.6) $\|{\mathbf{t}}_{c}(t)\|=\sqrt{-k_{h}(t)}=a_{c}\mskip
2.0mu\mathrm{dn}\mskip 2.0mu\tilde{u}\,$. The identity

$\mskip 2.0mu\mathrm{dn}(\tilde{u}+K)=\frac{\sqrt{1-m^{2}}}{\mskip
2.0mu\mathrm{dn}(\tilde{u})}$

implies

$\mskip 2.0mu\mathrm{dn}(\tilde{u})\cdot\mskip
2.0mu\mathrm{dn}(\tilde{u}+K)=\frac{b_{c}}{a_{c}}\,,\mskip
9.0mu\mbox{hence}\mskip 9.0mu\sqrt{k_{h}(\tilde{u})\cdot
k_{h}(\tilde{u}+K)}=a_{c}\,b_{c}\,.$

This confirms the claim. ∎

######  Theorem 8.

If an $N$-periodic billiard in an ellipse $e$ with an ellipse as caustic is
given with $N\equiv 0\pmod{2}$, then the distances
$r_{i}=\overline{Q_{i-1}P_{i}}$ and $l_{i}=\overline{P_{i}Q_{i}}$ satisfy

$\left.\begin{array}[]{rl}\mbox{for}\mskip 9.0muN=4n:&r_{i}\cdot
r_{i+n}=l_{i}\cdot l_{i+n}\\\\[2.84526pt] \mbox{for}\mskip
9.0muN=4n\\!+\\!2:&r_{i}\cdot l_{i+n}=l_{i}\cdot
r_{i+n+1}\end{array}\right\\}=k_{e}\,.$

###### Proof.

From the expressions for $r_{i}$ and $l_{i}$ in Lemma 5 follows, by virtue of
Lemma 7, for $N=4n$

$r_{i}\cdot
r_{i+n}=\frac{k_{e}}{a_{c}^{2}b_{c}^{2}}\,\|{\mathbf{t}}_{c}(t_{i-1}^{\prime})\|\,\|{\mathbf{t}}_{c}(t_{i})\|\,\|{\mathbf{t}}_{c}(t_{i+n-1}^{\prime})\|\,\|{\mathbf{t}}_{c}(t_{i+n})\|=k_{e}$

and the same result for $l_{i}\cdot l_{i+n}\,$. In the case $N=4n+2$ we obtain
similarly

$r_{i}\cdot l_{i+n}=l_{i}\cdot r_{i+n+1}=k_{e}\,,$

as stated. ∎

The following corollary is an immediate consequence of Theorem 8.

######  Corollary 9.

Let $s_{i}=\overline{P_{i}P_{i+1}}=l_{i}+r_{i+1}$ for $i=1,\dots,N$ be the
side lengths of the $N$-periodic billiard with even $N$ and
$s_{i}^{\prime}=\overline{P_{i}^{\prime}P_{i+1}^{\prime}}=r_{i+1}+l_{i+1}$
that of the conjugate billiard. Then,

$\begin{array}[]{rl}\mbox{for}\mskip
9.0muN=4n:&{\displaystyle\frac{s_{i+n}}{s_{i}}}={\displaystyle\frac{l_{i+n}}{r_{i+1}}}={\displaystyle\frac{r_{i+n+1}}{l_{i}}}\,,\\\\[7.11317pt]
\mbox{for}\mskip
9.0muN=4n\\!+\\!2:&{\displaystyle\frac{s_{i+n}}{s_{i-1}^{\prime}}}={\displaystyle\frac{l_{i+n}}{l_{i}}}={\displaystyle\frac{r_{i+n+1}}{r_{i}}}\,.\end{array}$

Finally we prove the invariance of k117 in (80, , Table 2).

######  Theorem 10.

Referring to the notation in Lemma 5, for even $N$ the products

$r_{1}r_{2}\dots r_{N}=l_{1}l_{2}\dots l_{N}=k_{e}^{N/2}$

are invariant against billiard motions. For $N\equiv 0\pmod{4}$ this is
already true for the products

$r_{1}r_{2}\dots r_{N/2}=l_{1}l_{2}\dots l_{N/2}=k_{e}^{N/4}.$

###### Proof.

For $N\equiv 0\pmod{4}$ the statements are a direct consequence of Theorem 8
and the central symmetry of the billiard which exchanges $r_{i}$ with
$r_{i+N/2}$ and $l_{i}$ with $l_{i+N/2}\,$.
In the remaining case $N=2n+2$ we note that by (3.2) $R(u):=r_{1}r_{2}\dots
r_{N}=l_{1}l_{2}\dots l_{N}$. Hence, by virtue of Theorem 8,

$R^{2}(u)=\prod_{i=1}^{N}\ (r_{i}\,l_{i+n})=k_{e}^{N},$

which yields the stated result. ∎

## References

  * (1) A. Akopyan, R. Schwartz, S. Tabachnikov, Billiards in ellipses revisited, Eur. J. Math. 2020, doi:10.1007/s40879-020-00426-9
  * (2) M. Bialy, S. Tabachnikov, Dan Reznik’s identities and more. Eur. J. Math. 2020, doi:10.1007/s40879-020-00428-7
  * (3) A.I. Bobenko, W.K. Schief, Y.B. Suris, J. Techter, On a discretization of confocal quadrics. A geometric approach to general parametrizations, Int. Math. Res. Not. IMRN 2020, no. 24, 10180–10230
  * (4) A. Chavez-Caliz, More about areas and centers of Poncelet polygons, Arnold Math. J. 7, 91–106 (2021).
  * (5) V. Dragović, M. Radnović, Integrable billiards and quadrics. Russian Math. Surveys 65/2, 319–379 (2010)
  * (6) V. Dragović, M. Radnović, Poncelet porisms and beyond. Integrable Billiards, Hyperelliptic Jacobians and Pencils of Quadrics. Birkhäuser/Springer Basel AG, Basel (2011)
  * (7) J.J. Duistermaat, Discrete Integrable Systems. QRT Maps and Elliptic Surfaces. Springer Monographs in Mathematics 304, Springer, New York (2010)
  * (8) G. Glaeser, H. Stachel, B. Odehnal, The Universe of Conics, From the ancient Greeks to 21st century developments, Springer Spectrum, Berlin, Heidelberg (2016)
  * (9) R. Hoppe, Elliptische Integrale und Funktionen nach Jacobi. http://www.dfcgen.de/wpapers/elliptic.pdf (2004–2015), accessed May 2021
  * (10) I. Izmestiev, S. Tabachnikov, Ivory’s Theorem revisited, J. Integrable Syst. 2/1, xyx006 (2017) R. Hoppe: Elliptische Integrale und Funktionen nach JACOBI M. Koecher, A. Krieg, Elliptische Funktionen und Modulformen, 2\. Aufl., Springer-Verlag, Berlin, Heidelberg (2007)
  * (11) D. Reznik, R. Garcia, J. Koiller, Eighty New Invariants of N-Periodics in the Elliptic Billiard, arXiv:2004.12497v11[math.DS] (2020)
  * (12) H. Stachel, Recalling Ivory’s Theorem, FME Transactions 47, No 2, 355–359 (2019)
  * (13) H. Stachel, The Geometry of Billiards in Ellipses and their Poncelet Grids. arXiv:2105.03362[math.MG] (2021).
  * (14) S. Tabachnikov, Geometry and Billiards, American Mathematical Society, Providence/Rhode Island (2005)


# Charge State-Dependent Symmetry Breaking of Atomic Defects in Transition
Metal Dichalcogenides

Feifei Xiang † 0000-0002-4930-5919 nanotech@surfaces Laboratory, Empa – Swiss
Federal Laboratories for Materials Science and Technology, Dübendorf 8600,
Switzerland  Lysander Huberich† nanotech@surfaces Laboratory, Empa – Swiss
Federal Laboratories for Materials Science and Technology, Dübendorf 8600,
Switzerland  Preston A. Vargas† Department of Materials Science and
Engineering, University of Florida, Gainesville, FL, 32611, USA  Riccardo
Torsi  Department of Materials Science and Engineering, The Pennsylvania State
University, University Park, PA 16082, USA  Jonas Allerbeck
0000-0002-3912-3265 nanotech@surfaces Laboratory, Empa – Swiss Federal
Laboratories for Materials Science and Technology, Dübendorf 8600, Switzerland
Anne Marie Z. Tan  Department of Materials Science and Engineering, University
of Florida, Gainesville, FL, 32611, USA  Institute of High Performance
Computing (IHPC), Agency for Science, Technology and Research (A*STAR),
Singapore 138632, Republic of Singapore  Chengye Dong  Two-Dimensional Crystal
Consortium, The Pennsylvania State University, University Park, PA 16802, USA
Pascal Ruffieux  0000-0001-5729-5354 nanotech@surfaces Laboratory, Empa –
Swiss Federal Laboratories for Materials Science and Technology, Dübendorf
8600, Switzerland  Roman Fasel  0000-0002-1553-6487 nanotech@surfaces
Laboratory, Empa – Swiss Federal Laboratories for Materials Science and
Technology, Dübendorf 8600, Switzerland  Oliver Gröning  nanotech@surfaces
Laboratory, Empa – Swiss Federal Laboratories for Materials Science and
Technology, Dübendorf 8600, Switzerland  Yu-Chuan Lin  Department of Materials
Science and Engineering, The Pennsylvania State University, University Park,
PA 16082, USA  Department of Materials Science and Engineering, National Yang
Ming Chiao Tung University, Hsinchu City 300, Taiwan  Richard G. Hennig
0000-0003-4933-7686 Department of Materials Science and Engineering,
University of Florida, Gainesville, FL, 32611, USA  Joshua A. Robinson
0000-0002-1513-7187 Department of Materials Science and Engineering, The
Pennsylvania State University, University Park, PA 16082, USA  Two-Dimensional
Crystal Consortium, The Pennsylvania State University, University Park, PA
16802, USA  Department of Chemistry and Department of Physics, The
Pennsylvania State University, University Park, PA, 16802, USA  Bruno Schuler
0000-0002-9641-0340<EMAIL_ADDRESS>nanotech@surfaces Laboratory, Empa –
Swiss Federal Laboratories for Materials Science and Technology, Dübendorf
8600, Switzerland

###### Abstract

The functionality of atomic quantum emitters is intrinsically linked to their
host lattice coordination. Structural distortions that spontaneously break the
lattice symmetry strongly impact their optical emission properties and spin-
photon interface. Here we report on the direct imaging of charge state-
dependent symmetry breaking of two prototypical atomic quantum emitters in
mono- and bilayer MoS2 by scanning tunneling microscopy (STM) and non-contact
atomic force microscopy (nc-AFM). By substrate chemical gating different
charge states of sulfur vacancies (Vac${}_{\text{S}}$) and substitutional
rhenium dopants (Re${}_{\text{Mo}}$) can be stabilized.
Vac${}_{\text{S}}^{-1}$ as well as Re${}_{\text{Mo}}^{0}$ and
Re${}_{\text{Mo}}^{-1}$ exhibit local lattice distortions and symmetry-broken
defect orbitals attributed to a Jahn-Teller effect (JTE) and pseudo-JTE,
respectively. By mapping the electronic and geometric structure of single
point defects, we disentangle the effects of spatial averaging, charge
multistability, configurational dynamics, and external perturbations that
often mask the presence of local symmetry breaking.

Mastering the controlled introduction of defects and impurities in
semiconductors has proven pivotal to their transformative technological
success. As miniaturization reaches atomic length scales and quantum
applications emerge, stochastic distributions make way for single impurity
engineering 1; 2. This need for precise defect manipulation extends to the
realm of artificial atom qubits, a cornerstone of quantum sensing and quantum
communication 3; 4; 3. Atomic quantum emitters based on defects in solids,
rely on spin-selective optical decay pathways enabling the high-fidelity spin
initialization and readout even at room temperature 5. Two-dimensional (2D)
materials such as hexagonal boron nitride or transition metal dichalcogenides
(TMDs) emerged as a game-changing platform to host such quantum emitters
because they lack surface dangling bonds that degrade spin coherence or photon
distinguishability 6; 7. Furthermore, 2D materials exhibit increased
extraction efficiency 8, integrate seamlessly with quantum optical devices 9;
10, and support electrical fine-tuning 11; 12 and spatial engineering 10; 13
of defects and their associated electronic states. Regardless of the host
material, the local defect symmetry governs the level structure and optical
selection rules that ultimately provide the spin-photon interface. Spontaneous
symmetry breaking from Jahn-Teller distortions and strain fields can therefore
dramatically change the radiative recombination rates, emission wavelength and
polarization 14; 15. Jahn-Teller systems, in particular, can strongly
influence spin-orbit coupling that gives rise to the inter-system crossing,
which yields the spin contrast 16; 17; 15. Therefore, understanding the
mechanisms that break the coordination symmetry of single atom qubits as a
function of their charge state and external factors such as strain are
decisive to tailor their functionality.

Although it is well known that TMDs and other transition metal compounds
exhibit enhanced susceptibility towards Jahn-Teller effect (JTE), detecting
JTE experimentally is challenging, because configurational dynamics and
spatial averaging can obscure the minute local distortions 18. Additionally,
these distortions are often affected by external perturbation, which can
complicate the interpretation of experimental results 18. Therefore, even
though the JTE and its extensions 19 has been discovered almost a century ago,
its abundance and significance in different materials systems has only
gradually been revealed and is still mostly inferred theoretically by ab
initio methods 20; 18. While a number of TMD point defects including chalcogen
vacancies 21; 22 and substitutional transition metal dopants 23 have been
predicted to undergo spontaneous symmetry breaking, it has only rarely been
verified experimentally using aberration-corrected transmission electron
microscopy 24, and no direct experimental observation of symmetry-broken
electronic states has been reported to date.

Here we present direct experimental and theoretical evidence for a Jahn-Teller
driven symmetry breaking of chalcogen vacancies 25 and rhenium-based quantum
emitters 26 in MoS2. Scanning tunneling microscopy (STM) orbital imaging 27,
CO-tip noncontact atomic force microscopy (nc-AFM) measurements 28, and
density functional theory (DFT) reveal the symmetry-broken defect orbitals and
the structural distortion for negatively charged sulfur vacancy
(Vac${}_{\text{S}}^{-1}$) and rhenium substituting for molybdenum in its
neutral and negative charge state (Re${}_{\text{Mo}}^{0}$,
Re${}_{\text{Mo}}^{-1}$). We reveal that Vac${}_{\text{S}}^{-1}$ coexist in
both a symmetric and distorted geometry, which may explain why the most
studied TMD defect has evaded the spectroscopic detection of spontaneous
symmetry breaking so far. Charge state tristability of Re dopants is achieved
by chemical gating via the underlying substrate. We find larger domains where
the structural distortions are aligned, which may indicate a non-negligible
strain field in the epitaxially-grown MoS2 samples. Surprisingly, the
comparatively broad defect resonances of some Re impurities exhibit a
continuous transformation of the measured defect orbital, suggesting a
configurational continuum as a consequence of the flat potential energy
surface experienced by the Re impurity.

## Jahn-Teller Driven Symmetry-Breaking of Vac${}_{\text{S}}^{-1}$ in MoS2

Chalcogen vacancies are the most prevalently discussed point defect in TMDs
due to their low formation energy 21, frequent detection in transmission
electron microscopy 29; 30, and attributed sub-band gap emission in optical
spectroscopy 31; 32. Despite their disputed presence at ambient conditions due
to their high reactivity 33, recent strategies involving deliberate vacancy
generation by annealing or ion bombardment and protection by vacuum or inert
capping layers made it possible to directly measure unpassivated vacancies in
TMDs 34; 25; 32.

Fig. 1a shows a STM image of monolayer MoS2 grown by metal-organic chemical
vapor deposition (MOCVD) on quasi-freestanding epitaxial graphene on
6H-silicon carbide (0001) (QFEG). S vacancies (Vac${}_{\text{S}}$) were
induced by annealing the sample up to 900 °C in ultrahigh-vacuum 34. Apart
from common (unintentional) as-grown defects like O${}_{\text{S}}$,
CH${}_{\text{S}}$, Cr${}_{\text{Mo}}$, and W${}_{\text{Mo}}$ 35 (see Fig. S2)
that are also present at much lower annealing temperatures, we find the
anticipated Vac${}_{\text{S}}$ in both the upper and lower S layer as
previously identified 34; 36. All S vacancies are negatively charged because
the QFEG substrate sets the Fermi level relatively high up in the MoS2 band
gap, such that the lowest Vac${}_{\text{S}}$ in-gap state becomes occupied.
The negative charge can be inferred from the band bending-induced dark halo
around the defect at positive sample bias 37.

While some Vac${}_{\text{S}}^{-1}$ appear threefold symmetric, others
unexpectedly exhibit a ‘distorted’ appearance with only a two-fold symmetry at
negative sample bias (Fig. 1a,b), even though both type of defects appear in
nc-AFM as sulfur vacancies in the top sulfur layer (Fig. 1c). Also in scanning
tunneling spectroscopy (STS) measurements the electronic structure between the
symmetric and distorted vacancy species are distinct as shown in Fig. 1d. The
symmetric Vac${}_{\text{S}}^{-1}$ has two occupied in-gap defect states
(HOMO-1 and HOMO) and one LUMO resonance close to the Fermi energy (blue),
while the distorted Vac${}_{\text{S}}^{-1}$ features one HOMO resonance (with
vibronic side-bands) and two LUMO resonances (red and orange). The symmetric
Vac${}_{\text{S}}^{-1}$ preserves the underlying C3v lattice symmetry in its
frontier orbitals (see Fig. S4e,f), whereas the HOMO of the distorted
Vac${}_{\text{S}}^{-1}$ variant has a lower C2v symmetry (see Fig. 1e). Such a
behavior was not observed for Vac${}_{\text{S}}^{0}$ in WS2 34, which
indicates the decisive role of the defect charge state.

We attribute the reduced point group symmetry of Vac${}_{\text{S}}^{-1}$ to a
JTE as previously predicted 22. Upon charging, the degenerated S vacancy state
(Fig. 1f) becomes unstable and relaxes, thereby lifting the orbital degeneracy
(Fig. 1g). The distorted Vac${}_{\text{S}}^{-1}$ geometry (Fig. S5c-e) has a
lower total energy by 116 meV 22. The calculated density of states (DOS) map
of the frontier orbitals shown in Fig. 1i are in excellent agreement with the
experimental differential conductance (d$I$/d$V$) maps. However, the reason
for the co-existence of both the symmetric and distorted
Vac${}_{\text{S}}^{-}1$ is not entirely clear. While a doubly negatively
charged Vac${}_{\text{S}}^{-2}$ would have a symmetric ground state 22, it is
unlikely to be formed due to the large Coulomb energy associated with the
highly confined defect wave function, in-line with the calculated charge
transition level from -1 to -2, which is above the conduction band minimum 22.
Possibly, the system is trapped in a local minimum with a potential barrier
too high to relax to the global potential energy minimum.

Motivated by the surprising discovery of symmetry breaking in one of the best
studied defects in TMDs, we designed a defect system that we can stabilize in
different charge states by chemical gating via the substrate. Next we will
discuss the charge state tristability of rhenium dopants in MoS2 and their
charge state-dependent symmetry breaking. Rhenium substituted for molybdenum
(Re${}_{\text{Mo}}$) has a D3h point symmetry group and produces fundamentally
different in-gap states as compared to Vac${}_{\text{S}}$, which makes it an
ideal candidate to test the generality of our findings.

## Charge State Tristability of Re Dopants in MoS2

Chemical doping of 2D semiconductors has been a popular but challenging route
to control their conductivity due to the large defect ionization energies
associated with the tightly confined defect wave functions in 2D 38. Recently,
a scalable MOCVD process has been developed to introduce Re dopants in TMDs in
controllable concentrations from hundreds of ppm to percentage concentrations
39; 40. Substitutional Re dopants in WSe2 on QFEG are positively charged
(ionized) and feature a series of closely-spaced unoccupied defect states
below the conduction band minimum. Here we leverage the higher electron
affinity $\chi$ of monolayer MoS2 as compared to WSe2, to position the Fermi
level in-between these Re${}_{\text{Mo}}$ defect states. In MoS2/QFEG, we find
that most Re impurities are charge neutral while some remain positively
charged, which also applies to Re-doped bilayer MoS2 grown on QFEG. By using
epitaxial graphene (EG) without hydrogen intercalation as a substrate, which
has a substantially smaller work function than QFEG, we push the Fermi level
even higher41; 42. Counter-intuitive for an $n$-type dopant, this results in
the Re impurities becoming negatively charged, i.e. they accept an additional
electron. By substrate chemical gating, we establish a charge state
tristability of Re dopants in MoS2.

In Fig. 2b we present a large-scale STM topography of Re-doped (5%) MoS2 on
partially H-intercalated EG. In areas with H-intercalation (gray plateau in
Fig. 2b and structural model shown in Fig. 2c), we observe the Re dopants in
the positive (Re${}_{\text{Mo}}^{+1}$) or neutral (Re${}_{\text{Mo}}^{0}$)
charge state, whereas on EG substrates (blue canyon in Fig. 2b, and structural
model in Fig. 2c) all Re are negatively charged (Re${}_{\text{Mo}}^{-1}$). The
charge states can be discriminated by the characteristic upwards (anionic
impurity) or downwards (cationic impurity) band bending in STS. This
translates into a dark depression (Fig. 2a left) or bright protrusion (Fig. 2a
right) of the STM topography if measured at the conduction band edge. Most
importantly, Re${}_{\text{Mo}}^{+1}$ preserves the underlying D3h lattice
symmetry, as previously reported for Re${}_{\text{W}}^{+1}$ in WSe2 39, while
the STM topography of Re${}_{\text{Mo}}^{0}$ and Re${}_{\text{Mo}}^{-1}$
appear slightly distorted.

The pink d$I$/d$V$ spectrum in Fig. 2e and defect orbitals depicted in Fig. S6
for Re${}_{\text{Mo}}^{+1}$ closely resemble the Re${}_{\text{W}}^{+1}$ states
observed in WSe2 39, with several defect states just above the Fermi level.
Re${}_{\text{Mo}}^{0}$ (Fig. 2e blue, and Fig. S7) exhibits several unoccupied
defect states above Fermi, along with an unusually broad defect resonance
below it. Lastly, Re${}_{\text{Mo}}^{-1}$ shown in green in Fig. 2e and Fig.
S8, has only a broad defect resonance below the Fermi level. Based on the STS
spectra we derive a simplified schematic level diagram for the differently
charged Re${}_{\text{Mo}}$ impurities in Fig. 2d. The EG and QFEG act as a
Fermi sea where electrons can be donated to or withdrawn from the
Re${}_{\text{Mo}}$. Upon (de)charging, many-body effects renormalize the
energies of all defect states because of their pronounced localization. For
instance, the unoccupied defect states for Re${}_{\text{Mo}}^{-1}$ are pushed
higher in energy, whereas they become closely spaced at low energy for
Re${}_{\text{Mo}}^{+1}$.

## Symmetry Breaking of Re${}_{\text{Mo}}^{-1}$ and Re${}_{\text{Mo}}^{0}$

To probe whether the striking electronic asymmetry of Re${}_{\text{Mo}}^{0}$
and Re${}_{\text{Mo}}^{-1}$ results from a (measurable) lattice distortion, we
use CO-tip nc-AFM 28 as shown in Fig. 3d-f. Due to the high surface
sensitivity of CO-tip nc-AFM, only the surface S atoms can be resolved as
bright (repulsive) features. The Mo atom positions can be inferred from the
slightly more attractive interaction as compared to the hollow site,
establishing the full unit cell of 1H-MoS2 indicated as an inset in Fig. 3d-f.
Although the Re impurity itself cannot be directly resolved, its location
(marked by the colored triangle) can be deduced from correlation with the STM
topography (Fig. 3a-c). For Re${}_{\text{Mo}}^{0}$ and Re${}_{\text{Mo}}^{-1}$
we find that one of the neighboring S atoms appears brighter, indicating an
outwards relaxation from the S plane, which can be seen in the line-cut below
the respective nc-AFM image. Occasionally, the protruding S atom is found to
switch during nc-AFM measurements at close tip–sample distance (Fig. S9). For
Re${}_{\text{Mo}}^{+1}$ no such relaxation could be resolved. Apart from the
vertical relaxation, which is most pronounced for Re${}_{\text{Mo}}^{-1}$,
also a lateral relaxation of the two opposing S atoms is observed (Fig. 3g).
The vertical protrusion and lateral strain profile are in excellent agreement
with the simulated CO-tip nc-AFM image (Fig. 3h) based on the calculated
relaxed geometry shown in Fig. 3i. 43

We employ DFT with a 2D charge correction scheme 44 to investigate the atomic
and electronic structure of Re${}_{\text{Mo}}$ in MoS2, considering the
positive, neutral, and negative charge states. While earlier ab initio studies
failed to capture the symmetry breaking in the neutral defect state 45; 46;
47, our analysis demonstrates that enforcing integer electronic occupation
successfully reproduces a distorted ground state geometry for
Re${}_{\text{Mo}}^{0}$ and Re${}_{\text{Mo}}^{-1}$. In Fig. 4a-c, we present
the adiabatic potential energy surface (APES) for all three charge states. In
the positive charge state, the APES exhibits a single minimum at the high
symmetry position. However, for Re${}_{\text{Mo}}^{0}$ and
Re${}_{\text{Mo}}^{-1}$, three equivalent off-centered minima are observed at
distances of 13 pm and 25 pm, respectively, towards one of the neighboring S
atoms. Although the structural relaxation of Re away from the center is twice
as large for the negative charge state, the energy gain of approximately -38
meV is comparable for both neutral and negative charge states. The obtained ab
initio APES supports the experimentally observed symmetry-broken ground states
of Re${}_{\text{Mo}}^{0}$ and Re${}_{\text{Mo}}^{-1}$.

To gain insight into the driving force behind the observed symmetry breaking,
we examine the orbital-projected density of states (pDOS) of
Re${}_{\text{Mo}}^{-1}$ in both the symmetric and distorted ground state
geometries, as shown in Fig. 4d and e, respectively. In contrast to the
Vac${}_{\text{S}}^{-1}$ case, the symmetric configuration of
Re${}_{\text{Mo}}^{-1}$ does not occupy a degenerate orbital at the Fermi
energy. Instead, two electrons occupy the the non-degenerate $d_{z^{2}}$ HOMO
(red curve in Fig. 4d). The excited state spectrum reveals degeneracies with
dominant contributions from $d_{xz}$ and $d_{yz}$ (green solid and dashed
line) orbitals and states composed of $d_{xy}$ and $d_{x^{2}-y^{2}}$ (blue
solid and dashed line) orbitals. The structural distortion lowers the energy
of the HOMO by creating a hybridized $d_{z^{2}}-d_{x^{2}-y^{2}}$ orbital that
is doubly occupied and shifts down to the mid-gap region. A priori it is
unclear if the hybridized orbital results from occupying a $d_{z^{2}}$ state
that takes on $d_{x^{2}-y^{2}}$ orbital character, or vice versa, which would
either suggest that the mechanism is attributed to a pseudo-JTE or hidden-JTE,
respectively. To differentiate between a pseudo-JTE and a hidden-JTE, one can
look at the energy along the distortion direction in the APES (Fig. 4a). A
hidden-JTE would typically manifest as an energy barrier between the symmetric
and distorted geometry, with the initially occupied state experiencing an
energy increase that is eventually surpassed by the energy gain of the
previously unoccupied state 19. Consequently, the absence of such an energy
barrier confidently identifies the distortion mechanism in this system as
being governed by the pseudo-JTE.

Similarly, the distortion in the neutral charge state is also caused by the
pseudo-JTE (cf. Fig. S13). In the neutral charge state, the odd number of
electrons gives rise to a spin-polarized ground state, wherein a singly
occupied $d_{z^{2}}-d_{x^{2}-y^{2}}$ SOMO is established (Fig. 4f). This leads
to a net magnetic moment that is absent in the positive or negative charged
defects. The presence of a paramagnetic ground state in the neutral defect
with a singly occupied spin channel is not captured by previous ab initio
descriptions of Re${}_{\text{Mo}}^{0}$ in MoS2 46. However, experimental
studies using electron paramagnetic resonance have successfully detected Re
doping in natural MoS2 samples, providing further validation for our findings
48. Furthermore, the spatial distribution of the highest occupied
($\gamma^{\prime}$) and lowest unoccupied ($\beta^{\prime}$)
Re${}_{\text{Mo}}^{0}$ orbitals shown in Fig. 4g exhibit the same spatial
distribution in perfect agreement with the experimental d$I$/d$V$ images, as
shown in Fig. 5b (resonances C3 and D1). Therefore, we conclude that the
symmetry-broken paramagnetic ground state observed in the neutral
Re${}_{\text{Mo}}$ substitution reconciles our theoretical predictions with
experimental observations based on scanning probe microscopy and previous
measurements using electron paramagnetic resonance.

## Configurational Continuum and Directional Alignment

The unusually broad defect resonances on the order of 300 meV observed for
some Re${}_{\text{Mo}}^{0}$ and Re${}_{\text{Mo}}^{-1}$ states raises the
question about the broadening mechanism. For many TMD defects, the intrinsic
lifetime broadening primarily caused by the weakly interacting graphene
substrate, is only a few meV 34. We do not expect a qualitatively different
substrate interaction for Re dopants, as it is mainly influenced by the
distance to the substrate, and the measured defect state broadening does not
decrease for Re in bilayer MoS2. In low-temperature tunneling spectroscopy of
single defects in 2D semiconductors, significant broadening of electronic
resonances can occur due to strong electron-phonon coupling 49. However, in
most cases the confined electronic states couple strongly only to a small set
of localized modes with energies of a few tens of meV, given the weight of the
atoms involved 49. Therefore, it is typically possible to differentiate the
vibronic ground state resonance (zero-phonon line) and the first few vibronic
excited states with an intrinsic broadening of only 3-4 meV, as can be seen
for instance in the d$I$/d$V$ spectra of the sulfur vacancy in Fig. 1d. Hence,
vibronic broadening is unlikely the main reason for the large defect state
broadening.

Remarkably, we observe that the orbital image of the broad defect state
resonance undergoes a continuous evolution of contrast, depending on the
energy at which the defect state is probed, as shown in Fig. 5a,b. This
behavior is observed on both mono- and bilayer MoS2, although the direction of
the progression varies depending on the precise location within the sample
(cf. series Fig. 5b and Fig. S7c). This unusual behavior might be a
consequence of the comparably flat adiabatic potential energy surface
introduced by the pseudo-JTE 50; 18, establishing a configurational continuum
without significant energy barrier between the states (Fig. 4a-c). This
phenomena is common for Jahn-Teller systems and is typically referred to as
dynamic JTE 50; 18. Such a configurational continuum may explain the unusual
progression of the STS contrast within the defect resonance and its
significant broadening.

In most regions of the sample, the symmetry breaking direction is not a random
distribution. Rather, there are larger domains where all
Re${}_{\text{Mo}}^{0}$ or Re${}_{\text{Mo}}^{-1}$ relax in the same direction,
as illustrated in Fig. 5c-e. Directionally-aligned domains have been observed
not only in highly doped (5%) samples, but also in samples with comparably low
doping concentrations ($<0.1$%). Therefore, we suspect that residual strain in
the epitaxially grown TMD, rather than direct dopant-dopant interactions,
could be responsible for the aligned domains. The flat APES may render the
symmetry-broken Re particularly susceptible to external perturbations such as
strain. In turn, this susceptibility may serve as a local strain sensor.
Further studies are needed to establish a firm correlation between the
mesoscale strain tensor and local Jahn-Teller distortion of the individual Re
dopants.

## Conclusions

In summary, we have directly visualized the charge state-dependent symmetry
breaking at S vacancies and Re substituting for Mo in mono- and bilayer MoS2.
Vac${}_{\text{S}}^{-1}$ (in the top and bottom) coexists in both a threefold
symmetric and symmetry-broken state. By controlling the Fermi level alignment
through the substrate chemical potential, we are able to stabilize Re in three
different charge states. While Re${}_{\text{Mo}}$ in the positive charge state
retains the underlying D3h lattice symmetry, Re${}_{\text{Mo}}^{0}$ and
Re${}_{\text{Mo}}^{-1}$ exhibit an increasingly distorted geometry. We assign
the minute local lattice relaxation detected by CO-tip nc-AFM and significant
electronic reconfiguration of the defect states measured by STS to a pseudo-
JTE, corroborated by DFT calculations. A progressive defect orbital evolution
observed within the unusually broad defect resonances is attributed to a flat
potential energy surface, which leads to a continuum of the ground state
geometry configuration. Complementary measurements are needed to corroborate
the link between the intriguing observation of directionally aligned Re
domains and an inherent strain profile in the MoS2 layers. The charge-state
dependent JTE exhibited by chalcogen vacancies and Re dopants offers a
powerful means to precisely control the photo-physics of these prototypical
atomic quantum emitters.

## References

## References

  * (1) Philips, S. G. _et al._ Universal control of a six-qubit quantum processor in silicon.  _Nature_ 609, 919–924 (2022).
  * (2) Kiczynski, M. _et al._ Engineering topological states in atom-based semiconductor quantum dots.  _Nature_ 606, 694–699 (2022).
  * (3) Awschalom, D. D., Hanson, R., Wrachtrup, J. & Zhou, B. B.  Quantum technologies with optically interfaced solid-state spins.  _Nat. Photon_ 12, 516–527 (2018).
  * (4) Atatüre, M., Englund, D., Vamivakas, N., Lee, S.-Y. & Wrachtrup, J.  Material platforms for spin-based photonic quantum technologies.  _Nat. Rev. Mater._ 3, 38–51 (2018).
  * (5) Dolde, F. _et al._ Room-temperature entanglement between single defect spins in diamond.  _Nat. Phys._ 9, 139–143 (2013).
  * (6) Tran, T. T., Bray, K., Ford, M. J., Toth, M. & Aharonovich, I.  Quantum emission from hexagonal boron nitride monolayers.  _Nat. Nanotechnol._ 11, 37–41 (2016).
  * (7) Gottscholl, A. _et al._ Initialization and read-out of intrinsic spin defects in a van der waals crystal at room temperature.  _Nat. Mater._ 19, 540–545 (2020).
  * (8) Aharonovich, I., Englund, D. & Toth, M.  Solid-state single-photon emitters.  _Nat. Photon._ 10, 631–641 (2016).
  * (9) Wu, S. _et al._ Monolayer semiconductor nanocavity lasers with ultralow thresholds.  _Nature_ 520, 69–72 (2015).
  * (10) Liu, X. & Hersam, M. C.  2d materials for quantum information science.  _Nat. Rev. Mater._ 4, 669–684 (2019).
  * (11) Grosso, G. _et al._ Tunable and high-purity room temperature single-photon emission from atomic defects in hexagonal boron nitride.  _Nat. Commun._ 8, 1–8 (2017).
  * (12) Noh, G. _et al._ Stark tuning of single-photon emitters in hexagonal boron nitride.  _Nano Lett._ 18, 4710–4715 (2018).
  * (13) Klein, J. _et al._ Site-selectively generated photon emitters in monolayer mos2 via local helium ion irradiation.  _Nat. Commun._ 10, 2755 (2019).
  * (14) Ulbricht, R. _et al._ Jahn-teller-induced femtosecond electronic depolarization dynamics of the nitrogen-vacancy defect in diamond.  _Nat. Commun._ 7, 13510 (2016).
  * (15) Lee, Y. _et al._ Spin-defect qubits in two-dimensional transition metal dichalcogenides operating at telecom wavelengths.  _Nat. Commun._ 13, 7501 (2022).
  * (16) Thiering, G. & Gali, A.  Ab initio calculation of spin-orbit coupling for an nv center in diamond exhibiting dynamic jahn-teller effect.  _Phys. Rev. B_ 96, 081115 (2017).
  * (17) Streltsov, S. V. & Khomskii, D. I.  Jahn-Teller Effect and Spin-Orbit Coupling: Friends or Foes? _Phys. Rev. X_ 10, 031043 (2020).
  * (18) Bersuker, I. B.  Jahn–Teller and Pseudo-Jahn–Teller Effects: From Particular Features to General Tools in Exploring Molecular and Solid State Properties.  _Chem. Rev._ 121, 1463–1512 (2021).
  * (19) Bersuker, I.  Four modifications of the jahn–teller effects. the problem of observables: spin–orbit interaction, tunneling splitting, and orientational polarization of solids.  _Phys. Chem. Chem. Phys_ (2023).
  * (20) Gali, Á.  Ab initio theory of the nitrogen-vacancy center in diamond.  _Nanophotonics_ 8, 1907–1943 (2019).
  * (21) Komsa, H.-P. & Krasheninnikov, A. V.  Native defects in bulk and monolayer MoS2 from first principles.  _Phys. Rev. B_ 91, 125304 (2015).
  * (22) Tan, A. M. Z., Freysoldt, C. & Hennig, R. G.  Stability of charged sulfur vacancies in 2D and bulk MoS2 from plane-wave density functional theory with electrostatic corrections.  _Phys. Rev. Materials_ 4, 064004 (2020).
  * (23) Cheng, Y. C., Zhu, Z. Y., Mi, W. B., Guo, Z. B. & Schwingenschlögl, U.  Prediction of two-dimensional diluted magnetic semiconductors: Doped monolayer MoS2 systems.  _Phys. Rev. B_ 87, 100401 (2013).
  * (24) Bertoldo, F. _et al._ Intrinsic Defects in MoS2 Grown by Pulsed Laser Deposition: From Monolayers to Bilayers.  _ACS Nano_ 15, 2858–2868 (2021).
  * (25) Klein, J. _et al._ Site-selectively generated photon emitters in monolayer MoS2 via local helium ion irradiation.  _Nat. Commun._ 10, 2755 (2019).
  * (26) Gupta, S., Yang, J.-H. & Yakobson, B. I.  Two-Level Quantum Systems in Two-Dimensional Materials for Single Photon Emission.  _Nano Lett._ 19, 408–414 (2019).
  * (27) Repp, G., Meyer, G., Stojkovic, S. M., Gourdon, A. & Joachim, C.  Molecules on Insulating Films: Scanning-Tunneling Microscopy Imaging of Individual Molecular Orbitals.  _Phys. Rev. Lett._ 94, 026803 (2005).
  * (28) Gross, L., Mohn, F., Moll, N., Liljeroth, P. & Meyer, G.  The Chemical Structure of a Molecule Resolved by Atomic Force Microscopy.  _Science_ 325, 1110 (2009).
  * (29) Komsa, H.-P. _et al._ Two-Dimensional Transition Metal Dichalcogenides Under Electron Irradiation: Defect Production and Doping.  _Phys. Rev. Lett._ 109, 035503 (2012).
  * (30) Zhou, W. _et al._ Intrinsic Structural Defects in Monolayer Molybdenum Disulfide.  _Nano Lett._ 13, 2615–2622 (2013).
  * (31) Tongay, S. _et al._ Defects Activated Photoluminescence in Two-Dimensional Semiconductors: Interplay Between Bound, Charged, and Free Excitons.  _Sci. Rep._ 3, 2657 (2013).
  * (32) Mitterreiter, E. _et al._ The role of chalcogen vacancies for atomic defect emission in MoS2.  _Nat. Commun._ 12, 3822 (2021).
  * (33) Barja, S. _et al._ Identifying Substitutional Oxygen as a Prolific Point Defect in Monolayer Transition Metal Dichalcogenides with Experiment and Theory.  _Nat. Commun._ 10, 3382 (2019).
  * (34) Schuler, B. _et al._ Large Spin-Orbit Splitting of Deep In-Gap Defect States of Engineered Sulfur Vacancies in Monolayer WS2.  _Phys. Rev. Lett._ 123, 076801 (2019).
  * (35) Schuler, B. _et al._ How Substitutional Point Defects in Two-Dimensional WS2 Induce Charge Localization, Spin-Orbit Splitting, and Strain.  _ACS Nano_ 13, 10520–10534 (2019).
  * (36) Mitterreiter, E. _et al._ Atomistic Positioning of Defects in Helium Ion Treated Single-Layer MoS2.  _Nano Lett._ 20, 4437–4444 (2020).
  * (37) Aghajanian, M. _et al._ Resonant and bound states of charged defects in two-dimensional semiconductors.  _Phys. Rev. B_ 101, 081201(R) (2020).
  * (38) Zhu, G.-J., Xu, Y.-G., Gong, X.-G., Yang, J.-H. & Yakobson, B. I.  Dimensionality-Inhibited Chemical Doping in Two-Dimensional Semiconductors: The Phosphorene and MoS2 from Charge-Correction Method.  _Nano Lett._ 21, 6711–6717 (2021).
  * (39) Kozhakhmetov, A. _et al._ Scalable Substitutional Re-Doping and its Impact on the Optical and Electronic Properties of Tungsten Diselenide.  _Adv. Mater._ 32, 2005159 (2020).
  * (40) Torsi, R. _et al._ Dilute Rhenium Doping and its Impact on Intrinsic Defects in MoS2 (2023).  ArXiv:2302.00110 [cond-mat].
  * (41) Pan, Y. _et al._ WSe2 homojunctions and quantum dots created by patterned hydrogenation of epitaxial graphene substrates.  _2D Mater._ 6, 021001 (2019).
  * (42) Subramanian, S. _et al._ Tuning transport across MoS2/graphene interfaces via as-grown lateral heterostructures.  _npj 2D Mater. Appl._ 4, 1–5 (2020).
  * (43) Hapala, P. _et al._ Mechanism of high-resolution STM/AFM imaging with functionalized tips.  _Phys. Rev. B_ 90, 085421 (2014).
  * (44) Freysoldt, C. & Neugebauer, J.  First-principles calculations for charged defects at surfaces, interfaces, and two-dimensional materials in the presence of electric fields.  _Phys. Rev. B_ 97, 205425 (2018).
  * (45) Lin, Y.-C. _et al._ Properties of Individual Dopant Atoms in Single-Layer MoS2: Atomic Structure, Migration, and Enhanced Reactivity.  _Adv. Mater._ 26, 2857–2861 (2014).
  * (46) Zhao, P. _et al._ Electronic and magnetic properties of Re-doped single-layer MoS2: A DFT study.  _Comp. Mater. Sci._ 128, 287–293 (2017).
  * (47) Hallam, T. _et al._ Rhenium-doped MoS2 films.  _Appl. Phys. Lett._ 111, 203101 (2017).
  * (48) Stesmans, A., Schoenaers, B. & Afanas’ev, V. V.  Variations of paramagnetic defects and dopants in geo-MoS2 from diverse localities probed by ESR.  _J. Chem. Phys._ 152, 234702 (2020).
  * (49) Cochrane, K. A. _et al._ Spin-dependent vibronic response of a carbon radical ion in two-dimensional WS2.  _Nat. Commun._ 12, 7287 (2021).
  * (50) Bersuker, I. B.  The Jahn-Teller and pseudo Jahn-Teller effect in materials science.  _J. Phys.: Conf. Ser._ 833, 012001 (2017).

## Figures

Figure 1:  Symmetric and symmetry-broken (distorted) negatively charged sulfur
vacancies (Vac${}_{\text{S}}^{-1}$). a, STM topography of annealing-induced
Vac${}_{\text{S}}^{-1}$ (colored circles) in monolayer MoS2/QFEG/SiC. b, Zoom-
in STM topography of the white dashed box in (a) containing a distorted (top
right) and a symmetric (bottom left) Vac${}_{\text{S}}^{-1}$. c, Nc-AFM image
with a metallic tip of the same area as in (b) proving that both defects are
Vac${}_{\text{S}}$ in the upper S layer. d, d$I$/d$V$ spectra of the frontier
in-gap defect orbitals of the symmetric (blue) and distorted (red)
Vac${}_{\text{S, top}}^{-1}$, as well as the distorted Vac${}_{\text{S,
bottom}}^{-1}$ (orange) ($V_{\text{mod}}$ = 2 mV). e, Corresponding constant
height dI/dV maps of the distorted S top (upper row) and bottom (bottom row)
vacancies highest occupied (HOMO) and lowest unoccupied defect orbitals (LUMO,
LUMO+1), labelled in d ($V_{\text{mod}}$ = 20 mV). f,g, Projected density of
states (pDOS) of Vac${}_{\text{S}}^{-1}$ in MoS2 in the symmetrized (f) and
distorted (g, ground state) geometry, calculated in a 4$\times$4 supercell.
The pDOS onto the Mo $d_{xy}$ and $d_{x^{2}-y^{2}}$ orbitals are shown in red.
h,i, Constant height DOS map of the HOMO and LUMO(+1) orbitals for the
symmetrized (h) and diststorted (i) Vac${}_{\text{S, top}}^{-1}$ geometry
Figure 2:  Charge state tristability of Re dopants in MoS2 on EG and QFEG on
SiC. a, STM topographies of single Re dopants in the three different charge
states (Re${}_{\text{Mo}}^{-1}$, Re${}_{\text{Mo}}^{0}$, and
Re${}_{\text{Mo}}^{+1}$). b 3D representation of a large-scale STM image of
Re-doped monolayer MoS2 on domains of EG (center blue region) or domains of
QFEG substrate (outer gray region). At positive bias (here 0.6 V), Re dopants
are imaged as a dark blue pits on EG whereas they appear as red hills on MoS2
grown on QFEG, depending on their charge states. c, Atomic model of the
MoS2/EG and MoS2/QFEG interface. On EG only Re${}_{\text{Mo}}^{-1}$ are
observed, whereas on QFEG both Re${}_{\text{Mo}}^{0}$ and
Re${}_{\text{Mo}}^{+1}$ are found. d, Schematic electronic level diagram of
the MoS2/EG (left) and MoS2/QFEG interface (right) with the Re${}_{\text{Mo}}$
level occupation indicated. The work function of EG $\Phi_{\text{EG}}$ is
substantially smaller than the work function of QFEG $\Phi_{\text{QFEG}}$.
$\chi_{\text{MoS${}_{2}$}}$ denotes the electron affinity of MoS2. e, dI/dV
spectra of Re${}_{\text{Mo}}^{-1}$ (green), Re${}_{\text{Mo}}^{0}$ (blue), and
Re${}_{\text{Mo}}^{+1}$ (pink), pristine monolayer MoS2 on EG (grey) and on
QFEG (black), respectively. Bias modulation amplitude: $V_{\text{mod}}$ = 20
mV.  Figure 3:  Charge state-dependent Jahn-Teller distortions. a-c, STM
topography of Re${}_{\text{Mo}}^{-1}$, Re${}_{\text{Mo}}^{0}$, and
Re${}_{\text{Mo}}^{+1}$, respectively. d-f, CO-tip nc-AFM image of the same Re
impurity as shown in a-c. Top sulfur atoms are repulsive (bright) and
molybdenum atoms attractive (dark). The MoS2 unit cell is displayed at the
bottom left in each image (S: yellow, Mo: turquoise). The Re atoms are located
in the center of the dashed triangle, bound to the three bright S atoms on the
top surface. A $\Delta\text{f}$ line profile across a sulfur row (arrow in
image) is shown in the lower panel, highlighting a protruding S atom next to
the Re impurity for Re${}_{\text{Mo}}^{-1}$ and slightly weaker for
Re${}_{\text{Mo}}^{0}$. g,h, Measured (g) and simulated (h) CO-tip nc-AFM
image of Re${}_{\text{Mo}}^{-1}$. The apparent sulfur-sulfur distances in the
top S layer is indicated by the red-to-yellow lines. Measurement and
simulation show significant lateral compressive strain introduced by the Re
impurity at nearest-neighbor S atoms opposite to the protruding S atom. i,
Calculated geometry of Re${}_{\text{Mo}}^{-1}$ exhibiting a considerable
distortion. The red arrows indicate the atom relaxations as compared to the
pristine MoS2 lattice (gray ball and stick model). For clarity, the arrows are
ten times longer than the actual relaxation. Figure 4:  Adiabatic potential
energy surface (APES) and projected DOS (pDOS) of Re${}_{\text{Mo}}$ in
different charge states. a-c, Potential energy surface of the Re dopant for
the negative, neutral, and positive charge state respectively. The energies
correspond to the $\Gamma$ point energy in a $5\times 5$ supercell. The black
and white dots indicate the symmetric and distorted positions of Re,
respectively. The sampled area for the APES is marked by a red parallelogram
in each structural model inset. d,e, Projected density of states of
Re${}_{\text{Mo}}^{-1}$ in the constrained symmetric (d) and distorted ground
state geometry (e), exhibiting electronic reconfiguration resulting from the
pseudo-JTE. f, Spin-polarized Re${}_{\text{Mo}}^{0}$ band structure resulting
in a paramagnetic ground state. g, Constant height DOS maps of different
energy windows in the Re${}_{\text{Mo}}^{0}$ band structure, as indicated in
f. The $\alpha^{\prime}$ image reflects the STM contrast at low bias. The
$\beta^{\prime}$ and $\gamma^{\prime}$ image corresponds to the lowest
unoccupied and highest occupied Re${}_{\text{Mo}}^{0}$ orbital, respectively
(cf. Fig. 5b C3 and D1).  Figure 5:  Configurational continuum and
directionally aligned domains. a, dI/dV spectroscopy of Re${}_{\text{Mo}}^{0}$
on bilayer MoS2 ($V_{\text{mod}}$ = 10 mV). The spectroscopy locations on the
defect are indicated with a dot of the same color in the inset STM image, and
the pristine MoS2 (2ML) spectrum is shown in black. b, Constant height dI/dV
maps ($V_{\text{mod}}$ = 20 mV) at the sample voltages labelled in a. A
progression of the orbital image is observed within certain resonances (cf.
D1-D5) that might explain their large broadening. The resemblance of the
highest occupied (D1) and lowest unoccupied defect orbital (C3) is indicative
of a singly occupied defect state. c, STM topography of neutral Re dopants
(top layer) in bilayer MoS2 on QFEG. The orientation of the distorted STM
contrast along the crystal axes is indicated by colored arrows. d,e, STM
topography of Re${}_{\text{Mo}}^{0}$ in bilayer MoS2 in a different location
also showing a preferential orientation. Here, the D${}_{\text{3h}}$ symmetry
breaking is observed at positive and negative bias, however the respective
symmetry axis is canted by 60°.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Acknowledgements

We thank Adam Gali for valuable discussions. F.X. acknowledges the Walter
Benjamin Programme from Deutsche Forschungsgemeinschaft (DFG). L.H., J.A. and
B.S. appreciate funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (Grant agreement
No. 948243). R.T., Y.-C. L, and J.A.R. acknowledge funding from NEWLIMITS, a
center in nCORE as part of the Semiconductor Research Corporation (SRC)
program sponsored by NIST through award number 70NANB17H041 and the Department
of Energy (DOE) through award number DESC0010697. R.T., C.D., and J.A.R.
acknowledge funding from the 2D Crystal Consortium, National Science
Foundation Materials Innovation Platform, under cooperative agreement
DMR-1539916. For the purpose of Open Access, the author has applied a CC BY
public copyright license to any Author Accepted Manuscript version arising
from this submission.

## Author contributions

† These authors contribute equally to this work.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information The online version contains supplementary material
available at https://doi.org/

Correspondence and requests for materials should be addressed to Bruno
Schuler.