\section{Introduction}
Kondo insulators or Kondo semiconductors (KSs) are strongly correlated electron systems where a narrow (pseudo)gap opens at low temperatures in the electronic density of states (DOS) at the Fermi level~\cite{Aeppli92}.
This is a result of the hybridization of localized $f$ electrons with the conduction $p$ or $d$ electrons~\cite{Riseborough00}.
Even though known since decades~\cite{Aeppli92}, KSs remain at the forefront of condensed matter research.
Recently it has been suggested that KSs may be candidates for strongly correlated materials with topologically protected surface states~\cite{Dzero10,Dzero12,Alexandrov13}.
So far, experiments along these lines have focused on the cubic KS SmB$_6$~\cite{Cooley99,Kim14}. \\
Here CeRu$_4$Sn$_6$ comes into play.
It has a tetragonal, non-centrosymmetric crystal structure ($I\bar{4}2m$)~\cite{Poettgen97} and therefore can shed light on the role crystal symmetry plays.
It was first synthesized in 1992~\cite{Das92} and a development of a narrow gap was found in several physical properties~\cite{Strydom05,Bruening10,Paschen10}.
Measurements on single crystals revealed large anisotropy in resistivity, specific heat, and thermopower~\cite{Paschen10,Winkler12,Haenel14}.
Optical conductivity measurements and LDA+DMFT calculations attributed this to a gapping in the tetragonal plane while the out-of-plane direction stays weakly metallic~\cite{Guritanu13}. \\
Here we report on the thermal conductivity $\kappa$ of CeRu$_4$Sn$_6$ single crystals, both along and perpendicular to the out-of-plane direction $c$, in the temperature range of 80\,mK to 80\,K.

\section{Experimental}
\label{exp}
The polycrystalline starting material for the crystal growth was synthesized by melting first Ce and Ru, and subsequently the resulting Ce/Ru alloy with Sn in a horizontal water cooled copper boat using high-frequency heating.
The purity of the starting materials was 99.99\% for Ce (Ames Lab) and Ru, and 99.9999\% for Sn.
The single crystal growth was performed by the floating zone melting technique using optical heating in a four-mirror furnace of Crystal Corporation.
All steps were performed under Ar 99.9999\% protective atmosphere after purging several times with Ar.
Since CeRu$_4$Sn$_6$ melts incongruently we used a self-flux technique \cite{Paschen10,Prokofiev12}. \\
The low-temperature (80\,mK~-~2\,K) thermal conductivity was measured in a $^3$He/$^4$He dilution refrigerator with a home-built setup employing a two-thermometer-one-heater technique using RuO$_2$ thick-film resistors as heater and thermometers.
The measurements in the range 2\,K to 80\,K were carried out in a $^4$He bath cryostat with thermocouples to measure the temperature gradient.

\section{Physical properties}
\label{phys}
Figure\,\ref{fig:k} shows the thermal conductivity $\kappa$ as a function of temperature $T$ in the whole temperature range investigated, on both a linear (main panel) and a double-logarithmic (inset) scale.
$\kappa$ was measured along the out-of-plane direction $c$ (called $\kappa_c$) and within the tetragonal plane, perpendicular to $c$ ($\kappa_{\perp c}$).
We calculated the electronic part of the thermal conductivity $\kappa^{\rm{WF}}$ from the electrical resistivity $\rho$\,\cite{Winkler12} employing the Wiedemann-Franz law $\kappa^{\rm{WF}} = L_0 T / \rho $ with the Lorenz number $L_0=2.44 \times 10^{-8}$\,W$\Omega$K$^{-2}$.
$\kappa$ is essentially isotropic, in marked contrast to other physical properties such as electrical resistivity~\cite{Winkler12}, thermopower~\cite{Haenel14}, and optical conductivity~\cite{Guritanu13}. \\
At low temperatures $\kappa$ increases strongly with increasing temperature towards a maximum of roughly $\kappa_{\rm{max}} \approx$ 6.3\,Wm$^{-1}$K$^{-1}$ at $T_{\rm{max}} = $ 28\,K.

\begin{figure}[h]
\centering
\includegraphics[width=28pc]{Fig1-k+kWF-wlogInset.eps}
\caption{\label{fig:k}Thermal conductivity $\kappa$ of CeRu$_4$Sn$_6$ as a function of temperature $T$ measured parallel and perpendicular to the out-of-plane direction $c$. The electronic part of the thermal conductivity $\kappa^{\rm{WF}}$, as obtained from resistivity $\rho$~\cite{Winkler12} via the Wiedemann-Franz law, is plotted as solid lines. The inset shows the same data on a double-logarithmic scale.}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=28pc]{Fig2-k-ph-wCallaway-lin.eps}
\caption{\label{fig:kCallaway}Phononic part of the thermal conductivity $\kappa_{\rm{ph}}$ of CeRu$_4$Sn$_6$ along and perpendicular to the $c$ direction as a function of temperature $T$. Solid lines are Callaway fits to the data (see text).}
\end{figure}
$\kappa^{\rm{WF}}$ is much lower than the total $\kappa$ and thus $\kappa$ is clearly phonon-dominated as expected in KSs.
We estimate the phononic contribution $\kappa_{\rm{ph}}$ as $\kappa_{\rm{ph}} = \kappa - \kappa^{\rm{WF}}$.
The result is plotted in Fig.\,\ref{fig:kCallaway}, together with least-square fits to the Callaway model~\cite{Callawy59}
\begin{equation}
\kappa_{ph} = \frac{k_B^4 \cdot T^3}{2 \cdot \pi^2 \cdot v_g \cdot \hbar^3} \cdot \left[\int^{\frac{\theta_D}{T}}_0 {\tau \cdot \frac{x^4 \cdot e^x}{\left(e^x-1\right)^2} \cdot dx} \right] \quad ,
\label{eq:callaway}
\end{equation}
where $v_g$ is the sound velocity, $\theta_D$ is the Debye temperature, $\tau$ is the phonon relaxation rate, $x = \frac{\hbar \omega}{k_B T}$, and the other symbols have the usual meaning.
Different scattering mechanisms influence $\tau$ according to the Matthiesen's rule
\begin{equation}
\frac{1}{\tau} = \sum_i{\frac{1}{\tau_i}} \quad.
\label{eq:matt}
\end{equation}
The various scattering mechanisms have different dependencies on the frequency $\omega$ and subsequently lead to different temperature dependencies of $\kappa_{\rm{ph}}$.
We have taken phonon-electron, -boundary, -defect and Umklapp scattering processes into account.
Note that we have neglected N processes here which would lead to a correction term in Eqn.\,\ref{eq:callaway}.
However the influence of these processes is expected to be small in materials with complex crystal structure.
Table\,\ref{tab:fits} lists the results of the fits. \\
\begin{center}
\begin{table}[h]
\caption{\label{tab:fits}Fit parameters of the Callaway fits, where the scattering rates are phonon-electron scattering $\tau^{-1}_{\rm{pe}} = A \cdot \omega^2$, phonon-boundary scattering $\tau^{-1}_{\rm{b}} = B$, phonon-defect scattering $\tau^{-1}_{\rm{pd}} = C \cdot \omega^4 $, and Umklapp scattering $\tau^{-1}_{\rm{U}} = D \cdot \omega^2 \cdot T \cdot exp(\frac{\theta_D}{-3T})$. The values in the table are the constants $A$, $B$, $C$, and $D$.}
\centering
\begin{tabular}{lcccc}
\hline
Direction & $A$ & $B$ & $C$ & $D$\\
  & (s) & (s$^{-1}$) & (s$^3$) & (sK$^{-1}$)\\
\hline
$c$ & \quad $2.2 \times 10^7$ \quad & \quad $2.8 \times 10^6$ \quad  & \quad 7250 \quad & \quad $9.16 \times 10^4$ \\
$\perp c$ & \quad $2.5 \times 10^7$ \quad   & \quad $3.0 \times 10^6$ \quad  & \quad 6878  \quad  & \quad $9.06 \times 10^4$ \\
\hline
\end{tabular}
\end{table}
\end{center}
\begin{figure}[h]
\centering
\includegraphics[width=28pc]{Fig3-k-ph-wCallaway-log.eps}
\caption{\label{fig:kCallawayLog}Phonon thermal conductivity $\kappa_{\rm{ph}}$ of CeRu$_4$Sn$_6$ along and perpendicular to the $c$ direction as a function of temperature $T$ on a double-logarithmic scale. Solid red and blue lines are fits to the data according to the Callaway model~\cite{Callawy59}. Black lines indicate different power laws: $T^3$ (solid), $T^2$ (dashed), and linear-in-$T$ (dotted).}
\end{figure}
Figure\,\ref{fig:kCallawayLog} shows the same data as Fig.\,\ref{fig:kCallaway} on a double-logarithmic scale, which reveals the different scattering contributions most clearly.
Straight lines indicate different power laws.
The solid line shows $T^3$ behaviour typical for phonon-boundary scattering.
The dashed line indicates $T^2$ behaviour which is seen for phonon-electron scattering, and the dotted line is linear in $T$. \\
The data follow the dashed line very closely in the range of 200\,mK to 2\,K.
Thus the dominant scattering mechanism in this range is identified as phonon-electron scattering.
At lower temperatures $\kappa_{\rm{ph}}$ shows a $T^3$ behaviour.
The fit parameters of the Callaway fits give the boundary scattering rates $\tau_{\rm{b},\perp c}^{-1} = 2.8 \times 10^6\,\rm{s}^{-1}$ and $\tau_{\rm{b},c}^{-1} = 3.0 \times 10^6\,\rm{s}^{-1}$ for the measurement perpendicular and parallel to $c$, respectively.
With a Debye temperature $\theta_D=250$\,K~\cite{Bruening10} the sound velocity can be computed to be $v_g = 2316$\,m/s.
With this $v_g$ the scattering times correspond to the mean free paths $l_{\perp c} = 0.83$\,mm and $l_{c} = 0.77$\,mm.
This is close to the smallest dimension of the samples and we therefore conclude that the phonons are predominantly scattered off the sample boundaries in this temperature range. \\
The onset of phonon-electron scattering, with the scattering rate $\tau_{\rm{pe}}^{-1}$, at higher temperatures leads to a further reduction of the thermal conductivity. According to Pippard's theory~\cite{Pippard55}, $\tau_{\rm{pe}}^{-1}$ is proportional to the charge carrier concentration $n$.
From Hall effect measurements $n$ was determined to be isotropic at both high and low temperatures~\cite{WinklerTh13} in agreement with the isotropic behaviour of $\kappa$ in these temperature ranges. However, at intermediate temperatures, signatures of the Kondo gap opening are anisotropic ~\cite{Guritanu13}. This may explain the small anisotropy seen in the thermal conductivity around 10\,K (hump in $\kappa_{\perp c}$, see Fig.\,\ref{fig:kCallaway}). As the gap opening is more pronounced for the plane perpendicular to $c$, the effect of electrons freezing out and thus disappearing as scattering partners for the phonons is more pronounced for $\kappa_{\perp c}$, which leads to slightly enhanced $\kappa_{\perp c}$ values. Similar anisotropies have been found in the KSs CeNiSn and CeRhSb and have likewise been attributed to anisotropies in the energy gap~\cite{Paschen00,Sera97,Isikawa91}.

