path: root/stokes.tex

\documentclass[aps,reprint,longbibliography,floatfix]{revtex4-2}

\usepackage[utf8]{inputenc} % why not type "Stokes" with unicode?
\usepackage[T1]{fontenc} % vector fonts
\usepackage[
  colorlinks=true,
  urlcolor=purple,
  citecolor=purple,
  filecolor=purple,
  linkcolor=purple
]{hyperref} % ref and cite links with pretty colors
\usepackage{amsmath, amssymb, graphicx, xcolor} % standard packages
\usepackage{newtxtext,newtxmath} % Times for PR
\usepackage[subfolder]{gnuplottex} % need to compile separately for APS

\begin{document}

\title{Analytic continuation over complex landscapes}

\author{Jaron Kent-Dobias}
\author{Jorge Kurchan}

\affiliation{Laboratoire de Physique de l'Ecole Normale Supérieure, Paris, France}

\date\today

\begin{abstract} \setcitestyle{authoryear,round}
  In this paper we follow up the study of `complex complex landscapes'
  [\cite{Kent-Dobias_2021_Complex}], rugged landscapes of many complex
  variables.  Unlike real landscapes, there is no useful classification of
  saddles by index. Instead, the spectrum at critical points determines their
  tendency to trade topological numbers under analytic continuation of the
  theory. These trades, which occur at Stokes points, proliferate when the
  spectrum includes marginal directions and are exponentially suppressed
  otherwise. This gives a direct interpretation of the `threshold' energy---which
  in the real case separates saddles from minima---where the spectrum of
  typical critical points develops a gap. This leads to different consequences
  for the analytic continuation of real landscapes with different structures:
  the global minima of ``one step replica-symmetry broken'' landscapes lie
  beyond a threshold and are locally protected from Stokes points, whereas
  those of ``many step replica-symmetry broken'' lie at the threshold and
  Stokes points immediately proliferate.
\end{abstract}

\maketitle

Consider an action $\mathcal S_\lambda$ defined on the phase space $\Omega$ and
depending on parameters $\lambda$. In the context of statistical mechanics,
$\mathcal S_{\beta,J}=-\beta H_J$ for some hamiltonian $H_J$ with quenched
parameters $J$ at inverse temperature $\beta$. A typical calculation stems from
the partition function
\begin{equation} \label{eq:partition.function}
  Z(\lambda)=\int_\Omega ds\,e^{\mathcal S_\lambda(s)}.
\end{equation}
This integral is often dominated by its behavior near stationary points of the
action, and understanding these points is usually important to evaluate the
partition function.

Recent developments have found that stationary points of the action are
important for understanding another aspect of the partition function: its
analytic continuation. The integral \eqref{eq:partition.function} is first
interpreted as a contour in a larger complex phase space, then deformed into a
linear combination of specially constructed contours each enumerated by a
stationary point. Analytic continuation of parameters preserves this
decomposition except at nongeneric points where contours intersect.

We investigate the plausibility of analytic continuation in complex models
where the action has a macroscopic number of stationary points. Such actions
are common in studies of glasses, spin glass, machine learning, black holes,
\dots We find that the geometry of the landscape, and in particular the
relative position and spectrum of stationary points, is key.

Analytic continuation of partition functions is useful for many reasons. Some
physically motivated theories have actions whose partition function is formally
infinite, but can be defined by continuing from a set of parameters where it
converges. Other theories have oscillatory actions that lead to a severe sign
problem in estimating the partition function, which can be addressed by taking
advantage of a deformed phase space where the phase of the action varies slowly.

Unfortunately the study is not so relevant for low-dimensional `rugged'
landscapes, which are typically constructed from the limits of series or
integrals of analytic functions which are not themselves analytic
\cite{Cavagna_1999_Energy}.

\section{Integration by Lefschetz thimble}

We return to the partition function \eqref{eq:partition.function}.  If
the action can be continued to a holomorphic function on the Kähler
manifold $\tilde\Omega\supset\Omega$ and $\Omega$ is orientable in $\tilde\Omega$,
then \eqref{eq:partition.function} can be considered a contour integral. In
this case, the contour can be freely deformed without affecting the value of
the integral. Two properties of this deformed contour would be ideal. First,
that as $|s|\to\infty$ the real part of the action goes to $-\infty$, to ensure
the integral converges. Second, that the contours piecewise correspond to
surfaces of slowing vary phase of the action, so as to ameliorate sign
problems.

Remarkably, there is an elegant recipe for accomplishing both these criteria at
once, courtesy of Morse theory. For a more thorough review, see
\citet{Witten_2011_Analytic}. Consider a stationary point of the action. The
union of all gradient descent trajectories on the real part of the action that
begin at the stationary point is known as a \emph{Lefschetz thimble}. Since
each point on the Lefschetz thimble is found through descent from the
stationary point, the real part of the action is bounded from above by its
value at the stationary point. Likewise, we shall see that the imaginary part
of the action is constant on a thimble.

Morse theory provides a universal correspondence between contours and thimbles.
For any contour $\Omega$, there exists a linear combination of thimbles such
that the relative homology of the combination with respect to decent int he
action is equivalent to that of the contour.  If $\Sigma$ is the set of
stationary points of the action and for each $\sigma\in\Sigma$ the set
$\mathcal J_\sigma\subset\tilde\Omega$ is its thimble, then this gives
\begin{equation} \label{eq:thimble.integral}
  Z(\lambda)=\sum_{\sigma\in\Sigma}n_\sigma\int_{\mathcal J_\sigma}ds\,e^{\mathcal S_\lambda(s)}.
\end{equation}
Each of these integrals is very well-behaved: convergent asymptotic series
exist for their value about each critical point. The integer weights $n_\sigma$
are fixed by comparison with the initial contour. For a real action, all maxima
in $\Omega$ contribute in equal magnitude.

Under analytic continuation, the form of \eqref{eq:thimble.integral}
generically persists. When the relative homology of the thimbles is unchanged
by the continuation, the integer weights are likewise unchanged, and one can
therefore use the knowledge of these weights in one regime to compute the
partition function in the other. However, their relative homology can change,
and when this happens the integer weights can be traded between critical
points. These trades occur when two thimbles intersect, or alternatively when
one stationary point lies in the gradient descent of another. These places are
called \emph{Stokes points}, and the gradient descent trajectories that join
two stationary points are called \emph{Stokes lines}.

The prevalence (or not) of Stokes points in a given continuation, and whether
those that do appear affect the weights of critical points of interest, is a
concern for the analytic continuation of theories. If they do not occur or
occur order-one times, one could reasonably hope to perform such a procedure.
If they occur exponentially often, there is little hope of keeping track of the
resulting weights.

\section{Gradient descent dynamics}

The `dynamics' describing thimbles is defined by gradient descent on the real
part of the action.
\begin{equation} \label{eq:flow.coordinate.free}
  \dot s
  =-\operatorname{grad}\operatorname{Re}\mathcal S
  =-\left(\frac\partial{\partial s^*}\operatorname{Re}\mathcal S\right)^\sharp
  =-\frac12\frac{\partial\mathcal S^*}{\partial s^*}g^{-1}\frac\partial{\partial s},
\end{equation}
where $g$ is the metric and the holomorphicity of the action was used to set
$\partial^*\mathcal S=0$.

We will be dealing with actions where it is convenient to refer to coordinates
in a higher-dimensional embedding space. Let $z:\tilde\Omega\to\mathbb C^N$ be
an embedding of phase space into complex euclidean space. This gives
\begin{equation}\label{eq:flow.raw}
  \dot z
  =-\frac12\frac{\partial\mathcal S^*}{\partial z^*}(Dz)^* g^{-1}(Dz)^T\frac\partial{\partial z}
\end{equation}
where $Dz=\partial z/\partial s$ is the Jacobian of the embedding.
The embedding induces a metric on $\tilde\Omega$ by $g=(Dz)^\dagger Dz$.
Writing $\partial=\partial/\partial z$, this gives
\begin{equation} \label{eq:flow}
  \dot z=-\tfrac12(\partial\mathcal S)^\dagger(Dz)^*[(Dz)^\dagger(Dz)]^{-1}(Dz)^T
  =-\tfrac12(\partial \mathcal S)^\dagger P
\end{equation}
which is nothing but the projection of $(\partial\mathcal S)^*$ into the
tangent space of the manifold, with $P=(Dz)^*[(Dz)^\dagger(Dz)]^{-1}(Dz)^T$.
Note that $P$ is hermitian.

Gradient descent on $\operatorname{Re}\mathcal S$ is equivalent to Hamiltonian
dynamics with the Hamiltonian $\operatorname{Im}\mathcal S$. This is because
$(\tilde\Omega, g)$ is Kähler and therefore admits a symplectic structure, but
that the flow conserves $\operatorname{Im}\mathcal S$ can be shown using
\eqref{eq:flow} and the holomorphic property of $\mathcal S$:
\begin{equation}
  \begin{aligned}
    \frac d{dt}&\operatorname{Im}\mathcal S
    =\dot z\partial\operatorname{Im}\mathcal S+\dot z^*\partial^*\operatorname{Im}\mathcal S \\
    &=\frac i4\left(
      (\partial \mathcal S)^\dagger P\partial\mathcal S-(\partial\mathcal S)^TP^*(\partial\mathcal S)^*
    \right) \\
    &=\frac i4\left(
      (\partial\mathcal S)^\dagger P\partial\mathcal S-[(\partial\mathcal S)^\dagger P\partial\mathcal S]^*
    \right) \\
    &=\frac i4\left(
      \|\partial\mathcal S\|^2-(\|\partial\mathcal S\|^*)^2
    \right)=0.
  \end{aligned}
\end{equation}
As a result of this conservation law, surfaces of constant imaginary action
will be important when evaluting the possible endpoints of trajectories. A
consequence of this conservation is that the flow in the action takes a simple
form:
\begin{equation}
  \dot{\mathcal S}
  =\dot z\partial\mathcal S
  =-\frac12(\partial\mathcal S)^\dagger P\partial\mathcal S
  =-\frac12\|\partial\mathcal S\|^2.
\end{equation}
In the complex-$\mathcal S$ plane, dynamics is occurs along straight lines in
the negative real direction.

Let us consider the generic case, where the critical points of $\mathcal S$ have
distinct energies. What is the topology of the $C=\operatorname{Im}\mathcal S$ level
set? We shall argue its form by construction. Consider initially the situation
in the absence of any critical point. In this case the level set consists of a
single simply connected surface, locally diffeomorphic to $\mathbb R^{2D-1}$. Now, `place' a generic
(nondegenerate) critical point in the function at $z_0$. In the vicinity of the critical
point, the flow is locally
\begin{equation}
  \begin{aligned}
    \dot z
    &\simeq-\tfrac12(\partial\partial\mathcal S|_{z=z_0})^\dagger P(z-z_0)^*
  \end{aligned}
\end{equation}
The matrix $(\partial\partial\mathcal S)^\dagger P$ has a spectrum identical to that of
$(\partial\partial\mathcal S)^\dagger$ save marginal directions corresponding to the normals to
manifold. Assuming we are working in a diagonal basis, this becomes
\begin{equation}
  \dot z_i=-\tfrac12\lambda_i\Delta z_i^*+O(\Delta z^2,(\Delta z^*)^2)
\end{equation}
Breaking into real and imaginary parts gives
\begin{equation}
  \begin{aligned}
    \frac{d\Delta x_i}{dt}&=-\frac12\left(
      \operatorname{Re}\lambda_i\Delta x_i+\operatorname{Im}\beta\lambda_i\Delta y_i
    \right) \\
    \frac{d\Delta y_i}{dt}&=-\frac12\left(
      \operatorname{Im}\lambda_i\Delta x_i-\operatorname{Re}\beta\lambda_i\Delta y_i
    \right)
  \end{aligned}
\end{equation}
Therefore, in the complex plane defined by each eigenvector of
$(\partial\partial\mathcal S)^\dagger P$ there is a separatrix flow of the form in
Figure \ref{fig:local_flow}. The effect of these separatrices in each complex
direction of the tangent space $T_{z_0}M$ is to separate that space into four
quadrants: two disconnected pieces with greater imaginary part than the
critical point, and two with lesser imaginary part. This partitioning implies
that the level set of $\operatorname{Im}\mathcal S=C$ for
$C\neq\operatorname{Im}\mathcal S(z_0)$ is split into two disconnected pieces, one
lying in each of two quadrants corresponding with its value relative to that at
the critical point.

\begin{figure}
  \includegraphics[width=\columnwidth]{figs/local_flow.pdf}
  \caption{
    Gradient descent in the vicinity of a critical point, in the $z$--$z^*$
    plane for an eigenvector $z$ of $(\partial\partial\mathcal S)^\dagger P$. The flow
    lines are colored by the value of $\operatorname{Im}H$.
  } \label{fig:local_flow}
\end{figure}

Continuing to `insert' critical points whose imaginary energy differs from $C$,
one repeatedly partitions the space this way with each insertion. Therefore,
for the generic case with $\mathcal N$ critical points, with $C$ differing in
value from all critical points, the level set $\operatorname{Im}\mathcal S=C$ has
$\mathcal N+1$ connected components, each of which is simply connected,
diffeomorphic to $\mathbb R^{2D-1}$ and connects two sectors of infinity to
each other.

When $C$ is brought to the same value as the imaginary part of some critical
point, two of these disconnected surfaces pinch grow nearer and pinch together
at the critical point when $C=\operatorname{Im}\mathcal S$, as in the black lines of
Figure \ref{fig:local_flow}. The unstable manifold of the critical point, which
corresponds with the portion of this surface that flows away, produce the
Lefschetz thimble associated with that critical point.

Stokes lines are trajectories that approach distinct critical points as time
goes to $\pm\infty$. From the perspective of dynamics, these correspond to
\emph{heteroclinic orbits}. What are the conditions under which Stokes lines
appear? Because the dynamics conserves imaginary energy, two critical points
must have the same imaginary energy if they are to be connected by a Stokes
line. This is not a generic phenomena, but will happen often as one model
parameter is continuously varied. When two critical points do have the same
imaginary energy and $C$ is brought to that value, the level set
$C=\operatorname{Im}\mathcal S$ sees formally disconnected surfaces pinch together in
two places. We shall argue that generically, a Stokes line will exist whenever
the two critical points in question lie on the same connected piece of this
surface.

What are the ramifications of this for disordered Hamiltonians? When some
process brings two critical points to the same imaginary energy, whether a
Stokes line connects them depends on whether the points are separated from each
other by the separatrices of one or more intervening critical points.
Therefore, we expect that in regions where critical points with the same
energies tend to be nearby, Stokes lines will proliferate, while in regions
where critical points with the same energies tend to be distant compared to
those with different energies, Stokes lines will be rare.

\textcolor{teal}{
  Here we make a generic argument that, for high-dimensional landscapes with
  exponentially many critical points, the existence of exponentially many
  Stokes points depends on the spectrum of the Hessian $\partial\partial H$ of
  critical points. 
}

\section{p-spin spherical models}

The $p$-spin spherical models are statistical mechanics models defined by the
action $\mathcal S=-\beta H$ for the Hamiltonian
\begin{equation} \label{eq:p-spin.hamiltonian}
  H(x)=\sum_{p=2}^\infty\frac{a_p}{p!}\sum_{i_1\cdots i_p}^NJ_{i_1\cdots i_p}x_{i_1}\cdots x_{i_p}
\end{equation}
where the $x\in\mathbb R^N$ are constrained to lie on the sphere $x^2=N$.  The tensor
components $J$ are complex normally distributed with zero mean and variances
$\overline{|J|^2}=p!/2N^{p-1}$ and $\overline{J^2}=\kappa\overline{|J|^2}$, and
the numbers $a$ define the model. The pure real $p$-spin model has
$a_i=\delta_{ip}$ and $\kappa=1$.

The phase space manifold $\Omega=\{x\mid x^2=N, x\in\mathbb R^N\}$ has a
complex extension $\tilde\Omega=\{z\mid z^2=N, z\in\mathbb C^N\}$. The natural
extension of the hamiltonian \eqref{eq:p-spin.hamiltonian} to this complex manifold is
holomorphic. The normal to this manifold at any point $z\in\tilde\Omega$ is
always in the direction $z$.  The projection operator onto the tangent space of
this manifold is given by
\begin{equation}
  P=I-\frac{zz^\dagger}{|z|^2},
\end{equation}
where indeed $Pz=z-z|z|^2/|z|^2=0$ and $Pz'=z'$ for any $z'$ orthogonal to $z$.

To find critical points, we use the method of Lagrange multipliers. Introducing $\mu\in\mathbb C$, 

\subsection{2-spin}

The Hamiltonian of the $2$-spin model is defined for $z\in\mathbb C^N$ by
\begin{equation}
  H_0=\frac12z^TJz.
\end{equation}
$J$ is generically diagonalizable by a complex orthogonal matrix. In a diagonal basis, $J_{ij}=\lambda_i\delta_{ij}$. Then $\partial_i H=\lambda_iz_i$. We will henceforth assume to be working in this basis. To find the critical points, we use the method of Lagrange multipliers. Introducing $\epsilon$, the constrained Hamiltonian is
\begin{equation}
  H=H_0+\epsilon(N-z^2)
\end{equation}
As usual, $\epsilon$ is equivalent to the energy per spin at any critical point.
Critical points must satisfy
\begin{equation}
  0=\partial_iH=(\lambda_i-2\epsilon)z_i
\end{equation}
which is only possible for $z_i=0$ or $\epsilon=\frac12\lambda_i$. Generically the $\lambda_i$ will all differ, so this can only be satisfied for one $\lambda_i$ at a time, and to be a critical point all other $z_j$ must be zero. In the direction in question,
\begin{equation}
  \frac1N\frac12\lambda_iz_i^2=\epsilon=\frac12\lambda_i,
\end{equation}
whence $z_i=\pm\sqrt{N}$. Thus there are $2N$ critical points, each corresponding to $\pm$ the cardinal directions in the diagonalized basis.

Suppose that two critical points have the same imaginary energy; without loss
of generality, assume these are associated with the first and second
cardinal directions. Since the gradient is proportional to $z$, any components that are
zero at some time will be zero at all times. The dynamics for the components of
interest assuming all others are zero are
\begin{equation}
  \begin{aligned}
    \dot z_1
    &=-z_1^*\left(d_1^*-\frac{d_1^*z_1^*z_1+d_2^*z_2^*z_2}{|z_1|^2+|z_2|^2}\right) \\
    &=-(d_1-d_2)^*z_1^*\frac{|z_2|^2}{|z_1|^2+|z_2|^2}
  \end{aligned}
\end{equation}
and the same for $z_2$ with all indices swapped.  Since $\Delta=d_1-d_2$ is
real, if $z_1$ begins real it remains real, with the same for $z_2$. Since the
critical points are at real $z$, we make this restriction, and find
\begin{equation}
  \begin{aligned}
    \frac d{dt}(z_1^2+z_2^2)=0\\
    \frac d{dt}\frac{z_2}{z_1}=\Delta\frac{z_2}{z_1}
  \end{aligned}
\end{equation}
Therefore $z_2/z_1=e^{\Delta t}$, with $z_1^2+z_2^2=N$ as necessary.  Depending
on the sign of $\Delta$, $z$ flows from one critical point to the other over
infinite time. This is a Stokes line, and establishes that any two critical
points in the 2-spin model with the same imaginary energy will possess one.
These trajectories are plotted in Fig.~\ref{fig:two-spin}.

\begin{figure}
  \begin{gnuplot}[terminal=epslatex, terminaloptions={size 8.65cm,5.35cm}]
    set xlabel '$\Delta t$'
    set ylabel '$z(t) / \sqrt{N}$'

    plot 1 / sqrt(1 + exp(2 * x)) t '$z_1$', \
         1 / sqrt(1 + exp(- 2 * x)) t '$z_2$'
  \end{gnuplot}
  \caption{
    The Stokes line in the 2-spin model when the critical points associated
    with the first and second cardinal directions are brought to the same
    imaginary energy. $\Delta$ is proportional to the difference between the
    real energies of the first and the second critical point; when $\Delta >0$
    flow is from first to second, while when $\Delta < 0$ it is reversed.
  } \label{fig:two-spin}
\end{figure}

Since they sit at the corners of a simplex, the critical points of the 2-spin
model are all adjacent: no critical point is separated from another by the
separatrix of a third. This means that when the imaginary energies of two
critical points are brought to the same value, their surfaces of constant
imaginary energy join.

\subsection{Pure \textit{p}-spin}

\begin{equation}
  H_p=\frac1{p!}\sum_{i_1\cdots i_p}J_{i_1\cdots i_p}z_{i_1}\cdots z_{i_p}
\end{equation}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex, terminaloptions={size 8.65cm,5.35cm}]
    set parametric
    set hidden3d
    set isosamples 100,25
    set samples 100,100
    unset key
    set dummy u,r
    set urange [-pi:pi]
    set vrange [1:1.5]
    set cbrange [0:2]
    set xyplane 0

    set xlabel '$\operatorname{Re}\epsilon$'
    set ylabel '$\operatorname{Im}\epsilon$'
    set zlabel '$r$'
    set cblabel '$\frac\epsilon{\epsilon_{\mathrm{th}}}$'

    p = 4
    set palette defined (0 "blue", 0.99 "blue", 1.0 "white", 1.01 "red", 2 "red")
    set pm3d depthorder border linewidth 0.5

    s(r) = sqrt(0.75 * log(9 * r**4 / (1 + r**2 + r**4)) / (8 * r**4 - r**2 - 1))
    x(u, r) = cos(u) * s(r) * sqrt(1 + 5 * r**2 + 5 * r**4 + r**6)
    y(u, r) = sin(u) * s(r) * sqrt((r**2 - 1)**3)
    thres(u, r) = ((x(u,r) / (r**(p - 2) + 1))**2 + (y(u,r) / (r**(p - 2) - 1))**2) / ((p - 1) / (2 * p * r**(p - 2)))

    splot "++" using (x(u, r)):(y(u, r)):2:(thres(u, r)) with pm3d lc palette
  \end{gnuplot}
  \caption{
    The surface of extant states for the 4-spin model, that is, those for which
    the complexity is zero.
  }
