# On the energy of topological defect lattices

B. Berche${}^{\textsf{\footnotesize{\mbox{?}}}}$, S.
Fumeron${}^{\textsf{\footnotesize{\mbox{?}}}}$, F.
Moraes${}^{\textsf{\footnotesize{\mbox{?}}}}$

(Received January 28, 2020, in final form February 27, 2020)

###### Abstract

Since the logarithm function is the solution of Poisson’s equation in two
dimensions, it appears as the Coulomb interaction in two dimensions, the
interaction between Abrikosov flux lines in a type II superconductor, or
between line defects in elastic media, and so on. Lattices of lines
interacting logarithmically are, therefore, a subject of intense research due
to their manifold applications. The solution of the Poisson equation for such
lattices is known in the form of an infinite sum since the late 1990’s. In
this article we present an alternative analytical solution, in closed form, in
terms of the Jacobi theta function.

Key words: topological defect, cosmic string, flux line

###### Abstract

Îñêëüêè ëîãàðèôìíà ôóíêöÿ ðîçâ’ÿçêîì ðâíÿííÿ Ïóàñîíà ó äâîõ âèìðàõ, âîíà
ôàêòèíî êóëîíâñüêîþ âçàìîäþ ó äâîâèìðíîìó âèïàäêó, âçàìîäþ ìæ ëíÿìè ïîòîêó
Àáðèêîñîâà ó íàäïðîâäíèêó II òèïó, àáî ìæ ëíéíèìè äåôåêòàìè ó ïðóæíèõ
ñåðåäîâèùàõ, ò. ï. Îñü îìó ãðàòêè ëíé, ùî âçàìîäþòü ëîãàðèôìíî ïðåäìåòîì
íòåíñèâíèõ äîñëäæåíü çàâäÿêè ¿õ áàãàòîãðàííèì çàñòîñóâàííÿì. Ùå ç êíöÿ 1990-õ
ðîêâ ðîçâ’ÿçîê ðâíÿííÿ Ïóàñîíà äëÿ òàêèç ãðàòîê áóâ âäîìèé ó âèãëÿä íåñêíåíî¿
ñóìè. Ó äàíé ñòàòò ïðåäñòàâëåíî àëüòåðíàòèâíèé àíàëòèíèé ðîçâ’ÿçîê ó
çàìêíóòîìó âèãëÿä â òåðìíàõ òåòà-ôóíêö¿ ßêîá.

Ключов слова: òîïîëîãíèé äåôåêò, êîñìíà ñòðóíà, ëíÿ ïîòîêó

A paper dedicated to our friend Ihor Mryglod

on the occasion of his 60th birthday.

##  1 Introduction

As intriguing as topological defects might be, crystals made out of them
appear to be even more exotic, like the soliton lattice that forms in doped
polyacetylene [1], for instance. Examples of topological defect lattices
abound in Condensed Matter Physics where one might find lattices of parallel
screw dislocations in solids [2], vortex lattices in rotating superfluids [3]
and in Bose-Einstein condensates [4], as well as the much studied magnetic
flux lattices in type II superconductors [5]. Liquid crystals contribute with
lattices of disclinations in nematics [6], and with lattices of screw
dislocations in cholesterics (known as twist grain boundaries) [7].
Nevertheless, the most fascinating topological defect lattices are found in
the realm of chiral liquid crystals [8] where skyrmions [9], hopfions (3D
skyrmions) [10], merons (half-skyrmions) [11] and even knots [12] may form
regular arrays. Off this planet one might have magnetic flux tube lattices in
neutron stars [13] and crystals of cosmic strings or of cosmic domain walls,
which have been considered as possible candidates for solid dark matter models
[14, 15]. All this zoo of topological defects shares a common origin: phase
transitions involving break of symmetry. Not surprisingly thus, the Kibble-
Zurek mechanism [16, 17] of defect formation applies both to cosmic strings
and disclinations in nematic liquid crystals [18]. A common feature of most of
the aforementioned topological defect crystals is a logarithmic interaction of
line defects in the lattice, for large enough separation between them, so the
defect internal structure may be neglected. This leads to the problem of
performing infinite log sums, much tackled in the 1990’s.

An important step was done in the calculation of energies and forces between
particles interacting logarithmically by [19, 20] which much improved the
efficiency of computer simulations. The expressions were obtained in terms of
products of elementary trigonometric or hyperbolic functions. In this work we
move a step forward and obtain for the solution of Poisson equation a closed
form for the logarithmic sum involving Jacobi theta functions. These functions
are special functions of complex variables which appear in the theory of
elliptic functions which are ubiquitous in mathematical physics.

Logarithmic potential appears as a solution of the two-dimensional Poisson
equation, and thus, describes the interaction in a two-dimensional Coulomb
gas, so there is no mystery in the appearance of elliptic and Jacobi theta
functions in the vortices-driven Berezinskii-Kosterlitz-Thouless transition
[21, 22, 23, 24]. Indeed, applying the theory of conformal mappings which hold
at the critical point of second order phase transitions in two-dimensional
systems, these special functions enable us to obtain a closed expression for
the correlation functions in $XY$ models [25]. More generically, they appear
in the Schwarz-Christoffel mapping and related conformal mappings in the
complex plane [26, 27, 28, 29].

We consider an infinite lattice of parallel string-like defects in 3D or,
equivalently, point-like defects in 2D, interacting logarithmically. Our
interest is to find the energy and, consequently, the force on a test defect
due to its interaction with the lattice. Although previous results for finite
lattices with periodic boundary conditions have been reported [30, 31, 32,
33], to the best of our knowledge, this is the first time where a closed form
for the logarithmic sum is achieved for the infinite lattice.

###  1.1 Rectangular lattice

Let us consider a rectangular Bravais lattice in $\mathbb{R}^{2}$ generated by
the basis vectors $\vec{a}_{1}=a\hat{x}$ and $\vec{a}_{2}=b\hat{y}$ such that
a point of the lattice located at $\vec{R}_{mn}=m\vec{a}_{1}+n\vec{a}_{2}$ is
associated to the pair $(m,n)\in\mathbb{Z}^{2}$. To each point of the lattice
we associate a defect. It is our purpose to find the potential due to this
array of defects, assuming the superposition principle. That is, we want to
perform the sum

$\displaystyle
V(\vec{r})=\lambda\sum_{(m,n)\in\mathbb{Z}^{2}}\ln|\vec{r}-\vec{R}_{mn}|^{2},$
(1.1)

where $\vec{r}=x\hat{x}+y\hat{y}$ is the position of a test defect and
$\lambda$ is the “charge” of the logarithmic interaction. Obviously, the
function defined by equation (1.1) is a solution of the 2D Poisson equation