Above 2 K, $\kappa_{\rm{ph}}$ deviates from the $T^2$ law to lower values. Here scattering on point defects and U processes become the dominant effects. The values of the prefactors for these processes, as listed in Tab.\,\ref{tab:fits}, are similar to what was found in other systems~\cite{Yan13}. However, for a quantitative discussion it is necessary to know the phonon band structure which is not available in this system. In addition, other aspects that play an important role in CeRu$_4$Sn$_6$ and are neglected in a simple Callaway model may come into play at elevated temperatures.

\section{Conclusion}
\label{con}
In summary we have investigated single-crystalline CeRu$_4$Sn$_6$ by means of thermal conductivity measurements in the temperature range of 80\,mK to 80\,K.
The data are essentially isotropic in contrast to findings for other physical properties. The reason is that the thermal conductivity is dominated by its phonon contribution, even at the lowest temperature where the behaviour
can be well explained by boundary and phonon-electron scattering.
This gives further support to the previous conjecture that the strongly anisotropic electrical resistivity is mostly due to anisotropies in the charge carrier mobility, with the charge carrier concentration being essentially isotropic.

\section{Acknowledgement}
We gratefully acknowledge financial support from the Austrian Science Fund (FWF projects I 623-N16, W1243, and I 2535-N27) and from the European Research Council (ERC grant 227378).
VM acknowledges the FAPERJ (Nota 10).


\section{Introduction}

Deformable registration is a major challenge in medical image processing.
The result is a dense mapping showing pixel-wise non-linear correspondences between a pair of images that best aligns the input image
$ I $ onto the target image $ T $ by means of some similarity definition $ \mathcal{L} $.
Deformable registration is applied in the analysis of patient-specific temporal or anatomical changes, e.g. from pre-operative to post-operative state, or to show inter-patient variances \cite{Sotiras2013}.
Deformable registration is also performed in atlas-based segmentation, where an input image is matched onto a target image with known segmentation \cite{Cabezas2011}.

Existing registration methods can be separated into two categories.
The first category is based on non-learning methods which estimate a registration $ w $ by optimizing a cost function of the form
\begin{equation}
    \argmin_{w} \Big\{ \mathcal{L}(T, w \circ I) + \lambda \mathcal{R}(w) \Big\} ~ ,
    \label{eq:cost}
\end{equation}
where $ w \circ I $ denotes $ I $ warped by $ w $.
A common assumption of $ w $ is a displacement or velocity vector field $ \bm{u}(\bm{x}) $.
The final deformation results in $ \bm{\phi}(\bm{x}) = \bm{x} + \bm{u}(\bm{x}) $ which maps every pixel coordinate $ \bm{x} $ to other pixel coordinates.
The first term in (\ref{eq:cost}) is referred to as data term, which is typically chosen to be a pixel intensity error measure. Optimization of the data term alone is considered ill-posed.
The second term $ \mathcal{R} $, weighted by trade-off factor $ \lambda $,  is a regularizer that shapes the registration by any chosen prior, which helps solving the ill-posed problem.
Common regularization is done by enforcing smoothness onto the displacement vector field by penalizing first or higher order spatial derivatives of $ \bm{u} $ \cite{Werlberger2010}.
The result of the registration algorithm heavily depends on the cost function and therefore on the chosen prior of $ \mathcal{R} $.

The second category implicitly learns the regularization prior by training a convolutional network on a large database of domain-specific images.
Early approaches rely on ground truth registrations \cite{Sokooti2017}, which are hard to obtain especially in medical imaging.
More recent methods \cite{Balakrishnan2019} propose unsupervised registration using the spatial transformer function \cite{Jaderberg2015}.
However, these methods either only support small displacements or require segmentation maps of the image pairs during training to assist the convergence \cite{Hu2018}.
Additionally, the trained networks are limited to register images from the training domain (e.g. CT or MRI).

Inspired by the idea of deep image priors \cite{Lempitsky2018}, we subsequently propose our learning-free method for deformable medical image registration using the structure of an untrained convolutional network as regularization prior.

\section{Methods}

\begin{figure}[t]
    \centering
    \includegraphics[width=1\textwidth]{overview.pdf}
    \caption{Overview of our method. The randomly-initialized generator network $ f_\theta $ acts as parameterization of the registration field $ \bm{\phi} $. The parameters $ \theta $ are optimized for every image pair individually by gradient descent.}
    \label{fig:overview}
\end{figure}

Lempitsky et al. have recently shown that excellent performance of CNNs for inverse image problems, such as denoising, is not only based on their ability to learn image priors from data, but is also based on the structure of a convolutional image generator itself \cite{Lempitsky2018}.
They gave evidence that the structure of a network alone is sufficient to capture enough image statistics to provide state-of-the-art performance in inverse image tasks.

Leveraged by this idea, we reformulate the task of deformable image registration by using the structure of a convolutional network as regularizer (see Fig.~\ref{fig:overview}).
An image generator network $ \bm{u} = f_{\theta}(z) $ with randomly-initialized parameters $ \theta $ is interpreted as parameterization of the dense displacement field $ \bm{u} \in \mathbb{R}^{2 \times H \times W} $ from which the deformation $ \bm{\phi} = \bm{x} + \bm{u} $ between an input image $ I \in \mathbb{R}^{C \times H \times W} $ and a target image $ T \in \mathbb{R}^{C \times H \times W} $ can be obtained by adding to the identity warp $ \bm{x} $.
The input $ z \in \mathbb{R}^{C' \times H \times W} \sim \mathcal{N}(0, 0.1) $ has the same spatial dimensions as $ \bm{\phi} $ and is sampled from a random normal distribution in every iteration.
This leads to the following optimization problem
\begin{equation}
    \argmin_{\theta} \Big\{ \mathcal{L} \left( T, \left( \bm{x} + f_{\theta}(z) \right) \circ I \right) \Big\} ~ ,
    \label{eq:task}
\end{equation}
where $ \left( \bm{x} + f_{\theta}(z) \right) \circ I $ denotes the differentiable spatial transformer function \cite{Jaderberg2015}. Eq.
(\ref{eq:task}) is optimized for every image pair $ \left\{ I, T \right\} $ using the Adam gradient descent optimizer \cite{Kingma2014}.
As data term, we chose pixel-wise mean absolute error
$ \mathcal{L}(T, \bm{\phi} \circ I) = \vert (\bm{\phi} \circ I) - T \vert $.
The architecture of the image generator network $ f_{\theta} $ is chosen according to \cite{Lempitsky2018}.
It has an encoder-decoder structure with skip connections between the encoding and decoding part.
To begin the optimization from close to an identity warp, we initialize the parameters with $ \theta \sim \mathcal{N} (0, 0.01) $.

\section{Results \& Conclusion}

\begin{figure}[t]
    \centering
    \includegraphics[height=5cm]{results_img.pdf} ~
    \includegraphics[height=5cm]{results-crop.pdf}
    \caption{Results of our approach compared to a state-of-the-art method from the Insight ToolKit (ITK). Left: Two example MRI pairs from the data set. Right: Boxplots of the means of $ \det ( J_{\bm{\phi}} ) $ and SSIM between $ T $ and $ \bm{\phi} \circ I $.}
    \label{fig:results}
\end{figure}
\enlargethispage{1\baselineskip}

We demonstrate our approach on the task of 2D brain magnetic resonance imaging (MRI) registration.
The data used in this work contain 109 pairs of MRI scans from The Cancer Genome Atlas \cite{NCI} showing lower-grade gliomas.
We use the structural similarity index (SSIM) \cite{Wang2004} between $ \bm{\phi} \circ I $ and $ T $ and the mean of the determinants of Jacobians $ \det ( J_{\bm{\phi}} ) $ \cite{Ashburner2007} of the deformation as evaluation metrics.
The latter metric shows regularity of $ \bm{\phi} $.
We compare our method to state-of-the-art methods from the Insight ToolKit (ITK) registration framework by combining an initial affine registration and a subsequent deformable displacement field registration \cite{Avants2012}.
Results for exemplary image pairs and boxplots of results for all image pairs are shown in Fig.\,\ref{fig:results}.
Additional results including registration fields are shown in appendix \ref{sec:results_app}.

The results reveal that the structure of a convolutional network can act as regularization in deformable medical image registration with state-of-the-art performance.
This connects traditional non-learning methods and learning-based methods by using randomly-initialized convolutional networks as prior.

\midlacknowledgments{This research has received funding from the European Union as being part of the EFRE OPhonLas project.}


\section{Introduction}\label{intro}

Following \cite{nemethi}, we review some terminology of singularity theory. Let $f:(\C^3,0)\to (\C,0)$ be the germ of a complex analytic function and let
\[K=f^{-1}(0)\cap S_{\varepsilon}\]
be the \emph{link of the singularity}, where $S_{\varepsilon}$ is the sphere of radius $\varepsilon$ centered at the origin.
Milnor proved in \cite{milnor} that there exists $\varepsilon_0>0$ such that $\forall\, 0<\varepsilon<\varepsilon_0$ the map
\[f/|f|:S_{\varepsilon}\smallsetminus K\to S^1=\{z\in\C:\,|z|=1\}\]
is a smooth fibration. For any such $\varepsilon$ there exists $\delta_{\varepsilon}$ so that $\forall\, 0<\delta<\delta_{\varepsilon}$ the restriction
\[
f:B_{\varepsilon}\cap f^{-1}(\partial D_{\delta}) \to \partial D_{\delta}
\]
is a smooth fibration, whose diffeomorphism type does not depend on $\varepsilon$ and $\delta$. This is what we refer to as the \emph{Milnor fibration} of $f$. The fiber
\[F=F_{\varepsilon,\delta}=B_{\varepsilon}\cap f^{-1}(\delta)\]
is the \emph{Milnor fiber} of $f$ (we omit $\varepsilon$ and $\delta$ from the notation of $F$). If $f$ is the germ of an \emph{isolated} singularity, then we have a diffeomorphism $\partial F\simeq K$, but in the case of non-isolated singularity $K$ is not smooth (while $F$ and $\partial F$ are always smooth manifolds). The boundary of the Milnor fiber comes with an extra structure, as explained below. Recall that:

\begin{ndef}
A contact structure on a 3-manifold $M$ is a nowhere integrable planes distribution $\xi$. If there exists an embedded disk $D\hookrightarrow M$ such that $\xi$ agrees with $TD$ along the boundary $\partial D$, then the contact structure $\xi$ is said to be \emph{overtwisted}, otherwise it is called \emph{tight}. A tight structure $\xi$ is called \emph{universally tight} if its pullback to the universal cover $\widetilde{M}$ is tight. The tight structure $\xi$ is called \emph{virtually overtwisted} if its pullback to some finite cover is overtwisted.
\end{ndef}

\begin{nrem}
A consequence of the geometrization conjecture is that the fundamental group of any 3-manifold is \emph{residually finite} (i.e. any non trivial element is in the complement of a normal subgroup of finite index), and this implies that any tight contact structure is either universally tight or virtually overtwisted (see \cite{honda}).
\end{nrem}

Contact topology shows up in singularity theory in the following way: the Milnor fiber $F$ of a singularity comes with a Stein structure $J$ which makes it a Stein filling of its boundary $\partial F$ equipped with the contact structure
\[\xi=T\partial F\cap JT\partial F,\]
which is always tight (see \cite{eliashberg}). Therefore, from a complex germ $f:(\C^3,0)\to (\C,0)$ we obtain a contact 3-manifold $(\partial F,\xi)$ with a Stein filling $(F,J)$ of it. We have a dichotomy:

\begin{itemize}
\item the singularity is isolated.
In this case the structure $\xi$ is universally tight (see \cite{lekili}). Another work on this topic is \cite{dynamics}, where the authors show that the link of the hypersurface singularity
\[z^p+2xy=0\]
is $L(p,p-1)$ with its unique tight contact structure (universally tight). We will prove (see corollary \ref{L(p,p-1)}) that this is the only lens space arising as the link of an isolated hypersurface singularity.

\item The singularity is not isolated. By contrast with previous point, this is the only case where a virtually overtwisted contact structure can arise. A good source of examples is given by the \emph{Hirzebruch singularity}
\[z^2+xy^n=0\]
with $n>1$, for which the boundary of the associated Milnor fiber is $L(2n,1)$. This type of singularity is studied in \cite[section 6]{pichon}.
\end{itemize}

\noindent Given a pair of coprime integers $p>q>1$, we consider the continued fraction expansion
\[\frac{p}{q}=[a_1,a_2,\ldots, a_n]=a_1-\frac{1}{a_2-\frac{1}{\ddots -\frac{1}{a_n}}}\]
with $a_i\geq 2$ for every $i$. As a smooth oriented 3-manifold, the lens space $L(p,q)$ is the integral surgery on a chain of unknots with framings $-a_1,-a_2,\ldots, -a_n$ (figure \ref{chain}).

The goal of this article is to study those lens spaces $L(p,q)$ with a tight contact structure $\xi$ arising as the boundary of the Milnor fiber of a hypersurface singularity $f:(\C^3,0)\to (\C,0)$. Theorem \ref{mainthm} gives a partial answer to a question raised by \cite[open problems 24.4.2]{nemethi}.

\begin{thm}\label{mainthm}

Let $\xi_{vo}$ be a virtually overtwisted structure on $L(p,q)$. If we are in one of the cases below, then $(L(p,q),\xi_{vo})$ is not the boundary of the Milnor fiber of any complex hypersurface singularity:

\begin{itemize}

\item[a)] $p/q=[a_1,a_2,\ldots,a_n]$ and $a_i$ is odd for some $i$;

