††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Critical behaviors of cascading dynamics on multiplex two-dimensional
lattices

Jeehye Choi  Research Institute for Nanoscale Science and Technology, Chungbuk
National University, Cheongju, Chungbuk 28644, Korea  Byungjoon Min  Research
Institute for Nanoscale Science and Technology, Chungbuk National University,
Cheongju, Chungbuk 28644, Korea  Department of Physics, Chungbuk National
University, Cheongju, Chungbuk 28644, Korea  K.-I. Goh  Department of Physics,
Korea University, Seoul 02841, Korea

###### Abstract

We study the critical phenomena of viable clusters in multiplex two-
dimensional lattices using numerical simulations. We identify viable sites on
multiplex lattices using two cascading algorithms: the cascade of activations
($\mathrm{CA}$) and deactivations ($\mathrm{CD}$). We found that the giant
viable clusters identified by $\mathrm{CA}$ and $\mathrm{CD}$ processes
exhibit different critical behaviors. Specifically, the critical phenomena of
$\mathrm{CA}$ processes are consistent with the ordinary bond percolation on a
single layer but $\mathrm{CD}$ processes exhibit the critical behaviors
consistent with mutual percolation on multiplex lattices. In addition, we
computed the susceptibility of cascading dynamics by using the concept of
ghost field. Our results suggest that the $\mathrm{CA}$ and $\mathrm{CD}$
processes generate viable clusters in different ways.

##  I Introduction

The elements of a complex system often function properly only if multiple
resources are provided through multiple layers of networks Kullman1979 ;
Rinaldi2001 ; Buldyrev2010 ; Min2014resource . For example, computer-
controlled systems function properly only when the Internet and power grids
are supplied interdependently Rinaldi2001 ; Buldyrev2010  In addition, for a
city to function properly, resources such as water, gas, and electricity must
be supplied by separate channels Kullman1979 ; Min2014resource ; Min2015 . For
this reason, the mutual connectivity of multiple layers of networks has
received much attention for several years Buldyrev2010 ; Son2012 ; Baxter2012
; Min2014resource ; Kivela2014 ; Lee2015 ; Lee2018 . A mutually-connected
component, which is a central concept in mutual percolation, is defined as a
set of nodes that every node pair in the component has at least one path,
composed of nodes within the same cluster, in each and every layer of networks
Buldyrev2010 ; Baxter2012 . Pioneering works on the mutual percolation have
studied the robustness of interdependent networks by using the notion of a
mutually connected giant component and showed that an abrupt transition can
appear between percolating and non-percolating phases Buldyrev2010 ;
Baxter2012 ; Son2012 ; Min2014Network . These studies imply that a small
perturbation can cause an abrupt collapse of the entire system, posing a
potentially catastrophe to complex systems Buldyrev2010 ; Son2012 ;
Min2014resource ; Lee2016 ; Zhou2014 ; Baxter2014 ; Hackett2016

Since many real-world systems are embedded in low dimensions, percolation
problems on multi-layered two-dimensional lattices have been also of interest
Son2011 ; Li2012 ; Bashan2013 . Unlike multiplex networks, the mutually
connected giant cluster on two-dimensional lattices emerges continuously as
the probability of bond occupation increases Son2011 ; Li2012 . However, the
studies on the critical behaviors of the mutual percolation on two-dimensional
lattices yielded conflicting results Son2011 ; Berezin2013 ; Son2013 ;
Grassberger2015 . Son et al. demonstrated that the percolation transition for
interdependent diluted lattices belongs to a different universality class with
a larger order parameter exponent than that of ordinary percolation Son2011 ;
Son2013 ; Grassberger2015 . However, another study reported that the mutual
percolation on two-dimensional lattices belongs to the same universality class
of ordinary percolation Berezin2013 . In addition to mutual percolation,
several percolation problems in multi-layer lattices have been also studied
such as history dependent percolation Jovanovic1994 ; Hu2020 ; Li2020 , and
percolation in porous media Cai1990 ; Hansen1993 .

From a perspective different from topological aspect, there have been studies
on the viability of nodes as a consequence of cascading dynamics following the
iterative activations or deactivations of nodes Min2014resource ; Min2015 ;
Baxter2014 . A model for the viability deals with systems that demand more
than one type of vital resource to be produced and distributed by source nodes
in multiplex networks Min2014resource . Viable nodes in this model are
identified by using two different dynamical processes: cascade of activations
($\mathrm{CA}$) and deactivations ($\mathrm{CD}$). The $\mathrm{CA}$ process
starts with all the sites in the deactivated state, and finds mutually-
connected clusters through the diffusion of activations from the source nodes;
in contrast, the $\mathrm{CD}$ process finds mutually-connected clusters by
iteratively removing any sites that are not mutually-connected. On multiplex
networks, these processes can produce two different final configurations of
viable nodes, corresponding to different stable solutions of a single mean-
field equation Min2014resource ; Baxter2014 . Viability is strongly related to
mutual percolation on multiplex networks because the dynamical consequence of
$\mathrm{CD}$ becomes identical to mutual percolation in the limit when the
fraction of initial source nodes goes to be zero. While cascading dynamics on
multiplex networks has been studied Min2014resource ; Min2015 , it has not
been explored on low-dimensional systems, especially the critical behaviors of
cascading dynamics.

In this work, we study the critical phenomena of cascading dynamics on
multiplex two-dimensional lattices using extensive Monte-Carlo simulations. We
found that two cascading processes, $\mathrm{CA}$ and $\mathrm{CD}$, can lead
to different critical behaviors in viable clusters statistics. While
$\mathrm{CD}$ exhibits the same set of critical exponents reported for mutual
percolation on interdependent networks Son2011 , $\mathrm{CA}$ exhibits the
critical behaviors akin to the ordinary percolation on two-dimensional
lattices. In addition, we computed the susceptibility of cascading dynamics by
using the ghost field concept Kasteleyn1969 ; Reynolds1977  and confirmed that
the hyperscaling relation is satisfied.

##  II Model

Figure 1:  (a) A multiplex lattice consists of double-layered two-dimensional
diluted lattices. (b,c) Examples of the final configurations of (b)
$\mathrm{CA}$ and (c) $\mathrm{CD}$ processes are depicted. Viable, non-
viable, and source sites are respectively denoted by black filled circles,
black open circles, and grey filled circles.

Let us imagine an example for cascading processes on a multiplex lattice. A
site needs two resources, say $\mathcal{A}$ and $\mathcal{B}$, to be viable,
e.g., a city requires water and electricity for functioning. Each resource is
delivered along the corresponding layer of the lattice. Initially, the
resources are generated from a few source sites in each layer. The important
point is that only a viable site can convey resources to the neighboring sites
along connected bonds. In this setting, we consider the following two
processes for finding viable sites from the given lattice structures and
initial source configurations. i) Assume that all sites except the source
sites are initially in the non-viable state. Sites that satisfy the viable
condition change into the viable state. These viable sites are able to convey
resources to their neighbors. There will be sites that change to the viable
state in turn and we call this process the cascade of activation
($\mathrm{CA}$). ii) On the other hand, one may start at the beginning when
all sites are putatively in a viable state. In this case, sites that do not
satisfy the viable condition must change from the viable to non-viable state
iteratively until no further deactivation is necessary. This process is called
the cascade of deactivation ($\mathrm{CD}$).

The specific simulation rules are as follows. There are two layers of two-
dimensional lattices, say $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ with
$N=L\times L$, where $L$ denotes the size of the lattice. Each bond in the
lattice is occupied with probability $p$, and empty with probability $1-p$. On
top of multiplex diluted lattices, resources $\mathcal{A}$ and $\mathcal{B}$
are distributed at their source sites. In this study, we assume that the
sources of resources $\mathcal{A}$ and $\mathcal{B}$ are located at every site
on the boundary of lattices $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$,
respectively [Fig. 1(a)]. We use cylindrical boundary conditions. Then, each
site on a multiplex lattice is viable only if both resources $\mathcal{A}$ and
$\mathcal{B}$ are supplied. Moreover, only viable sites can deliver resources
to their neighbors connected by occupied bonds. Otherwise a site becomes non-
viable and cannot convey resources further. Note that resources $\mathcal{A}$
and $\mathcal{B}$ are respectively supplied only through the chain of viable
sites on $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$. Using the two processes,
$\mathrm{CA}$ and $\mathrm{CD}$, we find viable sites in the steady state. We
then identify the viable clusters as a set of viable sites for every pair of
sites in the cluster has at least one path composed of viable sites within
each and every layer of the lattice. We also define the fraction of the giant
viable cluster as $V_{\mathrm{CA}}$ for $\mathrm{CA}$ as shown in Fig. 1(b),
and $V_{\mathrm{CD}}$ for $\mathrm{CD}$ as shown in Fig. 1(c).

##  III Results

We conducted Monte Carlo simulations for $\mathrm{CA}$ and $\mathrm{CD}$
processes on multiplex lattices with various sizes $L$. We then identify the
configuration of the viable clusters and analyze their critical behaviors at
the steady state. The fractions $V$ of the giant viable cluster for the two
processes are in general different from each other as shown in Fig. 2(a), as
in the case of multiplex networks Min2014resource . The giant viable cluster
for $\mathrm{CD}$ process appears at a smaller $p$ compared with that for
$\mathrm{CA}$. In addition, the fraction $V_{\mathrm{CD}}$ of the giant viable
cluster for $\mathrm{CD}$ is generally larger than $V_{\mathrm{CA}}$.

The critical behaviors of percolation problems are typically characterized by
the divergence of the average cluster size $\chi$ and correlation length $\xi$
at the critical point $p_{c}$ in the thermodynamic limit $L\rightarrow\infty$.
Near the critical point, the order parameter $V$, susceptibility $\chi$, and
correlation length $\xi$ follow power laws with the critical exponents
$\beta$, $\gamma$, and $\nu$ as follows:

$\displaystyle V(p)\sim(p-p_{c})^{\beta},$ (1)
$\displaystyle\chi(p)\sim|p-p_{c}|^{-\gamma},$ (2)
$\displaystyle\xi(p)\sim|p-p_{c}|^{-\nu}.$ (3)

In order to obtain the critical exponents, we use the conventional scaling
ansatz for finite-size systems Stauffer2018 ,

$\displaystyle V(p,L)=L^{-\beta/\nu}f_{V}[(p-p_{c})L^{1/\nu}],$ (4)
$\displaystyle\chi(p,L)=L^{\gamma/\nu}f_{\chi}[(p-p_{c}(L))L^{1/\nu}].$ (5)

where $p_{c}(L)$ is the location of the peak of susceptibility.

Figure 2:  (a) Numerical results of the fraction $V$ for the $\mathrm{CA}$
(red) and $\mathrm{CD}$ (blue) processes as a function of the bond occupation
probability $p$ on multiplex lattices with $L=800$, averaged over $100$
different realizations are shown. (b) The relation $|p_{c}-p_{c}(L)|\sim
L^{-1/\nu}$ for $\mathrm{CA}$ (red) and $\mathrm{CD}$ (blue) is shown. The
points represent numerical simulation results averaged over $400$ samples.

While computing the fraction $V$ of the giant viable cluster is
straightforward, how to measure susceptibility $\chi$ is far from trivial. For
the ordinary percolation, the susceptibility $\chi$ is directly measured by
the average size of the finite clusters. In our model, however, the quantity
corresponding to the susceptibility is not readily provided by the size of the
average finite viable clusters. To obtain the susceptibility in $\mathrm{CA}$
and $\mathrm{CD}$, we use the concept of “ghost field” analogous to magnetic
susceptibility Kasteleyn1969 ; Reynolds1977 ; Fortuin1972 . We add one
additional ghost site that is connected to each site in the system with a
probability $H$. By applying a small amount of $H$ corresponding to the
external magnetic field in spin systems, we measure the change in order
parameter $V$ such that $\chi=\partial V/\partial H$. In practice, we choose
one site that does not belong to the giant viable cluster. Next, we assume
that this site has become part of the giant viable cluster through the
connection to the ghost site. We then identify a set of sites that newly
entered the giant viable cluster because of the presence of the ghost site. We
call the set of sites as “susceptible cluster” and measure the size of the
susceptible cluster as susceptibility. In the numerical simulations, we
attempted to select every site that does not belong to the giant viable
cluster and measured the distribution $\phi(s)$ that a site belongs to a
susceptible cluster with size $s$. The susceptibility $\chi$ is obtained as