$(\partial^{2}_{xx}+\partial^{2}_{yy})V=2\piup\lambda\sum_{(m,n)\in\mathbb{Z}^{2}}\delta\left(x-ma\right)\delta\left(y-nb\right)$
(1.2)

and, therefore, we will not be concerned with additive constants appearing in
the logarithmic sum. This is the essence of the regularization process that we
need to use since the “raw” sum in equation (1.1) naturally diverges.

Now, defining

$\displaystyle\varphi$ $\displaystyle=x+\mathrm{i}y,$
$\displaystyle\bar{\varphi}$ $\displaystyle=x-\mathrm{i}y,$ (1.3)
$\displaystyle\sigma$ $\displaystyle=ma+\mathrm{i}nb,$
$\displaystyle\bar{\sigma}$ $\displaystyle=ma-\mathrm{i}nb$ (1.4)

we write equation (1.1) as

$V(x,y)=\lambda\sum_{(m,n)\in\mathbb{Z}^{2}}\ln\left[(\varphi-\sigma)(\bar{\varphi}-\bar{\sigma})\right].$
(1.5)

As mentioned above, the sums in equations (1.1) and (1.5) diverge but can be
regularized by subtracting constant divergent terms as we will see below.

Choosing to first perform the sum over $n$ in equation (1.5), we have

$\displaystyle V(x,y)$ $\displaystyle=$
$\displaystyle\lambda\sum_{m=-\infty}^{\infty}\ln\left[(\varphi-
ma)(\bar{\varphi}-ma)\right]$ (1.6) $\displaystyle+$
$\displaystyle\lambda\sum_{m=-\infty}^{\infty}\sum_{n=1}^{\infty}\ln\left\\{[(\varphi-
ma)^{2}+n^{2}b^{2}][(\bar{\varphi}-ma)^{2}+n^{2}b^{2}]\right\\}.$

Note that, in changing the sum over $\mathbb{Z}^{-}$ into a sum over
$\mathbb{Z}^{+}$, we have $\sigma\rightarrow\bar{\sigma}$ such that
$(\varphi-\sigma)(\bar{\varphi}-\bar{\sigma})$ from the $\mathbb{Z}^{-}$ sum
becomes $(\varphi-\bar{\sigma})(\bar{\varphi}-\sigma)$.

Now, using the identity

$\prod_{n=1}^{\infty}\left(1+\frac{z}{n^{2}A+B}\right)=\sqrt{\frac{B}{B+z}}\frac{\sinh(\piup\sqrt{B+z}/\sqrt{A})}{\sinh(\piup\sqrt{B/A})}\,,$
(1.7)

equation (1.6) writes

$\displaystyle V(x,y)$
$\displaystyle=\lambda\sum_{(m,n)\in\mathbb{Z}^{2}}\ln\left(m^{2}a^{2}+n^{2}b^{2}\right)+\lambda\sum_{m=-\infty}^{\infty}\ln\left(m^{2}a^{2}\right)$
$\displaystyle+\lambda\sum_{m=-\infty}^{\infty}\ln\left\\{\frac{\sinh[\piup(\varphi-
ma)/b]}{\sinh(\piup
ma/b)}\right\\}+\lambda\sum_{m=-\infty}^{\infty}\ln\left\\{\frac{\sinh[\piup(\bar{\varphi}-ma)/b]}{\sinh(\piup
ma/b)}\right\\},$ (1.8)

which results in

$\displaystyle V(x,y)$
$\displaystyle=\lambda\sum_{(m,n)\in\mathbb{Z}^{2}}\ln\left(m^{2}a^{2}+n^{2}b^{2}\right)+2\lambda\sum_{m=1}^{\infty}\ln\left(m^{2}a^{2}\right)+\lambda\ln\left(b^{2}/\piup^{2}\right)$
$\displaystyle+\lambda\ln[\sinh(\piup\varphi/b)\sinh(\piup\bar{\varphi}/b)]$
$\displaystyle+\lambda\sum_{m=1}^{\infty}\ln\left\\{\frac{\sinh[\piup(\varphi-
ma)/b]\sinh[\piup(\varphi+ma)/b]}{\sinh^{2}(\piup ma/b)}\right\\}$
$\displaystyle+\lambda\sum_{m=1}^{\infty}\ln\left\\{\frac{\sinh[\piup(\bar{\varphi}-ma)/b]\sinh[\piup(\bar{\varphi}+ma)/b]}{\sinh^{2}(\piup
ma/b)}\right\\}.$ (1.9)

The first three terms in equation (1.9) just add up to an infinite constant
and can be removed from the potential since it will still be a solution of
equation (1.2).

The remaining terms can be evaluated by the use of the identity

$\prod_{m=1}^{\infty}\left[\cos^{2}z+\sin^{2}z\coth^{2}(m\piup\chi)\right]=\csc
z\
\frac{\vartheta_{1}(z,\mathrm{e}^{-\piup\chi})}{\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup\chi})}\,,$
(1.10)

in terms of the Jacobi theta function $\vartheta_{1}$ and its first derivative
w.r.t. $z$. As warned by Abramowitz and Stegun [34], there is a bewildering
variety of notations for the theta functions. The one we use here, given by
[34] and [35], is such that they have the following Fourier representation:

$\displaystyle\vartheta_{1}(z,\mathrm{e}^{-\piup\chi})=2\sum_{n=0}^{\infty}(-1)^{n}\mathrm{e}^{-\piup\chi(n+1/2)^{2}}\sin\left[(2n+1)z\right],$
(1.11)

where $z$ and $\chi$ are complex numbers. With the help of the above
relations, equation (1.9) and, therefore, the regularized equation (1.1) takes
a surprisingly simple form

$V(x,y)=\lambda\ln\left[\frac{-|\vartheta_{1}(\mathrm{i}\piup\varphi/b,\mathrm{e}^{-\piup
a/b})|^{2}}{\sin(\mathrm{i}\piup\varphi/b)\sin(\mathrm{i}\piup\bar{\varphi}/b)}\right]-2\lambda\ln\left[\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup
a/b})\right].$ (1.12)

This results in the following expression in terms of the coordinates $x$ and
$y$:

$V(x,y)=\lambda\ln(x^{2}+y^{2})+\lambda\ln\left[\frac{|\vartheta_{1}(\frac{\piup}{b}(\mathrm{i}x-y),\mathrm{e}^{-\piup
a/b})|^{2}}{\cosh^{2}(\piup x/b)-\cos^{2}(\piup
y/b)}\right]-2\lambda\ln\left[\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup
a/b})\right].$ (1.13)