\item[b)] $p/q=[2x_1,2x_2]$;

\item[c)] $p/q=[2x_1,2x_2,\ldots,2x_n]$, with $x_i>1$ for every $i$ ($n\geq 3$) and either:
\begin{itemize}
\item[i)] $q^2\not\equiv 1 \mod p$;
\item[ii)] $q^2\equiv 1 \mod p$ and $n$ is even.
\end{itemize}

\end{itemize}
\end{thm}

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.35]{chain.pdf}
\caption{Surgery link producing $L(p,q)$.}
\label{chain}
\end{figure}

The first step in proving this theorem is to characterize those contact structures $\xi$ which can appear in the context of hypersurface singularity: the fact that $c_1(\xi)$ vanishes imposes certain conditions on the coefficients in the continued fraction expansion of $p/q$, which allow us to prove part $(a)$ of theorem \ref{mainthm}. Section \ref{secvan} deals with these numerical restrictions, adapted to the language of contact geometry. To prove parts $(b)$ and $(c)$ we need to look closely at the topology of the Milnor fibration and analyze its monodromy. In order to derive our statements we study the integral orthogonal group of the intersection form of the Milnor fiber, imposing a further restriction coming from a theorem of \cite{acampo}. This is explained in section \ref{mainproof}. The article ends with an open question which focuses on the limits of theorem \ref{mainthm}: where are these techniques failing?

\paragraph{Acknowledgments} The author wishes to thank András Némethi for raising this problem and for the help in trying to tackle it. Thanks are due also to András Stipsicz for the hospitality at the Renyi Institute of Mathematics in Budapest, where this work originated. Lastly, a lot of support and key insights came from Paolo Lisca, who pointed out the results from \cite{gerstein} and \cite{c1planar}. This article will be part of the PhD thesis of the author.

\section{Vanishing of the rotation numbers}\label{secvan}

The goal of this section is to prove theorem \ref{thmc1=0}, which is a special case of \cite[corollary 1.5]{c1planar} (this result applies since, by \cite{schonenberger}, every contact structure on $L(p,q)$ is planar). We prove it using elementary techniques that do not involve the Heegaard-Floer contact invariant. Theorem \ref{thmc1=0} will be the starting point in the proof of theorem \ref{mainthm}.

\begin{thm}\label{thmc1=0}
Let $L$ be a Legendrian linear chain of unknots in the standard contact $S^3$ and let $(L(p,q),\xi)$ be the contact 3-manifold obtained by Legendrian surgery on $L$. If $c_1(\xi)=0$, then $\rot (L_i)=0$ for every component $L_i$ of the link $L$.
\end{thm}

To specify a tight contact structure on $L(p,q)$ we put the link of figure \ref{chain} into Legendrian position with respect to the standard tight contact structure of $S^3$, in order to form a linear chain of Legendrian unknots, see for example figure \ref{eg}. We do this in such a way that the Thurston-Bennequin number of the $i^{\mbox{th}}$ component is $a_i-1$. When such a Legendrian representative has been chosen, we get a tight structure $\xi$ on $L(p,q)$ by performing $(-1)$-contact surgery on it.

\begin{figure}[h!]
\centering
\begin{subfigure}[t]{.33\textwidth}
  \centering
  \includegraphics[scale=0.3]{eg1}
  \caption{$L(28,15)$}
\end{subfigure}%
\begin{subfigure}[t]{.33\textwidth}
  \centering
 \includegraphics[scale=0.3]{eg2}
 \caption{$L(34,7)$}
\end{subfigure}
\begin{subfigure}[t]{.33\textwidth}
  \centering
  \includegraphics[scale=0.3]{eg3}
  \caption{$L(17,7)$}
\end{subfigure}%
\caption{Examples of Legendrian links.}
\label{eg}
\end{figure}

\noindent The rotation number of an \emph{oriented} Legendrian knot in $(S^3,\xi_{st})$ can be computed in the front projection by the formula
\[\rot(K)=\frac{1}{2}(c_D-c_U),\]
where $c_D$ and $c_U$ are respectively the number of down and up cusps. Since, as theorem \ref{thmc1=0} suggests, we are interested only in those Legendrian realizations with rotation number zero, we can avoid specifying an orientation.

To prove theorem \ref{thmc1=0} we need some notation that makes the computation easier, and a few results more.
Let $Q$ be the linking matrix of the link of figure \ref{chain}, which represents the intersection form of the plumbed 4-manifold associated to the linear graph
\begin{center}
\begin{tikzpicture}
        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-a_1$] (1)                  {};
        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-a_2$] (2) [right=of 1] {}
        edge [-]               (1);
        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-a_3$] (3) [right=of 2] {}
        edge [-]               (2);
        \node[shape=circle,fill=black,inner sep=1.5pt] (4) [right=of 3] {} ;
                 \node at ($(3)!.5!(4)$) {\ldots};

        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-a_n$] (5) [right=of 4] {}
        edge [-]               (4);
\end{tikzpicture}
\end{center}
written in the ordered basis given by the meridian of each curve. We know that
\[H^2(L(p,q))\ni c_1(\xi)=0 \Leftrightarrow \PD(c_1(\xi))=0\in H_1(L(p,q)).\]
By looking at the linear plumbing graph we can express
\[\PD(c_1(\xi))=\sum_{i=1}^nr_i\mu_i,\]
where $r_i$ and $\mu_i$ are respectively the rotation number $\rot(L_i)$ and the meridian of the $i^{th}$ component (compare with \cite{ozbagci}). The $\mu_i$'s are related by a set of equations coming from $Q$:
\[
\begin{cases}
-a_1\mu_1+\mu_2=0	\\
\mu_1-a_2\mu_2+\mu_3=0 \\
\vdots \\
\mu_j-a_{j+1}\mu_{j+1}+\mu_{j+2}=0 \\
\vdots \\
\mu_{n-1}-a_n\mu_n=0.
\end{cases}
\]
By choosing $\mu_1\in H_1(L(p,q))$ as a generator we can consider it to be $1\in \Z/p\Z\simeq H_1(L(p,q))$ and therefore identify all the other $\mu_i$'s with numbers (which are well-defined modulo $p$). We derive relations which hold over $\Z$, but this will be enough to us thanks to lemma \ref{lemma0}.
Therefore we get recursive expressions of the form
\[
\begin{cases}
\mu_1=1\\
\mu_2=a_1\\
\mu_i=a_{i-1}\mu_{i-1}-\mu_{i-2}.
\end{cases}
\]
Then $c_1(\xi)$ is 0 exactly when
\begin{equation}\label{chernrotations}
\sum_{i=1}^n r_i\mu_i\equiv 0\mbox{ modulo } \det(Q).
\end{equation}
Define
\[
\begin{cases}
\Delta[-1]=0;\\
\Delta[0]=1;\\
\Delta[i]=-a_i\Delta[i-1]-\Delta[i-2]
\end{cases}
\]
and note that
\begin{itemize}
\item $\det(Q)=\Delta[n]$;
\item $\sign(\Delta[i])=(-1)^i$, hence $\Delta[i]=(-1)^i|\Delta[i]|$;
\item $\mu_i> \mu_{i-1}$;
\item $|r_i|\leq a_i-2$ and $r_i\equiv a_i \mod 2$;
\item $\Delta[i]=(-1)^i\mu_{i+1}$ (proved by induction), therefore $|\Delta[i]|>|\Delta[i-1]|$.
\end{itemize}

\begin{lem}\label{lemma0}
Equation \ref{chernrotations} is satisfied in $\Z/\det(Q)\Z$ if and only if it is satisfied in $\Z$.
\end{lem}