$\chi(p,L)=\frac{\sum_{s}s\phi(s,p)}{\sum_{s}\phi(s,p)},$ (6)

where the sum is over all cluster sizes excluding the giant viable cluster.

We estimated the viability threshold and critical exponent $\nu$ for both
$\mathrm{CA}$ and $\mathrm{CD}$ using a finite size scaling ansatz. We first
obtained the threshold as $p_{c}^{\mathrm{CA}}\simeq 0.6743$ and
$p_{c}^{\mathrm{CD}}\simeq 0.5761$. As in the case of $\mathrm{CA}$ and
$\mathrm{CD}$ in multiplex networks Min2014resource , the giant viable
clusters in $\mathrm{CA}$ processes emerges at a larger $p$ compared to
$\mathrm{CD}$. Next, we obtained the exponent of the correlation length $\nu$
from the relation $|p_{c}(L)-p_{c}|\sim L^{-1/\nu}$ as shown in Fig. 2(b).
Here, $p_{c}(L)$ is a threshold with size $L$ and is estimated by the value of
$p$ at the maximum value of susceptibility $\chi$. The fitted values of $\nu$
for $\mathrm{CA}$ and $\mathrm{CD}$ are respectively $1.342(3)$ and
$1.197(7)$.

Figure 3:  (a) Log-log plot of $V_{\mathrm{CA}}L^{\beta/\nu}$ at
$p_{c}^{\mathrm{CA}}$ with respect to $L$ for various $p$. (b) Data collapse
of the scaled order parameter $V_{\mathrm{CA}}L^{\beta/\nu}$ vs.
$(p-p_{c})L^{1/\nu}$ with $\beta=0.137(1)$ and $\nu=1.342(3)$. (c) Log-log
plot of $\chi_{\mathrm{CA}}$ at $p_{c}(L)$ with respect to $L$, showing the
relation $\chi\sim L^{\gamma/\nu}$. (d) Data collapse of the scaled
susceptibility $\chi_{\mathrm{CA}}L^{-\gamma/\nu}$ vs. $(p-p_{c})L^{1/\nu}$
with $\gamma=2.41(1)$ and $\nu=1.342(3)$. Numerical results of $\mathrm{CA}$
were obtained for various lattice sizes $L$ from 50 to 800, averaged over 400
realizations.  Figure 4:  (a) Log-log plot of $V_{\mathrm{CD}}L^{\beta/\nu}$
at $p_{c}^{\mathrm{CD}}$ with respect to $L$ for various $p$. (b) Data
collapse of the scaled order parameter $V_{\mathrm{CD}}L^{\beta/\nu}$ vs.
$(p-p_{c})L^{1/\nu}$ with $\beta=0.163(2)$ and $\nu=1.197(7)$. (c) Log-log
plot of $\chi_{\mathrm{CD}}$ at $p_{c}(L)$ with respect to $L$, showing the
relation $\chi\sim L^{\gamma/\nu}$. (d) Data collapse of the scaled
susceptibility $\chi_{\mathrm{CD}}L^{-\gamma/\nu}$ vs. $(p-p_{c})L^{1/\nu}$
with $\gamma=2.07(1)$ and $\nu=1.197(7)$. Numerical results of $\mathrm{CD}$
were obtained for various lattice sizes $L$ from 50 to 800, averaged over 400
realizations.

We then obtained $\beta/\nu$ and $\gamma/\nu$ using the scalings respectively
$V\sim L^{-\beta/\nu}$ and $\chi\sim L^{\gamma/\nu}$. The results of critical
behaviors for $\mathrm{CA}$ processes are shown in Fig. 3. Figure 3(a) shows
the scaling of $V\sim L^{-\beta/\nu}$ with various values of $p$. The scaling
shows the estimates of the exponents $\beta=0.137(1)$ and $\nu=1.342(3)$. Data
collapse curves of $VL^{\beta/\nu}$ with respect to $(p-p_{c})L^{1/\nu}$ in
Fig. 3(b) confirm the estimated values of the critical exponents. Figure 3(c)
shows the relationship $\chi\sim L^{\gamma/\nu}$ with $\gamma/\nu=1.794(6)$
which is confirmed by the finite-size-scaled data collapse of $\chi
L^{-\gamma/\nu}$ with respect to $(p-p_{c})L^{1/\nu}$, as shown in Fig. 3(d).
We found that the critical exponents for $\mathrm{CA}$ processes are
consistent within error bars for ordinary bond percolation on two-dimensional
lattices.

For $\mathrm{CD}$ processes, we also obtain $\beta$ and $\gamma$ using scaling
$V\sim L^{-\beta/\nu}$ and $\chi\sim L^{\gamma/\nu}$ as shown in Fig. 4. The
estimated values of the critical exponents for the $\mathrm{CD}$ processes are
$\beta=0.163(2)$ and $\nu=1.197(7)$ [Fig. 4(a)]. In addition, figure 4(c)
exhibits the relationship $\chi\sim L^{\gamma/\nu}$ with $\gamma/\nu=1.73(1)$.
We also confirm the critical exponents by using the data collapse of
$VL^{\beta/\nu}$ with respect to $(p-p_{c})L^{1/\nu}$ [Fig. 4(b)] and $\chi
L^{-\gamma/\nu}$ with respect to $(p-p_{c})L^{1/\nu}$ [Fig. 4(d)]. The values
of the exponents $\beta$, $\gamma$, and $\nu$ are distinct from those for
$\mathrm{CA}$. In addition, the critical behaviors of $\mathrm{CD}$ are
consistent within error bars of the mutual percolation on two-dimensional
lattices reported in Son2011 ; Grassberger2015 .

Figure 5: Probability distributions $\phi(s)$ that a node belongs to a
susceptible cluster with size $s$ for (a) $\mathrm{CA}$ and (b) $\mathrm{CD}$.
Data collapses of $\phi(s)s^{\tau-1}$ vs. $s/L^{d_{f}}$ for (c) $\mathrm{CA}$
with $\tau=2.06(3)$ and $d_{f}=1.8976$ and for (d) $\mathrm{CD}$ with
$\tau=2.07(1)$ and $d_{f}=1.864$. The fractal dimension $d_{f}$ is estimated
by $d_{f}=d-\beta/\nu$. Monte-carlo simulations were performed at the critical
point $p=0.6743$ for $\mathrm{CA}$ and $p=0.5761$ for $\mathrm{CD}$. Numerical
results were obtained for various lattice sizes $L$ from 50 to 800, averaged
over 400 realizations.

We also obtain the size distribution $\phi(s)$ that a site belongs to
susceptible clusters with size $s$. At $p_{c}$, the size distribution follows
the power-law form $\phi(s)\sim s^{1-\tau}$ in the thermodynamic limit. We
apply a scaling ansatz to the distribution

$\displaystyle\phi(s)=s^{1-\tau}g(s/L^{d_{f}})$ (7)

where $g$ is a scaling function and $d_{f}$ is the fractal dimension given by
$d_{f}=d-\beta/\nu$ where $d$ is the spatial dimension of lattices. Figures
5(a,b) show the power-law decay of $\phi(s)$ at $p_{c}$ for (a) $\mathrm{CA}$
and (b) $\mathrm{CD}$, respectively. We found the estimate of critical
exponents as $\tau=2.06(3)$ for $\mathrm{CA}$ and $\tau=2.07(1)$ for
$\mathrm{CD}$. We confirm the scaling ansatz using data collapses, as shown in
Figs. 5(c,d).

The obtained values of the threshold and critical exponents are presented in
Table 1. We observed that the critical exponents for $\mathrm{CA}$ and
$\mathrm{CD}$ are distinct. Moreover, all critical exponents of $\mathrm{CA}$
are consistent within error bars with ordinary bond percolation on two-
dimensional lattices. However, the critical exponents of $\mathrm{CD}$ are
consistent with those of the mutual percolation on two-dimensional lattices,
which shows a discrepancy with $\mathrm{CA}$. In addition, the sets of
exponents for both $\mathrm{CA}$ and $\mathrm{CD}$ satisfy the hyperscaling
relations $d\nu=2\beta+\gamma$ and $\tau=d/d_{f}+1$.

Table 1: Critical points and exponents | $\mathrm{CA}$ | Ordinary | $\mathrm{CD}$ | Mutual
---|---|---|---|---
Percolation | Percolation Grassberger2015
$p_{c}$ | 0.6743(1) | 1/2 | 0.5761(1) | 0.576132(5)
$\nu$ | 1.342(3) | 4/3 | 1.197(7) | 1.200(7)
$\beta$ | 0.137(1) | 5/36 | 0.163(2) | 0.163(2)
$\gamma$ | 2.41(1) | 43/18 | 2.07(2) | $\cdot$
$\tau$ | 2.06(3) | 187/91 | 2.07(1) | $\cdot$
$1/\nu$ | 0.745(2) | 3/4 | 0.836(5) | 0.833
$\beta/\nu$ | 0.1024(5) | 5/48 | 0.136(1) | 0.1358
$\gamma/\nu$ | 1.794(6) | 43/24 | 1.73(1) | 1.728
Figure 6: Joint probability distributions $J(s,s_{\textrm{dbc}})$ where $s$ is
the size of susceptible clusters and $s_{\textrm{dbc}}$ is that of double-bond
clusters are shown for (a) $\mathrm{CA}$ and (b) $\mathrm{CD}$. Joint
probability distributions $J(s,s_{\textrm{mcc}})$ where $s_{\textrm{mcc}}$ is
the size of mutually connected clusters are shown for (c) $\mathrm{CA}$ and
(d) $\mathrm{CD}$. Numerical results are calculated at the critical points
$p=0.6743$ for $\mathrm{CA}$ and $p=0.5761$ for $\mathrm{CD}$ with $L=100$,
averaged over $5000$ realizations.

To understand the differences between the underlying mechanisms forming
susceptible clusters, we examine the joint probability distribution between
the sizes of susceptible clusters $s$ and that of double-bond clusters
$s_{\textrm{dbc}}$, and that of mutually connected clusters
$s_{\textrm{mcc}}$. A double bond refers to a connection in which a bond
between two adjacent sites exists in both layers, and the double-bond cluster
stands for a set of sites that are connected solely via double bonds. The
mutually connected component is defined as a set of sites that every pair of
sites in the cluster has at least one path within each and every layer of the
lattice.

Figures 6(a,b) show the joint probability distribution $J(s,s_{\textrm{dbc}})$
for (a) $\mathrm{CA}$ and (b) $\mathrm{CD}$ to examine the effect of double
bonds in forming susceptible clusters. The joint probability distributions are
computed at the critical point, $p=0.6743$ for $\mathrm{CA}$ and $p=0.5761$
for $\mathrm{CD}$. We found that the size of the double-bond clusters are
strongly correlated to that of the susceptible clusters for $\mathrm{CA}$.
This means that the double bonds play a dominant role in the $\mathrm{CA}$
process. For $\mathrm{CD}$ process, however, small susceptible clusters are
formed by double bonds, while large susceptible clusters are not. We found
that large-sized susceptible clusters in $\mathrm{CD}$ are constructed by
mutually connected clusters as shown in Fig. 6(d). On the other hand, for the
$\mathrm{CA}$ process there is already a mutually connected giant component
near the critical point as shown in Fig. 6(c). Therefore, mutually connected
clusters cannot contribute crucially the emergence of the large-sized
susceptible clusters in $\mathrm{CA}$ process. In conclusion, we found that
large-sized susceptible clusters are mainly constructed by double bonds in the
$\mathrm{CA}$ process whereas they are formed by mutual connections in the
$\mathrm{CD}$ process.