Equation (1.13) must have logarithmic singularities at the defect sites since
it is a compact version of equation (1.1). Nevertheless, it seems to have
extra singularities at $(x,y)=(0,nb)$. Since both $\cos(\piup y/b)$ and the
$\vartheta$ function are periodic in $y$ with periodicity $b$, it suffices to
examine this question near $(x,y)=(0,0)$. A closer look at equation (1.12)
indicates that there is no extra singularity there since
$\vartheta_{1}(0,\mathrm{e}^{-\piup\chi})=0$ and

$\lim_{z\rightarrow 0}\frac{\vartheta_{1}(z,\mathrm{e}^{-\piup\chi})}{\sin
z}=\lim_{z\rightarrow
0}\frac{\vartheta_{1}^{\prime}(z,\mathrm{e}^{-\piup\chi})}{\cos
z}=\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup\chi})=2\sum_{n=0}^{\infty}(-1)^{n}\mathrm{e}^{-\piup\chi(n+1/2)^{2}}(2n+1)\,,$
(1.14)

where the values of $\vartheta_{1}(0,q)$ and $\vartheta_{1}^{\prime}(0,q)$
were obtained from equation (1.11). Hence,
$V(0,0)=\lambda\ln(x^{2}+y^{2})|_{(x,y)=(0,0)}$ as it should. Likewise, by
changing the origin in equation (1.13) to $(0,nb)$, it follows that
$V(0,nb)=\lambda\ln[x^{2}+(y-nb)^{2}]|_{(x,y)=(0,nb)}$.

In equation (1.14) we see that, since $\chi$ is real,
$\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup\chi})$ is also real and,
therefore,
$[\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup\chi})]^{2}=|\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup\chi})|^{2}$.
This way, we rewrite equation (1.13) as

$V(x,y)=\lambda\ln\left[\frac{x^{2}+y^{2}}{\cosh^{2}(\piup x/b)-\cos^{2}(\piup
y/b)}\right]+\lambda\ln\left|\frac{\vartheta_{1}(\frac{\piup}{b}(\mathrm{i}x-y),\mathrm{e}^{-\piup
a/b})}{\vartheta_{1}^{\prime}(0,\mathrm{e}^{-\piup a/b})}\right|^{2}.$ (1.15)

In figure 1 we present a plot of this function for a generic rectangular
lattice and for the particular case of a square lattice.

Figure 1: (Colour online) Three-dimensional plot of the “potential function”
$V(x,y)$ for the rectangular and square lattices, respectively.

###  1.2 Triangular lattice

For the triangular lattice we consider $\vec{a}_{1}=a\hat{x}$ and
$\vec{a}_{2}=a\cos(\piup/3)\hat{x}+a\sin(\piup/3)\hat{y}$ such that a point of
the lattice located at $\vec{R}_{mn}=m\vec{a}_{1}+n\vec{a}_{2}$ will lead to

$|\vec{r}-\vec{R}_{mn}|^{2}=(\varphi-\eta)(\bar{\varphi}-\bar{\eta})\,,$
(1.16)

where $\varphi$ is given by equation (1.3) and

$\eta=\left(m+\mathrm{e}^{\mathrm{i}\piup/3}n\right)a.$ (1.17)

Following the steps of the previous section,

$\displaystyle V(x,y)$
$\displaystyle=\lambda\sum_{(m,n)\in\mathbb{Z}^{2}}\ln\left[(\varphi-\eta)(\bar{\varphi}-\bar{\eta})\right]=\lambda\sum_{m=-\infty}^{\infty}\ln\left[(\varphi-
ma)(\bar{\varphi}-ma)\right]$
$\displaystyle+\lambda\sum_{m=-\infty}^{\infty}\sum_{n=1}^{\infty}\ln[(\varphi^{2}-2ma\varphi+m^{2}a^{2}+e^{-i\piup/3}n^{2}a^{2})$
$\displaystyle\times(\bar{\varphi}^{2}-2ma\bar{\varphi}+m^{2}a^{2}+\mathrm{e}^{\mathrm{i}\piup/3}n^{2}a^{2})].$
(1.18)

The analogue of equation (1.8) is then

$\displaystyle V(x,y)$
$\displaystyle=\frac{\lambda}{2}\sum_{(m,n)\in\mathbb{Z}^{2}}\ln\left(m^{2}a^{2}+\mathrm{e}^{\mathrm{i}\piup/3}n^{2}a^{2}\right)+\frac{\lambda}{2}\sum_{(m,n)\in\mathbb{Z}^{2}}\ln\left(m^{2}a^{2}+\mathrm{e}^{-\mathrm{i}\piup/3}n^{2}a^{2}\right)$
$\displaystyle+\lambda\sum_{m=-\infty}^{\infty}\ln\left\\{ma\frac{\sinh[\piup(\varphi-
ma)\mathrm{e}^{\mathrm{i}\piup/6}/a]}{\sinh(\piup
m\mathrm{e}^{\mathrm{i}\piup/6})}\right\\}$
$\displaystyle+\lambda\sum_{m=-\infty}^{\infty}\ln\left\\{ma\frac{\sinh[\piup(\bar{\varphi}-ma)\mathrm{e}^{-\mathrm{i}\piup/6}/a]}{\sinh(\piup
m\mathrm{e}^{-\mathrm{i}\piup/6})}\right\\}.$ (1.19)

After discarding the additive constants, the above expression becomes

$\displaystyle V(x,y)$
$\displaystyle=\lambda\ln[\sinh(\piup\varphi\mathrm{e}^{\mathrm{i}\piup/6}/a)\sinh(\piup\bar{\varphi}\mathrm{e}^{-\mathrm{i}\piup/6}/a)]$
$\displaystyle+\lambda\sum_{m=1}^{\infty}\ln\left\\{\frac{\sinh[\piup(\varphi-
ma)\mathrm{e}^{\mathrm{i}\piup/6}/a]\sinh[\piup(\varphi+ma)\mathrm{e}^{\mathrm{i}\piup/6}/a]}{\sinh^{2}(\piup
m\mathrm{e}^{\mathrm{i}\piup/6})}\right\\}$
$\displaystyle+\lambda\sum_{m=1}^{\infty}\ln\left\\{\frac{\sinh[\piup(\bar{\varphi}-ma)\mathrm{e}^{-\mathrm{i}\piup/6}/a]\sinh[\piup(\bar{\varphi}+ma)\mathrm{e}^{-\mathrm{i}\piup/6}/a]}{\sinh^{2}(\piup
m\mathrm{e}^{-\mathrm{i}\piup/6})}\right\\},$ (1.20)

in analogy with equation (1.9).

In terms of the coordinates $x$ and $y$, the final expression for the
regularized potential is then