\begin{prf}
We prove by induction on $n$ that
\begin{equation}\label{inequality}
|\det(Q)|>\abs[\Big]{\sum_{i=1}^n r_i\mu_i}.
\end{equation}
This will tell that equation \ref{chernrotations} can only be satisfied with the zero-multiple of $\det(Q)$, since clearly for any non-zero $m\in\Z$ we have
\[|m\cdot \det(Q)|>|\det(Q)|.\]
If $n=1$ then the two sides of the inequality \ref{inequality} are respectively $|-a_1|$ and $|r_1|$, so it is true.
The first interesting case is then $n=2$:
\[a_2a_1-1>r_1+r_2a_1?\]
We have:
\begin{align*}
r_1+r_2a_2 < &(a_1-2)+(a_2-2)a_1\\
= & a_1a_2-a_1-2\\
< & a_1a_2-1.\qquad \checkmark
\end{align*}
The general case now:
\begin{align*}
\abs[\Big]{\sum_{i=1}^nr_i\mu_i}= & \abs[\Big]{\sum_{i=1}^nr_i\Delta[i-1]}\\
<&\left(\abs[\Big]{\sum_{i=1}^{n-1}r_i\Delta[i-1]}\right)+|r_n\Delta[n-1]| && \mbox{(induction)} \\
<&|\Delta[n-1]|+|r_n\Delta[n-1]|\\
\leq &|\Delta[n-1]|+|(a_n-2)\Delta[n-1]|.
\end{align*}
There are two possibilities for the right-hand side, according to the parity of $n$.
\begin{itemize}
\item[1)] $n$ is even (hence $\Delta[n]>0$):
\begin{align*}
|\Delta[n-1]+|(a_n-2)\Delta[n-1] |=& -\Delta[n-1]-a_n\Delta[n-1]+2\Delta[n-1]\\
=&-a_n\Delta[n-1]+\Delta[n-1]\\
<& -a_n\Delta[n-1]-\Delta[n-2] \mbox{ (because $|\Delta[i]|>|\Delta[i-1]|$ and $n$ is even)}\\
=&\Delta[n]\\
=&|\det(Q)|.
\end{align*}
\item[2)] $n$ is odd (hence $\Delta[n]<0$):
\begin{align*}
|\Delta[n-1]+|(a_n-2)\Delta[n-1] |=& a_n\Delta[n-1]-\Delta[n-1]\\
<& a_n\Delta[n-1]+\Delta[n-2] \mbox{ (because $|\Delta[i]|>|\Delta[i-1]|$ and $n$ is odd)}\\
=&-\Delta[n]\\
=&|\det(Q)|.
\end{align*}
\end{itemize}
\end{prf}

\begin{lem}\label{lemma1}
\[a_n\mu_n-a_{n-1}\mu_{n-1}>0.\]
\end{lem}

\begin{prf} We prove it by induction on $n$. If $n=2$ then the formula is just
\[a_2\mu_2-a_1\mu_1=a_2a_1-a_1>0.\qquad\checkmark\]
In general, assuming the result true for $n-1$, we have
\begin{align*}
a_n\mu_n-a_{n-1}\mu_{n-1}=& a_n(a_{n-1}\mu_{n-1}-\mu_{n-2})-a_{n-1}\mu_{n-1}\\
=& a_{n-1}\mu_{n-1}(a_n-1)-a_n\mu_{n-2}&& \mbox{(induction)}\\
>& a_{n-2}\mu_{n-2}(a_n-1)-a_n\mu_{n-2}\\
=& (a_{n-2}a_n-a_{n-2}-a_n)\mu_{n-2}\\
>& 0.
\end{align*}
\end{prf}

\begin{lem}\label{lemma2} If $r_n\neq 0$, then we have
\[|r_n\mu_n|-|r_{n-1}\mu_{n-1}|>0.\]
\end{lem}

\begin{prf}
\begin{align*}
|r_n\mu_n|-|r_{n-1}\mu_{n-1}|>& |\mu_n|-(a_{n-1}-2)\mu_{n-1}\\
=& a_{n-1}\mu_{n-1}-\mu_{n-2}-a_{n-1}\mu_{n-1}+2\mu_{n-1}\\
=& 2\mu_{n-1}-\mu_{n-2}\\
>& 0.
\end{align*}
\end{prf}

\begin{lem}\label{lemma3}
If $r_n\neq 0$, then we have
\[|r_n\mu_n|-\abs[\Big]{\sum_{i=1}^{n-1}r_i\mu_i}>0.\]
\end{lem}

\begin{prf} We prove it by induction on $n$. If $n=2$ the formula is just
\[|r_2\mu_2|-|r_1\mu_1|=|r_2a_1|-|r_1|>a_1-(a_1-2)=2>0.\qquad\checkmark\]
Now we do the general case. Assume the inequality holds for $n-1$ and let $j\leq n-2$ be the biggest integer such that $r_j\neq 0$ (note that if $r_i=0,\;\forall i\leq n-2$, then by lemma \ref{lemma2} we would be done after the second line in the following computation).
We have
\begin{align*}
|r_n\mu_n|-\abs[\Big]{\sum_{i=1}^{n-1}r_i\mu_i} >& |\mu_n|-\abs[\Big]{\sum_{i=1}^{j-1}r_i\mu_i}-|r_{n-1}\mu_{n-1}| &&  \mbox{(induction)}\\
>& \mu_n-|r_j\mu_j|-|r_{n-1}\mu_{n-1}| \\
\geq & \mu_n-(a_j-2)\mu_j-(a_{n-1}-2)\mu_{n-1}\\
= & a_{n-1}\mu_{n-1}-\mu_{n-2}-(a_j-2)\mu_j-(a_{n-1}-2)\mu_{n-1} \\
= & 2\mu_{n-1}-\mu_{n-2}-(a_j-2)\mu_j	\\
> & \mu_{n-1}-(a_j-2)\mu_j\\
= & \mu_{n-1}+\mu_j+\mu_j-a_j\mu_j && (n-2\geq j\Rightarrow n-1\geq j+1)\\
>&\mu_{j+1}+\mu_{j-1}+\mu_{j-1}-a_j\mu_j && (\mu_{j+1}+\mu_{j-1}=a_{j}\mu_{j})\\
=&\mu_{j-1}\\
>&0.
\end{align*}
\end{prf}

\begin{lem} \label{lemma4}
\[\sum_{i=1}^{n}r_i\mu_i=0\Longrightarrow r_i=0,\;\forall i.\]
\end{lem}

\begin{prf}
\[\sum_{i=1}^{n}r_i\mu_i=0\Longrightarrow r_n\mu_n=-\sum_{i=1}^{n-1}r_i\mu_i\Longrightarrow |r_n\mu_n|=\abs[\Big]{\sum_{i=1}^{n-1}r_i\mu_i}.
\]
But if $r_n\neq 0$, then we should have a strict inequality by lemma \ref{lemma3}, hence $r_n=0$. By applying this repeatedly we get to
\[r_n=r_{n-1}=\ldots=r_1=0.\]
\end{prf}

\noindent We can finally give the following:

\begin{prf}[of theorem \ref{thmc1=0}]
By combining lemmas \ref{lemma0} and \ref{lemma4}, we have that
\[c_1(\xi)=0\Longleftrightarrow \sum_{i=1}^n r_i\mu_i\equiv 0\mod p \Longleftrightarrow \sum_{i=1}^n r_i\mu_i=0\Longleftrightarrow r_i=0\,\forall i,\]
where $r_i=\rot (L_i)$.
\end{prf}

\section{Proof of main theorem }\label{mainproof}

\begin{center}
Are the virtually overtwisted structures on lens spaces realizable as the boundary of the Milnor fiber of some complex hypersurface singularity?
\end{center}

\noindent In the universally tight case we know (see \cite{nemethippp} and \cite{bhupal}) that all the fillings come from algebraic geometry, and the article \cite{park} shows that complex surface singularities produce all the fillings of small Seifert 3-manifold equipped with the  canonical contact structure.

\vspace{0.5cm}

\noindent \textbf{Key fact:} suppose there is a polynomial function $f:(\C^3,0)\to (\C,0)$ such that
\[\partial(F,J)=(L(p,q),\xi),\]
where $(F,J)$ is the Milnor fiber of $f$ and $\xi$ is the contact structure on the boundary induced by complex tangencies, as explained in section \ref{intro}. The links of figure \ref{eg} represent more than contact 3-manifolds: the components $L_i$ of any of those links can be thought as the attaching circles of the 2-handles of the Stein 4-manifold $(F,J)$. The first Chern class $c_1(F,J)\in H^2(F;\Z)$ evaluates on each 2-handle as the correspondent rotation number, and $c_1(\xi)$ is the restriction of $c_1(F,J)$. But the tangent bundle of $F$ is stably trivial, hence $c_1(F,J)=0$, and also $c_1(\xi)=0$ on the boundary $L(p,q)$.

\begin{prf}[of theorem \ref{mainthm}a]
If $(L(p,q),\xi_{vo})$ is the boundary of the Milnor fiber of a complex hypersurface singularity, then $c_1(\xi_{vo})=0$ and, by theorem \ref{thmc1=0}, all $\rot(L_i)$ are zero. \\
\noindent Let $p/q=[a_1,a_2,\ldots,a_n]$ and remember that
\[a_i\equiv \rot(L_i)\mod 2.\]
Since all the rotation numbers are zero, the conclusion follows.
\end{prf}

Theorem \ref{thmc1=0} implies also the following corollary, which will be used later.

\begin{cor}\label{uniqueness}
If $(L(p,q),\xi)$ is the boundary of the Milnor fiber of a hypersurface singularity, then it has a unique Stein filling, which is the Milnor fiber itself.
\end{cor}

\begin{prf}
The fact that $\rot(L_i)=0$ for every $i$ implies, by \cite[theorem 1.3]{menke}, that from the chain of Legendrian unknots producing $(L(p,q),\xi)$ we can forget about those components with $\tb(L_i)\neq -1$ and look for Stein fillings of the 3-manifold $Y$ which is left. Then, all the Stein fillings of $(L(p,q),\xi)$ will be uniquely obtained by attaching the 2-handles (corresponding to the forgotten components) to the Stein fillings of $Y$. From the link diagram of $(L(p,q),\xi)$ we see that $Y$ is a connected sum of $(L(n_j,n_j-1),\xi_{st})$, where each of the prime factor corresponds to a string of $-2$ in the expansion of $p/q$.
By \cite[theorem 16.9]{steinandback}, $Y$ admits a unique Stein filling, because each factor $(L(n_j,n_j-1),\xi_{st})$ does. This concludes the proof.
\end{prf}

\vspace{0.5cm}

\noindent Another consequence of theorem \ref{thmc1=0} is:

\begin{cor}\label{L(p,p-1)}
Let $(L(p,q),\xi)$ be a lens space with a contact structure arising as the link of an isolated hypersurface singularity. Then $q=p-1$.
\end{cor}

\begin{prf}
From theorem \ref{thmc1=0} we have that all the rotation numbers are zero and \cite[theorem 2.1]{lekili} says that $\xi$ is universally tight. By the work of \cite{honda} we know that a universally tight structure on a lens space is the result of contact $(-1)$-surgery on a link where all the stabilizations appear on the same side, i.e. when the (absolute values of the) rotation numbers are maximal. Therefore every Legendrian knot must have Thurston-Bennequin number equal to $-1$:
\[-\frac{p}{q}=[-2,-2,\ldots,-2]\qquad\Rightarrow\qquad q=p-1.\]
\end{prf}

\subsection*{Proof of Theorem 1b}