##  IV Discussion

In this work, we study the critical phenomena of cascading dynamics on
multiplex lattices considering two different processes to identify viable
clusters. We found that viable clusters identified by the $\mathrm{CA}$ and
$\mathrm{CD}$ processes exhibits different critical behaviors. The
$\mathrm{CA}$ process exhibits critical phenomena consistent to ordinary bond
percolation but the $\mathrm{CD}$ process exhibits those to mutual percolation
on multiplex lattices. In addition, we explicitly calculate susceptibility by
using the ghost field. Our study shows that viability that requires multiple
connectivity can be maintained by diverse mechanisms in a multiplex low-
dimensional system. The cascading dynamics that we consider here is an example
of cooperative couplings between multiple layers. Similar cooperative
interactions can be realized in real-world systems in a variety of ways such
as cooperative infections Chen2013 ; Min2020 , social contagions Lee2014 ;
Min2018 , and interdependencies in multi-layered systems Buldyrev2010 ;
Pocock2012 . In this regard, our research shows a glimpse of the complexity
generated by cooperative interactions in multiplex low-dimensional systems.
Further studies may be required to examine the cascading dynamics on an
arbitrary number of layers and partially coupled layers.

###### Acknowledgements.

This work was supported by the National Research Foundation of Korea (NRF)
grants funded by the Korea government (MSIT) (No. 2020R1I1A3068803 (BM) and
No. 2020R1A2C2003669 (K-IG)).

## References

  * (1) D. E. Kullman, The utilities problem, Mathematics Magazine 52 (5) (1979) 299–302.
  * (2) S. Rinaldi, J. Peerenboom, T. Kelly, Identifying, understanding, and analyzing critical infrastructure interdependencies, IEEE Control Systems Magazine 21 (6) (2001) 11–25.  doi:10.1109/37.969131.
  * (3) S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, S. Havlin, Catastrophic cascade of failures in interdependent networks, Nature 464 (7291) (2010) 1025–1028.  doi:10.1038/nature08932.
  * (4) B. Min, K.-I. Goh, Multiple resource demands and viability in multiplex networks, Phys. Rev. E 89 (4) (2014) 040802(R).
  * (5) B. Min, S. Lee, K.-M. Lee, K.-I. Goh, Link overlap, viability, and mutual percolation in multiplex networks, Chaos, Solitons, and Fractals 72 (2015) 49\.
  * (6) S.-W. Son, G. Bizhani, C. Christensen, P. Grassberger, M. Paczuski, Percolation theory on interdependent networks based on epidemic spreading, Europhysics Letters 97 (1) (2012) 16006.  doi:10.1209/0295-5075/97/16006.
  * (7) G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, J. F. F. Mendes, Avalanche Collapse of Interdependent Networks, Phys. Rev. Lett. 109 (24) (2012) 248701\.  doi:10.1103/PhysRevLett.109.248701.
  * (8) M. Kivela, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, M. A. Porter, Multilayer networks, Journal of Complex Networks 2 (3) (2014) 203–271.  doi:10.1093/comnet/cnu016.
  * (9) K.-M. Lee, B. Min, K.-I. Goh, Towards real-world complexity: an introduction to multiplex networks, The European Physical Journal B 88 (48) (2015) 48.
  * (10) D. Lee, B. Kahng, Y. S. Cho, K. I. Goh, D. S. Lee, Recent advances of percolation theory in complex networks, Journal of the Korean Physical Society 73 (2) (2018) 152–164.
  * (11) B. Min, S. D. Yi, K.-M. Lee, K.-I. Goh, Network robustness of multiplex networks with interlayer degree correlations, Phys. Rev. E 89 (4) (2014) 042811\.  doi:10.1103/PhysRevE.89.042811.
  * (12) D. Lee, S. Choi, M. Stippinger, J. Kertész, B. Kahng, Hybrid phase transition into an absorbing state: Percolation and avalanches, Phys. Rev. E 93 (4) (2016) 042109.  doi:10.1103/PhysRevE.93.042109.
  * (13) D. Zhou, A. Bashan, R. Cohen, Y. Berezin, N. Shnerb, S. Havlin, Simultaneous first- and second-order percolation transitions in interdependent networks, Phys. Rev. E 90 (1) (2014) 012803.  doi:10.1103/PhysRevE.90.012803.
  * (14) G. J. Baxter, S. N. Dorogovtsev, J. F. F. Mendes, D. Cellai, Weak percolation on multiplex networks, Phys. Rev. E 89 (4) (2014) 042801.
  * (15) A. Hackett, D. Cellai, S. Gómez, A. Arenas, J. P. Gleeson, Bond percolation on multiplex networks, Phys. Rev. X 6 (2) (2016) 021002.
  * (16) S.-W. Son, P. Grassberger, M. Paczuski, Percolation Transitions Are Not Always Sharpened by Making Networks Interdependent, Phys. Rev. Lett. 107 (19) (2011) 195702.  doi:10.1103/PhysRevLett.107.195702.
  * (17) W. Li, A. Bashan, S. V. Buldyrev, H. E. Stanley, S. Havlin, Cascading failures in interdependent lattice networks: The critical role of the length of dependency links, Phys. Rev. Lett. 108 (22) (2012) 228702.
  * (18) A. Bashan, Y. Berezin, S. V. Buldyrev, S. Havlin, The extreme vulnerability of interdependent spatially embedded networks, Nat. Phys. 9 (10) (2013) 667–672.
  * (19) Y. Berezin, A. Bashan, S. Havlin, Comment on “Percolation Transitions Are Not Always Sharpened by Making Networks Interdependent”, Phys. Rev. Lett. 111 (18) (2013) 189601.  doi:10.1103/PhysRevLett.111.189601.
  * (20) S.-W. Son, P. Grassberger, M. Paczuski, Son, Grassberger, and Paczuski Reply:, Phys. Rev. Lett. 111 (18) (2013) 189602.  doi:10.1103/PhysRevLett.111.189602.
  * (21) P. Grassberger, Percolation transitions in the survival of interdependent agents on multiplex networks, catastrophic cascades, and solid-on-solid surface growth, Phys. Rev. E 91 (6) (2015) 062806.  doi:10.1103/PhysRevE.91.062806.
  * (22) B. Jovanović, S. V. Buldyrev, S. Havlin, H. E. Stanley, Punctuated equilibrium and “history-dependent” percolation, Phys. Rev. E 50 (1994) R2403–R2406.
  * (23) M. Hu, Y. Sun, D. Wang, J.-P. Lv, Y. Deng, History-dependent percolation in two dimensions, Phys. Rev. E 102 (2020) 052121.  doi:10.1103/PhysRevE.102.052121.
  * (24) M. Li, L. Lü, Y. Deng, M.-B. Hu, H. Wang, M. Medo, H. E. Stanley, History-dependent percolation on multiplex networks, National Science Review 7 (8) (2020) 1296–1305.
  * (25) Z.-X. Cai, S. D. Mahanti, S. A. Solin, T. J. Pinnavaia, Dual percolation threshold in two-dimensional microporous media, Phys. Rev. B 42 (1990) 6636–6641.
  * (26) A. Hansen, E. L. Hinrichsen, D. Stauffer, Percolation in layered media — A conductivity approach, Transport in Porous Media 11 (1) (1993) 45–52.
  * (27) P. Kasteleyn, C. Fortuin, Phase transitions in lattice systems with random local properties, Physical Society of Japan Journal Supplement 26 (1969) 11.
  * (28) P. J. Reynolds, H. E. Stanley, W. Klein, Ghost fields, pair connectedness, and scaling: exact results in one-dimensional percolation, Journal of Physics A: Mathematical and General 10 (11) (1977) L203.
  * (29) D. Stauffer, A. Aharony, Introduction to percolation theory, Taylor & Francis, 2018\.
  * (30) C. M. Fortuin, On the random-cluster model. III. The simple random-cluster model, Physica 59 (4) (1972) 545–570.  doi:10.1016/0031-8914(72)90087-0.
  * (31) L. Chen, F. Ghanbarnejad, W. Cai, P. Grassberger, Outbreaks of coinfections: The critical role of cooperativity, Europhysics Letters 104 (5) (2013) 50001.  doi:10.1209/0295-5075/104/50001.
  * (32) B. Min, C. Castellano, Message-passing theory for cooperative epidemics, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2) (2020) 023131.  doi:10.1063/1.5140813.
  * (33) K.-M. Lee, C. D. Brummitt, K.-I. Goh, Threshold cascades with response heterogeneity in multiplex networks, Phys. Rev. E 90 (2014) 062816.  doi:10.1103/PhysRevE.90.062816.
  * (34) B. Min, M. San Miguel, Competition and dual users in complex contagion processes, Scientific Reports 8 (1) (2018) 14580.  doi:10.1038/s41598-018-32643-4.
  * (35) M. J. O. Pocock, D. M. Evans, J. Memmott, The robustness and restoration of a network of ecological networks, Science 335 (6071) (2012) 973–977.  doi:10.1126/science.1214915.


# Spherical curves whose curvature
depends on distance to a great circle
Dedicated to Professor Óscar Garay, in memoriam

Ildefonso Castro  Departamento de Matemáticas
Universidad de Jaén
23071 Jaén, Spain and IMAG, Instituto de Matemáticas de la Universidad de
Granada<EMAIL_ADDRESS>,  Ildefonso Castro-Infantes  Departamento de
Matemáticas
Universidad de Alicante
Alicante, Spain<EMAIL_ADDRESS>and  Jesús Castro-Infantes
Departamento de Geometría y Topología
Universidad de Granada
18071 Granada, Spain<EMAIL_ADDRESS>
###### Abstract.

Motivated by a problem posed by David A. Singer in 1999 and by the elastic
spherical curves, we study the spherical curves whose curvature is expressed
in terms of the distance to a great circle (or from a point). By introducing
the notion of spherical angular momentum, we provide new characterizations of
some well known curves, like the mentioned elastic curves, spherical
catenaries, loxodromic-type spherical curves, the Viviani’s curve, and the
spherical Archimedean spirals curves. Furthermore, we show that they may be
obtained as critical points of some energy curvature functionals. We also find
out several new families of spherical curves whose intrinsic equations are
expressed in terms of elementary functions or Jacobi elliptic functions, and
we are able to get arc length parametrizations of them.

###### 2010 Mathematics Subject Classification:

Primary 53A04; Secondary 74H05

##  1\. Introduction

The plane curves are uniquely determined up to rigid motion by its intrinsic
equation giving its curvature $\kappa$ as a function of its arc-length.
However, such curves are impossible to find explicitly in practice in most
cases, due to the difficulty in solving the three quadratures appearing in the
integration process. In [S99], David A. Singer considered a different sort of
problem: Can a plane curve be determined if its curvature is given in terms of
its position?

Probably, the most interesting solved problem in this setting corresponds to
the Euler elastic curves, whose curvature is proportional to one of the
coordinate functions, e.g. $\kappa(x,y)=c\,y$. Motivated by the above question
and by the classical elasticae, the authors studied in [CCI16] the plane
curves whose curvature depends on the distance to a line (say the $x$-axis,
and so $\kappa=\kappa(y)$) and in [CCIs17] the plane curves whose curvature
depends on the distance from a point (say the origin, and so
$\kappa=\kappa(r)$, $r=\sqrt{x^{2}+y^{2}}$) requiring in both cases the
computation of three quadratures too. But the simple case $\kappa(r)=r$, where
elliptic integrals appear, illustrated that the fact that the corresponding
differential equation is integrable by quadratures does not mean that it is
easy to perform the integrations. In [S99], only the very pleasant special
case of the classical Bernoulli lemniscate, $r^{2}=3\cos 2\theta$ in polar
coordinates, was solved explicitly, where the corresponding elliptic integral
becomes elementary.