$\displaystyle
V(x,y)=\lambda\ln\left\\{\frac{\vartheta_{1}\left[\mathrm{i}\frac{\piup}{a}\left(x+\mathrm{i}y\right)\mathrm{e}^{\mathrm{i}\piup/6},-\mathrm{i}\mathrm{e}^{-\sqrt{3}\piup/2}\right]\cdot\vartheta_{1}\left[\mathrm{i}\frac{\piup}{a}\left(x-\mathrm{i}y\right)\mathrm{e}^{-\mathrm{i}\piup/6},\mathrm{i}\mathrm{e}^{-\sqrt{3}\piup/2}\right]}{\cosh^{2}\left[\frac{\piup}{2a}\left(y-\sqrt{3}x\right)\right]-\cos^{2}\left[\frac{\piup}{2a}\left(x+\sqrt{3}y\right)\right]}\right\\}.$
(1.21)

A graphic representation of this function can be seen in figure 2.

Figure 2: (Colour online) Three-dimensional plot of the “potential function”
$V(x,y)$ for the triangular lattice.

Due to the linearity of equation (1.2), the above result can also be obtained
from the superposition of the potentials of two rectangular lattices displaced
relatively to each other in such a way as to form the triangular lattice (see
figure 3).

Figure 3: (Colour online) Representation of the triangular lattice as two
rectangular lattices shifted with respect to each other.

##  2 Conclusion

In this paper, we performed infinite logarithmic sums, with proper
regularization, to determine the interaction energy of rectangular and
triangular lattices of line defects having their axes along the $z$-direction.
By adjusting the defect strength $\lambda$, along with parameters governing
the geometry of a cell (namely $a$, $b$), one has the possibility to perform
defect engineering, that is tailoring material properties from controlled
defect arrays [36, 37].

Particles moving inside a lattice of topological defects may be highly
sensitive to initial conditions and hence the dynamics of these particles is
likely to lead to exponential divergence of initially closed trajectories. The
motion of fast electrons in a silicon crystal endowed with periodically
distributed atomic strings is known to be chaotic [38] and deserves a separate
treatment involving the statistical tools of dynamic hamiltonian systems. This
will be the object of a next study.

## Acknowledgements

F.M. is thankful for the financial support and warm hospitality of the
Statistical Physics Group at Université de Lorraine. This work was partially
funded by INCT nanocarbono, CNPq, CAPES and FACEPE (Brazilian agencies) and
PICS CNRS (France). We also thank Eudes Gomes for helping with figure 3.

## References

  * [1] Chen J., Chung T.C., Moraes F., Heeger A., Solid State Commun., 1985, 53, No. 9, 757–763,
doi:10.1016/0038-1098(85)90214-5

  * [2] Kosevich A., Low Temp. Phys., 2004, 30, No. 3, 247–252, doi:10.1063/1.1645186
  * [3] Lounasmaa O.V., Thuneberg E., PNAS, 1999, 96, No. 14, 7760–7767, doi:10.1073/pnas.96.14.7760
  * [4] Adhikari S., J. Phys.: Condens. Matter, 2019, 31, No. 27, 275401, doi:10.1088/1361-648x/ab14c5
  * [5] Abrikosov A.A., Sov. Phys. JETP, 1957, 5, 1174–1182.
  * [6] Murray B.S., Pelcovits R.A., Rosenblatt C., Phys. Rev. E, 2014, 90, No. 5, 052501,
doi:10.1103/physreve.90.052501

  * [7] Renn S.R., Lubensky T.C., Phys. Rev. A, 1988, 38, No. 4, 2132, doi:10.1103/physreva.38.2132
  * [8] Bahr C., Kitzerow H.S., Chirality in Liquid Crystals, Springer, New York, 2001.
  * [9] Fukuda J.i., Žumer S., Nat. Commun., 2011, 2, 246, doi:10.1038/ncomms1250
  * [10] Ackerman P.J., van de Lagemaat J., Smalyukh I.I., Nat. Commun., 2015, 6, 6012, doi:10.1038/ncomms7012
  * [11] Duzgun A., Selinger J.V., Saxena A., Phys. Rev. E, 2018, 97, No. 6, 062706, doi:10.1103/physreve.97.062706
  * [12] Tai J.S.B., Smalyukh I.I., Science, 2019, 365, No. 6460, 1449–1453, doi:10.1126/science.aay1638
  * [13] Mazur D., Heyl J.S., Phys. Rev. D, 2015, 91, No. 6, 065019, doi:10.1103/physrevd.91.065019
  * [14] Bucher M., Spergel D., Phys. Rev. D, 1999, 60, No. 4, 043505, doi:10.1103/physrevd.60.043505
  * [15] Battye R.A., Chachoua E., Moss A., Phys. Rev. D, 2006, 73, No. 12, 123528, doi:10.1103/physrevd.73.123528
  * [16] Kibble T.W., J. Phys. A: Math. Gen., 1976, 9, No. 8, 1387, doi:10.1088/0305-4470/9/8/029
  * [17] Zurek W.H., Nature, 1985, 317, No. 6037, 505, doi:10.1038/317505a0
  * [18] Bowick M., Chandar L., Schiff E., Srivastava A., Science, 1994, 263, No. 5149, 943–945,
doi:10.1126/science.263.5149.943

  * [19] Grønbech-Jensen N., Int. J. Mod. Phys. C, 1996, 7, No. 06, 873–881, doi:10.1142/s0129183196000727
  * [20] Grønbech-Jensen N., Comput. Phys. Commun., 1999, 119, No. 2–3, 115–121,
doi:10.1016/S0010-4655(99)00200-3.

  * [21] Berezinskii V.L., Sov. Phys. JETP, 1971, 32, 493.
  * [22] Kosterlitz J.M., Thouless D.J., J. Phys. C: Solid State Phys., 1972, 5, No. 11, L124,
doi:10.1088/0022-3719/5/11/002

  * [23] Kosterlitz J.M., J. Phys. C: Solid State Phys., 1974, 7, No. 6, 1046, doi:10.1088/0022-3719/7/6/005
  * [24] Kosterlitz J.M., J. Phys.: Condens. Matter, 2016, 28, No. 48, 481001, doi:10.1088/0953-8984/28/48/481001
  * [25] Berche B., J. Phys. A: Math. Gen., 2003, 36, No. 3, 585, doi:10.1088/0305-4470/36/3/301
  * [26] Lavrentiev M., Chabat B., Méthodes de la Théorie des Fonctions d’une Variable Complexe, Mir, Moscou, 1972, (in French).
  * [27] Di Francesco P., Mathieu P., Sénéchal D., Conformal Field Theory, Springer, New York, 1997\.
  * [28] Talapov A., Andreichenko V., Dotsenko V., Shchur L., Int. J. Mod. Phys. C, 1993, 4, 787–804,