Corollary \ref{uniqueness} says that if $(L(p,q),\xi_{vo})$ arises as $\partial (F,J)$, then the Stein filling $F$ is uniquely determined: topologically it is given by the plumbing of spheres according to the expansion of $p/q$. The monodromy $\varphi$ of the Milnor fibration induces, in cohomology, a homomorphism
\[\varphi^*:H^*(F;\Z)\to H^*(F;\Z)\]
such that the alternating sum of the traces is zero, by \cite[theorem 1]{acampo}:
\[\trace(\varphi^0)-\trace(\varphi^1)+\trace(\varphi^2)-\trace(\varphi^3)+\trace(\varphi^4)=0.\]
In our case, the Stein fillings of those lens spaces with $c_1(\xi)=0$ are simply connected, hence $\varphi^1=\varphi^3=0$. Moreover, $\varphi^0:\Z\to\Z$ is the identity and $\varphi^4=0$. Therefore previous equation simply reads as:
\begin{equation}\label{trace}
1+\trace(\varphi^2)=0.
\end{equation}
The intersection form of $F$ must be preserved by the homomorphism $\varphi^2$ and we are therefore led to study its isometry group. If in this group there is no element whose trace is $-1$, then formula \ref{trace} cannot be satisfied, and we conclude that the polynomial function $f:(\C^3,0)\to (\C,0)$ with $(\partial F,\xi)=(L(p,q),\xi_{vo})$ does not exist.

\begin{prf}[of theorem \ref{mainthm}b]
Our goal is to prove that there is no $f:(\C^3,0)\to (\C,0)$, with non-isolated singularity, whose Milnor fiber has boundary $L(p,q)$. Note that
\[x_1x_2>1,\]
otherwise the induced contact structure is universally tight.

Assume by contradiction that such $f$ exists. Then, by corollary \ref{uniqueness}, we know that the Milnor fiber $F$ has negative-definite intersection form
\[-M=
\begin{bmatrix}
   -2x_1 & 1 \\
   1 & -2x_2 \\
\end{bmatrix}.
\]
By equality \ref{trace}, we must have $\trace(\varphi^2)=-1$. The morphism $\varphi^2:H^2(F;\Z)\to H^2(F;\Z)$ is induced by a diffeomorphism which preserves the intersection form and therefore it is represented by an integral matrix $A$ with
\[\begin{cases}
|\det(A)|=1 \\
\trace(A)=-1
\end{cases}
\]
and such that
\[A(-M)A^T=-M.\]
We show now that such matrix cannot exist. We change sign to work with a positive definite matrix:
\[AMA^T=M\Longrightarrow
 \begin{bmatrix}
   a_1 & a_2 \\
   a_3 & a_4 \\
\end{bmatrix}
\begin{bmatrix}
   2x_1 & -1 \\
   -1 & 2x_2 \\
\end{bmatrix}
\begin{bmatrix}
   a_1 & a_3 \\
   a_2 & a_4 \\
\end{bmatrix} =
\begin{bmatrix}
   2x_1 & -1 \\
   -1 & 2x_2 \\
\end{bmatrix}. \]
We get equations:
\[\begin{cases}
2x_1a_1^2-2a_1a_2+2x_2a_2^2=2x_1 \\
2x_1a_3^2-2a_3a_4+2x_2a_4^2=2x_2
\end{cases}
\]
that can be rewritten as

\begin{subequations}
\begin{empheq}[left=\empheqlbrace]{align}
  & (2x_1-1)a_1^2+(a_1-a_2)^2+(2x_2-1)a_2^2=2x_1 \label{star1}
  \\
  & (2x_1-1)a_3^2+(a_3-a_4)^2+(2x_2-1)a_4^2=2x_2. \label{star2}
\end{empheq}
\end{subequations}

\noindent From equations \ref{star1} and \ref{star2} it follows that $a_1^2\leq 1$ and $a_4^2\leq 1$. But since $\trace(A)=a_1+a_4=-1$, we have that either
\[ \begin{cases}
a_1=0 \\
a_4=-1
\end{cases}
\qquad \mbox{or}\qquad
\begin{cases}
a_1=-1 \\
a_4=0.
\end{cases}\]
We do the first case $(a_1=0,\,a_4=-1)$, the other one is the same. From $AMA^T=M$ we also get
\[
 \begin{bmatrix}
   0 & a_2
\end{bmatrix}
\begin{bmatrix}
   2x_1 & -1 \\
   -1 & 2x_2 \\
\end{bmatrix}
\begin{bmatrix}
   a_3 \\
   -1  \\
\end{bmatrix} = -1, \]
which gives the equation $a_2(a_3+2x_2)=1$. Hence $a_2=a_3+2x_2=1$ (or both $-1$, but the conclusion is the same).
From \ref{star1} we have $x_2a_2^2=x_1$, which gives $x_2=x_1$. Then
\[\pm 1=\det(A)=\det
\begin{bmatrix}
   0 & \pm 1 \\
   a_3 & -1  \\
\end{bmatrix}
\Longrightarrow a_3=\pm 1.\]
We are in the case where $a_3+2x_2=1$, so either $x_1=x_2=0$ or $x_1=x_2=1$, which are both contradicting the condition $x_1x_2>1$.
\end{prf}

\subsection*{Proof of Theorem 1c}

Remember that in order to have $c_1(\xi)=0$ we need all the rotation numbers to be zero and in particular all the coefficients in the expansion to be even. Let $-M$ be the negative definite intersection lattice associated to the linear plumbing of spheres
\begin{center}
\begin{tikzpicture}
        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-2x_1$] (1)                  {};
        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-2x_2$] (2) [right=of 1] {}
        edge [-]               (1);
       \node[shape=circle,fill=black,inner sep=1.5pt,label=$-2x_3$] (3) [right=of 2] {}
        edge [-]               (2);
        \node[shape=circle,fill=black,inner sep=1.5pt] (4) [right=of 3] {} ;
                         \node at ($(3)!.5!(4)$) {\ldots};

        \node[shape=circle,fill=black,inner sep=1.5pt,label=$-2x_n$] (5) [right=of 4] {}
        edge [-]               (4);
\end{tikzpicture}
\end{center}

We are looking for a matrix $A$ representing the monodromy of a Milnor fibration on the second cohomology group, that respect the intersection form of the Milnor fiber (i.e. $A(-M)A^T=M$) and whose trace is $-1$.

What we need to understand is the integral orthogonal group $O_{\Z}(-M)$ of the negative definite lattice $(\Z^n,-M)$, which is isomorphic to $O_{\Z}(M)$. The latter is studied in the article \cite{gerstein}, where the following theorem is proved:

\begin{thm*}[\cite{gerstein}]
Let $M$ be the integer matrix
\[
\begin{bmatrix}
    2x_1       & -1  \\
 -1       & \ddots & \ddots \\
 & \ddots & \ddots & \ddots \\
  & & \ddots &\ddots &-1\\
  & & & -1 & 2x_n\\
\end{bmatrix}
\]
with $n\geq 2$ and $2x_i\geq 3\;\forall i$.
\begin{itemize}
\item[i)] If $x_i\neq x_{n+1-i}$ for some $i$, then $O_{\Z}(M)=\{\pm \id\}$.
\item[ii)] If $x_i=x_{n+1-i}$ for every $i=1,\ldots, n$, then $O_{\Z}(M)=\{\pm \id,\pm \rho\}$, where $\rho$ is the isometry that inverts the order of a basis.
\end{itemize}
\end{thm*}

Therefore we can derive:

\begin{prf}[of theorem \ref{mainthm}c]
The condition $q^2\not\equiv 1 \mod p$ can be rephrased in terms of the coefficients of the expansion by saying that $x_i\equiv x_{n+1-i}$ for some $i$ (see \cite[appendix]{orlik}). On the other hand, $x_i=x_{n+1-i}$ for all $i$ if and only if $q^2= 1 \mod p$.

In the first case, previous theorem tells that if there is an integer matrix $A$ with $AMA^T=M$, then $A=\pm	\id$. Since this does not have trace $-1$, there cannot be a hypersurface singularity whose Milnor fiber has boundary $(L(p,q),\xi_{vo})$, otherwise we would have a contradiction with A'Campo's formula \ref{trace}.

In the second case, again A'Campo's formula cannot be satisfied because a hypersurface singularity would produce a Milnor fibration with monodromy $A\in\{\pm \id,\pm \rho\}$:
\[\rho=
\begin{bmatrix}
     & & & 1  \\
  &   & 1  \\
  & \reflectbox{$\ddots$} \\
   1\\
\end{bmatrix}
\]
and, if $n$ is even, then $\trace(\rho)=0\neq -1$.
\end{prf}

\begin{que} There are cases which are not covered by theorem \ref{mainthm}: every time that in the continued fraction expansion of $-p/q$ appears a $-2$, the orthogonal group $O_{\Z}(M)$ is harder to understand. Nevertheless, in the easier case when $q^2\equiv 1\mod p$ we have a complete description of $O_{\Z}(M)$, but inside this group there is a matrix with trace -1 if the length of the expansion is odd. A simple case is for example
\[-\frac{p}{q}=-\frac{12}{7}=[-2,-4,-2].\]
Our techniques indeed do not exclude that the isometry
\[-\rho=
\begin{bmatrix}
    0 & 0 &-1  \\
  0 &-1  & 0  \\
   -1 &0  & 0  \\
\end{bmatrix}
\]
is the morphism induced by the monodromy of the Milnor fibration of a certain non-isolated hypersurface singularity producing $L(12,7)$. How can we deal with cases like this and exclude the existence of such a monodromy? Further works will hopefully clarify this problem and either find the polynomial function or rule out this possibility as well.
\end{que}

\bibliographystyle{alpha}


\section{Introduction \& Related Work}\label{sec:intro}

Motion planning is a core problem for robotic systems that must successfully navigate their surroundings. In uncertain and dynamic environments, the planning problem becomes especially challenging since a previously feasible solution can quickly become infeasible and must be recomputed.
Broadly speaking, motion planning approaches can be grouped into the categories of search, sampling, and optimization.

Search-based algorithms, such as A*~\cite{astar} guarantee completeness and optimality, but quickly succumb to the curse of dimensionality as discretization of the environment is necessary.
Sampling-based algorithms, including PRM~\cite{kavraki1996probabilistic} and RRT~\cite{kuffner2000rrt}, are probabilistically complete, generalize better to higher dimensional problems, and have been extended to provide (asymptotic) optimality~\cite{karaman2010incremental}. While these methods are well suited to solve complex problems in environments like 2D mazes, in practice they tend to be very computationally expensive for high dimensional systems and require additional post-processing to smooth the solutions.
Conversely, trajectory optimization algorithms initialize a trajectory from start to goal state and then optimize the path by minimizing some cost function~\cite{zucker2013chomp,schulman2014motion}. GPMP2~\cite{Mukadam-IJRR-18} represents smooth, continuous trajectories as samples from a Gaussian Process (GP) and solves the motion planning problem by performing probabilistic inference on a factor graph.  While these methods are much faster for high dimensional systems, they are sensitive to initialization and can get stuck in poor local optima. In practice, straight-line initializations have been shown to work well and random restarts are performed if that fails~\cite{zucker2013chomp,schulman2014motion}. Initializing with a feasible solution from a sampling-based planner is also an option~\cite{zucker2013chomp}, but this incurs extra computational burden.