In this paper, we pay our attention to the Singer’s problem version for curves
lying in a sphere:

> Can a spherical curve be determined when its curvature is given in terms of
> its position?

The geodesic curvature $\kappa$ of a spherical curve $\xi$ given as a function
of its arc length $s$ determines the curve (up to isometries of the sphere) by
integration of its Frenet equations. However, it is expectable that if the
curvature $\kappa$ of $\xi=(x,y,z)$ is given by a function of its position,
i.e. $\kappa=\kappa(x,y,z)$, the situation becomes quite complicated since the
general form of this problem is equivalent to solving the non linear
differential equation

$\left|\begin{array}[]{ccc}x(s)&y(s)&z(s)\\\
\dot{x}(s)&\dot{y}(s)&\dot{z}(s)\\\
\ddot{x}(s)&\ddot{y}(s)&\ddot{z}(s)\end{array}\right|=\kappa(x(s),y(s),z(s))$

with the constraints

$x(s)^{2}+y(s)^{2}+z(s)^{2}=1\quad\text{and}\quad\dot{x}(s)^{2}+\dot{y}(s)^{2}+\dot{z}(s)^{2}=1.$

The purpose of this article is the study of the aforementioned cases of the
classical Singer’s problem in the setting of spherical curves, considering
geodesics in the role of lines. Concretely, we consider a curve $\xi=(x,y,z)$
lying in the unit sphere $\mathbb{S}^{2}$ centred at the origin and write
$z=\sin\varphi$ ($\varphi$ being the latitude of $\xi$); we aim to control
those curves $\xi$ whose geodesic curvature $\kappa$ satisfies the condition
$\kappa\\!=\\!\kappa(\varphi)\\!\Leftrightarrow\\!\kappa\\!=\\!\kappa(z)$. We
point out that this condition includes both types of problems involving
curvature and distance, since $\varphi$ is the distance to the equator (the
great circle $\varphi\\!=\\!0\\!\Leftrightarrow\\!z\\!=\\!0$) and the
colatitude $\pi/2-\varphi$ is the distance to the North pole (the point
(0,0,1)).

As in the Euclidean case, it will be necessary again the computation of three
quadratures when $\kappa(z)$ is a continuous function (see Theorem 2.1) and
the key tool will be the notion of spherical angular momentum, which
completely determines a spherical curve (up to a family of distinguished
isometries) in relation with its relative position with respect to a fixed
geodesic.

In this way, we find out several interesting new families of spherical curves
whose intrinsic equations can be expressed in terms of elementary or Jacobi
elliptic functions. We also provide new characterizations of some well known
curves, like elastic-type spherical curves, spherical catenaries, loxodromic-
type spherical curves, the Viviani’s curve, and the spherical Archimedean
spirals curves. In addition, we show that they may be obtained as critical
points of some energy curvature functionals.

In a forthcoming paper [CCIs21] we have studied helicoidal minimal surfaces in
$\mathbb{S}^{3}$ by considering surfaces that are invariant under a helicoidal
motion in the $3$-sphere, that is, the composition of two independent
rotations in $\mathbb{S}^{3}$.

##  2\. Spherical curves such that $\kappa=\kappa(z)$ and the spherical
angular momentum

We introduce a smooth function associated to any spherical curve, which
completely determines it (up to a family of distinguished isometries) in
relation with its relative position with respect to a fixed geodesic.

Indeed, let $\xi=\xi(s)\colon I\subseteq\mathbb{R}\rightarrow\mathbb{S}^{2}$
be an immersed curve parametrized by the arc length, i.e.
$|\xi(s)|=|\dot{\xi}(s)|=1$, for any $s\in I$, where $I$ is some interval in
$\mathbb{R}$. Along the paper, $\dot{\,}$ will denote derivative with respect
to $s$ and $\langle\cdot,\cdot\rangle$ and $\times$ the Euclidean inner
product and the cross product in $\mathbb{R}^{3}$ respectively. Let
$T=\dot{\xi}$ be the unit tangent vector and $N=\xi\times\dot{\xi}$ the unit
normal vector of $\xi$. If $\nabla$ is the connection in $\mathbb{S}^{2}$, the
oriented geodesic curvature $\kappa$ of $\xi$ is given by the Frenet equation
$\nabla_{T}T=\kappa N$, which implies

(2.1) $\ddot{\xi}=-\xi+\kappa N,\quad\dot{N}=-\kappa\,\dot{\xi}$

and so $\kappa=\det(\xi,\dot{\xi},\ddot{\xi})$.

We are interested in the geometric condition that the curvature of $\xi$
depends on the distance to a geodesic of $\mathbb{S}^{2}$. If ${\bf
e}\in\mathbb{R}^{3}$ is a unit length vector, then $\langle\xi,{\bf e}\rangle$
is the signed distance to the orthogonal plane to $\bf e$ passing through the
origin. Without restriction, we consider ${\bf e}:=(0,0,1)$ and write
$\xi=(x,y,z)$ with $x^{2}+y^{2}+z^{2}=1$. So we can pay our attention to study
the condition $\kappa=\kappa(z)$ since $z=\langle\xi,{\bf e}\rangle$
represents the signed distance to the great circle
$\mathbb{S}^{2}\cap\\{z=0\\}$. Concretely, we use geographical coordinates in
$\mathbb{S}^{2}$ and write

$\xi=(\cos\varphi\cos\lambda,\cos\varphi\sin\lambda,\sin\varphi),\
-\pi/2\leq\varphi\leq\pi/2,\,-\pi<\lambda\leq\pi.$

Then it is interesting to notice that the latitude $\varphi$ is the signed
distance to the equator $\mathbb{S}^{2}\cap\\{z=0\\}\equiv\varphi=0$ and, in
addition, the colatitude $\pi/2-\varphi$ gives the distance to the North pole
$(0,0,1)$.

At a given point $\xi(s)$ on the curve, we introduce the spherical angular
momentum (with respect to the $z$-axis) $\mathcal{K}(s)$ as the (signed)
volume of the parallelepiped spanned by the position $\xi(s)$, the unit
tangent $T(s)$ and the vector ${\bf e}:=(0,0,1)$. Concretely, we define

(2.2) $\mathcal{K}(s):=-\det(\xi(s),T(s),{\bf e})=-\langle N(s),{\bf
e}\rangle=\dot{x}(s)y(s)-x(s)\dot{y}(s).$

In physical terms, as a consequence of Noether’s Theorem (cf. [A78]),
$\mathcal{K}$ may be described as the angular momentum of a particle of unit
mass with unit speed and spherical trajectory $\xi(s)$. We point out that
$\mathcal{K}$ assumes values in $[-1,1]$ and it is well defined, up to the
sign, depending on the orientation of the normal to $\xi$. In geographical
coordinates, $\mathcal{K}$ is given by

(2.3) $\mathcal{K}=-\dot{\lambda}\cos^{2}\varphi.$

The unit-speed condition on $\xi$ implies that
$\dot{\varphi}^{2}+\dot{\lambda}^{2}\cos^{2}\varphi=1$ and, assuming $\varphi$
is non constant and using (2.3), we deduce that

(2.4)
$ds=\frac{d\varphi}{\sqrt{1-\dot{\lambda}^{2}\cos^{2}\varphi}}=\frac{\cos\varphi\,d\varphi}{\sqrt{\cos^{2}\varphi-\mathcal{K}^{2}}}=\frac{dz}{\sqrt{1-z^{2}-\mathcal{K}^{2}}}$

and

(2.5)
$d\lambda=-\frac{\mathcal{K}ds}{\cos^{2}\varphi}=\frac{\mathcal{K}ds}{z^{2}-1}.$

Hence, given $\mathcal{K}=\mathcal{K}(z)$ as an explicit function, looking at
(2.4) and (2.5), one may attempt to compute $z(s)$ (and so $\varphi(s)$) and
$\lambda(s)$ in three steps: integrate (2.4) to get $s=s(z)$, invert to get
$z=z(s)$ and integrate (2.5) to get $\lambda=\lambda(s)$. We remark that the
integration constants appearing in (2.4) and (2.5) simply mean a translation
of the arc parameter and a rotation around the $z$-axis respectively.

In addition, using (2.1) and (2.2), we have that
$\dot{\mathcal{K}}=-\langle\dot{N},{\bf e}\rangle=\kappa\langle\dot{\xi},{\bf
e}\rangle=\kappa\dot{z}$ and, if we take into account the assumption
$\kappa=\kappa(z)$ (being $z$ non constant), we finally arrive at

(2.6) $d\mathcal{K}=\kappa(z)dz,$

that is, $\mathcal{K}(z)$ can be interpreted as an anti-derivative of
$\kappa(z)$.

As a summary, we can determine by quadratures in a constructive explicit way
the spherical curves such that $\kappa=\kappa(z)$, in the spirit of [S99,
Theorem 3.1].

######  Theorem 2.1.

Let $\kappa=\kappa(z)$ be a continuous function. Then the problem of
determining locally a spherical curve whose curvature is $\kappa(z)$ —$z$
representing the (non constant) signed distance to the great circle
$z\\!=\\!0$— with spherical angular momentum $\mathcal{K}(z)$ satisfying
(2.6), is solvable by quadratures considering the unit speed curve
$\\!\xi(s)\\!=\\!(x(s),y(s),z(s))$, with
$x(s)=\cos\varphi(s)\cos\lambda(s),\,y(s)=\cos\varphi(s)\sin\lambda(s),\,z(s)=\sin\varphi(s)$,
where $\varphi(s)$ and $\lambda(s)$ are obtained through (2.4) and (2.5) after
inverting $s=s(z)$. Such a curve is uniquely determined by $\mathcal{K}(z)$ up
to a rotation around the $z$-axis (and a translation of the arc parameter
$s$).

In other words, any spherical curve
$\xi=(x,y,z):I\subseteq\mathbb{R}\rightarrow\mathbb{S}^{2}$, with $z$ non-
constant, is uniquely determined by its spherical angular momentum
$\mathcal{K}$ as a function of its coordinate $z$, that is, by
$\mathcal{K}=\mathcal{K}(z)$. The uniqueness is modulo rotations around the
$z$-axis. Moreover, the curvature of $\xi$ is given by
$\kappa(z)=\mathcal{K}^{\prime}(z)$.

######  Remark 2.2.

If we prescribe a continuous function $\kappa=\kappa(z)$ as curvature, the
proof of Theorem 2.1 clearly implies the computation of three quadratures,
following the sequence:

  1. (i)

A one-parameter family of anti-derivatives of $\kappa(z)$:

$\int\\!\kappa(z)dz=\mathcal{K}(z).$

  2. (ii)

Arc-length parameter $s$ of $\xi=(x,y,z)$ in terms of $z$, defined —up to
translations of the parameter— by the integral:

$s=s(z)=\int\\!\frac{dz}{\sqrt{1-z^{2}-\mathcal{K}(z)^{2}}},$

where $-\sqrt{1-z^{2}}<\mathcal{K}(z)<\sqrt{1-z^{2}}$, and inverting $s=s(z)$
to get $z=z(s)$ and so the latitude of $\xi$ is $\varphi(s)=\arcsin z(s)$.

  3. (iii)

Longitude of $\xi=(\cos\varphi\cos\lambda,\cos\varphi\sin\lambda,\sin\varphi)$
in terms of $s$, defined —up to a rotation around the $z$-axis— by the
integral:

$\lambda(s)=\int\\!\frac{\mathcal{K}(z(s))}{z(s)^{2}-1}ds,$

where $|z(s)|<1$.

We note that we get a one-parameter family of spherical curves satisfying
$\kappa=\kappa(z)$ according to the spherical angular momentum
$\mathcal{K}(z)$ chosen in (i) and verifying $\mathcal{K}(z)^{2}+z^{2}<1$. It
will distinguish geometrically the curves inside a same family by their
relative position with respect to the equator (or the $z$-axis).

We now show two illustrative examples applying steps (i)-(iii) of the
algorithm described in Remark 2.2:

######  Example 2.3 ($\kappa\\!\equiv\\!0$).

Then $\mathcal{K}\\!\equiv\\!c\in\mathbb{R}$, and
$s\\!=\\!\int\frac{dz}{\sqrt{1-c^{2}-z^{2}}}=\arcsin\frac{z}{\sqrt{1-c^{2}}}$,
with $|c|<1$. Therefore $z(s)=\sqrt{1-c^{2}}\sin s$. This gives that
$\lambda(s)=-\arctan(c\tan s)$ and finally $\xi(s)=(\cos s,-c\sin
s,\sqrt{1-c^{2}}\sin s)$. It corresponds to the great circle
$\mathbb{S}^{2}\cap\\{\sqrt{1-c^{2}}\,y+c\,z=0\\}$. Up to rotations around the
$z$-axis, they provide arbitrary great circles in $\mathbb{S}^{2}$, except the
equator. As a consequence of Theorem 2.1, we deduce that the great circle
$\mathbb{S}^{2}\cap\\{\sqrt{1-c^{2}}\,y+c\,z=0\\}$ is the only spherical curve
(up to rotations around the $z$-axis) with constant spherical angular momentum
$\mathcal{K}\\!\equiv\\!c$ (see Figure 1).

Figure 1. Great circles $\mathbb{S}^{2}\cap\\{\sqrt{1-c^{2}}\,y+c\,z=0\\}$:
$\mathcal{K}\equiv c\in(-1,1)$.

For example, taking $c=0$, we get the meridian $\mathbb{S}^{2}\cap\\{y=0\\}$
and hence the meridians are the only spherical curves with null spherical
angular momentum. We notice that the limiting cases $c=\pm 1$ lead to the
equator $\mathbb{S}^{2}\cap\\{z=0\\}$.

######  Example 2.4 ($\kappa\equiv k_{0}>0$).

Now $\mathcal{K}(z)=k_{0}z+c$, $c\in\mathbb{R}$. In this case, it is not
difficult to get that

$z(s)=\frac{1}{1+k_{0}^{2}}\left(\sqrt{1-c^{2}+k_{0}^{2}}\,\sin\left(\sqrt{1+k_{0}^{2}}\,s\right)-c\,k_{0}\right)$

with $|c|<\sqrt{1+k_{0}^{2}}$. But the expression of $\lambda$ is far from
trivial and depends on the values of $c$. After a long computation, we deduce:

  * •

If $|c|\neq k_{0}$:
$\lambda(s)=\arctan\left(\frac{\sqrt{1-c^{2}+k_{0}^{2}}+(1-ck_{0}+k_{0}^{2})\tan({\frac{1}{2}\sqrt{1+k_{0}^{2}}s})}{(k_{0}-c)\sqrt{1+k_{0}^{2}}}\right)+$

$+\arctan\left(\frac{\sqrt{1-c^{2}+k_{0}^{2}}+(1+ck_{0}+k_{0}^{2})\tan({\frac{1}{2}\sqrt{1+k_{0}^{2}}s})}{(k_{0}+c)\sqrt{1+k_{0}^{2}}}\right).$

  * •

If $c=k_{0}$:
$\lambda(s)=\arctan\left(\frac{1-(1+2k_{0}^{2})\tan(\frac{1}{2}\sqrt{1+k_{0}^{2}}s)}{2k_{0}\sqrt{1+k_{0}^{2}}}\right).$

  * •

If $c=-k_{0}$:
$\lambda(s)=\arctan\left(\frac{1+(1+2k_{0}^{2})\tan(\frac{1}{2}\sqrt{1+k_{0}^{2}}s)}{2k_{0}\sqrt{1+k_{0}^{2}}}\right).$

Of course, up to rotations around the $z$-axis, we get all the non-parallel
small circles of $\mathbb{S}^{2}$. The parameter $c$ distinguishes the
position of the circle with respect to the equator. If $0\leq|c|<1$ the
circles intersect the equator transversely; in particular, when $c=0$ we
obtain the orthogonal circles to the equator. If $c=\pm 1$, the circles are
tangent to the equator. Finally, if $1<|c|<\sqrt{1+k_{0}^{2}}$, the circles do
not intersect the equator (see Figure 2).

Figure 2.  Small circles: $\mathcal{K}(z)=k_{0}\,z+c,\ k_{0}>0$;
$0\leq|c|<1$ (left), $c=\pm 1$ (center), $1<|c|<\sqrt{1+k_{0}^{2}}$ (right).

For example, if $c=0$ we arrive (up to rotations around the $z$-axis) at
$\lambda(s)=\arctan\left(\cos(\sqrt{1+k_{0}^{2}}\,s)/k_{0}\right)$ and finally

$\xi(s)=\frac{1}{\sqrt{1+k_{0}^{2}}}\left(k_{0},\cos\left(\sqrt{1+k_{0}^{2}}\,s\right),\sin\left(\sqrt{1+k_{0}^{2}}\,s\right)\right),$

that corresponds to the small circle
$\mathbb{S}^{2}\cap\\{x=k_{0}/\sqrt{1+k_{0}^{2}}\\}$. This is the only
spherical curve, up to rotations around the $z$-axis, with spherical angular
momentum $\mathcal{K}(z)=k_{0}z$ (see Figure 2).

######  Remark 2.5.

The main difficulties one can find carrying on the strategy described in
Remark 2.2 (or in Theorem 2.1) to determine a spherical curve whose curvature
is $\kappa=\kappa(z)$ are the following:

  1. (1)

The integration of $s=s(z)$: Even in the case that $\mathcal{K}(z)$ were
polynomial, the corresponding integral is not necessarily elementary. For
example, when $\mathcal{K}(z)$ is a quadratic polynomial, it can be solved
using Jacobian elliptic functions (see e.g. [BF71]). This is equivalent to
$\kappa(z)$ be linear, i.e. $\kappa(z)=2az+b$, $a,b\in\mathbb{R}$. We will
study such spherical curves in Section 3.

  2. (2)

The previous integration gives us $s=s(z)$; it is not always possible to
obtain explicitly $z=z(s)$, what is necessary to determine the curve.

  3. (3)

Even knowing explicitly $z=z(s)$, the integration to get $\lambda(s)$ may be
impossible to perform using elementary or known functions.

Nevertheless, along the paper we will study different families where we are
successful with the procedure described in Remark 2.2 and we will recover some
known curves and find out new spherical curves characterized by their
spherical angular momentum.

##  3\. Elastic-type curves on the sphere

###  3.1. A new characterization and a generalization of elastic curves

A unit speed spherical curve $\xi$ is said to be an elastica under tension
$\sigma$ (see [LS84]) if its curvature $\kappa$ satisfies the differential
equation

(3.1) $2\ddot{\kappa}+\kappa^{3}+(2-\sigma)\kappa=0$

for some $\sigma\in\mathbb{R}$. They are the critical points of the elastic
energy functional

$\int_{\xi}(\kappa^{2}+\sigma)ds$

acting on spherical curves with suitable boundary conditions. If $\sigma=0$,
then the constraint on arc length is removed and $\xi$ is called a free
elastica.

A possible generalization of free elasticae was considered in [AGM06], where
the authors studied the elastic curves in $\mathbb{S}^{2}$ which are circular
at rest. They are called $\lambda$-elastic curves. These curves are critical
points of the functional

$\int_{\xi}(\kappa+\lambda)^{2}ds,\,\lambda\in\mathbb{R}$

and are characterized by the Euler-Lagrange equation

(3.2) $2\ddot{\kappa}+\kappa^{3}+(2-\lambda^{2})\kappa+2\lambda=0.$

It is obvious that the 0-elastic curves are the free elasticae. The main
result of this section deals with the spherical curves which are the critical
points of the bending energy for variations with constant length (including
both previous types of elastic curves) relating them with the case commented
in part (1) of Remark 2.5.

######  Theorem 3.1.

Let $\xi$ be a spherical curve whose curvature $\kappa$ satisfies

(3.3) $\kappa=2a\langle\xi,{\bf e}\rangle+b,\ a\neq 0,\,b\in\mathbb{R},$

for some ${\bf e}\in\mathbb{R}^{3},\ |{\bf e}|=1$. Then:

  1. (i)

the spherical angular momentum $\mathcal{K}$ of $\xi$ is given by

$\mathcal{K}=a\langle\xi,{\bf e}\rangle^{2}+b\langle\xi,{\bf e}\rangle+c,$

for some $c\in\mathbb{R}$.

  2. (ii)

$\xi$ is a critical point of the functional

$\int_{\xi}(\kappa^{2}-2b\kappa+b^{2}-4ac)\,ds$

and so $\kappa$ satisfies the corresponding Euler-Lagrange equation

(3.4) $2\ddot{\kappa}+\kappa^{3}+\left(2-(b^{2}-4ac)\right)\kappa-2b=0.$

If $b=0$, $\xi$ is an elastica under tension $\sigma=-4ac$; and if $c=0$,
$\xi$ is a $\lambda$-elastic curve with $\lambda=-b$.

Conversely, if $\xi$ is a critical point of the functional

${\mathcal{F}}_{\sigma}^{\lambda}(\xi):=\int_{\xi}\left((\kappa+\lambda)^{2}+\sigma\right)\,ds,\
\lambda,\sigma\in\mathbb{R},$

then the curvature of $\xi$ can be written as in (3.3).

######  Remark 3.2.

It is remarkable the similitude of condition (3.3) characterizing the
generalized elasticae considered in Theorem 3.1 with the geometric property
satisfied by the classical Euler elastic curves in the plane: their curvature
is proportional to one of the coordinate functions, say $\kappa(x,y)=2\lambda
y+\mu$, $\lambda\neq 0,\,\mu\in\mathbb{R}$ (see Section 3 in [CCI16] and
Section 1 in [S99]). Even something similar happens to spacelike and timelike
elastic curves in Lorentz-Minkowski plane (see [CCIs18a] and [CCIs18b]).

###### Proof of Theorem 3.1.

From (2.2), (2.1) and (3.3), we get that

$\frac{d}{ds}\left(-\mathcal{K}+a\langle\xi,{\bf
e}\rangle^{2}+b\langle\xi,{\bf e}\rangle\right)=0.$

This proves part (i).

We also have from (3.3), (2.1) and (2.2), that
$\dot{\kappa}=2a\langle\dot{\xi},{\bf e}\rangle$ and
$\ddot{\kappa}=2a(-\langle\xi,{\bf e}\rangle-\kappa\mathcal{K})$. Now we can
easily check that $\kappa$ given by (3.3) satisfies (3.4) since, after a
straightforward computation using (i) and putting $\langle\xi,{\bf
e}\rangle=(\kappa-b)/2a$, we arrive at (3.4). We observe that if $b=0$ then
$\xi$ satisfies (3.1) with $\sigma=-4ac$ and if $c=0$ then $\xi$ satisfies
(3.2) with $\lambda=-b$. Moreover, in [AGM03] it is shown that for a given
differentiable function $P(\kappa)$, the critical points of the functional
$\int_{\xi}P(\kappa)ds$ are characterized by the Euler-Lagrange equation

(3.5)
$(\kappa^{2}+1)P^{\prime}(\kappa)+\frac{d^{2}(P^{\prime}(\kappa))}{ds^{2}}=\kappa
P(\kappa).$