\end{figure}

\subsection{(2 + 4)-spin}

\begin{equation}
  H_2+H_4
\end{equation}

\section{Numerics}

To confirm the presence of Stokes lines under certain processes in the $p$-spin, we studied the problem numerically.


\begin{acknowledgments}
  MIT mathematicians have been no help
\end{acknowledgments}

\bibliography{stokes}

\appendix

\section{Geometry}

The surface $M\subset\mathbb C^N$ defined by $N=f(z)=z^2$ is an $N-1$ dimensional
\emph{Stein manifold}, a type of complex manifold defined by the level set of a
holomorphic function \cite{Forstneric_2017_Stein}. One can define a Hermitian
metric $h$ on $M$ by taking the restriction of the standard metric of $\mathbb
C^N$ to vectors tangent along $M$. For any smooth function $\phi:M\to\mathbb
R$, its gradient is a holomorphic vector field given by
\begin{equation}
  \operatorname{grad}\phi=\bar\partial^\sharp\phi
\end{equation}
Suppose $u:M\to\mathbb C^{N-1}$ is a coordinate system. Then
\begin{equation}
  \operatorname{grad}\phi=h^{\bar\beta\alpha}\bar\partial_{\bar\beta}\phi\frac{\partial}{\partial u^\alpha}
\end{equation}
Let $z=u^{-1}$.
\begin{equation}
  \frac\partial{\partial u^\alpha}=\frac{\partial z^i}{\partial u^\alpha}\frac\partial{\partial z^i}