\begin{figure}[!t]
    \centering
    \begin{subfigure}[t]{0.48\linewidth}
		\centering
        \includegraphics[trim={95 95 95 95},clip,width=1\linewidth]{switch_1}
        \caption{}
		\label{fig:switch1}
    \end{subfigure}
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.48\linewidth}
		\centering
        \includegraphics[trim={95 95 95 95},clip,width=1\linewidth]{switch_2}
        \caption{}
		\label{fig:switch2}
    \end{subfigure}
    \caption{\small \text{POSH}\xspace being used in a dynamic environment where the robot (gray circle) is tasked to reach the goal (green circle) while avoiding obstacles (gray squares). At some time step $i$ In (a) the robot has traveled along the green path when \text{POSH}\xspace returns a new best (lowest-cost and feasible) trajectory (red). This represents a homotopy switch from the previous best trajectory (black). It executes the first step of this new trajectory and then at the next time step $i+1$ in (b) \text{POSH}\xspace again returns a new best trajectory (red) that is in a different homotopy class.}
    \label{fig:dynamic_forest}
	\vspace{-3mm}
\end{figure}

Recently, hybrid approaches to planning have attempted to combine the complimentary benefits of different planning schemes. For example, BIT*~\cite{gammell2015batch} combines graph-based and sample-based planning by using a heuristic function to guide a series of random geometric graphs towards the goal state. Its successor, RABIT*~\cite{choudhury2016rabit} performs BIT* while running trajectory optimization when constructing edges between samples, thus reducing sampling complexity by decreasing the number of samples discarded.
GPMP-GRAPH~\cite{Huang-ICRA-17}, an extension of GPMP2~\cite{Dong-RSS-16}, constructs a graph of interconnected trajectory initializations, similar to the one shown in Fig.~\ref{fig:dynamic_forest}, on a factor graph that when optimized can simultaneously evaluate an exponential number of initializations. After optimization, any path through this graph provides a feasible trajectory. Thus, it has the ability to find solutions in multiple unique homotopy-classes all at once.

While these planning approaches are typically employed in static environments, a number of algorithms~\cite{Park-ITOMP} have been proposed to handle navigating dynamic environments. Lifelong Planning A*~\cite{Koenig-LPA}, D* lite~\cite{SunYK10} and comparable online planning algorithms outperform repeatedly replanning from scratch by incrementally finding shortest paths on a graph with changing edge costs. RRT-X~\cite{Otte2014RRTXRM} an extension of RRT is suited to environments with moving obstacles and provides comparable runtime performance. RAMP~\cite{Vannoy-RAMP} draws from evolutionary computation by maintaining a population of trajectories and evaluating their respective quality (fitness) with respect to a cost function and has been shown to discover and exploit new homotopies. Since solutions can go in and out of feasibility as the environment changes access to solutions in multiple homotopy classes is beneficial. However, there exists a gap in current research with respect to solving the online planning problem by leveraging trajectories across different homotopy classes.

In this work, we present a novel trajectory optimization algorithm, Planning Online by Switching Homotopies (\text{POSH}\xspace), to handle such scenarios. We build on GPMP-GRAPH~\cite{Huang-ICRA-17} where multiple trajectories are inter-connected and represented as a factor graph upon which probabilistic inference is performed to optimize the entire graph. However, instead of retaining the same optimized trajectory from the initial time step for execution as GPMP-GRAPH does, POSH maintains and updates the entire graph at every time step. Specifically, at any time step the graph is pruned to remove unreachable (in time) states and is then reoptimized considering changes in the environment to find the new optimal trajectory. This grants our algorithm the unique ability to dynamically switch between different homotopy classes as illustrated in Fig.~\ref{fig:dynamic_forest} and allows it to better contend with dynamic environments as demonstrated by our experiments.

\section{Background}\label{sec:back}

In this section we review GPMP-GRAPH, closely following the explanation of Huang et al.~\cite{Huang-ICRA-17}, as well as other concepts, such as homotopy classes and their application to motion planning, necessary to explain our approach.

\subsection{GPMP-GRAPH}
GPMP-GRAPH is a graph-based trajectory optimization algorithm  designed to address the local minima problem by evaluating an exponential number of initializations simultaneously and, as a consequence, finding solutions in multiple homotopy classes. The graph is constructed as a series of chains between the start and goal state, with connections between the chains. Optimization happens over the entire graph, yielding improvements in planning time over sequentially optimizing a single chain with different initializations~\cite{Huang-ICRA-17}.

GPMP-GRAPH treats the motion planning problem as probabilistic inference. Given desired events $e$, it finds the maximum a posteriori (MAP) trajectory $x^*$
$$x^* = \arg\max_x p(x|e) = \arg\max_x p(x)l(x;e)$$
where $p(x)$ is the prior distribution over trajectories and $l(x;e)$ is the likelihood.

Continuous-time trajectories are represented as samples from a vector-valued Gaussian process (GP) $x(t)\sim \mathcal{GP}(\mu(t), \mathcal{K}(t, t'))$ where $\mu(t)$ is a vector-valued mean function and $\mathcal{K}(t, t')$ is a matrix-valued covariance function.
Following the definition of GP, trajectories can be parameterized at discrete times $t = \{t_0, \dots, t_n\}$ to be jointly Gaussian distributed
$$x = [x_0, \dots, x_n]^T\sim\mathcal{N}(\mu, \mathcal{K})$$
where:
$$\mu = [\mu(t_0), \dots, \mu(t_n)], \qquad \mathcal{K} = [\mathcal{K}(t_i, t_j)]_{i, j,  0 \leq i,  j \leq n}$$
The values $x_0, \dots, x_n$ are referred to as \textit{support states} and correspond to $x(0),\dots, x(n)$. This defines a prior distribution over trajectories
$$p(x) \propto \exp\Bigg\{-\frac12||x - \mu||^2_\mathcal{K}\Bigg\}$$
In practice, this distribution can be constructed from robot dynamics, for example, modeled as a linear time-varying stochastic differential equation (LTV-SDE)
$$\dot{x}(t) = A(t)x(t) + u(t) + F(t)w(t)$$
where $A$, $F$ are system matrices, $u$ is a known control input, and
$w(t) \propto \mathcal{GP}(0, Q_c\delta(t - t'))$ is a white noise process with $Q_c$ being the power spectral density matrix of the system and $\delta$ being Dirac delta function.
The first- and second-order moments of the solution to this SDE yield the mean and covariance of the GP, respectively. The prior penalizes the trajectory from deviations from the mean to maintain smoothness in the solutions. The likelihood $l(x;e)$, handles costs like obstacle avoidance and any other motion constraints.

\begin{figure}
    \centering
    \includegraphics[width=0.48\textwidth]{figs/fgfull.pdf}
    \caption{\small Factor Graph consisting of three Gauss-Markov Chains as illustrated in~\cite{Huang-ICRA-17}. The white circles represent support states. Black squares are collision factors and black circles are GP priors factors. Interpolated collision factors are omitted for clarity.}
    \label{fig:factor_graph}
	\vspace{-3mm}
\end{figure}

The posterior distribution $p(x|e)$ can be represented as a \textit{factor graph}. A factor graph is a bipartite graph $G = \{\Theta, F, E\}$, where $\Theta$ is a set of variables, $F$ is a set of factors, and $E$ is the set of edges connecting $F$ to $\Theta$. Then the posterior can be factorized as
$$p(x|e) = \prod_{t_i}f^{gp}_i(x_i, x_{i + 1})f^{obs}_i(x_i)\prod_{\tau}f^{intp}_{i, \tau} (x_i, x_{i + 1})$$
where $f^{gp}_i$ is a binary factor connecting consecutive states, $f_i^{obs}$ is a unary collision likelihood factor, and $f^{intp}_{i, \tau}$ is a collision likelihood factor for a state at $\tau$ between consecutive support states.
This factor graph has a sparse structure which can be exploited for rapidly solving the optimization problem~\cite{Mukadam-IJRR-18}.

So far this is identical to GPMP2~\cite{Mukadam-IJRR-18}, however GPMP-GRAPH constructs such a factor graph with multiple initial trajectories (chains), each with its own set of support states. Adjacent chains have interconnected support states with a GP prior factor and interpolated factors between two states from two different chains, allowing the robot to switch between chains. An example graph construction with three chains and interconnections is shown in Fig.~\ref{fig:factor_graph}. The graph's variables (corresponding to trajectory waypoints) are optimized using the Levenberg-Marquardt algorithm to maximize their probability given the desired events, such as being collision free and conforming to motion constraints. Then, from the optimized graph we can find the trajectory with the lowest cost or the trajectory that conforms to any predefined homotopy constraints set by the user. For more details please refer to~\cite{Huang-ICRA-17}.

GPMP-GRAPH suffers in practical performance because it is unable to exploit the different homotopy classes it identifies when the environment changes. By constructing a factor graph with multiple chains and then allowing those chains to converge to low-cost solutions in different homotopy-classes, the algorithm has many feasible trajectories available to it. However, dynamic obstacles may cause an initially good trajectory to decrease in quality if the obstacle moves in the way of the trajectory. Only optimizing the chain that the robot is on at that time may not allow the robot to switch homotopies, leading to collisions with the moving obstacle or inefficient solutions.

\subsection{Homotopy Classes}

GPMP-GRAPH improves upon its predecessor GPMP2, by initially considering multiple homotopy classes when optimizing over different initial trajectories in parallel.
Two continuous functions over the real numbers $x_1, x_2:~\mathbb{R}~\rightarrow~\mathbb{R}$ are said to be homotopically equivalent if there exists a continuous function
$$h: \mathbb{R} \times [0, 1] \rightarrow \mathbb{R}$$
such that $h(t, 0) = x_1(t)$ and $h(t, 1) = x_2(t), t \in \mathbb{R}$.
In the context of motion planning, it is said that two trajectories (with a fixed start and end state) belong to the same homotopy class if one trajectory can be \emph{continuously deformed} into the other without intersecting an obstacle~\cite{Bhattacharya-2012}. This deformation is given by the homotopy function $h$.

We use the concept of \textit{h-signature} to identify the frequency at which the robot changes the homotopy class of its trajectory. An h-signature is a unique identifier for a homotopy class in an environment. Two trajectories have the same h-signature if and only if they belong to the same homotopy class~\cite{Grigoriev-98}. The h-signature for a trajectory is determined by extending a vertical ray upwards towards \textit{y = +$\infty$} from the center of each obstacle and assigning a letter $t_k$ to obstacle $k$. We iterate through each point in the trajectory $x$ from start to goal and check at every time step $t$ if the trajectory crosses a ray between $t$ and $t - 1$. If the trajectory crosses obstacle $k$ from left to right we append $t_k$ to the h-signature of the trajectory and if the trajectory crosses obstacle $k$ from right to left we append $\bar{t_k}$ to the trajectory. The final h-signature is produced by reducing this signature by removing all instances of $t_k\bar{t_k}$ or $\bar{t_k}t_k$ in the signature. \cite{Grigoriev-98} proves that this procedure yields a signature that uniquely corresponds to a homotopy class in the environment.

By identifying trajectories in different homotopy classes, \text{POSH}\xspace is finding solutions to the problem that do not deform into each other. This is useful because trajectory-optimization methods like GPMP-GRAPH typically will not deform a solution through an obstacle, preferring to get caught in a local minima in the current homotopy class. By designing \text{POSH}\xspace such that trajectories in multiple homotopy-classes are maintained, the planner will better be able to adapt to changes in the environment that render the current trajectory suboptimal.

\section{Planning Online by Switching Homotopies}
We begin by motivating \text{POSH}\xspace with our key insight that switching homotopy classes is essential to navigating changing environments. Then we present our algorithm and its implementation details.

\subsection{Motivation}
In Fig.~\ref{fig:narrow_passage} we see GPMP-GRAPH try to navigate through a 2-Dimensional environment with a moving obstacle. While at time step $1$ the robot seems to have found a high quality trajectory, changes in the environment start to degrade that trajectory. Namely, reoptimizing a single chain means that the robot will consider only its current homotopy-class, as deforming the trajectory through an obstacle would lead to an increase in trajectory cost. As a consequence, the trajectory is pushed further into a corner as the obstacle approaches, leading to multiple collisions in the trajectory. Intuitively, it seems reasonable to travel along the left of the moving obstacle if it is moving closer, however such a solution would be in a different homotopy class. For this reason, we should maintain and reoptimize trajectories in both the homotopy classes, the one to the right of the obstacle and the one to the left of the obstacle. Then, at every time step the robot can replan a solution to the goal state using the reoptimized chains of the factor graph. As the obstacle moves closer, the low-quality trajectory to the right of the obstacle will suffer a higher cost than a trajectory that moves to the left of the obstacle. Thus, the robot will now respond to this change in trajectory quality by changing the homotopy class mid-execution to avoid collision.

\subsection{Algorithm}

\begin{figure*}[ht]
    \centering
    \begin{subfigure}[t]{0.241\linewidth}
        \includegraphics[trim={50 80 50 100},clip,width=\linewidth]{figs/plotA.eps}
        \caption{}\label{fig:unoptimized}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={50 70 50 100},clip,width=\linewidth]{figs/plotB.eps}
        \caption{}\label{fig:optimized}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={50 70 50 100},clip,width=\linewidth]{figs/plotC.eps}
        \caption{}\label{fig:prune}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={50 70 50 100},clip,width=\linewidth]{figs/plotD.eps}
        \caption{}\label{fig:new_htpy}
    \end{subfigure}
    \caption{\small In (a), the gray circle and green circle represent the robot and goal state respectively. The graph is warm-started from its previous state. In~(b), the graph is reoptimized and then A* is used to find a low-cost trajectory (red) on the graph. In~(c), edges that will no longer be traversable (dashed blue lines) after executing the first control are pruned. Finally in~(d), the robot executes its first control as the environment changes. We see that \text{POSH}\xspace finds a trajectory in a new homotopy class (red) after reoptimizing.}
    \label{fig:toy_breakdown}
	\vspace{-2mm}