It is an exercise to check that putting
$P(\kappa)=\kappa^{2}-2b\kappa+b^{2}-4ac$ in (3.5) we obtain (3.4). This
finishes the proof of part (ii).

Multiplying (3.4) by $\dot{k}$ we obtain a first integral

$\dot{\kappa}^{2}+\kappa^{4}/4+\left(1-(b^{2}-4ac)/2\right)\kappa^{2}-2b\kappa=E,$

where $E$ is a real constant. After a long computation involving (3.3), part
(i), (2.2) and using that $\langle\dot{\xi},{\bf
e}\rangle^{2}=1-\langle\xi,{\bf e}\rangle^{2}-\langle N,{\bf e}\rangle^{2}$,
since $\\{\xi,\dot{\xi},N\\}$ is an orthonormal basis in $\mathbb{R}^{3}$, we
get that

(3.6) $E=4a^{2}-b^{2}-(b^{2}-4ac)^{2}/4.$

On the other hand, suppose now that $\xi$ is a critical point of
${\mathcal{F}}_{\sigma}^{\lambda}$. Taking into account (3.5) we deduce that
$\kappa$ verifies the differential equation

(3.7) $2\ddot{\kappa}+\kappa^{3}+(2-\lambda^{2}-\sigma)\kappa+2\lambda=0.$

Multiplying (3.7) by $\dot{\kappa}$ and integration allow us to deduce the
energy $E\in\mathbb{R}$ of $\xi$:

(3.8)
$E:=\dot{\kappa}^{2}+\frac{\kappa^{4}}{4}+\left(1-\frac{\lambda^{2}+\sigma}{2}\right)\kappa^{2}+2\lambda\kappa.$

We want to prove that $\kappa$ has an expression like in (3.3). For this
purpose, we first observe that if $\lambda=0$ then $\kappa\equiv 0$ is a
trivial solution to (3.7). For example, we can take ${\bf
e}\in\mathbb{R}^{3},\ |{\bf e}|=1$ the unit normal vector orthogonal to the
vectorial plane containing the corresponding great circle in $\mathbb{S}^{2}$.
Now we must look for $a\neq 0$ and $b\in\mathbb{R}$ satisfying (3.3).
Comparing (3.4) and (3.7), we take $b=-\lambda$ and observe that $\sigma$ must
satisfy $\sigma=-4ac$, with $a\neq 0$, $c\in\mathbb{R}$. Using (3.8), we have
that
$4E+4\lambda^{2}+(\lambda^{2}+\sigma)^{2}=4\dot{\kappa}^{2}+4(\kappa+\lambda)^{2}+\left(\kappa^{2}-(\lambda^{2}+\sigma)\right)^{2}>0$.

But looking at (3.6), $E$ must satisfy $E=4a^{2}-b^{2}-(b^{2}-4ac)^{2}/4$ and,
eliminating $c$, we finally arrive at
$0<4E+4\lambda^{2}+(\lambda^{2}+\sigma)^{2}=16a^{2}$, which allows us to
obtain the searched value for $a$. ∎

The study of (free) elastic curves on the sphere has been considered under
different approaches (see for example [AGM03], [BC94], [J95], [LS84] or
[S08]), paying special attention to the closed ones. The closed
$\lambda$-elastic spherical curves were studied in [AGM06]. All the mentioned
articles are based on the study of the differential equation for the geodesic
curvature of the spherical curve, being sometimes integrated directly in terms
of Jacobi elliptic functions.

In our approach of Theorem 3.1, we can choose ${\bf e}=(0,0,1)$ without loss
of generality and so $\langle\xi,{\bf e}\rangle=z$. In this way we arrive at
the conditions

(3.9) $\kappa(z)=2az+b,\,\mathcal{K}(z)=az^{2}+bz+c,\,a\neq
0,b,c\in\mathbb{R},$

in the notation of Theorem 2.1. Thus, using suitable coordinates in the
sphere, we can conclude the following uniqueness result for the spherical
elastic curves considered in literature.

######  Corollary 3.3.

  1. (a)

The elasticae under tension $\sigma$ are the only spherical curves (up to
rotations around $z$-axis) with spherical angular momentum
$\mathcal{K}(z)=az^{2}+c,\,a\neq 0,c\in\mathbb{R}$, with $\sigma=-4ac$. In
particular, the free elasticae are characterized by the spherical angular
momentum $\mathcal{K}(z)=az^{2},\,a\neq 0$.

  2. (b)

The $\lambda$-elastic curves are the only spherical curves (up to rotations
around $z$-axis) with spherical angular momentum
$\mathcal{K}(z)=az^{2}+bz,\,a\neq 0,b\in\mathbb{R}$, with $\lambda=-b$.

######  Remark 3.4.

Inspired by the Langer and Singer work on the Kirchhoff elastic rod [LS96], in
[AGM04] it is obtained by geometrical means the first integrals of the Euler-
Lagrange equations of curvature energy functionals $\int_{\xi}P(\kappa)ds$,
where $P$ is a smooth function and $\kappa$ denotes the curvature of the
spherical curve $\xi$ in $\mathbb{S}^{3}$. Assuming that
$P^{\prime\prime}(\kappa)\neq 0$, the critical points of the functional
$\int_{\xi}P(\kappa)ds$ are characterized by a couple of differential
equations naturally related to a system of cylindrical coordinates in the
three-sphere adapted to the curve $\xi$.

In the case that the curve $\xi$ lies in $\mathbb{S}^{2}$, from Section 3 of
[AGM04], using concretely equation (33) with $\theta=0$ ($b=0$), then we have
that $\cos^{2}\psi=P^{\prime}(k)^{2}/a^{2}$, $a\neq 0$, where
$\psi=\pi/2-\varphi$ is the colatitude of $\xi$.

As a consequence, if a spherical curve $\xi$ is a critical point of
$\int_{\xi}P(\kappa)ds$, then there exist geographical coordinates
$(\varphi,\lambda)$ adapted to $\xi$ such that

$\sin\varphi=\delta\,P^{\prime}(\kappa),\,\delta\neq 0.$

This result is consistent with Theorem 3.1, since if
$P(\kappa)=(\kappa+\lambda)^{2}+\sigma$, then
$\sin\varphi=z=\frac{\kappa-b}{2a}=\delta\,P^{\prime}(\kappa)$, taking
$\delta=\frac{1}{4a}$ and $\lambda=-b$.

If one follows the strategy described in Remark 2.2 to determine the
generalized elasticae satisfying (3.9), it is necessary to perform the
integral

(3.10) $s=s(z)=\int\frac{dz}{\sqrt{P(z)}},\ P(z)=1-z^{2}-(az^{2}+bz+c)^{2}.$

Since $P(z)$ is a fourth order polynomial, (3.10) can be solved in terms of
Jacobi elliptic functions once the nature and multiplicity of the roots of
$P(z)$ are determined. After inverting $s=s(z)$ to get $z=z(s)$ we would
arrive at the expressions of $\kappa=\kappa(s)=2az(s)+b$ compatible with the
ones given in [AGM03], [BC94], [J95], [LS84] or [S08] and [AGM06], at least
when $b=0$ or $c=0$. Besides helices and circles, the borderline, orbitlike
and wavelike elasticae appear and, in addition, a more general case according
to the expressions of the geodesic curvature in terms of the arc parameter
(see [S08]). In order to get the explicit expression of the generalized
elasticae, we also have to compute

(3.11) $\lambda(s)=\int\\!\frac{az(s)^{2}+bz(s)+c}{z(s)^{2}-1}\,ds,$

which, in general, leads to complicated elliptic integrals. In the next two
sections we are going to illustrate the aforementioned computations in two
interesting and attractive cases.

###  3.2. Seiffert’s spherical spirals

The Seiffert’s spirals are defined as those spherical curves obtained when one
moves along the surface of a sphere with constant speed, while maintaining a
constant angular velocity with respect to a fixed diameter (cf. [E00]). These
curves are given in cylindrical coordinates $(r,\theta,z)$, $r^{2}+z^{2}=1$,
by the parametric equations

(3.12) $r=\operatorname{sn}(s,p),\ \theta=p\,s,\ z=\operatorname{cn}(s,p),\
(p>0),$

where $p$ is a positive constant and $\operatorname{sn}$ and
$\operatorname{cn}$ are the elementary Jacobi elliptic functions (cf. [BF71]
for instance). Erdös provided in [E00] a derivation of the equations of this
curve, as well as an analysis of its properties, including conditions for
obtaining periodic orbits. When $p>1$, the spiral is located entirely in the
northern hemisphere.

Now we prove that these curves are elastic curves with positive tension
corresponding to the conditions $\kappa(z)=2az$, $\mathcal{K}(z)=az^{2}-a$,
i.e. $b=0$, $a+c=0$ in (3.9). Then $\sigma=4a^{2}>0$ and there is no
restriction if we consider $a>0$. So (3.10) can be written as

$\left(\frac{dz}{ds}\right)^{2}=(1-z^{2})\left(1-a^{2}+a^{2}z^{2}\right),\
\frac{a^{2}-1}{a^{2}}<z^{2}<1,$

which implies that $z(s)=\operatorname{cn}(s,a)$ (cf. [BF71] for instance) and
thus $r(s)=\operatorname{sn}(s,a)$. Using that $a+c=0$ in (3.11) we get that
$\lambda(s)=a\,s$ and so we arrive at the Seiffert’s spirals (see Figure 3).
As a summary:

######  Corollary 3.5.

The Seiffert’s spirals (3.12) are the only spherical curves (up to rotations
around $z$-axis) with spherical angular momentum
$\mathcal{K}(z)=pz^{2}-p,\,p>0$.

Figure 3. Seiffert’s spirals ($p\approx 0$, $0<p<1$, $p\approx 1$, $p>1$):
$\mathcal{K}(z)=pz^{2}-p$, $p>0$.

###  3.3. Borderline spherical elastic curves

We study in this section elastic curves with null energy producing elasticae
under positive tension $\sigma>0$. This corresponds to the conditions
$\kappa(z)=2az$, $\mathcal{K}(z)=az^{2}-1$, i.e. $b=0$, $c=-1$ in (3.9). Then
$E=0$ according to (3.6) and we can take $a>0$ in order to $\sigma=4a>0$. Now
(3.10) leads to

$s=s(z)=\int\frac{dz}{z\,\sqrt{2a-1-a^{2}z^{2}}},\ \,2a-1-a^{2}z^{2}>0,$

which implies that $a>1/2$. The above integral becomes elementary and, after
inverting $s=s(z)$ and up to a translation on the parameter $s$, we get

(3.13) $z(s)=\frac{2a-1}{a}\operatorname{sech}(\sqrt{2a-1}\,s).$

Looking at (3.11), we get that if $a=1$ then $\lambda(s)=s$ (and
$z(s)=\operatorname{sech}s$), and when $a>1/2$ with $a\neq 1$ we obtain

(3.14)
$\lambda(s)=s+\arctan\left(\frac{\sqrt{2a-1}}{1-a}\tanh(\sqrt{2a-1}s)\right).$

This family corresponds to the “borderline elasticae” described in [S08] which
are asymptotic to the equator. We show some pictures of them in Figure 4. In
conclusion, we deduce the following uniqueness result:

######  Corollary 3.6.

The borderline elasticae given by (3.13) and (3.14) are the only spherical
curves (up to rotations around $z$-axis) with spherical angular momentum
$\mathcal{K}(z)=az^{2}-1,\,a>1/2$.

We remark that when $a=1$ we recover the Seiffert’s spiral corresponding to
$p=1$, since $\operatorname{cn}(s,1)=\operatorname{sech}s$ (see [BF71]).