\end{equation}
\begin{equation}
  \bar\partial_{\bar\beta}\phi=\frac{\partial\bar z^i}{\partial\bar u^{\bar\beta}}\frac{\partial\phi}{\partial\bar z^i}
\end{equation}
\begin{equation}
  \operatorname{grad}\phi
  =\frac{\partial\bar z^{\bar\jmath}}{\partial\bar u^{\bar\beta}}h^{\bar\beta\alpha}\frac{\partial z^i}{\partial u^\alpha}\frac{\partial\phi}{\partial\bar z^{\bar\jmath}}\frac\partial{\partial z^i}
\end{equation}
At any point $z_0\in M$, $z_0$ is parallel to the normal to $M$. Without loss of generality, take $z_0=e^N$. In a neighborhood of $z_0$, the map $u^\alpha=z^\alpha$ is a coordinate system.
its inverse is $z^i=u^i$ for $1\leq i\leq N-1$ and
\begin{equation}
  z^N=\sqrt{N-u^2}.
\end{equation}
The Jacobian is
\begin{equation}
  \frac{\partial z^i}{\partial u^\alpha}=\delta^i_\alpha-\delta_N^{i}\frac{u^\alpha}{\sqrt{N-u^2}}
\end{equation}
and therefore the Hermitian metric induced by the map is
\begin{equation}
  h_{\alpha\bar\beta}=\frac{\partial\bar z^i}{\partial\bar u^\alpha}\frac{\partial z^{\bar\jmath}}{\partial u^{\bar\beta}}\delta_{i\bar\jmath}
  =\delta_{\alpha\bar\beta}+\frac{\bar u^{\alpha}u^{\bar\beta}}{|N-u^2|}