\end{figure*}

\text{POSH}\xspace is summarized by Algorithm 1. The initial obstacle positions are obtained to create the initial SDF as shown by line 1. In line 2, a factor graph is constructed given initial robot states, a signed distance field (SDF) of the environment, and parameters like smoothness prior (desired trade-off between smoothness and collision avoidance) and ratio of interconnections (ranges from 0 i.e. no connections to 1 i.e. every state is interconnected). Given a number of desired trajectories (chains) to be interconnected and a number of time steps, the locations of the initial support states are computed by constructing a series of ellipses with the major radius stretching from the start state to the goal state and varying minor radii as shown in Fig.~\ref{fig:nointerconnection}. Based on the desired ratio of interconnections, states in spatially nearby chains are then interconnected to allow switching between them (see Fig.~\ref{fig:halfinterconnection} and Fig.~\ref{fig:interconnection}). In line 4, given the initial graph and SDF, \text{POSH}\xspace employs the Levenberg-Marquardt algorithm to optimize the factor graph.

\begin{algorithm}[!t]
\caption{POSH}\label{alg:plosh}
\begin{algorithmic}[1]
\State $SDF\gets initial\_obstacle\_poses$
\State $factor\_graph \gets \{parameters, robot\_states, SDF\}$
\For{$t = 1, \dots, T$}
    \State $factor\_graph.optimize()$
    \State $x \gets A^*(factor\_graph)$
    \State $factor\_graph.prune\_unreachable()$
    \State $execute\_step(x)$
    \State $SDF.update(obstacle\_poses)$
\EndFor
\end{algorithmic}
\end{algorithm}

Up to this step, \text{POSH}\xspace has executed the same steps as GPMP-GRAPH would in solving the motion planning problem. Following this, GPMP-GRAPH simply executes the steps in the most feasible trajectory initially obtained and re-optimize this same trajectory at every time step until the goal configuration is reached. However, as discussed earlier this is insufficient to successfully plan in a stochastic, dynamic environment. POSH reoptimizes the entire graph with an updated SDF and re-calculates the most feasible trajectory. These steps, shown in lines 4-8 in Algorithm 1, are executed at every time step until the goal state is reached.

We use the example in Fig.~\ref{fig:toy_breakdown} to illustrates the process \text{POSH}\xspace follows. Fig.~\ref{fig:unoptimized} shows the graph before an iteration of optimization. The robot is the gray circle in the bottom left and the green circle in the top right is the goal state. The graph's current state represents the optimal graph configuration found given the SDF at the last time step. In Fig.~\ref{fig:optimized}, \text{POSH}\xspace re-optimizes the factor graph based on the updated SDF. Then, as shown in line 5 of Algorithm 1, A* is run on this updated, optimized factor graph to return the lowest cost trajectory (plotted in red). Because A* is complete and optimal, and the graph connects the start and goal states, it is guaranteed that the returned solution will be the lowest cost solution out of possible feasible trajectories encoded in the graph.

\begin{figure}[!t]
	\centering
    \begin{subfigure}[t]{0.32\linewidth}
        \includegraphics[trim={380 50 330 100},clip,width=\linewidth]{figs/nointer.eps}
        \caption{}\label{fig:nointerconnection}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.32\linewidth}
        \includegraphics[trim={380 50 330 100},clip,width=\linewidth]{figs/halfinter.eps}
        \caption{}\label{fig:halfinterconnection}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.32\linewidth}
        \includegraphics[trim={380 50 330 100},clip,width=\linewidth]{figs/fullinter.eps}
        \caption{}\label{fig:interconnection}
    \end{subfigure}%
    \caption{\small Possible graph initializations with~(a) no connection between chains,~(b) connections between chains at every other state, and~(c) fully connected chains. The robot is the gray circle and the goal state is the green circle.}
	\vspace{-3mm}
\end{figure}

This is followed by line 6, in which \text{POSH}\xspace runs depth first search (DFS) to prune away a subgraph of the factor graph to only maintain the reachable states. Given the current state as the root and the goal state, DFS determines which support states are unreachable in time and are therefore irrelevant. In Fig.~\ref{fig:prune}, the pruned away subgraph is shown with dashed blue lines. By executing the first step (line 7 of Algorithm 1) in the trajectory returned by A* in Fig.~\ref{fig:prune}, the robot arrives at its new position depicted in Fig.~\ref{fig:new_htpy}. At the next time step, the obstacle moves closer to the trajectory the robot was planning to take as Fig.~\ref{fig:new_htpy} shows. Then in line 8 the SDF is updated given changes in the environment. After running another iteration of optimization and A*, the robot now identifies a new trajectory in a different homotopy class that is lower cost and feasible. By doing so, \text{POSH}\xspace avoids the local minima problem that GPMP-GRAPH is unable to contend with as seen in Fig.~\ref{fig:narrow_passage}, resulting in a collision. On the other hand, \text{POSH}\xspace successfully adapts to the dynamic environment by switching homotopies. This process repeats until the robot has reached the goal.

\subsection{Implementation Details}

\text{POSH}\xspace is implemented by building on the GTSAM~\cite{frank2012factor} and GPMP2 library. The factor graph is optimized using the Levenberg-Marquardt algorithm initialized with $\lambda = 10^{-5}$. In order to successfully plan in a stochastic, dynamic environment, \text{POSH}\xspace must optimize the factor graph based on an updated SDF and determine the most feasible trajectory. Only the first step is executed as this trajectory might not be feasible for future configurations of the environment. To replan efficiently at the next time step, the relevant portions of the optimized factor graph must be retained so that the graph is not reconstructed from scratch. \text{POSH}\xspace accomplishes this by generating a spanning tree of the graph rooted at the first state in the most feasible trajectory, before navigating to it. The subgraph of the factor graph that does not belong to this spanning tree is pruned away as the states in this subgraph become unreachable in the next time step. Additionally, to replan quicker, the next optimization of the graph is warm-started with the state of the previously initialized factor graph. This results in Levenberg-Marquardt taking less time to optimize the factor graph, since the initialization of the graph was a local minima of the cost function at the last time step, assuming the environment did not drastically change between time steps.

\section{Results}\label{sec:results}

We benchmark \text{POSH}\xspace against GPMP2 and GPMP-GRAPH, two trajectory optimization methods within the same family on multiple 2D dynamic environments\footnote{A video of the experiments is available at \newline \indent \url{https://youtu.be/FtUI9VR3iWI}}. GPMP2 and GPMP-GRAPH are set up as a single chain and interconnected multi-chain graphs respectively. At each time step the GPMP2 chain is optimized and a step is taken on the optimized path. Then, the past state is pruned, the noisy state measurement is added, and the chain is reoptimized for the next time step. Since GPMP-GRAPH is an offline batch method, to keep comparisons fair, we will optimize the full multi-chain graph at the first time step to get an optimal sequence of states (a chain), while the other states are pruned away. From then on it is treated as a GPMP2 chain and is optimized accordingly at every time step. In contrast, \text{POSH}\xspace will continue to keep a pruned graph after every time step to be reoptimized.

These algorithms' performance are tested on two simulated 2D datasets, a dynamic forest and narrow passageway environment. In both environments, the holonomic robot, starting from some initial state, is tasked with reaching the goal while avoiding obstacles. The robot is also subjected to stochastic executions and noisy measurements i.e. the robot's position at every time step is perturbed to simulate execution noise and localization noise is injected into state measurements respectively.
Sections IV-A and IV-B will describe the narrow passageway and dynamic forest environments respectively as well as analyze the algorithms' performance in each. We record the following metrics averaged across multiple runs in each environment: (i) success rate to quantify the runs where the robot gets to the goal without any collisions, (ii) collision intensity to quantify the percentage of the trajectory spent in collisions when they occur, (iii) total distance traveled indicative of smoothness in the path, and (iv) total number of homotopy switches made from reoptimizing at each time step.

For all algorithms, the safety distance, $\epsilon$, and obstacle cost weight, $\sigma_{obs}$ (represents the trade-off between a trajectory smoothness and collision avoidance), are kept constant. The execution prior (initial state) and localization noise levels are all kept constant as well. To account for stochasticity, each algorithm's performance was averaged over 10 Monte-Carlo runs for each experiment. One narrow passageway environment and 13 dynamic forest environments were used to test the algorithms' performance. In the case of GPMP-GRAPH and \text{POSH}\xspace, hyperparameters like, the number of chains $N_I$, the ratio of interconnections between chains $R_I$, and GP prior noise applied to the interconnections between chains $Q_I$, are tuned to give good performance for the respective algorithm.

\subsection{Dynamic Narrow Passageway Benchmark}

\begin{table}[!t]
	\caption{Benchmark results on narrow passage dataset.}
	\label{table:NP-methodsResults}
	\begin{center}
		\begin{tabular}[c]{rccc}
			\toprule
			Metric & \text{POSH}\xspace & GPMP-GRAPH & GPMP2 \\
			\midrule
			Success Rate (\%) & \textbf{100} & 10 & 10\\
			Collision Intensity (\%) & \textbf{0} & 6.5 & 7.5\\
			Distance (m) & \textbf{26.7} & 31.6 & 33.8\\
			Homotopy Switches & \textbf{1.2} & 0.2 & 1.1\\
			\bottomrule
		\end{tabular}
	\end{center}
	\vspace{-3mm}
\end{table}

\begin{figure*}[ht]
    \centering
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={60 80 100 90},clip,width=\linewidth]{figs/GRAPHOpt1Passage.eps}
        \caption{$t=1$}\label{fig:t1}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={60 80 100 90},clip,width=\linewidth]{figs/GRAPHOpt4Passage.eps}
        \caption{$t=4$}\label{fig:t6}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={60 80 100 90},clip,width=\linewidth]{figs/GRAPHOpt5Passage.eps}
        \caption{$t=5$}\label{fig:t8}
    \end{subfigure}
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={60 80 100 90},clip,width=\linewidth]{figs/GRAPHOpt11Passage.eps}
        \caption{$t=10$}\label{fig:tn}
    \end{subfigure}%
    \caption{\small GPMP-GRAPH in (a) finds the lowest-cost trajectory (red) from the initial graph construction of 2 chains. It optimizes and follows this path at each time step, shown in green. In (b)-(c) the optimized path (red) continues to get pushed laterally due to the moving object, forcing the robot to collide with the obstacle. The final trajectory in (d) shows the robot colliding with the passageway as a result of being unable to successfully plan around the moving obstacle.}
    \label{fig:narrow_passage}