Figure 4. Borderline elastic curves: $\mathcal{K}(z)=az^{2}-1,\,a>1/2$ (left:
$1/2<a\leq 1$; right: $a\geq 1$).

##  4\. Loxodromic-type curves on the sphere

In this section we are interested in critical points of the functional

$\mathcal{F}^{\mu}(\xi):=\int_{\xi}\sqrt{\kappa^{2}+\mu^{2}}\,ds,\,\mu>0.$

This functional was considered in Section 6 of [AGM03] acting on the space of
immersed closed curves in $\mathbb{S}^{2}$, motivated by total
$\mathbb{R}^{3}$-curvature type functionals. We aim to connect the critical
points of $\mathcal{F}^{\mu}$ with the classical loxodromes curves in
$\mathbb{S}^{2}$ and others spherical curves with similar characteristics.

Taking into account Remark 3.4, we are devoted to study the spherical curves
with curvature

(4.1) $\kappa(z)=\frac{\mu\,z}{\sqrt{\delta^{2}-z^{2}}},\ \mu>0,\delta\neq 0.$

We will distinguish cases according to $0<\mu<1$, $\mu=1$ and $\mu>1$.

###  4.1. Case $0<\mu<1$: spherical loxodromes

The loxodromes are interesting curves in the sphere studied, among others, by
Pedro Nunes in 1537, Simon Stevin in 1608 or Maupertuis in 1744. They are also
known as rhumb lines because they make a constant angle $\alpha\in(0,\pi/2)$
with the meridians (cf. [F93]). Analytically, using geographical coordinates
$(\varphi,\lambda)$, they are defined by the equation

(4.2) $d\lambda=\cot\alpha\frac{d\varphi}{\cos\varphi}.$

The aim of this section is the study of the spherical curves satisfying

(4.3) $\kappa(\varphi)=a\tan\varphi,\ 0<a<1,$

or, equivalently,

(4.4) $\kappa(z)=\frac{az}{\sqrt{1-z^{2}}},\ 0<a<1.$

So they correspond in (4.1) to the election $\mu=a\in(0,1)$ and $\delta=1$.

The trivial solution of (4.4) is given by the equator $z=0$. We follow the
method described in Theorem 2.1 and Remark 2.2, considering the spherical
angular momentum

$\mathcal{K}(z)=-a\sqrt{1-z^{2}},\ 0<a<1.$

Then we have:

$s\\!=\\!\int\frac{dz}{\sqrt{1-z^{2}-(-a\sqrt{1-z^{2}})^{2}}}\\!=\\!\int\frac{dz}{\sqrt{(1-a^{2})(1-z^{2})}}\\!=\\!\frac{1}{\sqrt{1-a^{2}}}\arcsin
z,$

and so $z(s)=\sin(\sqrt{1-a^{2}}\,s)$ and therefore

(4.5) $\varphi(s)=\sqrt{1-a^{2}}\,s.$

On the other hand, from (4.5), we get:

(4.6)
$\lambda(s)\\!=\\!\int\frac{a\,ds}{\cos\varphi(s)}\\!=\\!\frac{a}{\sqrt{1\\!-\\!a^{2}}}\log\left(\sec(\sqrt{1\\!-\\!a^{2}}\,s)+\tan(\sqrt{1\\!-\\!a^{2}}\,s)\right),$

where $|s|<\dfrac{\pi}{2\sqrt{1-a^{2}}}$.

We deduce from (4.5) and (4.6) that

$d\lambda=\frac{a}{\sqrt{1-a^{2}}}\,\frac{d\varphi}{\cos\varphi}$

and, taking into account (4.2) and that $\mathcal{K}(\varphi)=-a\cos\varphi$,
we conclude the following characterization of the loxodromes.

######  Corollary 4.1.

The loxodromes $d\lambda=\cot\alpha\,d\varphi/\cos\varphi$ are the only
spherical curves (up to rotations around $z$-axis) with spherical angular
momentum $\mathcal{K}(\varphi)=-\cos\alpha\cos\varphi$, $\alpha\in(0,\pi/2)$.

From (4.3) and (4.5), we arrive at the intrinsic equation of the loxodromes
(see Figure 5), given by

$\kappa(s)=\cos\alpha\,\tan(\sin\alpha\,s),\,\alpha\in(0,\pi/2).$

Figure 5. Loxodromes: $\mathcal{K}(\varphi)=-\cos\alpha\cos\varphi$,
$\alpha\in(0,\pi/2)$.

###  4.2. Case $\mu=1$

We now want to determine the spherical curves whose curvature is given by

(4.7) $\kappa(z)=\frac{z}{\sqrt{a-z^{2}}},\ 0<a<1.$

Looking at (4.1), we are now considering $\mu=1$ and $\delta=\sqrt{a}$.

The trivial solution is given by the equator $z=0$. We follow the method
described in Theorem 2.1 or Remark 2.2 for the non trivial case, considering
the spherical angular momentum

$\mathcal{K}(z)=-\sqrt{a-z^{2}}.$

In this way, we get:

$s=s(z)=\int\frac{dz}{\sqrt{1-z^{2}-(-\sqrt{a-z^{2}})^{2}}}=\int\frac{dz}{\sqrt{1-a}}=\frac{z}{\sqrt{1-a}},$

and thus $z(s)=\sqrt{1-a}\,s$.

We write $a=\sin^{2}\alpha$ , $0<\alpha<\pi/2$, and abbreviate
$c_{\alpha}=\cos\alpha,s_{\alpha}=\sin\alpha$. Hence we have:

$\lambda(s)=\int\frac{\sqrt{s_{\alpha}^{2}-c_{\alpha}^{2}s^{2}}}{1-c_{\alpha}^{2}s^{2}}\,ds,$

which gives

$\begin{array}[]{c}\lambda(s)=\dfrac{1}{c_{\alpha}}\arctan\left(\dfrac{c_{\alpha}s}{\sqrt{s_{\alpha}^{2}-c_{\alpha}^{2}s^{2}}}\right)\\\
-\dfrac{1}{2}\arctan\left(\dfrac{c_{\alpha}s+s_{\alpha}^{2}}{c_{\alpha}\sqrt{s_{\alpha}^{2}-c_{\alpha}^{2}s^{2}}}\right)-\dfrac{1}{2}\arctan\left(\dfrac{c_{\alpha}s-s_{\alpha}^{2}}{c_{\alpha}\sqrt{s_{\alpha}^{2}-c_{\alpha}^{2}s^{2}}}\right),\end{array}$

where $|s|<\tan\alpha$.

We observe that $|z(s)|<\sin\alpha$. Using (4.7) and $z(s)=c_{\alpha}s$, we
get the intrinsic equation

$\kappa(s)=\dfrac{c_{\alpha}s}{\sqrt{s_{\alpha}^{2}-c_{\alpha}^{2}s^{2}}},\
|s|<\tan\alpha,$

of this family of spherical curves of loxodromic-type (see Figure 6)
characterized by the geometric angular momentum
$\mathcal{K}(z)=-\sqrt{\sin^{2}\alpha-z^{2}}$, $0<\alpha<\pi/2$.

Figure 6. Spherical curve with $\mathcal{K}(z)=-\sqrt{a-z^{2}}$, $0<a<1$.

###  4.3. Case $\mu>1$

Finally we wish to study the spherical curves whose curvature is given by

(4.8) $\kappa(z)=\frac{az}{\sqrt{1-az^{2}}},\,a>1.$

We observe that corresponds in (4.1) with the elections $\mu=\sqrt{a}$ and
$\delta=1/\sqrt{a}$. The trivial solution is given by the equator $z=0$.
Otherwise, we make use of the method proposed in Theorem 2.1 or Remark 2.2,
considering the spherical angular momentum

$\mathcal{K}(z)=-\sqrt{1-az^{2}}.$

In this way, we obtain:

$s=s(z)=\int\frac{dz}{\sqrt{1-z^{2}-(-\sqrt{1-az^{2}})^{2}}}=\int\frac{dz}{\sqrt{a-1}\,|z|}=\frac{\log|z|}{\sqrt{a-1}},$

and so $|z(s)|=e^{\sqrt{a-1}\,s}>0$.

If we write $a=\cosh^{2}\delta$, $\delta>0$, then $z(s)=\pm
e^{\sinh\delta\,s}$ and

$\begin{array}[]{c}\lambda(s)=\int\frac{\sqrt{1-\cosh^{2}\delta\,e^{2\sinh\delta\,s}}}{1-e^{2\sinh\delta\,s}}ds=\\\
\\\
\\!=\\!-\dfrac{\operatorname{arctanh}\left(\sqrt{1\\!-\\!\cosh^{2}\delta\,e^{2\sinh\delta\,s}}\right)}{\sinh\delta}\\!+\\!\arctan\left(\dfrac{\sqrt{1\\!-\\!\cosh^{2}\delta\,e^{2\sinh\delta\,s}}}{\sinh\delta}\right),\end{array}$

where $s<-\log\cosh\delta/\sinh\delta$.

Using (4.8), we get the intrinsic equation

$|\kappa(s)|=\dfrac{\cosh^{2}\delta\,e^{\sinh\delta\,s}}{\sqrt{1-\cosh^{2}\delta\,e^{2\sinh\delta\,s}}},\
s<-\log\cosh\delta/\sinh\delta,$

of this family of spherical curves of loxodromic-type (see Figure 7)
characterized by the geometric angular momentum
$\mathcal{K}(z)=-\sqrt{1-\cosh^{2}\delta z^{2}}$.

Figure 7. Spherical curve with $\mathcal{K}(z)=-\sqrt{1-az^{2}}$, $a>1$.

##  5\. Spherical catenaries

In this section we are interested in critical points of the functional

$\mathcal{F}(\xi):=\int_{\xi}\sqrt{\kappa}\,ds.$

The above functional $\int_{\xi}\kappa^{1/2}\,ds$ was considered in Section 5
of [AGM03] acting on the space of convex ($\kappa>0$) closed curves in
$\mathbb{S}^{2}$, motivated by the study of ($r=1/2$)-generalized elastic
curves in $\mathbb{S}^{2}$. We aim to connect the critical points of
$\mathcal{F}$ with the classical catenaries curves in $\mathbb{S}^{2}$.

Taking into account Remark 3.4, we are devoted to study the spherical curves
with curvature

(5.1) $\kappa(z)=\frac{\delta^{2}}{4z^{2}},\,\delta\neq 0.$

The spherical catenaries are the equilibrium lines of an inelastic flexible
homogeneous infinitely thin massive wire included in a sphere, placed in a
uniform gravitational field. Like any catenaries, their centres of gravity
have the minimal altitude among all the curves with given length passing by
two given points. They were studied by Bobillier in 1829 and by Gudermann in
1846 (cf. [F93]). See Figure 8.

Figure 8. Closed spherical catenaries.

Using cylindrical coordinates $(r,\theta,z)$ in $\mathbb{R}^{3}$, they can be
described analytically by the following first integral of the corresponding
ordinary differential equation:

(5.2) $(z-z_{0})\,r^{2}\,\frac{d\theta}{ds}={\rm constant}.$

We study in this section spherical curves satisfying the condition

(5.3) $\kappa(z)=a/z^{2},\ a>0.$

So they correspond in (5.1) to the election $\delta=2\sqrt{a}$. For any $a>0$,
it is easy to prove that there exists an unique angle
$\varphi_{0}\in(0,\pi/2)$ such that $a=\tan\varphi_{0}\sin^{2}\varphi_{0}$.
Thus the parallel $z=\sin\varphi_{0}$ is a constant solution to (5.3).

We now apply Theorem 2.1 and Remark 2.2 considering (5.3) and

$\mathcal{K}(z)=-\frac{a}{z},\ a>0.$

Then we have:

$s=s(z)=\int\frac{z\,dz}{\sqrt{z^{2}(1-z^{2})-a^{2}}},$