\end{equation}
The metric can be inverted explicitly:
\begin{equation}
  h^{\bar\beta\alpha}
  =\delta^{\bar\beta\alpha}-\frac{\bar u^{\bar\beta}u^\alpha}{|N-u^2|+|u|^2}.
\end{equation}
Putting these pieces together, we find
\begin{equation}
  \frac{\partial\bar z^{\bar\jmath}}{\partial\bar u^{\bar\beta}}h^{\bar\beta\alpha}\frac{\partial z^i}{\partial u^\alpha}
  =\delta^{\bar\jmath i}-\frac{z^{\bar\jmath}\bar z^i}{|z|^2}
\end{equation}
\begin{equation}
  \operatorname{grad}\phi
  =\left(\delta^{\bar\jmath i}-\frac{z^{\bar\jmath}\bar z^i}{|z|^2}\right)
  \frac{\partial\phi}{\partial\bar z^{\bar\jmath}}\frac\partial{\partial z^i}
\end{equation}

\section{Numerics}

To study Stokes lines numerically, we approximated them by parametric curves.
If $z_0$ and $z_1$ are two stationary points of the action with
$\operatorname{Re}\mathcal S(z_0)>\operatorname{Re}\mathcal S(z_1)$, then we
take the curve
\begin{equation}
  z(t)
  =(1-t)z_0+tz_1+(1-t)t\sum_{i=0}^mg_it^i
\end{equation}
where the $g$s are undetermined complex vectors. These are fixed by minimizing
a cost function, which has a global minimum only for Stokes lines. Defining
\begin{equation}
  \mathcal L(t)
  = 1-\frac{\operatorname{Re}[\dot z^*(z(t))\cdot z'(t)]}{|\dot z(z(t))||z'(t)|}
\end{equation}
this cost is given by
\begin{equation}
  \mathcal C=\int_0^1 dt\,\mathcal L(t)
\end{equation}
$\mathcal C$ has minimum of zero, which is reached only by functions $z(t)$
whose tangent is everywhere parallel to the direction $\dot z$ of the dynamics.
Therefore, functions that minimize $\mathcal C$ are time-reparameterized Stokes
lines.

We explicitly computed the gradient and Hessian of $\mathcal C$ with respect to
the parameter vectors $g$. Stokes lines are found or not between points by
using the Levenberg--Marquardt algorithm starting from $g^{(i)}=0$ for all $i$,
and approximating the cost integral by a finite sum.

\end{document}