doi:10.1142/s0129183193000653

  * [29] Chatelain C., Berche B., Phys. Rev. E, 1999, 60, 3853–3865, doi:10.1103/PhysRevE.60.3853.
  * [30] Glasser M., J. Math. Phys., 1974, 15, No. 2, 188–189, doi:https://doi.org/10.1063/1.1666619
  * [31] Stremler M.A., J. Math. Phys., 2004, 45, No. 9, 3584–3589, doi:10.1063/1.1777403
  * [32] Tyagi S., Phys. Rev. E, 2004, 70, No. 6, 066703, doi:10.1103/physreve.70.066703
  * [33] Tyagi S., Mol. Phys., 2006, 104, No. 3, 359–363, doi:10.1080/00268970500358735.
  * [34] Abramovitz M., Stegun I.A. (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards Applied Mathematics Series, Vol. 55, U.S. Government Printing Office, Washington, D.C., 1964.
  * [35] Olver F.W., NIST Handbook of Mathematical Functions Hardback and CD-ROM, Cambridge University Press, Cambridge, 2010.
  * [36] Fumeron S., Berche B., Moraes F., Santos F.A.N., Pereira E., Eur. Phys. J. B, 2017, 90, No. 5, 95,
doi:10.1140/epjb/e2017-70384-5.

  * [37] Santos F., Fumeron S., Berche B., Moraes F., Nanotechnology, 2016, 27, No. 13, 135302,
doi:10.1088/0957-4484/27/13/135302

  * [38] Akhiezer A., Truten’ V., Shul’ga N., Phys. Rep., 1991, 203, No. 5, 289–343,
doi:10.1016/0370-1573(91)90005-7

Ukrainian 0

Ñòîñîâíî åíåðã¿ ãðàòîê ç òîïîëîãíèìè äåôåêòàìè Á.
Áåðø${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ñ.
Ôóìåðîí${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ô.
Ìîðàåñ${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ Äèíàìêà ñèìåòðÿ, ëàáîðàòîðÿ ôçèêè òåîðåòèíî¿
õì¿, CNRS - óíâåðñèòåò Ëîððåí, UMR 7019, Âàíäóâð ë Íàíñ, Ôðàíöÿ
${}^{\textsf{\footnotesize 2}}$ Êàôåäðà ôçèêè, Ôåäåðàëüíèé óíâåðñèòåò
Ïåðíàìáóêó 52171–900 Ðåñô, Áðàçèëÿ


# Reflections in an Octagonal Mirror Maze

David Eppstein Department of Computer Science, University of California,
Irvine. Research supported in part by NSF grant CCF-2212129.

###### Abstract

Suppose we are given an environment consisting of axis-parallel and diagonal
line segments with integer endpoints, each of which may be reflective or non-
reflective, with integer endpoints, and an initial position for a light ray
passing through points of the integer grid. Then in time polynomial in the
number of segments and in the number of bits needed to specify the coordinates
of the input, we can determine the eventual fate of the reflected ray.

##  1 Introduction

There are many problems in graphics and visibility testing where it is of
interest to determine the path that would be taken by a ray of light, through
an environment that may contain mirrors. Figure 1 gives a simple example of a
problem of this type. Even for very restricted environments such as the one
depicted, where the starting point of the ray and all endpoints of mirrored
segments have integer coordinates and where the mirrors are all either axis-
aligned or at $45^{\circ}$ angles to the axes, the path of such a ray may be
very complicated, taking a number of reflections that may depend on the
geometry of the input and not merely on its combinatorial complexity. For
instance, a ray that bounces diagonally between two parallel mirrors on
opposite sides of a long thin rectangle will only exit the rectangle after a
number of bounces proportional to the aspect ratio of the rectangle, even
though such an environment consists of only two segments. Nevertheless, in
that simple two-mirror example, the eventual path of the ray is easy to
predict, without resorting to separately simulating each bounce. What about
for environments of more than two mirrors? Is it still easy to ray-trace
reflected rays in such environments?

Figure 1: The reflections of a light ray (red) among mirrors that are axis-
aligned or at $45^{\circ}$ angles to the axes (blue)

We formalize this problem as follows. The input environment is described as a
collection of line segments, with integer endpoints and either parallel to a
coordinate axis or at a $45^{\circ}$ angle to the axes. Each side of each line
segment may be marked as reflective or non-reflective. We are also given an
integer position for the start of a light ray, and an integer vector
describing the initial direction of the light ray. The restricted orientation
of the mirrors ensures that each reflection of the ray in one of the
reflective segments continues to have integer slope, on a line through
infinitely many points of the integer grid. If, after repeated reflections,
the ray eventually hits a non-reflective segment, the endpoint of a segment,
or its initial position and orientation, it stops; otherwise, it must
eventually escape the environment along an unobstructed infinite ray. The
output of the problem is the eventual fate of the ray: the point where it
stops, or the ray along which it escapes. Our main result is that we can
determine this outcome in polynomial time.

Let $n$ denote the number of segments of the input, and suppose that all of
the integers in the input specification (including the ones specifying the
initial vector of the traced ray) have magnitude at most $N$. For these
problem size parameters, it would be trivial to solve the problem in time
polynomial in $N$ – simply trace the ray one reflection at a time – but this
time bound is not polynomial, as it is exponentially larger than the input
size. Our time bound is weakly polynomial, but not strongly polynomial: it is
a polynomial of the number of bits required to specify the input, which is
$O(n\log N)$. For the purposes of polynomial-time complexity, it would be
equivalent to allow input coordinates that are rational numbers, rather than
integers, with numerators and denominators of magnitude at most $N$. Clearing
denominators would produce an integer input whose coordinates have magnitude
$N^{O(n)}$ (the product of the input numerators and denominators), still
requiring only a polynomial number of bits to specify, $O(n^{2}\log N)$.

The main idea of our algorithm is to transform the problem into one of
determining the iterated behavior of a certain type of discrete one-
dimensional dynamic system, which in a related recent paper [4] we called an
_iterated integer interval exchange transformation_. In turn, following that
paper, we can transform the iterated integer interval exchange transformation
problem into a previously-studied problem in computational topology, of
following paths along _normal curves_ on triangulated topological surfaces. To
solve this path-following problem on normal curves in triangulated surfaces,
we apply algorithms of Erickson and Nayyeri [5].

The general topic of visibility and ray-shooting with reflection is one with
much prior work, both heuristic as part of the computer graphics rendering
pipeline and with more rigorous bounds in computational geometry, for which
see, e.g., [3, 2, 9, 8, 1, 7, 6, 10]. However, this past work has a
combinatorial complexity that blows up with the number of reflections. In
contrast, our results give a polynomial time algorithm whose complexity is
independent of the number of reflections.

##  2 Iterated interval exchange transformations

In a recent paper of the author [4] we investigated a broad class of problems,
involving computing the $n$th iterate of a polynomial-time bijective function.
One motivation for this investigation was in ray-tracing problems like the one
studied here: if an environment consists only of mirrors, with no absorbing
barriers for light, then (modulo representational issues involving whether
reflections preserve the integer nature of a light ray) the mapping from each
reflected position and direction of a light ray to the next is just such a
polynomial-time bijection. Most of the problems considered in our recent paper
have high computational complexity. However, our paper also identified a
special class of bijections, the _integer interval exchange transformations_ ,
for which iterates can be computed in polynomial time. We will use the
resulting _iterated integer interval exchange transformation problem_ as a
subroutine in our algorithm for finding the result of a sequence of
reflections. In this section we summarize the definitions needed to apply this
problem, following our previous paper.

We define an integer interval exchange transformation to be a certain type of
piecewise-linear bijective mapping on a range of consecutive integers. It may
be defined by a partition of the range into subintervals, and by a permutation
of those subintervals. The transformation then translates each subinterval
(meaning that it acts on this interval by addition of the same value to each
integer in the interval), so that the translated subintervals again form a
partition of the same range, reordered into the given permutation. An example,
used in Fig. 2, is the transformation on $[0,15)$ that permutes the intervals
$a=[0,3],b=[4,5],c=[6],d=[7,14]$ into the permuted order $b,d,c,a$. This
permutation describes the function

$x\mapsto\begin{cases}x+11,&\text{for }x\in[0,3]\\\ x-4,&\text{for
}x\in[4,5]\\\ x+4,&\text{for }x\in[6]\\\ x-5,&\text{for }x\in[7,14].\\\
\end{cases}$

A transformation of this type, with $m$ intervals on the range $[0,M-1]$, can
be specified by $O(m\log M)$ bits of information, specifying the endpoints and
permuted position of each subinterval. The resulting integer function may be
evaluated on any integer $x$ in its range in time $O(m)$, by a sequential
search of the listed subintervals to find the one containing $x$, and a second
scan of the subintervals to determine which ones have permuted positions
before the one containing $x$ and contribute to the translation offset for
$x$. Even faster evaluations are possible if the intervals are stored in
sorted order with their translation offsets. The _iterated interval exchange
transformation problem_ takes as input an integer interval exchange
transformation $\mu$, represented in either of these ways, an integer $x$ in
its range, and another non-negative integer $k$. The goal is to compute
$\mu^{(k)}(x)$, the result of repeatedly replacing $x$ by its transformed
value, $k$ times.

Following a suggestion of Mark Bell, our paper [4] shows that, for every
integer interval exchange transformation, it is possible to find a
corresponding triangulated two-dimensional manifold, and a _normal curve_ on
the surface, such that tracing the normal curve for a specified number of
steps corresponds to evaluating the integer interval exchange transformation
(Figure 2). Here, a normal curve is a curve through the triangles of the
surface, avoiding triangle vertices and passing straight across each triangle
from edge to edge, without crossing itself. It can be specified by a system of
numbers on each edge counting the number of segments of the curve that cross
that edge; this specification must obey certain consistency constraints (the
numbers of crossings on the three edges of each triangle must obey the
triangle inequality and sum to an even number). This specification is
sufficient to reconstruct the curve itself, up to topological equivalence.

Figure 2: A normal curve (light blue) on a triangulated double torus (black
triangles and red vertices, glued from top to bottom and from left side to
right side with the pairing indicated by the letters). Traversing the normal
curve upwards from its central horizontal line, through the glued edges from
top to bottom, and continuing upwards back to the same central line, permutes
the branches of the curve according to the integer interval exchange
transformation that maps $[0,3]\mapsto[11,14]$, $[4,5]\mapsto[0,1]$, $6\mapsto
10$, and $[7,14]\mapsto[2,9]$. From [4].

In the transformation from our paper [4], the integers in the range $[0,M-1]$
of an integer interval exchange transformation are represented as the sequence
of $M$ crossings of a normal curve, along a central horizontal edge of a
triangulated surface. Each of these crossing points $x$ is connected by the
normal curve, along a path of exactly $s$ segments for a number $s$ determined
as part of the construction, to its image $\mu(x)$ according to the integer
interval exchange transformation. Following this construction, the iterated
interval exchange transformation problem can then be reduced to finding the
crossing point that is $sk$ steps ahead of $x$ along the normal curve. This
problem, of tracing paths for a given number of steps on a normal curve, has
been given a polynomial-time solution by Erickson and Nayyeri [5]. It follows
that the iterated interval exchange transformation problem can also be solved
in polynomial time. More precisely, the time is $O(m^{2}\log M)$, after an
initial step in which the input parameter $k$ is reduced modulo the total
number of steps in (a component of) the normal curve [4]. As an integer
division of a $\log k$-bit number by a $\log M$-bit number, this reduction
step can be performed in time $O(\log k\log M)$ using naive division
algorithms.

##  3 Partial integer interval exchange

Reflection in a mirror is a reversible transformation on systems of rays, but
absorption by a non-reflective surface is not: many different rays could be
absorbed at the same point. To mimic this non-reversibility in an integer
exchange problem, while still allowing the polynomial-time procedure from our
previous paper to apply, it is convenient to generalize the integer interval
exchange problem to allow transformations that are only partially defined, as
follows.

We define a _partial integer interval exchange transformation_ , for the range
$[0,M-1]$, to be a system of disjoint subintervals of this range, together
with a transformation that offsets each of these subintervals to another
system of disjoint subintervals (necessarily of equal lengths). For instance,
by omitting the subinterval $[6]$ from the previous example, we obtain a
partial integer exchange transformation that maps that maps
$[0,3]\mapsto[11,14]$, $[4,5]\mapsto[0,1]$, and $[7,14]\mapsto[2,9]$. The
_domain_ of the transformation is the union of the given subintervals, and the
_codomain_ is the union of their target subintervals. This example has domain
$[0,5]\cup[7,14]$ and codomain $[0,9]\cup[11,14]$.

######  Lemma 3.1.

If a partial integer interval exchange transformation is repeatedly applied to
an input $x$ that does not belong to the codomain, it eventually reaches a
transformed value that does not belong to the domain.

######  Proof 3.2.

Consider the directed graph that connects each value to its transformed image.
This graph has in-degree and out-degree at most one at each vertex, and has
finitely many vertices, so it consists of a disjoint union of directed paths
and directed cycles. An input $x$ that does not belong to the codomain has no
incoming edge, so it is the starting vertex of a path, and is eventually
transformed into the ending vertex of the same path, a value that does not
belong to the domain.

We define the _iterated partial integer interval exchange transformation
problem_ to be a computational task that takes as input the description of a
partial integer interval exchange transformation (as a system of subintervals
and their targets) and a value $x$ that does not belong to the codomain, and
that produces the corresponding value that does not belong to the domain,
according to Lemma 3.1.

######  Lemma 3.3.

The iterated partial integer interval exchange transformation problem can be
solved in time $O(m^{2}\log M+\log^{2}M)$, polynomial in the input
representation size.

######  Proof 3.4.

We transform the problem in polynomial time into an equivalent instance of the
(non-partial) iterated integer interval exchange transformation problem. Let
$I_{1},I_{2},\dots$ be the intervals of the given partial transformation $f$,
and let $f(I_{1})$ etc. denote their images after the transformation. Suppose
also that the given partial transformation operates on the range $[0,M-1]$ of
length $M$. Let $m$ denote the total number of elements in this range that are
missed by $f$: they are not in its domain. We define a new transformation
$\bar{f}$ that operates on the range $[0,Mm+M+m-1]$ of length $Mm+M+m$, as
follows:

  * •

For each subinterval $I_{i}$ in the given transformation, with image
$f(I_{i})$, we include in $\bar{f}$ the mapping $I_{i}\mapsto f(I_{i})$.

  * •

For each maximal subinterval $J_{i}$ of $[0,M-]\setminus\cup I_{i}$ (a
subinterval not in the domain of $f$) we include in $\bar{f}$ a mapping from
$J_{i}$ to a subinterval of $[M,M+m-1]$, so that the images of these
subintervals are disjoint and completely cover $[M,M+m-1]$.

  * •

We include in $\bar{f}$ the mapping $[M,Mm+M-1]\mapsto[M+m,Mm+M+m-1]$.
Iterating this mapping eventually transforms any value in $[M,M+m-1]$ to a
value in $[Mm+M,Mm+M+m-1]$, but it takes $M$ iterations to do so.

  * •

For each maximal subinterval $K_{i}$ of $[0,M-]\setminus\cup f(I_{i})$ (a
subinterval not in the codomain of $f$) we include in $\bar{f}$ a mapping from
a subinterval of $[Mm+M,Mm+M+m-1]$ to $K_{i}$, so that the preimages of these
mappings are disjoint and completely cover $[Mm+M,Mm+M+m-1]$.

For instance, for the example partial integer interval exchange transformation
$f$ given above, $M=15$ and $m=1$ (there is only one missing value from the
transformation), and we have $\bar{f}$ mapping $[0,3]\mapsto[11,14]$,
$[4,5]\mapsto[0,1]$, $[7,14]\mapsto[2,9]$; $[6]\mapsto[15]$;
$[15,29]\mapsto[16,30]$; and $[30]\mapsto[10]$.

Suppose we apply the algorithm to the iterated integer interval exchange
transformation problem, with transformation $\bar{f}$, on an input value $x$
that does not belong to the codomain, and that the output of this algorithm is
a value $z$. If we iterate $\bar{f}$ for a total of $M$ iterations, starting
with a value $x$ that does not belong to the codomain, it will reach a value
$y$ that does not belong to the domain in fewer than $M$ iterations, by Lemma
3.1. The next iteration will map $y$ into a value $y^{\prime}$ in the
subinterval $[M,M+m-1]$, and the subsequent (again fewer than $M$) iterations
will each add $m$ to this value $y^{\prime}$. We may therefore obtain
$y^{\prime}$ by $z$ as the unique value in the subinterval $[M,M+m-1]$ that is
congruent to $z$ modulo $m$. From $y^{\prime}$, we may obtain the desired
value $y$ as $\bar{f}^{-1}(y^{\prime})$.

The time bound is obtained by plugging in the number of pieces of the
resulting transformation, $O(m)$, the range of its values, $O(Mm)$, and the
number of iterations, $O(M)$, into the previous time bound for iterated
integer interval exchange transformations.

##  4 Converting reflection to partial integer interval exchange

The reason that we have restricted our attention to reflections in line
segments that are horizontal, vertical, and diagonal is that these kinds of
reflections preserve the integrality of the reflected rays. We formalize these
observations as follows.

######  Lemma 4.1.

If a ray whose direction is specified by a vector $(x,y)$ is reflected by a
sequence of horizontal, vertical, or diagonal mirrors, then the resulting
ray’s direction can be specified by one of the eight vectors $(\pm x,\pm y)$
or $(\pm y,\pm x)$.

Figure 3: Cases for reflection of a ray from a horizontal, vertical, or
diagonal mirror

######  Proof 4.2.

Vertical mirrors negate the first coordinate, horizontal mirrors negate the
second coordinate, and diagonal mirrors exchange the two coordinates (possibly
also negating both of them); see Fig. 3. The set of eight vectors of the lemma
are preserved by these operations, so no other direction is possible.

Along with these eight directions, it will also be important to keep track of
the left-right orientation of each ray; that is, whether an image that follows
a thickened copy of the ray has its orientation preserved or flipped. This
changes on each reflection, and cannot be determined only from the ray’s
direction: a ray that reflects perpendicularly from a mirror will have its
orientation changed, whereas a ray that returns from a corner-reflector (two
perpendicular mirrors) will not.

######  Lemma 4.3.

If a ray that passes through infinitely many points of the integer grid is
reflected by a sequence of horizontal, vertical, or diagonal mirrors, each
with integer endpoints, then the reflected ray again passes through infinitely
many points of the integer grid.

######  Proof 4.4.

When a ray is reflected by a mirror, it passes after the reflection through
the same sequence of grid points that the unreflected ray would pass through
in the reflection of the grid. But with mirrors of the type described by the
lemma, the reflection of the grid is an integer grid with the same points.

Figure 4: Subdivision of grid segments and diagonals (blue) by the system of
all lines through integer points in the direction of the vector $(5,3)$ (red).

######  Observation 4.5.

Let $\mathcal{L}$ be the system of all lines in the plane that pass through
integer points in the direction of a vector $(x,y)$, where $x$ and $y$ are
integers in lowest terms (their greatest common divisor is one). Then the
lines of $\mathcal{L}$ subdivide each vertical unit segment of the integer
grid into $x$ equal-length pieces, and each horizontal unit segment of the
grid into $y$ equal-length pieces. They subdivide each diagonal segment of the
grid, of length $\sqrt{n}$ with slope of the same sign as the slope of the
lines in $\mathcal{L}$, into $|x-y|$ equal-length pieces, and each diagonal
segment of the opposite slope into $|x+y|$ equal -length pieces. (See Fig. 4.)

Putting these observations together, we can transform any octagonal mirror
maze into an equivalent partial integer interval exchange transformation.

######  Lemma 4.6.

Suppose we are given an environment, described as a collection of line
segments, each side of which may be marked as reflective or non-reflective,
with integer endpoints, an integer position for the start of a light ray, and
an integer vector describing the initial direction of the light ray. Then in
polynomial time we can construct a partial integer interval exchange
transformation $f$ whose values correspond to points of the environment
(either the starting point or points along the segments of the environment)
and directions of a reflected light ray emanating from that point, such that:

  * •

If a value $v$ belongs to the domain of $f$, then the ray described by $v$ is
eventually reflected by the environment, with its first reflection at a
position and direction described by $f(v)$.

  * •

If a value $v$ does not belong to the domain of $f$, then the ray described by
$v$ hits an absorbing barrier before being reflected, or escapes to infinity.

######  Proof 4.7.

We surround the given environment with a non-reflective bounding box that will
catch all escaping rays. By Lemma 4.1 and Lemma 4.3, we need only consider
rays through integer points, in eight directions and two left-right
orientations. By 4.5, we need only consider a discrete evenly-spaced set of
possible reflection points along each segment of the environment, of a size
that can be described by an integer with polynomially many bits. (Recall that
we are not restricting the lengths of our grid segments or our initial
direction to have polynomial magnitude, only to be integers with a polynomial
number of bits.) We will create a partial integer interval exchange
transformation $f$ whose range is partitioned into subintervals, one for each
combination of an environment segment, a direction of the emanating ray, and a
left-right orientation, with the length of these subintervals obtained by
multiplying the length of the segment by the number of reflection points per
unit length given by 4.5. For each of these subintervals, we further partition
it into sub-subintervals, in polynomial time, according to the next object in
the environment that a ray in that direction would hit, as depicted in Fig. 5.
(This is simply a lower envelope of $n$ disjoint line segments, for a
projection direction determined by the ray direction.) When the next object to
be hit is reflecting, we map a sub-subintervals to the interval describing the
corresponding reflection points along that object. When, instead, it is
absorbing, we omit that sub-subinterval from the partial integer interval
exchange transformation.

Figure 5: Reflection points along segment $AB$ of a given environment,
partitioned into subintervals according to the next object in the reflected
path for incoming rays of slope $-\tfrac{1}{2}$.

##  5 The main result

Putting these components together, we have the following result:

######  Theorem 5.1.

Suppose we are given an environment, described as a collection of line
segments, each side of which may be marked as reflective or non-reflective,
with integer endpoints, an integer position for the start of a light ray, and
an integer vector describing the initial direction of the light ray. Then in
time polynomial in the number of segments and in the number of bits needed to
specify the integers of the input, we can determine whether the reflected ray
is eventually absorbed or escapes to infinity. If it is absorbed, we can
determine where it is absorbed, what direction it comes from when it is
absorbed, and how many bounces it makes before this happens. If it escapes to
infinity, we can determine its eventual escape path, and how many bounces it
takes before reaching this path. The time bound for these algorithms is
$O(n^{2}\log N+\log^{2}N)$.

######  Proof 5.2.

We convert the input to an equivalent partial integer interval exchange
transformation according to Lemma 4.6, and then apply the polynomial-time
algorithm for the iterated partial integer interval exchange transformation
problem of Lemma 3.3.

An input of size $n$ and coordinate magnitude $N$ can be converted into a
partial integer interval exchange transformation whose number of subintervals
is $O(n)$ (each comes from a trapezoid in four trapezoidal decompositions with
sides parallel to one of the directions of the reflected rays) and whose range
is $O(N^{2})$ (combining the magnitude of the environment coordinates with the
number of reflection points along each grid segment). For inputs of this size,
the time to apply an algorithm for the iterated partial integer interval
exchange transformation is $O(n^{2}\log N+\log^{2}N)$.

It would be of interest to determine to what extent this algorithm can be
generalized. Can we determine the geometric length of the reflected path of a
light ray, and not just its number of bounces? Are there other systems of
slopes, beyond the axis-aligned and diagonal slopes, for which similar
algorithms can work? How does the complexity of the algorithm depend on the
system of slopes? For slopes that do not preserve the rationality of reflected
rays, what can be said about the computational complexity of the problem?

## References

  * [1] Mridul Aanjaneya, Arijit Bishnu, and Sudebkumar Prasant Pal.  Directly visible pairs and illumination by reflections in orthogonal polygons.  In Sylvain Petitjean, editor, Collection of abstracts of the 24th European Workshop on Computational Geometry, pages 241–244. INRIA-LORIA, March 18–20 2008.  URL: https://physbam.stanford.edu/~aanjneya/mridul_files/papers/vis.pdf.
  * [2] Boris Aronov, Alan R. Davis, Tamal K. Dey, Sudebkumar P. Pal, and D. Chithra Prasad.  Visibility with multiple reflections.  Discrete & Computational Geometry, 20(1):61–78, 1998.  doi:10.1007/PL00009378.
  * [3] Boris Aronov, Alan R. Davis, Tamal K. Dey, Sudebkumar P. Pal, and D. Chithra Prasad.  Visibility with one reflection.  Discrete & Computational Geometry, 19(4):553–574, 1998.  doi:10.1007/PL00009368.
  * [4] David Eppstein.  The complexity of iterated reversible computation.  Electronic preprint arxiv:2112.11607, 2021.
  * [5] Jeff Erickson and Amir Nayyeri.  Tracing compressed curves in triangulated surfaces.  Discrete & Computational Geometry, 49(4):823–863, 2013.  doi:10.1007/s00454-013-9515-z.
  * [6] Subir Kumar Ghosh, Partha P. Goswami, Anil Maheshwari, Subhas C. Nandy, Sudebkumar Prasant Pal, and Swami Sarvattomananda.  Algorithms for computing diffuse reflection paths in polygons.  The Visual Computer, 28(12):1229–1237, 2012.  doi:10.1007/s00371-011-0670-z.
  * [7] Salma Sadat Mahdavi, Ali Mohades, and Bahram Kouhestani.  Computing $k$-link visibility polygons in environments with a reflective edge.  In Proceedings of the 23rd Annual Canadian Conference on Computational Geometry, Toronto, Ontario, Canada, August 10-12, 2011, 2011.  URL: https://www.cccg.ca/proceedings/2011/papers/paper63.pdf.
  * [8] Joseph O’Rourke and Octavia Petrovici.  Narrowing light rays with mirrors.  In Proceedings of the 13th Canadian Conference on Computational Geometry, University of Waterloo, Ontario, Canada, August 13-15, 2001, pages 137–140, 2001.  URL: https://www.cccg.ca/proceedings/2001/orourke-13443.ps.gz.
  * [9] D. Chithra Prasad, Sudebkumar P. Pal, and Tamal K. Dey.  Visibility with multiple diffuse reflections.  Computational Geometry: Theory and Applications, 10(3):187–196, 1998.  doi:10.1016/S0925-7721(97)00021-7.
  * [10] Arash Vaezi and Mohammad Ghodsi.  Visibility extension via mirror-edges to cover invisible segments.  Theoretical Computer Science, 789:22–33, 2019.  doi:10.1016/j.tcs.2019.02.011.