2.1 Definition of pensive billiards

We begin with a planar domain $D\subset \mathbb {R}^2$ whose boundary $\gamma = \partial D$ is a smooth closed curve. The phase space of a pensive billiard system is given by

$$ \begin{align*} M=\{(q,v) ~|~q\in \gamma \text{ and } v\in T_q\mathbb{R}^2 \text{ an inward pointing unit vector}\}. \end{align*} $$

It represents the location and direction of motion of a billiard ball right after impact with the boundary. The classical billiard ball map $\mathsf {CB}\colon M\to M$ sends $(q_1,v_1)\mapsto (q_2,v_2)$ , as obtained by traveling the table from $q_1$ in direction $v_1$ until one hits $\gamma $ at $q_2$ , and then reflecting the velocity vector $v_1$ about $T_{q_2}\gamma $ to obtain $v_2$ according to the usual reflection law.

There are several ways to introduce coordinates in M. Parametrize $\gamma $ by arc length as $\gamma (s)$ for $s\in [0,2L]$ , and, given $v\in T_{\gamma (s)}\mathbb {R}^2$ , denote by $\theta \in (0,\pi )$ the angle between v and $\gamma '(s)$ . Thus M is a cylinder with coordinates $(s,\theta )$ . Alternatively, let $p=\mathrm {proj} \,v|_{T_q\gamma }$ be the projection of v to the tangent space to $\gamma $ at the point q. (Note that $p:=\cos \theta $ .) Then M can be identified with the (co)tangent bundle of unit balls,

$$ \begin{align*} M=T^{(*)}_{< 1}\gamma=\{(q, p)~|~q\in \gamma, \ -1 <p<1\} \end{align*} $$

(the Euclidean structure in $\mathbb {R}^2$ allows one to identify tangent and cotangent spaces). Recall that M is a symplectic manifold: Let $\lambda =p\,dq$ and $\omega =d\lambda $ be the Liouville one-form and the corresponding standard symplectic structure $\omega $ on $T^*\gamma $ .

Now we introduce pensive billiards. Fix an arbitrary smooth delay function $\ell (p)$ for $p\in (-1,1)$ (alternatively, sometimes it is convenient to rewrite this function as $\tilde \ell (\theta ):=\ell (\cos \theta )$ ). This function will stand for the length of the path along the boundary (in the arc-length parametrization) that a billiard ball will spend between hitting the boundary (at the angle $\theta $ where $p=\cos \theta $ ) and reflecting from it.

Figure 2 Pensive Billiard.

Alternatively, one can formulate this map in $(s, p)$ coordinates as

$$ \begin{align*} \mathsf{PB}:(s_1, p_1)\mapsto (s_2+\ell(p_2), p_2)\in T_{< 1}\gamma, \end{align*} $$

resembling the classical billiard which contemplated for time $\ell (p)$ before reflecting, hence the name.

In what follows, we will denote the components of pensive and classical billiard maps by

$$ \begin{align*} (S(s,p), P(s,p)) := \mathsf{PB}(s,p) \quad \text{and} \quad (S_{\mathsf{cl}}(s,p), P_{\mathsf{cl}}(s,p)) := \mathsf{CB}(s,p). \end{align*} $$

When it is more convenient to work in the $(s, \theta )$ coordinates, we will write, analogously,

$$ \begin{align*} (S(s,\theta), \Theta(s,\theta)) := \mathsf{PB}(s,\theta) \quad \text{and} \quad (S_{\mathsf{cl}}(s,\theta), \Theta_{\mathsf{cl}}(s,\theta)) := \mathsf{CB}(s,\theta). \end{align*} $$

While the domain D does not have to be convex for the pensive billiard map to be well defined, there are some additional good properties when this is the case, as we discuss below. Without convexity, the corresponding billiard maps (both classical and pensive) might not even be continuous.

Proof. First note that the pensive billiard map $\mathsf {PB}$ is a composition of the classical billiard $\mathsf {CB}$ and a shift map $\mathsf {Sh}: (s, \theta )\mapsto (s+\tilde \ell (\theta ), \theta )$ . Then the invariance for $\mathsf {PB}$ immediately follows from that for the classical billiard:

$$ \begin{align*}\mathsf{PB}^*\omega=\mathsf{Sh}^*(\mathsf{CB}^*\omega)=\mathsf{Sh}^*\omega=\sin\theta\,d\theta\wedge d(s+\tilde \ell(\theta))= \sin\theta\,d\theta\wedge ds=\omega, \end{align*} $$

as required. For completeness, we recall the proof of the corresponding lemma for the classical billiard map (see, e.g., [Reference Tabachnikov29]):

Figure 3 Puck Billiard. For a puck billiard, generating function (2.3) involves potential $V(p)=\widetilde {V}(\theta )$ defined by equation (2.1) and equal to the length of the geodesic’s part on the side surface of the puck.

Figure 4 Pensive Billiard arising as a limit of point vortex motion.

Lemma 2.4. The area form $\omega =\sin \theta \, d\theta \wedge dt$ on $M=T^*_{< 1}\gamma $ is $\mathsf {CB}$ -invariant (i.e., the classical billiard map is a symplectomorphism).

Indeed, let $f(s_1,s_2)=\|\gamma (s_1)- \gamma (s_2)\|$ , where $\|\cdot \|$ is the Euclidean distance function in $\mathbb {R}^2$ . Then, $\frac {\partial f}{\partial s_2}$ is the projection of $\nabla f$ onto the tangent space to the curve $\gamma $ at $\gamma (s_2)$ , so $\frac {\partial f}{\partial s_2}=\cos \theta _2$ , since $\nabla f$ is a unit vector making angle $\theta _2$ with (the positive tangent to) $\gamma $ . Similarly, $\frac {\partial f}{\partial s_1}=-\cos \theta _1$ , so

$$ \begin{align*} df=\frac{\partial f}{\partial s_1} ds_1 +\frac{\partial f}{\partial s_2}ds_2=-\cos\theta_1 \, d s_1 +\cos \theta_2 \, d s_2. \end{align*} $$

This can be understood as follows: on $M\times M$ the one-form in the right-hand side becomes a complete differential $df$ upon restriction to the graph of the billiard map $\mathsf {CB}$ . Now, by taking one more differential we obtain

$$ \begin{align*} 0=d^2f=d\Big{(} \frac{\partial f}{\partial s_1} ds_1 +\frac{\partial f}{\partial s_2}ds_2\Big{)}=\sin\theta_1 \,d\theta_1 \wedge ds_1 -\sin \theta_2\, d\theta_2 \wedge ds_2, \end{align*} $$

which completes the proof of both the lemma and theorem.

2.2 Examples of delay functions

Here are three examples of pensive billiard systems:

A) Satellite billiard (see also billiards with slipping in [Reference Fomenko, Vedyushkina and Zav’yalov12]). Imagine that a satellite moving with constant velocity on the orbit receives and sends back a signal according to the billiard law, but requires some time to process the signal. Assuming that the processing takes the same amount of time independent of the angle of incidence, we arrive at the definition of pensive billiard with constant delay function, $\ell (p)\equiv \mathrm {const}$ . The classical billiard corresponds to $\ell (p)\equiv 0$ .

B) Puck billiard (also called a coin billiard, see [Reference Bialy, Fierobe, Glutsyuk, Levi, Plakhov and Tabachnikov4]). Consider the geodesic flow on the surface of a cylinder of height h with base $\gamma =\partial D$ , glued on the top and on the bottom to two copies of D. On the top and the bottom of the cylinder the motion is along straight lines, while on the cylindrical part it goes along geodesics of the cylinder. Identifying top and bottom copies of D, such a motion is described by the pensive billiard, where the delay function $\tilde \ell (\theta )$ is the base of the right triangle of height h and angle $\theta $ opposite to the base, that is, $\tilde \ell (\theta )=h\cot \theta $ , or equivalently, $\ell (p)=hp/\sqrt {1-p^2}$ for $p=\cos \theta $ .

C) Vortex billiard. This pensive billiard arises naturally in the motion of point vortex dipoles in a bounded two-dimensional domain D. Here we only mention that it corresponds to the delay function $\ell (p)=L(1 -p/\sqrt {1+p^2})$ , where the length of $\gamma =\partial D$ is equal to $2L$ . We discuss the vortex billiard in detail in Section 4.

2.3 Variational principle

Here we prove that the pensive billiard with any delay function is subject to a variational principle.