which implies that $a<1/2$ and $1-\sqrt{1-4a^{2}}<2z^{2}<1+\sqrt{1-4a^{2}}$,
and it is not difficult to get:

(5.4) $z(s)=\sqrt{\frac{1+\sqrt{1-4a^{2}}\sin 2s}{2}}.$

In addition, we have that

(5.5) $d\lambda=\frac{a}{z(1-z^{2})}\,ds$

Looking at (5.2), taking into account that $r^{2}+z^{2}=1$ and
$\theta=\lambda$, we deduce from (5.5) that we get a spherical catenary (with
$z_{0}=0$ and constant $a\in(0,1/2)$). However, the explicit computation of
$\lambda$ in terms of the arc parameter $s$ requires elliptic integrals of the
first and third kind. As a summary, we have proved the following uniqueness
result.

######  Corollary 5.1.

The spherical catenaries (5.5) are the only spherical curves (up to rotations
around the $z$-axis) with spherical angular momentum $\mathcal{K}(z)=-a/z$
(and curvature $\kappa(z)=a/z^{2}$), $0<a<1/2$.

Combining (5.3) and (5.4) we have that the intrinsic equation of the spherical
catenaries is given by

$\kappa(s)=\frac{2a}{1+\sqrt{1-4a^{2}}\sin 2s},\ 0<a<1/2.$

##  6\. New and classical spherical curves

The purpose of this section is to find out new curves $\xi=(x,y,z)$ in
$\mathbb{S}^{2}$, expressed in terms of elementary functions or in terms of
Jacobi elliptic functions, prescribing their curvature as a function of the
distance from the equator in such a way we can avoid the difficulties
described in Remark 2.5. In addition, we provide uniqueness results for some
well known spherical curves in terms of the spherical angular momentum
introduced in Section 2.

###  6.1. Spherical curves such that $\kappa(\varphi)=p\cos
2\varphi/\cos\varphi,\,0<p<1$

The purpose of this section is to find out new curves in $\mathbb{S}^{2}$
expressed in terms of Jacobi elliptic functions prescribing in a suitable way
their curvature in terms of their latitude. Concretely, we aim to study the
spherical curves whose curvature is given by

(6.1) $\kappa(\varphi)=\frac{p\cos 2\varphi}{\cos\varphi},\ 0<p<1.$

Recalling that $z=\sin\varphi$, (6.1) is equivalent to

(6.2) $\kappa(z)=\frac{p(1-2z^{2})}{\sqrt{1-z^{2}}},\ 0<p<1.$

We follow the strategy proposed by Theorem 2.1 or Remark 2.2, considering the
spherical angular momentum

$\mathcal{K}(z)=p\,z\sqrt{1-z^{2}}.$

Then we get:

$s=\int\frac{dz}{\sqrt{(1-z^{2})(1-p^{2}z^{2})}}=\int\frac{d\varphi}{\sqrt{1-p^{2}\sin^{2}\varphi}}=F(\varphi,p)=F(\arcsin
z,p),$

where $F(\cdot,p)$ denotes the elliptic integral of first class of modulus $p$
(see e.g. [BF71]). Hence $\varphi(s)=\operatorname{am}(s,p)$ and
$z(s)=\operatorname{sn}(s,p)$, where $\operatorname{am}(\cdot,p)$ is the
Jacobi amplitude and $\operatorname{sn}(\cdot,p)$ is the Jacobi sine amplitude
of modulus $p$ (see e.g. [BF71]). In addition:

$\lambda(s)=-p\int\frac{\operatorname{sn}(s,p)}{\operatorname{cn}(s,p)}ds$

where $\operatorname{cn}(\cdot,p)$ is the Jacobi cosine amplitude of modulus
$p$. Using formula 316.01 of [BF71], we finally arrive at the following
expression for the longitude:

$\lambda(s)=-\frac{p}{2p^{\prime}}\log\left(\frac{\operatorname{dn}(s,p)+p^{\prime}}{\operatorname{dn}(s,p)-p^{\prime}}\right),$

where $\operatorname{dn}(\cdot,p)$ is the Jacobi delta amplitude of modulus
$p$ and $p^{\prime}=\sqrt{1-p^{2}}$ is the complementary modulus.

Using (6.2) and that $z(s)=\operatorname{sn}(s,p)$, joint to formula 124.02 of
[BF71], we get the intrinsic equation

$\kappa(s)=p\left(2\operatorname{cn}(s,p)-1/\operatorname{cn}(s,p)\right),\
0<p<1,$

of the only spherical curves (up to rotations around $z$-axis) with spherical
angular momentum $\mathcal{K}(z)=p\,z\sqrt{1-z^{2}}$ or, equivalently,
$\mathcal{K}(\varphi)=(p/2)\sin 2\varphi$, $0<p<1$.

These curves are embedded and closed since $\xi(s+4K(p))=\xi(s)$, where $K(p)$
is the complete elliptic integral of first class of modulus $p$ (see Figure
9).

Figure 9. Spherical curves with
$\mathcal{K}(z)\\!=\\!p\,z\sqrt{1\\!-\\!z^{2}}$, $0\\!<\\!p\\!<\\!1$.

###  6.2. Viviani’s curve and spherical Archimedean spirals

The Viviani’s curve is the intersection between a sphere of radius $R$ and a
cylinder of revolution with diameter $R$ such that a generatrix passes by the
centre of the sphere; so this curve is at the same time spherical and
cylindrical. We can obtain a Viviani’s curve by sticking the tip of a compass
inside a cylinder of revolution and tracing on this cylinder a “circle” with
radius equal to the cylinder diameter. It was studied by Vincenzo Viviani in
1692 (cf. [F93]). In geographical coordinates of $\mathbb{S}^{2}$, the
Viviani’s curve can be simply described as $\varphi=\lambda$ (see Figure 10).

We study in this section spherical curves satisfying the condition

(6.3) $\kappa(z)=\frac{z(3-z^{2})}{(2-z^{2})^{3/2}},$

applying Theorem 2.1 and Remark 2.2, considering the spherical angular
momentum

$\mathcal{K}(z)=\frac{z^{2}-1}{\sqrt{2-z^{2}}}.$

Then we have:

(6.4) $s=s(z)=\int\sqrt{\frac{2-z^{2}}{1-z^{2}}}\,dz=E(\arcsin z,1/2)$

which involves an elliptic integral $E$ of second kind (see e.g. [BF71]). In
addition, we get:

(6.5) $d\lambda=\frac{ds}{\sqrt{2-z^{2}}}$

Using that $z=\sin\varphi$, (6.4) and (6.5), we get easily that
$d\lambda=d\varphi$. Hence we have proved the following characterization of
the Viviani’s curve.

######  Corollary 6.1.

The Viviani’s curve $\varphi=\lambda$ is the only spherical curve (up to
rotations around the $z$-axis) with spherical angular momentum
$\mathcal{K}(\varphi)=-\cos^{2}\varphi/\sqrt{1+\cos^{2}\varphi}$.

Figure 10. Viviani’s curve.

The spherical Archimedean spiral curves are natural generalizations of the
Viviani’s curve, since they are described in geographical coordinates by
$\varphi=n\lambda$, $n>0$. A spherical Archimedean spiral is algebraic if and
only if $n\in\mathbb{Q}$. They were studied by Guido Grandi in 1728, also
called clelias. They are the loci of a point $P$ on a meridian of a sphere
rotating at constant speed $\omega$ around the polar axis, the point $P$ also
moving at constant speed $n\,\omega$ along this meridian (see Figure 11).
Therefore, physically, we obtain a clelia when peeling an orange or when
rewinding regularly a spherical wool ball.

Figure 11. Spherical Archimedean spiral curves: $\varphi=n\lambda$, $n=1/3$
and $n=3$.

A similar argument to the used in the preceding section, considering now

(6.6) $\kappa(z)=\frac{z(2n^{2}+1-z^{2})}{(n^{2}+1-z^{2})^{3/2}}$

give us the following uniqueness result, whose proof we will omit.

######  Corollary 6.2.

The spherical spiral curves $\varphi=n\lambda$, $n>0$, are the only spherical
curves (up to rotations around the $z$-axis) with spherical angular momentum
$\mathcal{K}(\varphi)=-\cos^{2}\varphi/\sqrt{n^{2}+\cos^{2}\varphi}$.

### Acknowledgements

The authors were partially supported by State Research Agency and European
Regional Development Fund via the grant no. PID 2020-117868GB-I00 supported by
MCIN/AEI/10.13039/ 501100011033. Second named author also by MICINN/FEDER
project PGC2018-097046-B-I00 and Fundación Séneca project 19901/GERM/15,
Spain. Third named author was supported by Predoc grant no. FPU16/03096
supported by MCIN/AEI/10.13039/501100011033.

We would also like to thank Álvaro Pámpano for helpful discussion about the
topic, specially for comments concerning Remark 3.4.

## References

  * [A78] V.I. Arnold.  Mathematical methods of classical mechanics. Springer Verlag, 1978.
  * [AGM03] J. Arroyo, O.J. Garay and J.J. Mencía.  Closed generalized elastic curves in $\mathbb{S}^{2}(1)$. J. Geom. Phys. 48 (2003), 339–353.
  * [AGM04] J. Arroyo, O.J. Garay and J.J. Mencía.  Extremals of curvature energy actions on spherical closed curves J. Geom. Phys. 51 (2004), 101–125.
  * [AGM06] J. Arroyo, O.J. Garay and J.J. Mencía.  Elastic circles in 2-spheres. J. Phys. A: Math. Gen. 39 (2006), 2307–2324.
  * [BC94] G. Brunnett and P.E. Crouch.  Elastic curves on the sphere Adv. Comp. Math. 2 (1994), 23–40.
  * [BF71] P.F. Byrd & M.D. Friedman.  Handbook of elliptic integrals for engineers and physicists. 2nd edition. Springer Verlag, 1971.
  * [CCI16] I. Castro and I. Castro-Infantes.  Plane curves with curvature depending on distance to a line. Diff. Geom. Appl. 44 (2016), 77–97.
  * [CCIs17] I. Castro, I. Castro-Infantes and J. Castro-Infantes.  New plane curves with curvature depending on distance from the origin. Mediterr. J. Math. (2017) 14: 108.
  * [CCIs18a] I. Castro, I. Castro-Infantes and J. Castro-Infantes.  Curves in Lorentz-Minkowski plane: elasticae, catenaries and grim-reapers. Open Math. 16 (2018), 747–766.
  * [CCIs18b] I. Castro, I. Castro-Infantes and J. Castro-Infantes. Curves in Lorentz- Minkowski plane with curvature depending on their position. Open Math. 18 (2020), 749–770.
  * [CCIs21] I. Castro, I. Castro-Infantes and J. Castro-Infantes, Helicoidal minimal surfaces in the 3-sphere, Preprint (2021).
  * [E00] P. Erdös.  Spiraling the Earth with C. G. J. Jacobi. Amer. J. Phys. 68 (2000), 888–895.
  * [F93] R. Ferréol.  Encyclopédie des formes mathématiques remarquables. www. mathcurve. com
  * [J95] V. Jurdjevic.  Non-Euclidean elastica. Amer. J. Math. 117 (1995), 93–124.
  * [LS84] J. Langer and D. Singer.  The total squared curvature of closed curves. J. Diff. Geom. 20 (1984), 1–22.
  * [LS96] J. Langer and D. Singer.  Lagrangian aspects of the Kirchhoff elastic rod. SIAM Rev. 38 (1996), 1–17.
  * [S99] D. Singer.  Curves whose curvature depends on distance from the origin. Amer. Math. Monthly 106 (1999), 835–841.
  * [S08] D. Singer.  Lectures on elastic curves and rods. Curvature and Variational Modelling in Physics and Biophysics, AIP Conf. Proc. vol. 1002 (2008), 3–32.