\end{figure*}

The dynamic narrow passageway consists of an obstacle moving back-and-forth with constant velocity in a 1D motion (left to right then right to left) in front of a narrow passageway that contains two entrances as shown in Fig.~\ref{fig:narrow_passage}. As this obstacle moves in front of the passageway it will alternately block one of the entrances. This environment is used to test the algorithms' online planning capabilities in a controlled setting. As we argued previously the ability to switch homotopies will become important in such a scenario.

The results are summarized in Table~\ref{table:NP-methodsResults}. \text{POSH}\xspace is able to successfully plan 100\% in narrow passage environment, while executing on average 1-2 homotopy switches to explore the both passage entrances. As a result, it also minimizes distance traveled due to this early, prudent homotopy switch. On the other hand, GPMP2 and GPMP-GRAPH are largely unsuccessful at solving this problem as they fail to successfully find collision free paths. While at first, the optimal path from start to goal is obvious for GPMP2, the moving obstacle quickly eliminates this as a viable path, forcing the algorithms' to improvise or fail. They get stuck in local minima as seen in Fig.~\ref{fig:narrow_passage} where the object's lateral motion continues to push the trajectory to the side until finding a collision free path is no longer possible under the given temporal and physical constraints. GPMP-GRAPH's initial solution that spans multiple chains also succumbs to the same issue since it is reduced to a single chain graph after the first time step. On the other hand, Fig.~\ref{fig:toy_breakdown} shows that \text{POSH}\xspace contends with this local minima problem by utilizing a multi-chain factor graph to explore a path in the other homotopy class once it becomes more viable than the previously desirable homotopy.
It is important to note that the homotopy switches made by GPMP-GRAPH and GPMP2 are a result of the dynamic environment changing the optimized path's \textit{h-signature} rather than an intended algorithmic or strategic homotopy switch.

\subsection{Dynamic Forest Benchmark}

\begin{figure*}[ht]
    \centering
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={80 80 50 90},clip,width=\linewidth]{figs/rf_t1.eps}
        \caption{$t = 1$}\label{fig:rf_t1}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={80 80 50 90},clip,width=\linewidth]{figs/rf_t2.eps}
        \caption{$t = 6$}\label{fig:rf_t2}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={80 80 50 90},clip,width=\linewidth]{figs/rf_t3.eps}
        \caption{$t = 10$}\label{fig:rf_t3}
    \end{subfigure}%
	\hfill\vline\hfill
    \begin{subfigure}[t]{0.24\linewidth}
        \includegraphics[trim={80 80 50 90},clip,width=\linewidth]{figs/rf_t4.eps}
        \caption{$t = 21$}\label{fig:rf_t4}
    \end{subfigure}%
    \caption{\small \text{POSH}\xspace successfully returns a collision free trajectory utilizing 4 chains. At each consecutive time-step $t$, the robot (gray circle) prunes away unreachable portion of the graph and reoptimizes the remaining graph until it reaches the goal state (green circle). The path the robot has traveled so far is shown green and the lowest-cost trajectory going forward is shown in red.}
    \label{fig:dynamic_forest_timelapse}
	\vspace{-5mm}
\end{figure*}

The dynamic forest environment consists of numerous objects exhibiting random 2D motion. This environment is simulated given a number of obstacles, their initial state and a random bounded acceleration applied at each time step. The dynamic forest is shown in Fig.~\ref{fig:dynamic_forest}. It is employed to benchmark the ability of the algorithms to contend with a harder replanning problem.

The dynamic forest benchmark results are shown in Table~\ref{table:DF-methodsResults}. Our algorithm generates paths that are more successful and when it fails it spends less time on average in collision. Again, this can be attributed to \text{POSH}\xspace constantly replanning and reoptimizing the multi-chain factor graph instead of sticking with a single chain factor graph, allowing for online homotopy switching. This process can be seen in Fig.~\ref{fig:dynamic_forest_timelapse} and Fig.~\ref{fig:dynamic_forest}. Since there are many more paths for \text{POSH}\xspace to consider, it returns a feasible path more often than the two other algorithms.
In Fig.~\ref{fig:dynamic_forest}, the black line represents the trajectory considered best at the previous time step while the red line represents the trajectory considered best at the current time step. Fig.~\ref{fig:switch1} signifies a homotopy switch as the red and black lines cannot be continuously deformed into each other without colliding with obstacles. Fig.~\ref{fig:switch2} shows the robot executing the first step of the red path from Fig.~\ref{fig:switch1} while again considering another homotopy switch.

Note that GPMP2's on average shorter distance traveled can be attributed to the straight line initialization that characterizes the region around the local minima, but this also results in a greater collision intensity and a much lower success rate.

\begin{table}[!t]
	\caption{Benchmark results on dynamic forest dataset.}
	\label{table:DF-methodsResults}
	\begin{center}
	\resizebox{\columnwidth}{!}{
		\begin{tabular}[c]{rccc}
			\toprule
			Metric & \text{POSH}\xspace & GPMP-GRAPH & GPMP2 \\
			\midrule
			Success Rate (\%) & \textbf{29.23} & 11.92 & 8.46\\
			Collision Intensity (\%) & \textbf{1.9} & 5.1 & 8.3\\
			Distance (m) & 59.05 & 66.52 & \textbf{53.41}\\
			Homotopy Switches & \textbf{5.15} & 4.23 & 3.54\\
			Avg. Comp. Time (s) at $t=0$ & 0.221 & 0.020 & \textbf{0.0197}\\
			Avg. Comp. Time (s) at any $t>0$ & 0.077 & 0.029 & \textbf{0.028}\\
			\bottomrule
		\end{tabular}
	}
	\end{center}
	\vspace{-2mm}
\end{table}

\begin{table}[!t]
	\caption{Performance of \text{POSH}\xspace with varying number of chains on dynamic forest dataset.}
	\label{table:DF_numPathsResults}
	\begin{center}
		\begin{tabular}[c]{rccc}
		\toprule
		Metric & $N_I$=2 & $N_I$=4 & $N_I$=6 \\
		\midrule
		Success Rate (\%) & 22.31 & \textbf{29.23} & 18.85\\
		Collision Intensity (\%) & \textbf{1.9} & \textbf{1.9} & 2.2\\
		Distance (m) & 59.1 & \textbf{59.0} & 60.6\\
		Homotopy Switches & \textbf{5.23} & 5.15 & 4.38\\
		\bottomrule
		\end{tabular}
	\end{center}
	\vspace{-3mm}
\end{table}

We also analyze the performance of \text{POSH}\xspace on factor graphs with different numbers of chains, $N_I=2$, $N_I=4$, and $N_I=6$. The results in Table~\ref{table:DF_numPathsResults} indicate that there is a sweet spot in the middle (at $N_I=4$) to get the best performance. An example with four chains is visualized in Fig.~\ref{fig:dynamic_forest_timelapse}. The two chain factor graph will have fallen into two of the four homotopy classes that the four chain factor graph will have, thus exhibiting similar metrics. The six chain factor graph exhibits a reduction in performance due to chains collapsing into redundant homotopy classes as a result of $Q_I$ and the number of time steps making the optimization problem too constrained. It is also important to consider the algorithmic computation time when deciding on the number of chains to include in the factor graph as increasing number of chains increases computation time. Ideally the computation time should be smaller than the time length of the trajectory between two time steps ($\delta t$) for \text{POSH}\xspace to be applicable online, which was indeed the case in our experiments. Table~\ref{table:DF-methodsResults} shows the average computation times for the first iteration at $t=0$ and for any other iteration at $t>0$, where $\delta t = 0.5$s.

\section{Discussion}\label{sec:conc}

In this paper, we introduced \text{POSH}\xspace, a novel online graph-based trajectory optimization planning algorithm. It builds on recent inference based trajectory optimization techniques to add capability to plan online and dynamically switch homotopies. We found that \text{POSH}\xspace can leverage knowledge of distinct trajectories across different homotopy classes to dynamically select a low cost solution in light of new information about the environment. In dynamic environments with moving obstacles, previous trajectory-optimization approaches are collision-prone as they become stuck in a local minima. Rather than committing to a low-quality trajectory in a local minima, \text{POSH}\xspace successfully generates smooth and collision-free trajectories by maintaining and optimizing multiple interconnected solutions to the planning problem.

Our work is currently limited to a finite time horizon setting i.e. the number of discretized time steps we are planning over is fixed. For long range tasks in dynamic environments, occasionally a receding horizon formulation is more beneficial to avoid wasting computation on latter portions of the plan far ahead in the future. Extending our approach to such settings can further improve its capability and applicability. \text{POSH}\xspace (and GPMP-GRAPH) can sometimes suffer from different trajectories collapsing into the same homotopy class. To contend with this problem it will be beneficial to further investigate principled ways of integrating sampling strategies for new support states such that \text{POSH}\xspace not only prunes unreachable states, but also add new states at every time step to allow for diversity and exploration.

\vspace{-1mm}
\bibliographystyle{IEEEtran}