Let $\gamma (s)$ be a point on the boundary $\gamma =\partial D$ and $A\in D$ a point inside the domain. Recall that the cosine of the incidence angle at $\gamma (s)$ is $\cos \theta =p_A(s) = \frac {(\gamma (s)-A) \cdot \gamma '(s)}{\|\gamma (s) - A\|}.$

Theorem 2.5. Given two points $A, B \in D$ , the pensive billiard trajectory that goes from A to B after hitting boundary once must hit $\gamma $ at point $\gamma (s)$ , where s is a critical point of the function

$$ \begin{align*} f_{A,B}(s) = \|A -\gamma(s)\| + V(p_A(s)) + \|\gamma(s + \ell(p_A(s))) - B\|, \end{align*} $$

where a “potential” V is defined by the integral

(2.1) $$ \begin{align} V(p) = \int_0^{p}q\,\ell'({q})\,d{q}. \end{align} $$

Proof. For convenience, define $d_A(s) = \|\gamma (s) - A\|$ , and observe that

$$ \begin{align*} \frac{d}{ds}d_A(s) = \nabla_Q\|Q - A\||_{Q = \gamma(s)}\cdot \gamma'(s) = \frac{(\gamma(s)-A)}{\|\gamma(s) - A\|}\cdot \gamma'(s)= p_A(s). \end{align*} $$

We then compute

$$ \begin{align*} \frac{d}{ds}f_{A,B}(s) &= \frac{d}{ds}d_A(s) + \frac{d}{ds}V(p_A(s)) + \frac{d}{ds}d_B(s + \ell(p_A(s))) \\ &= p_A(s) + V'(p_A(s))p_A'(s) + p_B(s + \ell(p_A(s)))(1 + \ell'(p_A(s))p_A'(s))\\ &= p_A(s) + p_A(s)\ell'(p_A(s))p_A'(s) + p_B(s + \ell(p_A(s)))(1 + \ell'(p_A(s))p_A'(s))\\ &= (p_A(s) + p_B(s + \ell(p_A(s))))(1 + \ell'(p_A(s))p_A'(s)). \end{align*} $$

Finally, recall that pensive billiard trajectories are characterized by $p_A(s) + p_B(s + \ell (p_A(s))) = 0$ , equivalent to the equality of the angles of incidence and reflection.

Remark 2.6. The puck billiard is given by $\ell (p) = h \cot (\arccos p) = h\frac {p}{\sqrt {1 - p^2}}$ . A direct computation gives that in this case potential (2.1) is

$$ \begin{align*} V(p) = \int_0^{p}q\,\ell'(q)\,dq = h\int_0^{p}\frac{q}{(1 - q^2)^{3/2}}dq= \frac{h}{\sqrt{1 - p^2}} = \frac{h}{\sin(\arccos p)}. \end{align*} $$

In this case, potential V has a natural interpretation: it is equal to the length of “curvilinear hypotenuse,” as shown in Figure 3. This interpretation also extends to the case of generalized puck billiards (see Section 3).

Remark 2.7. Geometrically, the potential

$$ \begin{align*}V(p) = \int_0^{p}q\,\ell'(q)\,dq = p\,\ell(p)- \int_0^{p} \ell(q)\,dq \end{align*} $$

is equal to the area above the graph of the delay function $\ell (q)$ inside the rectangle $[0,p]\times [0,\ell (p)]$ . Such an area naturally arises in the context of generating functions for twist maps, see Section 2.5.

2.4 Generating function

It turns out that a large class of pensive billiards admits explicit generating functions. To specify this class, we introduce the following assumption.

By the implicit function theorem, this condition implies that in a neighborhood U of $(s_0, S_0) \in \gamma \times \gamma $ , for any choice of points $(s, S) \in U$ , one could smoothly choose a direction $p(s,S)$ such that $\mathsf {PB}(s, p(s,S)) = (S, P)$ . (Note that global transitivity is equivalent to the twist condition, as discussed below.)

For classical billiards, this property is a consequence of convexity of the domain. For general pensive billiards, however, this assumption is not satisfied automatically. An example of a situation in which this assumption fails dramatically is a pensive billiard on the unit disk with delay function $\tilde {\ell }(\theta ) = -2\theta $ . One can easily observe that in this case pensive billiard map is the identity (see equation (2.8) below), and consequently the billiard is not transitive.

Note that the example above is nongeneric even within the class of linear delays: one can see that for any other delay function of the form $\tilde {\ell }(\theta ) = C\theta $ with $C \neq -2$ , pensive billiard on a disk is transitive.

Now we can introduce a generating function for transitive pensive billiards.

Definition 2.9. Consider a pensive billiard on $\gamma $ with smooth delay function $\ell $ . Assume that it is locally transitive at $(s_0,p_0)$ . Given $s, S \in (0, L)$ sufficiently close to $s_0, S_0$ denote by $p^* = p^*(s, S)$ a direction in which the pensive billiard ball needs to leave from $\gamma (s)$ to arrive at $\gamma (S)$ , that is, a solution of

(2.2) $$ \begin{align} S_{\mathsf{cl}}(s, p^*) + \ell(P_{\mathsf{cl}}(s, p^*)) = S, \end{align} $$

the existence of which is guaranteed by transitivity. Recall the potential

$$ \begin{align*} V(p) := \int_{0}^{p} q \,\ell'(q) dq, \end{align*} $$

and define a map $H: \mathbb {R}^2 \to \mathbb {R}$ as

(2.3) $$ \begin{align} \begin{aligned} H(s, S) &:= H_{\mathsf{cl}}(s, S_{\mathsf{cl}}(s, p^*)) + V(P_{\mathsf{cl}}(s,p^*)), \end{aligned} \end{align} $$

where $H_{\mathsf {cl}}(s,S) := ||\gamma (s) - \gamma (S)||$ is a standard generating function for classical billiards. Map H is called a generating function of a pensive billiard with a delay function $\ell $ .

This definition is motivated by the proposition (well known for classical billiards; see [Reference Tabachnikov29]):

Proposition 2.10. Any locally transitive pensive billiard possesses a generating function $ H(s, S)$ given by (2.3). Namely, if $(S, P) = \mathsf {PB}(s, p)$ then

(2.4) $$ \begin{align} \frac{\partial H }{\partial s}(s, S) = -p \quad \text{and} \quad \frac{\partial H }{\partial S}(s, S) = P. \end{align} $$

Proof. Proposition follows from the corresponding one for classical billiards. Recall from [Reference Tabachnikov29] that for $H_{\mathsf {cl}}$ we have that if $\mathsf {CB}(s, p) = (S_{\mathsf {cl}}, P_{\mathsf {cl}})$ , then

(2.5) $$ \begin{align} \frac{\partial H_{\mathsf{cl}} }{\partial s}(s, S_{\mathsf{cl}}) = -p \quad \text{and} \quad \frac{\partial H_{\mathsf{cl}} }{\partial S_{\mathsf{cl}}}(s, S_{\mathsf{cl}}) = P_{\mathsf{cl}}. \end{align} $$

In addition, observe that differentiating (2.2) with respect to s and S, we obtain

(2.6) $$ \begin{align} \begin{aligned} \frac{\partial S_{\mathsf{cl}}}{\partial s} + \frac{\partial S_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial s} + \ell'(P_{\mathrm{cl}})\left(\frac{\partial P_{\mathsf{cl}}}{\partial s} + \frac{\partial P_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial s}\right) &= 0 \quad \text{and}\\ \frac{\partial S_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial S} + \ell'(P_{\mathrm{cl}}) \frac{\partial P_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial S} &= 1. \end{aligned} \end{align} $$

Combining (2.5) and (2.6), we have

$$ \begin{align*} \begin{aligned} \frac{\partial H}{\partial s}(s, S) &= \frac{\partial }{\partial s}H_{\mathsf{cl}}(s, S_{\mathsf{cl}}(s, p^*)) + \frac{\partial }{\partial s}\int_{0}^{P_{\mathsf{cl}}(s,p^*)} q \,\ell'(q) dq\\ &= -p + P_{\mathsf{cl}}\bigg[ \frac{\partial S_{\mathsf{cl}}}{\partial s} + \frac{\partial S_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial s}\bigg] + P_{\mathsf{cl}} \ell'(P_{\mathrm{cl}})\bigg[ \frac{\partial P_{\mathsf{cl}}}{\partial s} + \frac{\partial P_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial s}\bigg]\\ &= -p + P_{\mathsf{cl}}\bigg[ \frac{\partial S_{\mathsf{cl}}}{\partial s} + \frac{\partial S_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial s} + \ell'(P_{\mathrm{cl}})\bigg(\frac{\partial P_{\mathsf{cl}}}{\partial s} + \frac{\partial P_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial s}\bigg)\bigg]= -p. \end{aligned} \end{align*} $$

Analogously,

$$ \begin{align*}\kern-28pt \begin{aligned} \frac{\partial H}{\partial S}(s, S) &= \frac{\partial }{\partial S}H_{\mathsf{cl}}(s, S_{\mathsf{cl}}(s, p^*)) + \frac{\partial }{\partial S}\int_{0}^{P_{\mathsf{cl}}(s,p^*)} q \,\ell'(q) dq\\ &= P_{\mathsf{cl}}\frac{\partial S_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial S} + P_{\mathsf{cl}} \ell'(P_{\mathrm{ cl}})\frac{\partial P_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial S}= P_{\mathsf{cl}}\bigg[\frac{\partial S_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial S} + \ell'(P_{\mathrm{cl}})\frac{\partial P_{\mathsf{cl}}}{\partial p}\frac{\partial p^*}{\partial S}\bigg]= P_{\mathsf{cl}}= P. \end{aligned}\\[-45pt] \end{align*} $$

Proposition above allows for the description of periodic pensive billiard trajectories as critical points of some action. Indeed, for N points $s_1, ..., s_N$ on $\gamma $ , suppose that there exist $p_1, ..., p_N$ such that a billiard is locally transitive at $(s_i, p_i)$ and $S(s_i, p_i) = s_{i + 1}$ . Then locally around $(s_i, s_{i+1})$ we can define generating functions $H_i$ . With that, define (with the convention $s_{N + 1} = s_1$ )

$$ \begin{align*} H^N(s_1, ..., s_N) := \sum_{i = 1}^N H_i(s_i, s_{i+1}). \end{align*} $$

Corollary 2.11. A sequence of points $\{\gamma (s_1), \gamma (s_2), \dots , \gamma (s_N)\}$ constitutes a periodic trajectory if and only if

$$ \begin{align*} \frac{\partial H^N}{\partial s_i} = 0 \end{align*} $$

for $i = 1, ..., N$ . In this case, cosines of angles of incidence at points $\gamma (s_i)$ are given by

$$ \begin{align*} p_i = -\frac{\partial H_i}{\partial s_i}(s_i, s_{i +1}). \end{align*} $$

2.5 Twist property

The next natural question is under which conditions pensive billiard is a twist map.

Definition 2.12 (see, e.g., [Reference Katok and Hasselblatt17, Definition 9.3.1]).

A diffeomorphism $f:S^1 \times (0,1)\rightarrow S^1 \times (0,1)$ is called a right (respectively, left) twist map, if its lift $F = (F_1, F_2)$ to the universal cover satisfies

$$ \begin{align*} \partial_2 F_1> 0 \quad (\textit{respectively,}\ \partial_2 F_1 < 0). \end{align*} $$

Recall that by definition for the pensive billiard we have $ S = S_{\mathsf {cl}} + \tilde {\ell }(\Theta _{\mathsf {cl}})$ so that

$$ \begin{align*} \frac{\partial S}{\partial \theta} = \frac{\partial S_{\mathsf{cl}}}{\partial \theta} + {\tilde{\ell}}'(\Theta_{\mathsf{cl}})\frac{\partial \Theta_{\mathsf{cl}}}{\partial \theta}. \end{align*} $$

For the classical billiards in smooth convex domains we have (see [Reference Katok, Strelcyn, Ledrappier and Przytycki18, Part V, Theorem 4.2])

$$ \begin{align*} \nabla \mathsf{CB}(s, \theta) = \begin{pmatrix} \frac{\partial S_{\mathsf{cl}}}{\partial s} & \frac{\partial S_{\mathsf{cl}}}{\partial \theta}\\[1ex] \frac{\partial \Theta_{\mathsf{cl}}}{\partial s} & \frac{\partial \Theta_{\mathsf{cl}}}{\partial \theta} \end{pmatrix} = \begin{pmatrix} \frac{\kappa(s)d(s, S_{\mathsf{cl}}) - \sin\theta}{\sin \Theta_{\mathsf{cl}}} & \frac{d(s, S_{\mathsf{cl}})}{\sin \Theta_{\mathsf{cl}}}\\[1ex] \frac{\kappa(S_{\mathsf{cl}})\kappa(s)d(s, S_{\mathsf{cl}}) - \kappa(S_{\mathsf{cl}})\sin\theta - \kappa(s)\sin \Theta_{\mathsf{cl}}}{\sin \Theta_{\mathsf{cl}}} & \frac{\kappa(S_{\mathsf{cl}})d(s, S_{\mathsf{cl}}) - \sin \Theta_{\mathsf{cl}}}{\sin \Theta_{\mathsf{cl}}} \end{pmatrix}, \end{align*} $$

where $d(s, S_{\mathsf {cl}})$ is the Euclidean distance between points $\gamma (s)$ and $\gamma (S_{\mathsf {cl}})$ , and $\kappa (s)$ is the curvature of $\gamma $ at $\gamma (s)$ . Consequently, the right (left) twist condition for the pensive billiard map becomes

$$ \begin{align*} 0 \overset{(>)}{<} \frac{\partial S}{\partial \theta} = \frac{\partial S_{\mathsf{cl}}}{\partial \theta} + \tilde{\ell}'(\Theta_{\mathsf{cl}})\frac{\partial \Theta_{\mathsf{cl}}}{\partial \theta} = \frac{d(s, S_{\mathsf{cl}})}{\sin \Theta_{\mathsf{cl}}} + \tilde{\ell}'(\Theta_{\mathsf{cl}})\frac{\kappa(S_{\mathsf{cl}})d(s, S_{\mathsf{cl}}) - \sin \Theta_{\mathsf{cl}}}{\sin \Theta_{\mathsf{cl}}}, \end{align*} $$

or equivalently, it can be reformulated as follows.

Proposition 2.13. If for all $(s,\theta )\in S^1 \times (0,\pi )$

(2.7) $$ \begin{align} d(s, S_{\mathsf{cl}}) + \tilde{\ell}'(\Theta_{\mathsf{cl}})(\kappa(S_{\mathsf{cl}})d(s, S_{\mathsf{cl}}) - \sin \Theta_{\mathsf{cl}}) \overset{(<)}{>} 0, \end{align} $$

then the pensive billiard map is a right (left) twist map.

We remark that the twist property is by no means automatic. The following counterexample shows that for any given nonconstant delay function $\ell $ , the pensive billiard on a sufficiently thin ellipse is not a twist map.

Example 2.14. Let D be an ellipse with axis lengths $\varepsilon $ and $1/\varepsilon $ , where $\varepsilon $ is sufficiently small (see Figure 5). Let $\theta _0$ be such that $\tilde {\ell }'(\theta _0) \neq 0$ . We claim that $\partial S /\partial \theta $ changes sign and hence the corresponding pensive billiard is not a twist.

Figure 5 Pensive billiard map is not a twist map on a sufficiently thin ellipse. In the figure, $\tilde {\ell }'(\theta _0)> 0$ .

Indeed, compute $\partial S /\partial \theta $ at points $(s_{\mathrm {min/ max}}, \theta _{\mathrm {min/ max}})$ from which the classical billiard hits at angle $\theta _0$ the ellipse vertices, that is, the points $\gamma (a)$ and $\gamma (b)$ of minimal and maximal curvature, respectively. From the proposition above (as well as a straightforward computation), these derivatives are given by

$$ \begin{align*} \begin{aligned} \frac{\partial S}{\partial \theta}\big(s_{\mathrm{min}}, \theta_{\mathrm{min}}\big) &= \frac{d(s_{\mathrm{min}}, a) + \tilde{\ell}'(\theta_0)(\kappa(a)d(s_{\mathrm{min}}, a) - \sin \theta_0)}{\sin\theta_0} \\ &= \frac{1}{\sin\theta_0}\left(\frac{2\varepsilon\sin\theta_0}{\sin^2\theta_0 + \varepsilon^4\cos^2\theta_0} + \tilde{\ell}'(\theta_0)\bigg(\varepsilon^3\frac{2\varepsilon\sin\theta_0}{\sin^2\theta_0 + \varepsilon^4\cos^2\theta_0} -\sin\theta_0\bigg)\right)\\ &= -\tilde{\ell}'(\theta_0) + O(\varepsilon), \end{aligned} \end{align*} $$

and

$$ \begin{align*} \begin{aligned} \frac{\partial S}{\partial \theta}\big(s_{\mathrm{max}}, \theta_{\mathrm{max}}\big) &= \frac{d(s_{\mathrm{max}}, b) + \tilde{\ell}'(\theta_0)(\kappa(b)d(s_{\mathrm{max}}, b) - \sin \theta_0)}{\sin\theta_0} \\ &= \frac{1}{\sin\theta_0}\left(\frac{2\varepsilon^3\sin\theta_0}{\cos^2\theta_0 + \varepsilon^4\sin^2\theta_0} + \tilde{\ell}'(\theta_0)\bigg(\frac{1}{\varepsilon^3}\frac{2\varepsilon^3\sin\theta_0}{\cos^2\theta_0 + \varepsilon^4\sin^2\theta_0} -\sin\theta_0\bigg)\right)\\ &=\tilde{ \ell}'(\theta_0)\bigg(\frac{2}{\cos^2\theta_0} - 1\bigg) + O(\varepsilon). \end{aligned} \end{align*} $$

Consequently, for sufficiently small $\varepsilon $ , $\partial S /\partial \theta $ changes sign, so the twist condition is not satisfied.

This calculation was done with the assumption of a fixed delay function and varying domain, so it does not immediately apply to vortex billiards (for which the delay function is predetermined by the domain’s perimeter). However, the same counterexample works in the case of vortex billiard as well, since for any fixed $\theta _0$ , the derivative $\tilde {\ell }'(\theta _0) = L \frac {\sin \theta _0}{(1 + \cos ^2\theta _0)^{3/2}} \sim \varepsilon ^{-1}$ still dominates $O(\varepsilon )$ term for small $\varepsilon $ .

Example 2.15. For a disk of radius R the condition is simplified to

$$ \begin{align*} 2R\sin \Theta_{\mathsf{cl}} + \tilde{\ell}'(\Theta_{\mathsf{cl}}) \left(\frac{1}{R}\cdot 2R\sin\Theta_{\mathsf{cl}} - \sin\Theta_{\mathsf{cl}}\right) \overset{(<)}{>} 0 \quad\Longleftrightarrow\quad \tilde{\ell}' \overset{(<)}{>} - 2R. \end{align*} $$

For the puck billiard’s delay function $\tilde {\ell }(\theta ) = h \cot \theta $ , the right twist condition is never satisfied unless $h = 0$ , while the left twist condition becomes $ h> 2R. $ Heuristically, this means that for sufficiently tall puck, as the angle of incidence grows, the decrease of the shift along the boundary dominates the increase of the arc-length coordinate due to the classical billiard motion.

These examples motivate the following theorem, which says that a large class of pensive billiards on domains with controlled curvature are twist maps:

Proof. From a comparison argument in ordinary differential equations that can be found, for example, in [Reference Blaschke5, Part 4, §24], condition $1/R \leq \kappa \leq 1/r$ implies that $\partial D$ is contained in the region between two circles of radii r and R tangent to it at $S_{\mathsf {cl}}$ (see illustration in Figure 6).

Figure 6 A curve with curvature bounded between $\frac {1}{R}$ and $\frac {1}{r}$ lies between circles of radii r and R tangent to it.

This, in turn, implies that

$$ \begin{align*} 2r\sin\Theta_{\mathsf{cl}} \leq d(s, S_{\mathsf{cl}}) \leq 2R\sin\Theta_{\mathsf{cl}}\quad \Rightarrow \quad 2\frac{r}{R}\sin\Theta_{\mathsf{cl}} \leq \kappa(S_{\mathsf{cl}}) d(s, S_{\mathsf{cl}}) \leq 2\frac{R}{r}\sin\Theta_{\mathsf{cl}}. \end{align*} $$

Consequently, if $\tilde {\ell }'> -\frac {2r}{2({R}/{r})-1}$

$$ \begin{align*} d(s, S_{\mathsf{cl}}) + \tilde{\ell}'(\Theta_{\mathsf{cl}})(\kappa(S_{\mathsf{cl}})d(s, S_{\mathsf{cl}}) - \sin \Theta_{\mathsf{cl}})> 2r\sin\Theta_{\mathsf{cl}} - \frac{2r}{2\frac{R}{r}-1}\cdot(2\frac{R}{r}-1)\sin\Theta_{\mathsf{cl}} = 0, \end{align*} $$

so (2.7) implies that $\mathsf {PB}$ is a right twist map. Analogously, if $\tilde {\ell }' < -\frac {2R}{2\frac {r}{R}-1}$

$$ \begin{align*} d(s, S_{\mathsf{cl}}) + \tilde{\ell}'(\Theta_{\mathsf{cl}})(\kappa(S_{\mathsf{cl}})d(s, S_{\mathsf{cl}}) - \sin \Theta_{\mathsf{cl}}) < 2R\sin\Theta_{\mathsf{cl}} - \frac{2R}{2\frac{r}{R}-1}\cdot(2\frac{r}{R}-1)\sin\Theta_{\mathsf{cl}} = 0, \end{align*} $$

so $\mathsf {PB}$ is a left twist map.

Now we verify these conditions on $\tilde {\ell }'$ for the puck and vortex billiards.

For the puck billiard with delay function $\tilde {\ell }(\theta ) = h\cot \theta $ we have $ \tilde {\ell }'(\theta ) = -\frac {h}{\sin ^2\theta } \in (-\infty , -h), $ so if $h> \frac {2R}{2({r}/{R})-1}$ , the condition of Theorem 2.16 on $\tilde {\ell }'$ is satisfied, and the puck billiard is a left twist map. For vortex billiard one has $\tilde {\ell }(\theta ) = L\big (1 - \frac {\cos \theta }{\sqrt {1 + \cos ^2\theta }}\big )$ , which implies that $\tilde {\ell }'(\theta ) = L\frac {\sin \theta }{(1 + \cos ^2\theta )^{3/2}}> 0 > -\frac {2r}{2({R}/{r})-1}$ and, therefore, the vortex billiard map is a right twist map.

Many other implications of the twist property can be found in [Reference Katok and Hasselblatt17, Chapter 9.3, 13]. As a sample of such applications we specialize the Birkhoff theorem on periodic orbits to pensive billiards.

2.6 Examples of dynamical behaviors

We now give three examples where dynamical properties of the pensive billiard system can generally be derived from the classical results.

2.6.1 Existence of periodic orbits

Recall that a Birkhoff periodic orbit of type $(p,q)$ is a periodic orbit with period q and rotation number $p/q$ , along with some additional properties – see Definition 9.3.6 in [Reference Katok and Hasselblatt17]. We have

Theorem 2.18. Let D be a smooth convex billiard table with the boundary curvature between $1/R$ and $1/r$ for some $0 < R/2 < r < R$ . Then

• • there exists $h_0 = h_0(r,R)$ such that for any $h> h_0$ , and any $p/q \in \mathbb {Q}$ with p and q relatively prime, puck billiard of height h in D has at least two Birkhoff periodic orbits of type $(p,q)$ .

• • for any $p/q \in \big (\frac {1}{2} - \frac {1}{2\sqrt {2}} ,\frac {3}{2} +\frac {1}{2\sqrt {2}} \big ) \cap \mathbb {Q}$ with p and q relatively prime, the vortex billiard in D has at least two Birkhoff periodic orbits of type $(p,q)$ .

Proof. In view of Theorem 2.16, the proof is an application of the Birkhoff theorem: such orbits exist on twist intervals. Hence one needs to find the corresponding twist intervals for the puck and vortex billiards. For the puck billiard its twist interval is

$$ \begin{align*}(\lim_{\theta\rightarrow \pi}S, \lim_{\theta\rightarrow 0}S) = (1 + \frac{1}{2L}\lim_{\theta\rightarrow \pi}\tilde{\ell}(\theta), \frac{1}{2L}\lim_{\theta\rightarrow 0}\tilde{\ell}(\theta)) = (-\infty, \infty), \end{align*} $$

while for the vortex billiard its twist interval is

$$ \begin{align*} (\lim_{\theta\rightarrow 0}S, \lim_{\theta\rightarrow \pi}S) = (\frac{1}{2L}\lim_{\theta\rightarrow 0}\tilde{\ell}(\theta), 1 + \frac{1}{2L}\lim_{\theta\rightarrow \pi}\tilde{\ell}(\theta)) = ( \frac{1}{2}(1 - \frac{1}{\sqrt{2}}), 1 + \frac{1}{2}(1 + \frac{1}{\sqrt{2}})), \end{align*} $$

and the statement follows.

2.6.2 Caustics on the disk

As the simplest example, let us analyze the case of $D = B_1(0)$ . The simplification comes from the fact that the angle $\theta $ at which the point leaves the boundary is the same as the angle at which it hits next. Consequently, the phase space of any pensive billiard in the disk splits into invariant sets $M_\theta = \gamma \times \{\theta \}$ , and the restriction of $\mathsf {PB}$ to each $M_\theta $ is a circle rotation by angle

(2.8) $$ \begin{align} \tilde{\theta} = 2\theta + \tilde{\ell}(\theta). \end{align} $$

Hence any result for a standard billiard in a disk can be adapted to the pensive billiard case. The induced transformation on this invariant circle is a rotation (depending also on the delay function $\tilde {\ell }(\theta )$ ) through a fixed angle. As a consequence, we have, for example:

Figure 7 Left: plot of realization of a pensive billiard on a polygon as an interval exchange transformation. Right: interpretation of the plot. The red vector corresponds to the red point on the plot; the blue vector is its image under $\mathsf {PB}$ . Here $\tilde {\ell }(\theta )= \cot \theta $ and $\theta _0=\pi /9$ .

Proposition 2.19. Any nondegenerate pensive billiard in the disk has concentric circle caustics.

We will return to this proposition in the context of vortex billiards in Section 4.4.

3.1 Length of geodesics on the side of a cylinder

It turns out, the potential V defined by equation 2.1 has a natural meaning of the length of “curvilinear hypotenuse” in the case of a generalized puck billiard. Indeed, normalize the height of the cylinder: $h=1$ . Assume that a geodesic $\{\lambda (t)\}$ enters the side N with the unit speed $v_0 = p\partial _s + \sqrt {1 - p^2}\partial _y$ , where, as usual, $p=\cos \theta $ , and it exits at the bottom after time T.

Theorem 3.2. For a cylinder with metric $ f(y)\,ds^2 + dy^2$ geodesic segment $\{\lambda (t)\}$ with initial velocity $v_0 = p\partial _s + \sqrt {1 - p^2}\partial _y$ has length

$$ \begin{align*} \mathcal{V}(p) \big|^T_0 := \int_0^1 \frac{1}{\sqrt{1 - \frac{p^2}{f(y)}}} dy , \end{align*} $$

and it has horizontal shift (along s-coordinate)

(3.1) $$ \begin{align} \ell (p)\big|^T_0 = \int_0^1 \frac{p}{f(y)}\frac{1}{\sqrt{1 - \frac{p^2}{f(y)}}}dy. \end{align} $$

In particular, $\mathcal {V}(p)$ is the potential for the delay function $\ell (p)$ : $\mathcal {V}(p) = \int _0^{p}q\,\ell '({q})\,d{q}$ .

Proof. We observe that the corresponding Lagrangian $L(s, y, \dot {s}, \dot {y}, t) = \sqrt {f(y)\dot {s}^2 + \dot {y}^2}$ is invariant under shifts in variables t and s. Hence from Noether’s theorem $E = f(y)\dot {s}^2 + \dot {y}^2$ and $P = f(y)\dot {s}$ are conserved along the trajectory. Consequently, its length is given by

$$ \begin{align*} \mathcal{V}(p)\big|^T_0 = \int_0^T |\dot{\lambda}|dt = \int_0^1 \sqrt{E} \bigg(\frac{dy}{dt}\bigg)^{-1}dy = \int_0^1 \frac{\sqrt{E}}{\sqrt{E - \frac{P^2}{f(y)}}} dy. \end{align*} $$

Similarly, the horizontal shift is given by

$$ \begin{align*} \ell(p)\big|^T_0 = \int_0^T \frac{ds}{dt} dt = \int_0^1 \frac{ds}{dt} \bigg(\frac{dy}{dt}\bigg)^{-1} dy = \int_0^1 \frac{P}{f(y)}\frac{1}{\sqrt{E - \frac{P^2}{f(y)}}}dy. \end{align*} $$

Evaluating E and P at time $t = 0$ we get $ E = 1$ and $ P = p$ , so that (3.1) holds.

To verify the relation between $\mathcal {V}$ and $\ell $ we compute

$$ \begin{align*} \int_0^p q\ell'(q) \,dq &= p \ell(p) - \int_0^p \ell(q)\, dq \\ &= \int_0^1\left( \frac{p^2}{f(y)}\frac{1}{\sqrt{1 - \frac{p^2}{f(y)}}} + \sqrt{1 - \frac{p^2}{f(y)}}\right) dy =\int_0^1 \frac{1}{\sqrt{1 - \frac{p^2}{f(y)}}} dy = \mathcal{V}(p), \end{align*} $$

so $\mathcal V$ is the potential for the delay function $\ell $ .

Corollary 3.3. Suppose a pensive billiard with delay function $\ell $ is given. If there exists a function $f: [0,1] \to [1, \infty )$ satisfying the boundary conditions $f(0)= f(1) $ and the relation

$$ \begin{align*} \ell(p) = \int_0^1 \frac{p}{f(y)}\frac{1}{\sqrt{1 - \frac{p^2}{f(y)}}}dy\, \end{align*} $$

then trajectories of the pensive billiard in D naturally correspond to geodesics on the closed cylinder $(D\times \{0\}) \cup (\partial D \times [0,1]) \cup (D\times \{1\})$ with the flat metric on the top and bottom and the metric on the sides given by $f(y)ds^2 + dy^2$ .

3.2 From the puck billiards to (almost) classical ones

Assume that the function $f(y)$ for $y\in [0, h]$ has the symmetry property: $f(y)= f(h-y)$ . Then instead of defining a delay function, one can consider the classical billiard on the “half-puck,” simply reflecting geodesic according to the standard billiard law in the curve $y=h/2$ . Such a geodesic will reappear at the top copy of D with the same delay $\ell (p)$ , as in the regular puck billiard.

This suggests the following consideration: extend the flat domain D beyond its boundary $\gamma =\partial D$ to a new domain $D_{h/2}$ by a band of width $h/2$ with a special metric $f(y)ds^2 + dy^2$ for $y\in [0, h/2]$ and consider the classical billiard in this, generally speaking, nonflat domain $D_{h/2}$ . Note that geodesics in $D_{h/2}$ refract at the points of $\gamma $ , where the flat metric of D changes to the nonflat cylindrical metric.

This point of view can be particularly interesting for small h. In particular, it would be interesting to study persistence of caustics or integrability in this enlargement of D for special domains.

4.1 Helmholtz–Kirchhoff point vortex system

We recall the classical derivation of the point vortex system. On simply connected planar domains (see [Reference Drivas and Elgindi8, §2.2] and [Reference Grotta-Ragazzo, Gustafsson and Koiller14] for a discussion of the general case), the two-dimensional Euler equation can be written as a closed evolution for the vorticity measure $\omega =\nabla ^\perp \cdot u$ :

$$ \begin{align*} \partial_t \omega + u\cdot\nabla\omega = 0, \\ u = \nabla^{\perp} \Delta^{-1}\omega, \end{align*} $$

which manifests in transport of vorticity $\omega $ by a velocity field u generated from $\omega $ through Biot-Savart law. In the point vortex model, we assume that the initial vorticity is given by a finite linear combination of Dirac $\delta $ -functions $\omega = \sum _{i=1}^{N} \Gamma _i\delta _{\mathbf {z}_i(t)}$ located at points $\mathbf {z}_i = (x_i, y_i)$ in a domain $D \subset \mathbb {R}^2$ and carrying circulations $\Gamma _i\in \mathbb {R}$ . Formally plugging $\omega $ into the Biot–Savart law to obtain velocity, we find

$$ \begin{align*} u(\cdot) = \nabla^{\perp}\Delta^{-1}\omega(\cdot) = \nabla^{\perp}\int_{D}G_D(\cdot ,\mathbf{w})\sum_{i=1}^{N} \Gamma_i\delta_{\mathbf{z}_i(t)}(\mathbf{w}) d\mathbf{w} = \nabla^{\perp}\sum_{i=1}^{N} \Gamma_iG_D({\mathbf{z}_i(t)}, \cdot), \end{align*} $$

where $G_D$ is Green’s function for the Dirichlet Laplacian on D. Consequently, dynamics of positions of point vortices obey the following system:

$$ \begin{align*} \dot{\mathbf{z}}_i(t) = \nabla^{\perp}\sum_{j=1}^{N} \Gamma_jG_D(\mathbf{z}_i(t),\mathbf{z}_j(t)), \quad \ i=1, \dots,N. \end{align*} $$

Note that this calculation is formal, as the $i = j$ term in the preceding sum is infinite. However, the contribution of self-interaction can be rigorously computed by observing that the singularity is radially symmetric, and the quantity

$$ \begin{align*} R_D({\mathbf{z}}) := \lim_{\mathbf{w} \to \mathbf{z}} \Big[ G_D(\mathbf{z}, \mathbf{w}) + \frac{\log(|\mathbf{z}- \mathbf{w}|)}{2 \pi}\Big] \end{align*} $$

is bounded (see, e.g., [Reference Marchioro and Pulvirenti22], [Reference Drivas, Glukhovskiy and Khesin9]). Consequently for any radially symmetric desingularization of the vortex, logarithmic singularity does not contribute to the motion, and the equations can be replaced by the regularized version:

$$ \begin{align*} \dot{\mathbf{z}}_i(t) = \nabla^{\perp}\sum_{j\neq i} \Gamma_jG_D(\mathbf{z}_i(t),\mathbf{z}_j(t)) + \nabla^{\perp}\Gamma_iR_D(\mathbf{z}_i), \quad \ i=1, \dots,N. \end{align*} $$

Kirchhoff observed that this system can be written in the canonical Hamiltonian form. Namely, the dynamics of N point vortices is governed by the following Hamiltonian system:

$$ \begin{align*} \Gamma_i \dot{x_i} = \frac{\partial H}{\partial y_i},\qquad \Gamma_i \dot{y_i} = -\frac{\partial H}{\partial x_i}, \qquad \text{for each} \ i=1, \dots,N, \end{align*} $$

where the Hamiltonian function $H:D^N\setminus \mathbf {D} \to \mathbb {R}$ is

$$ \begin{align*}H(\mathbf{z}_1,...,\mathbf{z}_N) = \sum_{1\leq i<j\leq N}\Gamma_i\Gamma_j G_D (\mathbf{z}_i, \mathbf{z}_j) + \frac{1}{2}\sum_{i=1}^{N}\Gamma_i^2R_D(\mathbf{z}_i) \end{align*} $$

(where $ \mathbf {D}$ consists of all the diagonals $\mathbf {z}_i=\mathbf {z}_j$ ), and the symplectic structure on $D^N \subset \mathbb {R}^{2N}$ is given by $\sum _{i=1}^N \Gamma _i\, dx_i\wedge dy_i$ .

We note that the motion of even a small number of point vortices in the presence of boundary is far from being completely analyzed. The complete analysis for one vortex in general domains was carried out by Gustafsson in [Reference Gustafsson15]. Details of motion of several vortices in special domains can be found in [Reference Kimura20, Reference Yang30, Reference Modin and Viviani24, Reference Khesin and Wang19] and many other works. In particular, in the special case of a half-plane and two equal vortices, $\Gamma _1=\Gamma _2$ , an interesting phenomenon of leapfrogging occurs [Reference Love21, Reference Mavroyiakoumou and Berkshire23, Reference Newton25, Reference Khesin and Wang19].

4.2 Derivation of vortex billiard system

The vortex billiard is derived as a limit of zero separation of a vortex dipole evolving as Kirchhoff point vortices. Specifically, consider two vortices with circulations $\varepsilon \Gamma $ and $-\varepsilon \Gamma $ for $\Gamma \in \mathbb {R}$ , initially located at

$$ \begin{align*} \mathbf{z}^\varepsilon_\pm := ({x}_0 \pm \varepsilon\cos\theta, {y}_0 \pm \varepsilon \sin \theta ), \end{align*} $$

for some $(x_0,y_0)\in D$ , provided $2\varepsilon < \mathrm {dist}((x_0,y_0),\partial D)$ . We separate the motion into three regions treated in the following steps.

Step 1: (far from boundary) We claim that away from the boundary singular dipole moves in a straight line, with the relative distance between vortices being constant to first order. To be more precise, let $\delta = \varepsilon ^h$ for some $1/2<h<1$ , and let

$$ \begin{align*} D^{in}_\delta = \{\mathbf{z} \in D \ |\ \text{dist}(\mathbf{z}, \partial D)> \delta\}. \end{align*} $$

We observe that restricted to $ D^{in}_\delta \times D^{in}_\delta $ , Green’s and Robin’s functions satisfy

$$ \begin{align*} \begin{aligned} \nabla G(\mathbf{z}, \mathbf{w})&= -\frac{1}{2 \pi} \frac{\mathbf{z}-\mathbf{w}}{\|\mathbf{z}-\mathbf{w}\|^2} + o\bigg(\frac{1}{\|\mathbf{z}-\mathbf{w}\|}\bigg),\qquad \nabla R(\mathbf{z}) = o\bigg(\frac{1}{\|\mathbf{z}-\mathbf{w}\|}\bigg), \end{aligned} \end{align*} $$

as $\|\mathbf {z}-\mathbf {w}\|\to 0$ . Consequently, to the leading order the equations of motions are given by those in the plane, where the fact that the dipole moves along a straight line with constant speed is evident. See, for instance, [Reference Boatto and Koiller6, Reference Drivas, Glukhovskiy and Khesin9] for further details.

Step 2: (near the boundary). We claim that as the dipole approaches boundary, the interaction is captured by that in the half-plane, which will be described in §4.3. Specifically, we have

Proposition 4.1. If dipole traveling at speed v hits the boundary at angle $\theta $ , it splits into two vortices traveling along the boundary in the opposite directions with speeds

$$ \begin{align*} v_{\pm} = v\left(\sqrt{1 + \cos^2\theta} \mp \cos{\theta}\right). \end{align*} $$

Remark 4.2 (Silver ratio).

In the limiting case of zero angle, that is, $\cos \theta =1$ , normalizing $v=1$ , the speeds $v_{\pm }$ of the separated vortices are given by the silver ratio $\chi = 1+ \sqrt {2}$ , as

$$ \begin{align*} v_+ = \frac{1}{\chi}, \qquad v_- = \chi. \end{align*} $$

Alternatively the silver ratio can be observed in the distances to the boundary after splitting, as by Theorem 4.3 the speed of a single vortex is inversely proportional to the distance to the boundary.

Figure 9 Dipole hitting the boundary “at zero angle” splits into vortices traveling along the boundary at speeds (respectively, distances to the boundary) proportional (respectively, inversely proportional) to the silver ratio and inverse silver ratio.

Proof for Proposition 4.1 in the half-plane will be given in the next section, while in the meantime, we provide an argument for the validity of the half-plane approximation. First, by the discussion above, the vortex dipole Hamiltonian in $ D^{b}_\delta := D \setminus D^{in}_\delta $ is given by its expression on the whole plane, up to an $o(\varepsilon )$ correction. For $x\in M$ lying within $\delta $ -neighborhood of the boundary, the Laplacian can be expressed as

$$ \begin{align*} \Delta = \delta^2(\Delta_{\tilde{z},\tilde{s}} + \delta^2 R_\Delta), \end{align*} $$

where $\Delta _{\tilde {z},\tilde {s}} = \partial ^2_{\tilde {z}} + \partial ^2_{\tilde {s}}$ is the Laplacian in the coordinates

$$ \begin{align*} \begin{aligned} \tilde{z} &:= \delta\text{ dist}(x, \partial D)\\ \tilde{s} &:= \delta\text{ proj}_{\partial D} x, \end{aligned} \end{align*} $$

and $R_\Delta $ is a bounded correction, see [Reference Drivas, Iyer and Nguyen10]. Consequently, Green’s functions of $\Delta $ and $\Delta _{\tilde {z}, \tilde {s}}$ in $D^{b}_\delta \cap \{\tilde {s} < \delta \}$ coincide up to $O(\delta ^2) = o(\varepsilon )$ correction. Note that $\Delta _{\tilde {z}, \tilde {s}}$ is essentially the (rescaled) half-plane Laplacian, hence the conclusion.

Step 3: (along the boundary) As vortices leave the region $ D^{b}_\delta \cap \{\tilde {s} < \delta \}$ , interaction between them (equivalently, the Green’s function term in the Hamiltonian) becomes of lower order compared to the interactions with the boundary (Robin’s function term). Consequently, to the leading order the dynamics is given by the following theorem, proved in [Reference Flucher and Gustafsson11, Corollary 14]:

Theorem 4.3 [Reference Flucher and Gustafsson11].

Let $ D\subset \mathbb {R}^2$ be simply connected. A vortex with circulation $\Gamma $ at a small distance $\varepsilon $ from $\partial D$ stays within a distance $\varepsilon + O(\varepsilon ^2)$ , and travels along the boundary at speed

$$ \begin{align*} v = -\frac{\Gamma}{4\pi \varepsilon} + O(\varepsilon). \end{align*} $$

Since in our renormalization of dipole $\Gamma \sim \varepsilon $ , to the first order the speed of vortices stays as described in Proposition 4.1. In particular, the boundary curvature does not contribute to the vortex’ speed.

Finally, after traveling along the boundary in opposite directions, vortices meet again the distance $2L\frac {v_+}{v_+ + v_-} = L\frac {\sqrt {1 + \cos ^2\theta } - \cos \theta }{\sqrt {1 + \cos ^2\theta }}$ away from the point of splitting in the direction of motion of the positive vortex. The argument in Step 2 can be run backward to conclude that vortices will merge into a dipole and enter $ D^{in}_\delta $ , and then Steps 1-3 repeat.

Figure 10 Cartoon of time evolution of a vortex dipole hitting the boundary.

This defines the vortex billiard as the limit of the dipole motion. Next we describe this limit as a pensive billiard with a specific delay function.

4.3 Dipoles on the half-plane: Fission-fusion rules

Dipoles in a half-plane can have different types of motion, depending on the strength of their interaction as compared to the interaction with their mirrors (or equivalently, “with the boundary”).

In the half-plane, Green’s function has the following form:

$$ \begin{align*} G_{\mathbb{R}^2_+}(\mathbf{z}, \mathbf{w}) = -\frac{1}{2\pi}\log{|\mathbf{z} - \mathbf{w}|} + \frac{1}{2\pi}\log{|\mathbf{z} - \overline{\mathbf{w}}|}, \end{align*} $$

where $\overline {\mathbf {w}}$ is the mirror image of $\mathbf {w}$ . Consequently,

$$ \begin{align*} R_{\mathbb{R}^2_+}(\mathbf{z}) := \lim_{\mathbf{w} \to \mathbf{z}} \big[ G_{\mathbb{R}^2_+}(\mathbf{z}, \mathbf{w}) + \frac{\log(|\mathbf{z}- \mathbf{w}|)}{2 \pi}\big] = \lim_{\mathbf{w} \to \mathbf{z}} \big[\frac{1}{2\pi}\log{|\mathbf{z} - \overline{\mathbf{w}}|}\big] = \frac{1}{2\pi}\log{|\mathbf{z} - \overline{\mathbf{z}}|}. \end{align*} $$

Combining two preceding equations, Hamiltonian for two vortices in a half-plane is written as

$$ \begin{align*} H(\mathbf{z}_1,\mathbf{z}_2) = \frac{\Gamma_1\Gamma_2}{2\pi} \log{\frac{|\mathbf{z}_1 - \overline{\mathbf{z}}_2|}{|\mathbf{z}_1 - \mathbf{z}_2|}} + \frac{\Gamma_1^2}{4\pi}\log{|\mathbf{z}_1 - \overline{\mathbf{z}}_1|} +\frac{\Gamma_2^2}{4\pi}\log{|\mathbf{z}_2 - \overline{\mathbf{z}}_2|}. \end{align*} $$

For a dipole with $\Gamma _1 = -\Gamma _2 = \Gamma $ , we denote for convenience

$$ \begin{align*} \begin{aligned} x_1 &= x_+ \quad x_2 = x_-,\\ y_1 &= y_+ \quad y_2 = y_-, \end{aligned} \end{align*} $$

to simplify the expression for the Hamiltonian. Namely, introducing relative coordinates

$$ \begin{align*} (x_{\mathsf{rel}}, y_{\mathsf{abs}}) := \left(x_+ - x_-, \frac{y_+ + y_-}{2}\right), \end{align*} $$

the Hamiltonian becomes

$$ \begin{align*} H(x_{\mathsf{rel}},y_{\mathsf{abs}}) = \frac{-\Gamma^2}{4\pi}\log\bigg[\frac{1}{\mu^2 + x_{\mathsf{rel}}^2} + \frac{1}{4y_{\mathsf{abs}}^2 - \mu^2}\bigg] \end{align*} $$

where

(4.1) $$ \begin{align} \mu := y_+ - y_- \end{align} $$

is a conserved “momentum” corresponding to the translational invariance.

Figure 11 contains different types of motion for a pair of vortices with the opposite circulations in a nutshell, as the interaction parameter $\mu $ changes.

Figure 11 Types of motion for two vortices of opposite circulations in the half-plane (bottom of each panel corresponds to the boundary). Top: As the parameter $\mu $ defined in (4.1) increases through the silver ratio, the dipole stops being formed and vortices start passing each other. Bottom: As the parameter $\mu $ increases beyond $\mu ^* = 1 + 2\phi + 2\sqrt {1 + 2\phi }\approx 8.35$ , the reverse motion ceases and vortices pass each other smoothly, see [Reference Khesin and Wang19].

We are interested in the strongest interaction when vortices form a dipole, which is depicted in the upper-left panel.

4.3.1 Fission Rule

If a vortex dipole “hits” the boundary, it splits into positive and negative vortices, traveling with different speed along the boundary in the opposite directions.

Proposition 4.4. In the half-plane, consider a dipole $\mathbf {z}^\varepsilon _+, \mathbf {z}^\varepsilon _-$ with circulations $\varepsilon \Gamma :=\varepsilon \Gamma _+ = - \varepsilon \Gamma _-$ , initially located at $(x_0 \pm \varepsilon \sin \theta , y_0 \mp \varepsilon \cos \theta )$ . Set $t_* =4\pi y_0/(\Gamma \sin \theta )$ to be the “time of hitting the boundary.” Then the motion of the dipole vortices converges to the following trajectories $\widetilde {\mathbf {z}}^\varepsilon _+(t), \widetilde {\mathbf {z}}^\varepsilon _-(t)$ in the sense that for any $t \geq 0$ , $|\mathbf {z}^\varepsilon _\pm (t)- \widetilde {\mathbf {z}}^\varepsilon _\pm (t)|\to 0$ as $\varepsilon \rightarrow 0$ : Prior to hitting the boundary (for $t< t_*$ ):

• • $\widetilde {\mathbf {z}}^\varepsilon _i(t) = (x_0 \pm \varepsilon \sin \theta , y_0 \mp \varepsilon \cos \theta ) + \frac {\Gamma }{4\pi }t(-\cos \theta , -\sin \theta )$ .

After hitting the boundary (for $t>t_*$ ):

• • $\widetilde {\mathbf {z}}^\varepsilon _i(t) = \Big ((t-t_*)\frac {\pm \Gamma }{4\pi \left (\sqrt {1 + \cos ^2\theta } \mp \cos \theta \right )}, \varepsilon \big (\sqrt {1 + \cos ^2\theta } \mp \cos \theta \big )\Big )$ .

Proof. From the conservation of Hamiltonian, we can express the equation of the trajectory in $(y_{\mathsf {abs}}, x_{\mathsf {rel}})$ –variables:

$$ \begin{align*} \frac{1}{\varepsilon^2\cos^2\theta + (x_{\mathsf{rel}}/2)^2} + \frac{1}{y_{\mathsf{abs}}^2 - \varepsilon^2\cos^2\theta} = \frac{1}{\varepsilon^2} + \frac{1}{y_0^2 - \varepsilon^2\cos^2{\theta}}. \end{align*} $$

Let $ 3/4< h< 1$ be fixed. Expressing $y_{\mathsf {abs}}$ as a function of $x_{\mathsf {rel}}$ , we see that for any $x_{\mathsf {rel}}> \varepsilon ^h$

$$ \begin{align*} \begin{aligned} y_{\mathsf{abs}} &= \sqrt{\frac{1}{\frac{1}{\varepsilon^2} + \frac{1}{y_0^2 - \varepsilon^2\cos^2{\theta}} - \frac{1}{\varepsilon^2\cos^2\theta + (x_{\mathsf{rel}}/2)^2}} + \varepsilon^2\cos^2\theta} = \sqrt{\varepsilon^2(1 + \cos^2\theta) + O(\varepsilon^{4-2h})}\\ & = \varepsilon\sqrt{1 + \cos^2\theta} + O(\varepsilon^{2 - h}), \end{aligned} \end{align*} $$

where implicit constant in big-O notation may depend on $\theta , y_{\mathsf {abs}},$ and h. Together with the conservation of momentum $\mu = y_+ - y_- = -2\varepsilon \cos \theta $ , we get that in $\{x_{\mathsf {rel}}> \varepsilon ^h\}$

$$ \begin{align*} y_\pm = y_{\mathsf{abs}} \pm \frac{\mu}{2} = \varepsilon(\sqrt{1 + \cos^2\theta} \mp \cos\theta) + O(\varepsilon^{2 - h}). \end{align*} $$

Analogously, for any $y_{\mathsf {abs}}> \varepsilon ^h$ , we have that $ x_{\mathsf {rel}} = 2\varepsilon \sin \theta + O(\varepsilon ^{2 - h}). $ This establishes the shape of the limiting trajectory.

To obtain the velocity, we also fix $0< h' < 1-h$ and notice that for sufficiently small $\varepsilon $ the whole trajectory lies in the region (see Figure 12)

$$ \begin{align*} \begin{aligned} A_x\cup A_O \cup A_y :=& \big\{|x_{\mathsf{rel}} - 2\varepsilon\sin(\theta)|< \varepsilon^{2-h-h'}, y_{\mathsf{abs}}> \varepsilon^h\big\}\cup\big\{0<x_{\mathsf{rel}} < \varepsilon^h,0<y_{\mathsf{abs}} < \varepsilon^h\big\} \cup \\ &\quad\quad\big\{x_{\mathsf{rel}}>\varepsilon^h,|y_{\mathsf{abs}} -\varepsilon\sqrt{1 + \cos^2\theta}| <\varepsilon^{2-h-h'}\big\}. \end{aligned} \end{align*} $$

Figure 12 Trajectory of a dipole in $(x_{\mathsf {rel}}, y_{\mathsf {abs}})$ coordinates for small $\varepsilon $ .

Recalling that evolution of $x_{\mathsf {rel}}, y_{\mathsf {abs}}$ reads

$$ \begin{align*} \dot{x}_{\mathsf{rel}} =\frac{\Gamma}{4\pi}\frac{\varepsilon}{\frac{1}{4\varepsilon^2} + \frac{1}{y_0^2 - 4\varepsilon^2\cos^2{\theta}}}\frac{8y_{\mathsf{abs}}}{(4y_{\mathsf{abs}}^2 - 4\varepsilon^2\cos^2\theta)^2}, \qquad \dot{y}_{\mathsf{abs}} =\frac{\Gamma}{4\pi}\frac{-\varepsilon}{\frac{1}{4\varepsilon^2} + \frac{1}{y_0^2 - 4\varepsilon^2\cos^2{\theta}}}\frac{2x_{\mathsf{rel}}}{(x_{\mathsf{rel}}^2 + 4\varepsilon^2\cos^2\theta)^2}. \end{align*} $$

we compute that inside $A_x, A_O$ , and $A_y$ we have, respectively,

$$ \begin{align*} \begin{cases} \dot{x}_{\mathsf{rel}} = O(\varepsilon^{3-3h}) \\ \dot{y}_{\mathsf{abs}} = -\frac{\Gamma}{4\pi}\sin\theta + O(\varepsilon^{1-h-h'}) \end{cases}\!\!\!\!\!\!,\begin{cases} \dot{x}_{\mathsf{rel}} \gtrsim \varepsilon^{3-3h} \\ \dot{y}_{\mathsf{abs}} < 0 \end{cases}\!\!\!\!\!\!,\mathrm{~and ~} \begin{cases} \dot{x}_{\mathsf{rel}} = \frac{\Gamma}{2\pi}\sqrt{1 + \cos^2\theta} + O(\varepsilon^{1-h-h'}) \\ \dot{y}_{\mathsf{abs}} = O(\varepsilon^{3-3h}) \end{cases}\!\!\!\!\!\!. \end{align*} $$

From the first system we conclude that for $t < t_* - O(\varepsilon ^{1-h-h'})$ the trajectory stays in $A_x$ and satisfies $ y_{\mathsf {abs}}(t) = \hat {y} - \frac {\Gamma }{4\pi }t\sin \theta + O(\varepsilon ^{1-h-h'}). $

From the second system, we conclude that dipole spends at most time $\varepsilon ^{4h-3}$ in $A_O$ before entering $A_y$ at time $t_{**} = t_* + O(\varepsilon ^{4h-3})$ .

From the third system we conclude that for $t> t_{**}$ the trajectory stays in $A_y$ and satisfies $ x_{\mathsf {rel}}(t) = \frac {\Gamma }{2\pi }(t - t_{**})\sqrt {1 + \cos ^2\theta } + O(\varepsilon ^{1-h-h'}). $ With that, $x_{\mathrm {abs}}(t)$ can be recovered by integrating its evolution equation. Then, individual positions of the vortices can be recovered from absolute and relative coordinates algebraically. Sending $\varepsilon \to 0$ concludes the proof.

4.3.2 Fusion rule

The following proposition is a counterpart of the previous one, except that we consider merging instead of splitting. First we note that a single vortex of strength $\Gamma $ in the half-plane $y>0$ starting at the point $(x,y)$ moves along the boundary, keeping constant height y at the constant speed $v:= |{\Gamma }/{4\pi {y}}|$ in the direction depending on the sign of $\Gamma $ .

Proposition 4.5. Consider two point vortices $\mathbf {z}^\varepsilon _+, \mathbf {z}^\varepsilon _-$ with circulations $\varepsilon \Gamma :=\varepsilon \Gamma _+ = - \varepsilon \Gamma _-$ in the half-plane, initially located at $({x}_+, \varepsilon \hat {y}_+)$ , $({x}_-, \varepsilon \hat {y}_-)$ (with ${x}_+ < {x}_-$ ). Set $t_* = (x_- - x_+)/(\frac {\Gamma }{4\pi \hat {y}_+} + \frac {\Gamma }{4\pi \hat {y}_-})$ to be the “time of merging.” As $\varepsilon \rightarrow 0$ , the motion converges to the following trajectories $\widetilde {\mathbf {z}}^\varepsilon _+(t), \widetilde {\mathbf {z}}^\varepsilon _-(t)$ in the sense that for any $t \geq 0$ , $|\mathbf {z}^\varepsilon _\pm (t)- \widetilde {\mathbf {z}}^\varepsilon _\pm (t)|\to 0$ as $\varepsilon \rightarrow 0$

• • Before merging, for $t<t_*$ , the vortices move along the boundary:

  $$ \begin{align*}\widetilde{\mathbf{z}}^\varepsilon_\pm(t) = \Big(x_i \pm t\frac{\Gamma}{4\pi\hat{y}_\pm}, \varepsilon\hat{y}_\pm\Big).\end{align*} $$

• • If $\chi ^{-2} < \frac {\hat {y}_+}{\hat {y}_-} <\chi ^2 $ where $\chi $ is the silver ratio, for $t> t_*$ the vortices merge into a dipole:

  $$ \begin{align*}\widetilde{\mathbf{z}}^\varepsilon_i(t) = (x_* \mp\varepsilon\sqrt{\hat{y}_+\hat{y}_-} \sin\theta, \pm \varepsilon\sqrt{\hat{y}_+\hat{y}_-}\cos\theta) + \frac{\Gamma}{4\pi\sqrt{\hat{y}_+\hat{y}_-}}t(-\cos\theta, -\sin\theta),\end{align*} $$
  where $x_* = \frac {x_+\hat {y}_+ + x_-\hat {y}_-}{\hat {y}_+ +\hat {y}_-}$ and $\theta = \tan ^{-1}\sqrt {\frac {(\hat {y}_+ - \hat {y}_-)^2}{-\hat {y}_+^2 + 6 \hat {y}_+\hat {y}_- - \hat {y}_-^2}}$ .

• • Otherwise, for $t> t_*,$ the vortices pass each other:

  $$ \begin{align*}\widetilde{\mathbf{z}}^\varepsilon_\pm(t) = \Big(x_i \pm t\frac{\Gamma}{4\pi\hat{y}_\pm}, \varepsilon\hat{y}_\pm\Big).\end{align*} $$

Proof. We only prove the convergence at the level of trajectories; upgrading to all-time convergence is analogous to the proof of the “fission” Proposition 4.4. Again from conservation of Hamiltonian,

$$ \begin{align*} \frac{1}{\mu^2 + x_{\mathsf{rel}}^2} + \frac{1}{4y_{\mathsf{abs}}^2 - \mu^2} = \exp\left({-\frac{4\pi H}{\varepsilon^2\Gamma^2}}\right), \end{align*} $$

where $\mu := y_+ - y_-$ is the conserved “momentum” corresponding to the translational invariance. If $\exp \left ({-\frac {4\pi H}{\varepsilon ^2\Gamma ^2}}\right ) < {1}/{\mu ^2}$ , this curve in the $(x_{\mathsf {rel}}, y_{\mathsf {abs}})$ -plane has a vertical asymptote, which implies that as vortices approach each other they merge into a dipole. Otherwise, $y_{\mathsf {abs}}$ is bounded, so the vortices pass each other. Evaluating all quantities at the initial data, we observe that

$$ \begin{align*} \exp\left({-\frac{4\pi H}{\varepsilon^2\Gamma^2}}\right) = \frac{1}{4\varepsilon^2\hat{y}_+\hat{y}_+} + O(1) \quad\text{and}\quad \frac{1}{\mu^2} = \frac{1}{\varepsilon^2(\hat{y}_+ - \hat{y}_-)^2}. \end{align*} $$

Consequently, for small $\varepsilon $ , condition $\exp \left ({-\frac {4\pi H}{\varepsilon ^2\Gamma ^2}}\right )\leq \frac {1}{\mu ^2}$ will be satisfied if and only if

$$ \begin{align*} (y_+ - y_-)^2 \leq 4y_+ y_- \Leftrightarrow 3 -2\sqrt{2} < \frac{y_+}{y_-} < 3 + 2\sqrt{2}. \end{align*} $$

Note that, in terms of the silver ratio, $\chi ^2= 3 + 2\sqrt {2} $ and $ \chi ^{-2} = 3 - 2\sqrt {2} $ . In this case, as $y_{\mathsf {abs}} \rightarrow \infty $ (corresponding to the limit $t \rightarrow \infty $ ) we have $ \exp \left ({\frac {4\pi H}{\varepsilon ^2\Gamma ^2}}\right ) \xrightarrow {t\to \infty } \mu ^2 + x_{\mathsf {rel}}^2, $ so

$$ \begin{align*} \tan\theta^\varepsilon(t) = \frac{\mu}{x_{\mathsf{rel}}(t)} \xrightarrow{t\to \infty} \frac{\mu}{\sqrt{\exp\left({\frac{4\pi H}{\varepsilon^2\Gamma^2}}\right) - \mu^2}} \xrightarrow{\varepsilon\to 0} \sqrt{\frac{(\hat{y}_+ - \hat{y}_-)^2}{-\hat{y}_+^2 + 6 \hat{y}_+\hat{y}_- - \hat{y}_-^2}} = \tan\theta. \end{align*} $$

Finally, the width of the merged dipole can again be found from conservation of energy:

$$ \begin{align*} \sqrt{\mu^2 + x_{\mathsf{rel}}^2} \xrightarrow{t\to \infty} \exp\left({\frac{2\pi H}{\varepsilon\Gamma}}\right) \xrightarrow{\varepsilon\to 0} \sqrt{4\varepsilon^2\hat{y}_+\hat{y}_-}= 2\varepsilon\sqrt{\hat{y}_+\hat{y}_-}. \end{align*} $$

The speed is computed from the fact that it is inversely proportional to the width: $ v = \frac {\Gamma }{4\pi \sqrt {\hat {y}_+\hat {y}_-}}. $

4.3.4 Vortex pairs, leapfrogging and the metallic means

It turns out that both the golden and silver ratios also manifest themselves in the motion of vortex pairs in the half-plane, that is, point vortices of the same sign and strength. We present it here to demonstrate ubiquitous appearance of the metallic means. Figure 13 contains different types of motions for a vortex pair for different values of $\mu $ .

Figure 13 Types of motion of two vortices of equal circulations in the half-plane. Top: As the parameter $\mu $ defined in (4.1) exceeds the golden ratio, the vortices stop reversing. Bottom: As the parameter $\mu $ exceeds the silver ratio, periodic leapfrogging motion ceases and the vortices start passing each other only once. See [Reference Khesin and Wang19].

Since vortex pairs behave differently at infinity, they need a different normalization. We use the fact that they necessarily happen to be at the same vertical at least once. Let $0<y_1^*\le y_2^*$ be their y-coordinates when the vortices are on the same vertical. We again normalize them so that the distance between the lowest vortex ${y}_1^*$ and its mirror image $\bar y_1^*$ equals 1: $|{y}_1^*-\bar y_1^*|=1$ . Then

Theorem 4.9. The motion for the vortex pair changes as follows depending on $|\mu |:=y_2^*- y_1^*$ :

• • for $0\le |\mu |<\phi $ the vortices leapfrog with periodic reversing of the top one;

• • for $|\mu |=\phi $ the upper vortex has an instantaneous stop, and its trajectory has a cusp;

• • for $\phi \le |\mu |<\chi $ the vortices have a periodic leapfrogging motion without reversing;

• • for $\chi \le |\mu |$ the vortices pass each other only once.

Proof. The golden ratio part was proved in [Reference Khesin and Wang19] in the normalization using the cross-ratio involving four points, $y_2^*, y_1^*, \bar y_1^*, $ and $\bar y_2^*$ . It is equivalent to the one above upon fixing $|{y}_1^*-\bar y_1^*|=1$ .

The silver ratio part follows from the paper [Reference Love21]. Indeed, Love proved that a periodic leapfrogging motion with or without reversing exists under the condition $3 - 2\sqrt {2} < {y}^*_2/{y}^*_1 < 3 + 2\sqrt {2} $ . As we have shown above, it is equivalent to $|\mu |<\chi $ .

There arises a natural question of whether other metallic means also appear in the problem of point vortices or in two-dimensional hydrodynamics in general.

4.4 Numerical experiments with vortex dipoles

In this section, we present a few numerical examples of point vortex billiards. One set of examples is for a family of domains where convexity is eventually lost. Another example comprises multiple dipole billiards on the same disk, resulting in a chaotic motion.

Let us start with the case of a single dipole on disk. Recall that the motion of two vortices in the disk is integrable (due to the rotational and time-translation symmetries of the system, the angular momentum and Hamiltonian serve as two integrals of motion in involution.), see, for example, [Reference Geshev and Chernykh13, Reference Yang30]. Both the motion of a nonsingular vortex dipole pair on the disk and the vortex billiard as its limit have caustics; see Figure 14 for a demonstration of Proposition 2.19 for such caustics. In particular, the corresponding trajectories can be found explicitly, rather than numerically. Bifurcation diagrams as interaction parameters vary have been studied in [Reference Ryabov and Shadrin26, Reference Ryabov, Sokolov and Palshin27, Reference Sokolov and Ryabov28].

Figure 14 The caustic circle for one dipole on the disk. Three different initial conditions, with decreasing separations left to right. Trajectories of vortices with positive and negative circulations are shown, respectively, in red and blue.

Even far from the singular limit, the motion resembles a billiard system (Figure 15, left panel). On the disk, interaction with image vortices sitting at “inverted” points outside the circle completely characterizes the interaction with the boundary (Figure 15, right panel).

Figure 15 Left: Vortex dipole on disk with large inter-vortex distance. Right: trajectories of the dipole together with its image charges.

4.4.1 Single dipole on Neumann ovals

We present here a numerical study of the point vortex system on Neumann ovals, a family of domains which interpolates from a disk domain to a two-disk domain through hourglass-type domains. All of these can be defined as conformal images of unit disk under the two-parameter ( $a,\lambda $ ) mapping

$$ \begin{align*} z=F_{a,\lambda}(Z) := \frac{aZ}{1 - \lambda^2 Z^2}. \end{align*} $$

The parameter a can be chosen to fix the area of the domain to be $\pi $ by taking $a(\lambda ) = \frac {1-\lambda ^4}{(1+\lambda ^4)^{1/2}}$ . The parameter $\lambda $ controls the shape of the oval: $\lambda =0$ is the disk, while for $\lambda \in (0,1)$ is an interpolating family of hourglass domains which limit as $\lambda \to 1$ to two disks.

The practical utility of this family of domains is that one can relate explicit Green’s function on the disk domain (obtained by the method of inversion), to Green’s function on the Neumann oval domain via the (inverse) conformal map. Specifically, identifying $\mathbb {R}^2$ with $\mathbb {C}$ and letting f be the inverse map of $F_\lambda $ , Green’s function on the Neumann oval is

$$ \begin{align*} G(z,w) = G_{\mathrm{disk}}(f(z),f(w)), \qquad z,w\in \mathrm{Neumann\ oval}. \end{align*} $$

The Hamiltonian for the vortex system is based on this adjusted Green’s function. We exploit this structure to numerically integrate two point vortices on Neumann oval domains, see Figure 16. For more details and interesting simulations of point vortices on these domains, see [Reference Ashbee, Esler and McDonald1] and [Reference Ashbee2, Appendix A.2.2].

Figure 16 Dipole on Neumann Oval domains with $\lambda =0,0.1,0.3,0.4,0.5,0.7$ , respectively.

We remark that, although the domains are conformally equivalent to the disk, the two-vortex system is not necessarily integrable on the Neumann oval despite being so on the disk. The reason is that, while the Green’s function is conformally invariant, the Robin function is not and hence the corresponding Hamiltonian is not conformally invariant. It would be interesting to study ergodic properties of pensive billiards on this family of domains, and whether or not caustics survive for small values of $\lambda $ .

4.4.2 Multiple interacting dipoles

There is an important difference between the vortex billiard system as compared to other pensive billiards (such as the puck billiard). Namely, if several classical billiards are superimposed on the same table, their interaction is almost surely absent. On the contrary, it is fully generic that different vortex billiards on the same domain will interact, although only on the boundary of the domain, during those time intervals when the dipole is split into two monopole vortices traveling separately along the boundary. Since this boundary is one-dimensional, collisions are guaranteed and can easily happen between vortices of different pairs. That is, vortices generically exchange partners and reappear inside the domain in different vortex pairs. See Figure 17 for a demonstration of this behavior.

Figure 17 Two dipole pairs on the disk. Trajectories of vortices with positive circulation are shown in red, while for vortices with negative circulation in blue. Initially, two dipoles are formed by the solid and dashed vortices. The separation in their initial positions is decreasing left to right.

For this reason, to consider more than one vortex dipole pair, one needs to derive the correct limiting dynamics from the original Kirchhoff point vortex system. However, it turns out that there is only one extra constraint in addition to the above study: if two vortices of the same sign meet when traveling along the boundary, they simply pass by each other in the limit of initial separations taken to zero. We prove this in greater generality here, allowing the circulations of the vortices in the pair to be not equal or opposite. The result is that there are no mergers of vortices of the same sign in the half-plane:

Proposition 4.10. Consider two point vortices of circulations $\varepsilon \Gamma _1$ and $\varepsilon \Gamma _2$ with $\Gamma _1 \neq -\Gamma _2$ in the half-plane, starting from initial positions $(x_1, \varepsilon \hat {y}_1)$ , $(x_2, \varepsilon \hat {y}_2)$ . Assume in addition that $x_1 \neq x_2$ . Then, there exists $\varepsilon _0> 0$ such that for all $0 < \varepsilon < \varepsilon _0$ vortices pass one another.

Proof. Introduce relative and absolute coordinates:

$$ \begin{align*} x_{\mathsf{rel}} = x_1-x_2, \quad y_{\mathsf{rel}} = y_1 - y_2,\quad x_{\mathsf{abs}}= \frac{\Gamma_1 x_1 + \Gamma_2 x_2}{\Gamma_1 + \Gamma_2},\quad y_{\mathsf{abs}}= \frac{\Gamma_1 y_1 + \Gamma_2 y_2}{\Gamma_1 + \Gamma_2}, \end{align*} $$

and denote

$$ \begin{align*} \hat{y}_{\mathsf{rel}} = \hat{y}_1 - \hat{y}_2,\quad \hat{y}_{\mathsf{abs}}= \frac{\Gamma_1 \hat{y}_1 + \Gamma_2 \hat{y}_2}{\Gamma_1 + \Gamma_2}. \end{align*} $$

Conservation of linear momentum $y_{\mathsf {abs}}$ and energy H yields the equation of trajectory in $(x_{\mathsf {rel}}, y_{\mathsf {rel}})$ coordinates

$$ \begin{align*} &\left(\hat{y}_{\mathsf{abs}} + \frac{\Gamma_2}{\Gamma_1 + \Gamma_2}\frac{y_{\mathsf{rel}}}{\varepsilon}\right)^{\Gamma_1^2} \left(\hat{y}_{\mathsf{abs}} - \frac{\Gamma_1}{\Gamma_1 + \Gamma_2}\frac{y_{\mathsf{rel}}}{\varepsilon}\right)^{\Gamma_2^2}\bigg[\frac{x_{\mathsf{rel}}^2 + \varepsilon^2(2\hat{y}_{\mathsf{abs}} - \frac{\Gamma_1 - \Gamma_2}{\Gamma_1 + \Gamma_2}\frac{y_{\mathsf{rel}}}{\varepsilon})^2}{x_{\mathsf{rel}}^2 + y_{\mathsf{rel}}^2}\bigg]^{\Gamma_1\Gamma_2} = \\ &\qquad \left(\hat{y}_{\mathsf{abs}} + \frac{\Gamma_2}{\Gamma_1 + \Gamma_2}\hat{y}_{\mathsf{rel}}\right)^{\Gamma_1^2} \left(\hat{y}_{\mathsf{abs}} - \frac{\Gamma_1}{\Gamma_1 + \Gamma_2}\hat{y}_{\mathsf{rel}}\right)^{\Gamma_2^2}\bigg[\frac{\hat{x}_{\mathsf{rel}}^2 + \varepsilon^2(2\hat{y}_{\mathsf{abs}} - \frac{\Gamma_1 - \Gamma_2}{\Gamma_1 + \Gamma_2}\hat{y}_{\mathsf{rel}})^2}{\hat{x}_{\mathsf{rel}}^2 +\varepsilon^2 \hat{y}_{\mathsf{rel}}^2}\bigg]^{\Gamma_1\Gamma_2}. \end{align*} $$

Rearranging, the above equality gives

$$ \begin{align*} x_{\mathsf{rel}}^2 &= \tfrac{y_{\mathsf{rel}}^2 \left(\hat{y}_{\mathsf{abs}} + \frac{k}{1+k}\hat{y}_{\mathsf{rel}}\right)^{1/k} \left(\hat{y}_{\mathsf{abs}} - \frac{1}{1+k}\hat{y}_{\mathsf{rel}}\right)^{k} \left(\frac{\hat{x}_r^2 + \varepsilon^2(2\hat{y}_{\mathsf{abs}} - \frac{1-k}{1+k}\hat{y}_{\mathsf{rel}})^2}{\hat{x}_{\mathsf{rel}}^2 + \varepsilon^2\hat{y}_{\mathsf{rel}}^2}\right) - \varepsilon^2(2\hat{y}_{\mathsf{abs}} - \frac{1 - k}{1 + k}\frac{y_{\mathsf{rel}}}{\varepsilon})^2\left(\hat{y}_{\mathsf{abs}} + \frac{k}{1+k}\frac{y_{\mathsf{rel}}}{\varepsilon}\right)^{1/k} \left(\hat{y}_{\mathsf{abs}} - \frac{1}{1+k}\frac{y_{\mathsf{rel}}}{\varepsilon}\right)^{k}}{ \left(\hat{y}_{\mathsf{abs}} + \frac{k}{1+k}\frac{y_{\mathsf{rel}}}{\varepsilon}\right)^{1/k} \left(\hat{y}_{\mathsf{abs}} - \frac{1}{1+k}\frac{y_{\mathsf{rel}}}{\varepsilon}\right)^{k} - \left(\hat{y}_{\mathsf{abs}} + \frac{k}{1+k}\hat{y}_{\mathsf{rel}}\right)^{1/k} \left(\hat{y}_{\mathsf{abs}} - \frac{1}{1+k}\hat{y}_{\mathsf{rel}}\right)^{k} \left(\frac{\hat{x}_{\mathsf{rel}}^2 + \varepsilon^2(2\hat{y}_{\mathsf{abs}} - \frac{1-k}{1+k}\hat{y}_{\mathsf{rel}})^2}{\hat{x}_{\mathsf{rel}}^2 + \varepsilon^2\hat{y}_{\mathsf{rel}}^2}\right) }\\ &=:\frac{g_{\varepsilon}(y_{\mathsf{rel}})}{f_{\varepsilon}(y_{\mathsf{rel}})}, \quad \text{where} \quad k := \frac{\Gamma_2}{\Gamma_1}. \end{align*} $$

We aim to show that, for sufficiently small $\varepsilon $ , $x_{\mathsf {rel}}$ is unbounded along the trajectory. This would prove that after the two vortices meet, they do not become bound forever. For that, it is sufficient to establish the existence of $y_{\mathsf {rel}}^{\text {crit}}$ for which

$$ \begin{align*} f_\varepsilon(y_{\mathsf{rel}}^{\text{crit}}) = 0, \quad f^{\prime}_\varepsilon(y_{\mathsf{rel}}) \neq 0, \quad \text{and } g_\varepsilon(y_{\mathsf{rel}}^{\text{crit}}) \neq 0. \end{align*} $$

To this end, consider

$$ \begin{align*} y_{\mathsf{rel}}^* =\arg\max_{y \in (-\frac{k}{1+k} \hat{y}_{\mathsf{abs}},\frac{1}{1+k} \hat{y}_{\mathsf{abs}})} \left(\hat{y}_{\mathsf{abs}} + \frac{k}{1+k}y\right)^{1/k} \left(\hat{y}_{\mathsf{abs}} - \frac{1}{1+k}y\right)^{k}, \end{align*} $$

and observe that for sufficiently small $\varepsilon>0$ , $f_{\varepsilon }(y_{\mathsf {rel}}^*)> 0 $ and $ f_{\varepsilon }(y_{\mathsf {rel}}) < 0$ . Hence, for sufficiently small $\varepsilon $ such $y_{\mathsf {rel}}^{\text {crit}}$ satisfying $f_{\varepsilon }(y_{\mathsf {rel}}^{\text {crit}}) = 0$ exists. The other two conditions follow from a direct verification.

Corollary 4.11. Only vortices of equal strength and opposite sign can exchange their pairs. They form a multiple dipole billiard system satisfying the fission–fusion rules. In other options the vortices of different pairs do not interact in the limit of zero initial separation.

As a consequence of Proposition 4.10, the system for multiple dipoles, strictly speaking, is no longer a billiard. Indeed, one must keep track of the relative speeds (or, equivalently, distances to the boundary) of vortices of opposite signs to determine the angle for their reentering the domain if they ever meet on the boundary. Unlike the case of one dipole, such an angle is no longer conserved, since a vortex might meet a stranger vortex, of different speed, before its original partner, while it travels along the boundary. One can check that the new “fusion rule” for vortices meeting along the boundary is:

1. 1. if $\Gamma _1 = -\Gamma _2$ and the ratio of their boundary velocities satisfies

   $$ \begin{align*} \chi^{-2}\leq \bigg|\frac{v_1}{v_2}\bigg| \leq \chi^2 \end{align*} $$
   where $\chi $ is the silver ratio, then the monopoles form a dipole which enters the interior at angle $\theta = \theta ({v_1}/{v_2})$ and moves at speed $\sqrt {-v_1v_2}$ . Here $\theta ({v_1}/{v_2})$ is the explicit function derived in Proposition 4.5:
   $$ \begin{align*} \theta(r) = \tan^{-1}\sqrt{\frac{(r - 1/r)^2}{-r^2 + 6 r - 1}}. \end{align*} $$

2. 2. otherwise, monopoles continue traveling along the boundary without interaction.

It is interesting that the addition of just one extra dipole of equal strength seems to yield a completely chaotic dynamical system. See Figure 17 for numerical simulations on a disk domain.

5.1 Definition of pensive outer billiards

Recall the definition of an outer billiard, following [Reference Tabachnikov29]. Given a smooth strictly convex curve $\gamma $ in the plane, the outer billiard is the following map from the exterior of $\gamma $ into itself. Let X be a point outside $\gamma $ . There exist two tangent lines to $\gamma $ passing through X. Choose one of them, say, the right one from the viewpoint of X, and reflect X in the tangency point P to obtain a new point $Y=\mathsf {OB}(X)$ (see left panel of Figure 18). The map $\mathsf {OB}$ is called the outer billiard map for the curve $\gamma $ .

Note that the definition above is a manifestation of the principle that notions related to the length in classical billiards are replaced with the area for outer billiards. As we show below, on the plane this definition is equivalent to the following one:

For $ a(r)\equiv 0$ or $\theta (r)\equiv 0$ one obtains the standard outer billiard $\mathsf {OB}$ . The equivalence of the two definitions on the plane follows from the following proposition.

Figure 18 Left: outer billiard. Right: pensive outer billiard.

Before proving this statement we introduce coordinates in which outer billiards are convenient to work with. Parametrize $\gamma $ by the angle $\alpha $ that the tangent makes with a given (say horizontal) direction. For each point $X \in \mathbb {R}^2\setminus D$ , let right and left tangents from X to $\gamma $ have lengths $r, \bar {r}$ and meet $\gamma $ at points $\gamma (\alpha ), \gamma (\bar {\alpha })$ . In this way, $(\alpha ,r)$ and $(\bar {\alpha }, \bar {r})$ form two coordinate systems on $\mathbb {R}^2\setminus D$ .

Proof. Indeed, these are essentially polar coordinates on the plane. This also follows from a direct computation using the formulas

$$ \begin{align*} (x,y) = \gamma(\alpha) + r\frac{\gamma'(\alpha)}{|\gamma'(\alpha)|}, \quad \text{where} \quad \frac{\gamma'(\alpha)}{|\gamma'(\alpha)|} = (\cos\alpha,\sin\alpha), \end{align*} $$

and the analogous ones for $(\bar {\alpha }, \bar {r})$ .

Now the proof of Proposition 5.3 is immediate.

Proof. This follows from the fact that the area swept by a tangent segment of length r traveling along $\gamma $ as $\alpha $ changes from $\alpha _1$ to $\alpha _2$ is

$$ \begin{align*} \int_{0}^{r}\int_{\alpha_1}^{\alpha_2}rdrd\alpha = \int_{0}^{r}rdr\int_{\alpha_1}^{\alpha_2}d\alpha =\frac{r^2}{2}(\alpha_2 - \alpha_1).\\[-45pt] \end{align*} $$

The main property of a pensive outer billiard, “inherited” from the standard outer billiard is given by the following

Proof. By definition of the outer billiard map, it sends a point with coordinates $(\alpha ,r)$ to the point with coordinates $(\bar {\alpha }, \bar {r}) = (\alpha , r)$ . It is then immediate that $\mathsf {OB}$ is area-preserving:

$$ \begin{align*} \mathsf{OB}^*\mu = \mathsf{OB}^*(\bar{r}d\bar{r}\wedge d\bar{\alpha}) = rdr\wedge d\alpha = \mu \end{align*} $$

We then note that by Proposition 5.3, pensive outer billiard map $\mathsf {PO}$ is a composition of $\mathsf {OB}$ and a shift map $\mathsf {Sh}$ , given in the $(\alpha , r)$ coordinates by $(\alpha ,r) \mapsto (\alpha + 2a(r)/r^2, r)$ . The latter map is also area-preserving:

$$ \begin{align*} \mathsf{Sh}^*\mu = \mathsf{Sh}^*(rdr\wedge d\alpha) = rdr\wedge d(\alpha+ 2a(r)/r^2) = rdr\wedge d\alpha = \mu.\\[-37pt] \end{align*} $$

5.2 Duality of pensive billiards on the sphere

Similarly to classical and pensive billiards, the definition of the pensive outer billiard can be generalized to nonflat domains. The following consideration of the pensive outer billiard on a spherical domain justifies Definition 5.1 with the swept area, as opposed to Definition 5.2 valid only in the flat case: on the sphere these two definitions are not equivalent!

Theorem 5.6. The pensive billiard with delay function $\tilde {\ell }(\theta )$ and pensive outer billiard with delay function $a(\theta )$ are projectively dual on the unit sphere, provided that

(5.1) $$ \begin{align} a(\theta) = \tilde{\ell}(\theta)(1 - \cos\theta). \end{align} $$

Figure 19 Duality of pensive billiards and pensive outer billiards on the sphere. Pensive billard map with respect to red curve $\gamma $ sends great circle x to great circle y. Pensive outer billiard map with respect to dual curve $\gamma _*$ sends point X dual to x to point Y dual to y. Moreover, functions $\tilde {\ell }$ and a are related by Equation (5.1)

This theorem mimics the one for the classical and outer billiards on the sphere, see [Reference Tabachnikov29]. The projective duality on the sphere interchanges points (poles) and the corresponding oriented great circles (equators), analogs of lines on the sphere. An oriented tangent line (great circle) to a curve $\gamma $ is sent to a point of the dual curve $\gamma _*$ , which by definition consists of poles of the tangent great circles to $\gamma $ . An immediate corollary of this duality is that the pensive outer billiard on the sphere preserves the spherical areas.

Compared to the classical case, we need the following additional geometric result, which can be regarded as a generalization of the Archimedes theorem (as we discuss below).

Proposition 5.7. Let $\sigma (s)$ be any curve in $S^2 \subset \mathbb {R}^3$ (not necessarily parametrized by arclength). The area swept by the segment of angular length $\theta $ tangent to $\sigma $ at $\sigma (s)$ as s varies between $s_1$ and $s_2$ is

$$ \begin{align*} (1 - \cos\theta)\int_{s_1}^{s_2}\frac{|\sigma"(s) \cdot (\sigma(s) \times \sigma'(s)))|}{|\sigma'(s)|^2}ds. \end{align*} $$

Proof. The region swept by the tangent segment of length $\theta $ is parametrized by

$$ \begin{align*} r(s,\beta) = \cos\beta\sigma(s) + \sin\beta\frac{\sigma'(s)}{|\sigma'(s)|}, \qquad (s, \beta) \in [s_1, s_2]\times[0, \theta]. \end{align*} $$

The area is consequently computed as

$$ \begin{align*} \begin{aligned} &\int_0^\theta\int_{s_1}^{s_2}\bigg|\frac{\partial r}{\partial s}\times \frac{\partial r}{\partial \beta}\bigg| \,ds\,d\beta = \int_0^\theta\int_{s_1}^{s_2}\bigg|r \cdot \bigg(\frac{\partial r}{\partial s}\times \frac{\partial r}{\partial \beta}\bigg)\bigg| \,ds\,d\beta \\ &= \int_0^\theta\int_{s_1}^{s_2}\sin\beta\frac{|\sigma"(s) \cdot (\sigma(s) \times\sigma'(s))|}{|\sigma'(s)|^2}\,ds\,d\beta =(1 -\cos \theta)\int_{s_1}^{s_2}\frac{|\sigma"(s) \cdot (\sigma(s) \times \sigma'(s)))|}{|\sigma'(s)|^2}ds, \end{aligned} \end{align*} $$

where we used that $\sigma (s)\cdot \sigma '(s) = 0$ together with the symmetries of the triple product.

Proof of Theorem.

The (pensive) billiard map acts on oriented curves: an incidence ray is mapped to the reflected one. For the classical billiard, at the moment of reflection in $\gamma $ the incidence ray a, the tangent p to $\gamma $ and the reflected ray b pass through the same point and make equal angles. On the sphere this translates to the corresponding three poles $A, P,$ and B, lying on the same line (great circle) and making distances equal to the angle of incidence, $|AP|=|{PB}| = \theta $ . The latter is the definition of the outer billiard.

For the pensive billiards, before reflecting the point slides along $\gamma $ for the distance $\tilde {\ell }(\theta )$ depending on the incidence angle $\theta $ and then reflected from the tangent $\tilde p$ . We shall see that the area swept by the corresponding segment moving along the dual curve is equal to

$$ \begin{align*} \tilde{\ell}(\theta)(1 - \cos\theta). \end{align*} $$

Given an arclength parametrization $\gamma (s)$ of $\gamma $ , dual curve $\gamma _*$ can be parametrized by

$$ \begin{align*} \gamma_*(s) = \gamma(s)\times \gamma'(s). \end{align*} $$

Hence, the statement will follow from Proposition 5.7 applied to $\sigma = \gamma _*$ if we prove that

$$ \begin{align*} \frac{|\gamma_*"(s) \cdot (\gamma_*(s) \times \gamma_*'(s)))|}{|\gamma_*'(s)|^2} = 1. \end{align*} $$

This is indeed the case. Expanding $\gamma "(s)$ in the orthonormal basis $(\gamma (s), \gamma '(s), \gamma (s)\times \gamma '(s))$ we find

$$ \begin{align*} \gamma"(s) = -\gamma(s) + (\gamma"(s)\cdot(\gamma(s)\times\gamma'(s)))\gamma(s)\times\gamma'(s). \end{align*} $$

Hence,

$$ \begin{align*} |\gamma_*'(s)|^2 = |\gamma(s)\times\gamma"(s)|^2| = |(\gamma"(s)\cdot(\gamma(s)\times\gamma'(s))|^2. \end{align*} $$

On the other hand,

$$ \begin{align*} \begin{aligned} |\gamma_*"(s) \cdot (\gamma_*(s) \times \gamma_*'(s)))| &= |\gamma_*'(s) \cdot (\gamma_*(s) \times \gamma_*"(s)))| = |(\gamma(s)\times\gamma"(s)) \cdot (\gamma_*(s) \times \gamma_*"(s)))| =\\ &=|\gamma"(s)\cdot(\gamma(s)\times\gamma'(s))| |\gamma'(s) \cdot (\gamma_*(s) \times \gamma_*"(s)))|\\ &=|\gamma"(s)\cdot(\gamma(s)\times\gamma'(s))| |\gamma_*"(s) \cdot \gamma(s)|\\ &=|\gamma"(s)\cdot(\gamma(s)\times\gamma'(s))| |(\gamma \times \gamma"'(s) + \gamma'(s) \times \gamma"(s)) \cdot \gamma(s)|\\ &=|\gamma"(s)\cdot(\gamma(s)\times\gamma'(s))||(\gamma'(s) \times \gamma"(s)) \cdot \gamma(s)|\\ &=|\gamma"(s)\cdot(\gamma(s)\times\gamma'(s))|^2. \end{aligned} \end{align*} $$

This concludes our proof.

Note that Proposition 5.7 includes a theorem due to Archimedes as a special case:

Proof of Corollary.

It is enough to prove the theorem for a region R bounded by two parallel planes $\{z = h_1\}$ and $\{z = h_2\}$ , where we can assume $0<h_2 < h_1$ . We note that this region is given by swiping a tangent segment along the circle lying in the first plane; one can see that the length $\theta $ of the segment needs to satisfy $h_1(1 - \cos \theta ) = h_1 - h_2$ . Applying lemma to this segment, note that by symmetry of the triple product, the integrand is given by the projection of $\gamma (s)$ onto $\gamma '(s)\times \gamma "(s)$ , which is here just a projection of $\gamma (s)$ onto the vertical direction, that is, $h_1$ .

Consequently, area of R is given by

$$ \begin{align*} (1 - \cos\theta)\int_{s_1}^{s_2}\frac{|\sigma"(s) \cdot (\sigma(s) \times \sigma'(s)))|}{|\sigma'(s)|^2}ds = \frac{h_1 - h_2}{h_1}\int_{0}^{2\pi}h_1 ds = 2\pi(h_1 - h_2), \end{align*} $$

which is the same as the area of the projection of R onto cylinder.

The above duality of two types of pensive billiards allows one to make hydrodynamical meaning for pensive outer billiards on the sphere: they are dual to systems of vortex dipoles in spherical regions in the limit of zero separation.

We note that, since many facts in dynamics often turn out to be simpler to prove for outer billiards than for classical ones, it would be interesting to extend other discussed results, including the twist property, generating functions, the golden and silver ratios, periodic trajectories, etc. to the domain of pensive outer billiards. One can consider also pensive outer billiards on other surfaces with constant curvature, for example, the hyperbolic plane. Various features, such as area preservation of the billiard map, can be extended to that case.

In conclusion, there are many interesting questions worth investigating for the pensive billiard system which we have not studied. For example,

• • The pensive billiard map has a variational formulation but, in general, is not a twist map. What can be said (using Morse theory) about the periodic orbits in this case?

• • We prove very little about the existence (and nonexistence) of caustics. Does KAM theory applies, and is there a version of Lazutkin’s theorem in this setting? What can be said toward Mather’s nonexistence criterion? First results in that direction for puck billiards were recently obtained in [Reference Barbieri and Clarke3].

• • The billiard system inside ellipses is integrable. For which combinations of the delay functions and shapes of billiard tables will the pensive billiards be integrable?

• • What is the multidimensional version of pensive billiards?