\Large\sf Discrete sine-Gordon Dynamics on networks\bigskip
Denys \textsc{Dutykh
In this study we consider the sine-Gordon equation formulated on domains which are not locally homeomorphic to any subset of the Euclidean space. More precisely, we formulate the discrete dynamics on trees and graphs. Each edge is assumed to be a 1D uniform lattice with end points identified with graph vertices. A special treatment is needed at the junctions in order to couple 1D lattices into a global communicating network. Our approach is based on considering the local conservation properties. Some preliminary numerical results are shown on a simple graph containing four loops. These results show the performance of the scheme in non-trivial realistic conditions. \bigskip Keywords: sine-Gordon; networks; kink solution
SG_Trees.tex
\RequirePackage{fix-cm} % Fix LaTeX2e bugs.
\documentclass[a4paper, oneside, reqno, dvips, 12pt]{article}
\usepackage{fixltx2e}   % Fix LaTeX2e bugs.

%\usepackage[english]{babel}
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}

\usepackage{eucal}
\usepackage{xspace}
\usepackage{amsgen}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{upgreek}
\usepackage{amsfonts}
\usepackage{MnSymbol}
\usepackage{mathtools}
\usepackage[nice]{nicefrac}
\usepackage{mathrsfs, dsfont}
\usepackage{esint}

\usepackage{a4wide}
\footskip=15mm \textwidth=150mm \textheight=220mm

\usepackage[numbers]{natbib}

%%% Definition of some colors:
\usepackage[table]{xcolor}
\usepackage{dcolumn}  % Align table columns on decimal point

\definecolor{oneblue}{rgb}{0.0, 0.0, 0.85}

\usepackage{psfrag}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{morefloats}
\usepackage{indentfirst}

\usepackage{acronym}
\usepackage{microtype}
\usepackage[labelfont={bf,sf},%
            labelsep=period,%
            justification=raggedright]{caption}

\usepackage[usenames, dvipsnames, pdf]{pstricks}
\usepackage{epsfig}
\usepackage{pst-grad} % For gradients
\usepackage{pst-plot} % For axes

%%% Hyperref links:
%\usepackage[colorlinks,
%            urlcolor=oneblue,
%            linkcolor=oneblue,
%            citecolor=oneblue,
%            bookmarksopen=false,
%            pagebackref]{hyperref}
%\hypersetup{pdfpagemode=UseNone}

\usepackage{titletoc}
\usepackage{titlesec}
\titleformat{\section}{\large\bfseries\scshape}{\thesection.}{0.5em}{}

%%% Other definitions:
\vfuzz2pt % Don't report over-full v-boxes if over-edge is small
\hfuzz2pt % Don't report over-full h-boxes if over-edge is small

\numberwithin{equation}{section}

\newtheorem{remark}{Remark}

\renewcommand{\baselinestretch}{1.05}

\acrodef{sg}[sG]{sine-Gordon}
\acrodef{nswe}[NSWE]{Nonlinear Shallow Water Equations}

\newcommand{\up}[1]{$\,^{\mathrm{\small\textsf{#1}}}$} % \up command from French babel

%%% symbols
\newcommand{\depth}{d}
\newcommand{\upi}{{\pi}}
\newcommand{\s}{{\sf s}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\I}{\mathbb{I}}
\newcommand{\M}{\mathbb{M}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Ee}{\mathbb{E}}
\newcommand{\D}{\mathcal{D}}
\newcommand{\E}{\mathcal{E}}
\newcommand{\depthcm}{\dbar}
\newcommand{\ud}{\mathrm{d}}
\newcommand{\ui}{\mathrm{i}}
\newcommand{\ue}{\mathrm{e}}
\newcommand{\Tt}{\mathbb{T}}
\newcommand{\T}{\mathcal{T}}
\newcommand{\Y}{\mathcal{Y}}
\newcommand{\Id}{\mathbb{I}}
\newcommand{\In}{\mathcal{I}}
\newcommand{\Ll}{\mathscr{L}}
\newcommand{\Nn}{\mathscr{N}}
\renewcommand{\u}{\mathbf{u}}
\newcommand{\Mo}{\mathcal{M}}
\renewcommand{\P}{\mathcal{P}}
\newcommand{\eps}{\varepsilon}
\renewcommand{\O}{\mathcal{O}}
\renewcommand{\H}{\mathcal{H}}
\renewcommand{\L}{\mathcal{L}}
\renewcommand{\S}{\mathcal{S}}
\newcommand{\n}{\boldsymbol{n}}
\newcommand{\x}{\boldsymbol{x}}
\renewcommand{\alpha}{\upalpha}
\newcommand{\Mat}{\mathrm{Mat}}
\newcommand{\z}{{\boldsymbol{z}}}
\newcommand{\J}{{\boldsymbol{J}}}
\renewcommand{\a}{{\boldsymbol{a}}}
\renewcommand{\b}{{\boldsymbol{b}}}
\renewcommand{\v}{{\boldsymbol{v}}}
\newcommand{\btau}{\boldsymbol{\tau}}
\newcommand{\vd}{{\boldsymbol{\delta}}}

% Comparisons:
\renewcommand{\geq}{\geqslant}
\renewcommand{\leq}{\leqslant}

% Operators
\newcommand{\wt}[1]{\widetilde{#1}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}

%%% math operators
\newcommand{\Div}{\grad\scal}
\newcommand{\scal}{\boldsymbol{\cdot}}
\newcommand{\grad}{\boldsymbol{\nabla}}
\newcommand{\diag}{\mathop{\mathrm{diag}}}
\newcommand{\sech}{\mathop{\mathrm{sech}}}
\newcommand{\sign}{\mathop{\mathrm{sign}}}
\newcommand{\pd}[2]{\frac{\partial\, #1}{\partial\/ #2}}
\newcommand{\od}[2]{\frac{\mathrm{d}\,#1}{\mathrm{d}\/#2}}

%%% rational numbers
\newcommand{\half}{{\textstyle{1\over2}}}
\newcommand{\third}{{\textstyle{1\over3}}}
\newcommand{\eighth}{{\textstyle{1\over8}}}
\newcommand{\twothird}{{\textstyle{2\over3}}}
\newcommand{\halfi}{{\textstyle{1\over2\ui}}}
\newcommand{\threehalf}{{\textstyle{3\over2}}}
\newcommand{\threeforth}{{\textstyle{3\over4}}}
\newcommand{\sixteenth}{{\textstyle{1\over16}}}
\newcommand{\threeeighth}{{\textstyle{3\over8}}}

\newcommand{\ie}{\emph{i.e.}~}
\newcommand{\eg}{\emph{e.g.}~}
\newcommand{\etc}{\emph{etc.}~}
\newcommand{\etal}{\emph{et al.}~}

\title{\Large\sf Discrete sine-Gordon Dynamics on networks\bigskip}

\author{
        Denys \textsc{Dutykh} \\
        {\small CNRS--LAMA UMR 5127} \\
        {\small Universit\'e Savoie Mont Blanc} \\
        {\small 73376 Le Bourget-du-Lac France} \and
        Jean-Guy \textsc{Caputo} \\
        {\small Laboratoire de Math\'ematiques} \\
        {\small INSA de Rouen} \\
        {\small 76801 Saint-Etienne du Rouvray France}}

\date{\bigskip\today}

\begin{document}

\maketitle

\begin{abstract}
In this study we consider the sine-Gordon equation formulated on domains which are not locally homeomorphic to any subset of the Euclidean space. More precisely, we formulate the discrete dynamics on trees and graphs. Each edge is assumed to be a 1D uniform lattice with end points identified with graph vertices. A special treatment is needed at the junctions in order to couple 1D lattices into a global communicating network. Our approach is based on considering the local conservation properties. Some preliminary numerical results are shown on a simple graph containing four loops. These results show the performance of the scheme in non-trivial realistic conditions.

\bigskip
\textbf{Keywords:} sine-Gordon; networks; kink solution
\end{abstract}

\newpage
\tableofcontents
\clearpage

\section{Introduction}

Currently there is a growing demand for modelling and understanding of various flow problems on networks. A generic network can be described as follows. Imagine a (usually finite) set of points or simply \emph{vertices} immersed in an Euclidean space $\Ee^2$ or $\Ee^3$ (depending on the application in hands). Some of the points are connected by 1D segments (or more generally curves, which are homeomorphic to segments), which are called the \emph{edges}. Mathematically the networks are formalized using the graph theory \cite{Gould2012}. However, there is an important subtle difference with the graph theory. Namely, in some applications the geometry of edges (\eg their length, shape, thickness) can matter, while in the graph theory the only relevant information is the fact that two points are connected by an edge. Such sensitive applications include, for example, the blood flow modelling \cite{Vassilevskii2011}. Thus, a network combines in a single data structure the corresponding geometrical and topological information on vertices and edges. The flow is usually modelled with Partial Differential Equations (PDEs), which come to replace the methods from discrete mathematics or Ordinary Differential Equations (ODEs). For a general recent review of this topic we refer to \cite{Bressan2014}. This field continues to attract more and more researchers from modelling, analysis, numerics, optimization and control theory \cite{Zuazua2013}.

One of the main difficulties of formulating evolution problems described by PDEs on networks (\ie graphs) consists in the fact that these objects are non-manifolds. Nowadays we can say that the formulation of Hamiltonian mechanics on manifolds does not pose any serious technical difficulties \cite{Arnold1997}. However, the modeling of various processes on networks (such as electric circuits, blood arteries, water-pipe supply) necessitates the generalization of the classical mechanics to non-manifolds like graphs and trees. Consider, for example, a $Y$- or $T$--junction. This domain is a \emph{semi-algebraic set}, but not a \emph{manifold}. The difficulty comes from the branching point whose neighbourhood is not \emph{homeomorphic} to any Euclidean space $\Ee^k$. It goes without saying that any composition of $Y$--junctions into, for example, a complex tree will not be a manifold either.

To our knowledge, the Hamiltonian problems on non-manifolds have not been systematically studied. A notorious exception is the work \cite{Bibikov2009}, where the wave scattering in the Klein--Gordon equation was investigated on a domain consisting of three semi-infinite straight lines having one common point. We can also mention the publication \cite{Bona2008} where the BBM equation was considered on a tree. However, their study does not rely on any particular variational structure of the governing equation.

In the present study we will consider the celebrated \acf{sg} equation which is Hamiltonian and integrable PDE \cite{Takhtadzhyan1974, Dashen1974}. However, the integrability of the \acs{sg} equation is not compulsory for our purposes. In the developments presented below we will use the Hamiltonian and Lagrangian structures to determine the relevant conservative quantities and construct an appropriate symplectic discretization. In the present study, we will consider the discrete dynamics of \acs{sg} on lattices composed into graphs. The transition rules between adjacent lattices at junction points follow from the discretization of local conservation laws. This approach was already used in \cite{Bibikov2009, Caputo2014}. What we hope is that the limit as the lattice parameter $\Delta\to 0$ will provide us with the continuous version of the Hamiltonian mechanics on non-manifolds.

The present study is organized as follows. In the following Section~\ref{sec:sg} we present some basic facts on the \acs{sg} equation. In Section~\ref{sec:discr} the discrete version of the \acs{sg} equation is constructed. Then, the formulation of the (discrete) \acs{sg} equation on graphs (non-manifolds) is presented in Section~\ref{sec:graphs}. Some preliminary numerical results on graphs are shown in Section~\ref{sec:num}. Finally, the main conclusions of this study are summarized in Section~\ref{sec:concl}.

\section{Continuous sine-Gordon equation}\label{sec:sg}

Consider the real space-time coordinates $(x,t) \in \R\times\R_+$. Then, the most common version of the \acf{sg} equation reads \cite{Scott2003}
\begin{equation}\label{eq:sg}
  u_{tt}\ -\ u_{xx}\ +\ \sin u \ =\ 0,
\end{equation}
where the subscripts $(\cdot)_{t}$, $(\cdot)_{x}$ denote the derivatives with respect to the time $t$ and space $x$ coordinates. In order to obtain a well-posed boundary value problem, equation \eqref{eq:sg} is completed by periodic or homogeneous Neumann boundary conditions. The linear part $\Box^2 u := u_{tt}\ -\ u_{xx}$ is the D'Alembertian or Laplacian in Minkowski space $\M^2$. The \acs{sg} equation is known to be Lorentz invariant and an integrable PDE \cite{Faddeev1987}.

The name for this equation comes presumably from M.~\textsc{Kruskal} by analogy with the linear Klein--Gordon equation:
\begin{equation}\label{eq:kg}
  u_{tt}\ -\ u_{xx}\ +\ u \ =\ 0,
\end{equation}
which is the linearization of \eqref{eq:sg}. The Klein--Gordon equation \eqref{eq:sg} in its turn should be rather called the Klein--Fock equation, see for more details \cite{Fock1926}.

Changing to the characteristic (\emph{light cone}) coordinates $(\xi, \zeta)$, the \acs{sg} equation can be written in another equivalent form
\begin{equation*}
  u_{\xi\zeta}\ =\ \sin u, \quad \xi\ :=\ \frac{x\ +\ t}{2}, \quad \zeta\ :=\ \frac{x\ -\ t}{2}.
\end{equation*}
However, in the present study we will privelege the inital formulation \eqref{eq:sg}.

\subsection{Variational structure}

The \acs{sg} equation can be derived as the Euler--Lagrange equation of the following Lagrangian density
\begin{equation*}
  \L_{\mathrm{sG}}\ :=\ \half\bigl(u_t^2\ -\ u_x^2\bigr)\ -\ 1\ +\ \cos u.
\end{equation*}
By using the Taylor expansion for the potential energy term, the Lagrangian density $\L_{sg}$ can be rewritten as
\begin{equation*}
  \L_{\mathrm{sG}}\ =\ \underbrace{\half\bigl(u_t^2\ -\ u_x^2\bigr)\ -\ \half u^2}_{\L_{\mathrm{KG}}}\ +\ \sum_{k=2}^{\infty}(-1)^{k}\frac{u^{2k}}{(2k)!}\ =\ \L_{\mathrm{KG}} \ +\ \sum_{k=2}^{\infty}(-1)^{k}\frac{u^{2k}}{(2k)!}.
\end{equation*}
The last form of the Lagrangian density is particularly suitable for the complexification of the \acs{sg} equation which derives from the following Lagrangian
\begin{equation}\label{eq:Lc}
  \L^{c}_{\mathrm{sG}}\ =\ \half\bigl(u_t\, u^*_t\ -\ u_x\, u_x^*\bigr)\ -\ \half u\, u^*\ +\ \sum_{k=2}^{\infty}(-1)^{k}\frac{(u\, u^*)^k}{(2k)!},
\end{equation}
where $u^*(x,t)$ is the complex conjugate of $u(x,t)$ (in the complex-valued version of the \acs{sg} equation).

Moreover, the \acs{sg} equation possesses also the Hamiltonian formulation
\begin{equation}\label{eq:ham}
  \z_t\ =\ \J\vd_\z\H, \qquad \J\ :=\ 
  \begin{pmatrix}
     0 & 1 \\
    -1 & 0
   \end{pmatrix},
\end{equation}
$\z := (u,v)$, $\vd_\z := (\delta_u, \delta_v)$ is the variational grafient, $\J$ is the symplectic operator and, finally, the Hamiltonian functional $\H$ is defined as
\begin{equation*}
  \H[\z]\ :=\ \int_{-\infty}^{+\infty}\Bigl[\,\half v^2\ +\ \half u_x^2\ +\ 1\ -\ \cos u\,\Bigr]\, \ud x.
\end{equation*}
Equations \eqref{eq:ham} can be rewritten in component-wise form for the sake of clearness
\begin{eqnarray*}
  u_t &=& \frac{\delta\H}{\delta v} = v, \\
  v_t &=& -\frac{\delta\H}{\delta u} = u_{xx}\ -\ \sin u.
\end{eqnarray*}
Consequently, the \acs{sg} equation can be viewed as a Hamiltonian system with phase space $(u,v)$ and the symplectic form
\begin{equation}\label{eq:sympl}
  \omega\ :=\ \int_{-\infty}^{+\infty} \ud u \wedge \ud v\, \ud x.
\end{equation}

The Hamiltonian $\H[u,v] \equiv \E(u)$ is a conserved quantity and it has the sense of the physical energy
\begin{equation*}
  \E(u)\ :=\ \int_{-\infty}^{+\infty}\Bigl[\,\half u_t^2\ +\ \half u_x^2\ +\ 1\ -\ \cos u\,\Bigr]\, \ud x.
\end{equation*}

\begin{remark}
One can remark also that there is another important functional conserved for the equation \eqref{eq:sg} which can be associated to the total momentum
\begin{equation*}
  \Mo[\z]\ :=\ \int_{-\infty}^{+\infty} vu_x\, \ud x.
\end{equation*}
The conservation of $\Mo[\z]$ can be readily checked by computing $\od{\Mo}{t}$.
\end{remark}

\subsection{Exact solutions}

The \acs{sg} equation has the following kink-type solution:
\begin{equation}\label{eq:kink0}
  u(x,t)\ =\ 4\arctan \ue^{\gamma(x + x_0 - ct)}, \qquad x_0\in\R,
\end{equation}
where $c\in[0,1)$ is the kink celerity and $\gamma$ is the so-called Lorentz factor
\begin{equation}\label{eq:kink}
  \gamma^2\ :=\ \frac{1}{1 - c^2}.
\end{equation}
The energy of the kink can be computed analytically as well $\E = 8\gamma$. In our scaling the speed of light is equal to $1$. The kink solution realizes a smooth transition between two ground states $0\ \leadsto\ 2\pi$. More generally, the kinks can link in the phase space two neighbouring states $2\pi k\ \leadsto\ 2\pi(k \pm 1)$.

There exist many other exact analytical solutions to the \acs{sg} equation \cite{Scott2003}. The most important one, besides the kink, is perhaps the breather which is out of scope of the present study.

\begin{figure}
  \centering
  \includegraphics[width=0.69\textwidth]{figs/Kink2.eps}
  \caption{\small\em The analytical kink solution \eqref{eq:kink0} to the sine-Gordon equation \eqref{eq:sg} for two values of the speed parameter $c = 0.5$ and $c = 0.95$.}
  \label{fig:kink}
\end{figure}

\section{Discrete sine-Gordon equation}\label{sec:discr}

In order to propose a discrete version of the \acs{sg} equation we will follow the variational framework. We showed above that the \acs{sg} equation is a Hamiltonian PDE. A natural way to convert it into a discrete dynamical system is to employ a symplectic discretization \cite{Leimkuhler2004}. The workflow is determined by the method of lines:
\begin{itemize}
  \item Discretize the Hamiltonian functional $\H[\z]$ in space on a lattice to obtain a system of coupled Hamiltonian ODEs.
  \item Discretization in time of the system of Hamiltonian ODEs using a symplectic scheme.
\end{itemize}
This programme will be realized below by following the gret lines of \cite{Leimkuhler2004}.

Consider a uniform lattice $\{x_j = k\Delta x\, |\, j=1,\ldots,n\}$, $\Delta x > 0$. The values of the \acs{sg} solution $u(x,t)$ at lattice points will be denoted by $u_j \approx u(x_j)$. For the moment we will consider only interior nodes. The junction points (end points of the lattice) will be discussed below. After the discretization, the phase space becomes finite dimensional, since $\{\z_j\}_{j=1}^{n} = \{(u_j, v_j)\in\R^2\}_{j=1}^{n}\in\R^{2n}$. The discrete symplectic form on this space becomes
\begin{equation*}
  \omega_n\ =\ \sum_{j=1}^{n}\ud u_j\wedge \ud v_j\, \Delta x,
\end{equation*}
which is a straightforward discretization of \eqref{eq:sympl}. The Hamiltonian functional will be approximated with the \emph{rectangular rule} as the following sum
\begin{equation*}
  \H_n[\{u_j,v_j\}]\ =\ \sum_{j=1}^{n}\biggl[\half v_j^2\ +\ \half\Bigl(\frac{u_j - u_{j-1}}{\Delta x}\Bigr)^2\ +\ 1\ -\ \cos u_j \biggr]\,\Delta x.
\end{equation*}
The system of Hamiltonian ODEs follows automatically
\begin{equation*}
  \od{\z_j}{t}\ =\ \J_n\cdot\grad_{\z_j}\H_n[\z], \qquad
  \J_n = \begin{pmatrix}
           \boldsymbol{0} & \I_n \\
           -\I_n & \boldsymbol{0}
         \end{pmatrix}, \qquad i=1,\ldots,n.
\end{equation*}
After computing the derivatives, the last semi-discrete system becomes
\begin{eqnarray}\label{eq:semid1}
  \od{u_j}{t} &=& v_j, \\
  \od{v_j}{t} &=& \frac{w_{j+1}\ -\ w_j}{\Delta x}\ -\ \sin u_j, \label{eq:semid2}
\end{eqnarray}
where $w_{j+1}(t) := \frac{u_{j+1} - u_j}{\Delta x}$ and $w_j(t) := \frac{u_j - u_{j-1}}{\Delta x}$. It can be shown \cite{Leimkuhler2004} that the semi-discrete scheme \eqref{eq:semid1}, \eqref{eq:semid2} satisfies a local energy conservation law
\begin{equation}\label{eq:ed}
  \od{}{t}\;\Bigl[\half v_j^2\ +\ \half w_j^2\ +\ 1\ -\ \cos u_j\Bigr]\ +\ \frac{F_{j+\frac12}\ -\ F_{j-\frac12}}{\Delta x}\ =\ 0,
\end{equation}
where the quantity in brackets $\E_j := \half v_j^2\ +\ \half w_j^2\ +\ 1\ -\ \cos u_j$ is the (semi-)discrete energy and $F_{j+\frac12} := -v_jw_{j+1}$, $F_{j-\frac12} := -v_{j-1}w_j$ are the energy fluxes.

In order to obtain a fully discrete scheme, the system of ODEs \eqref{eq:semid1}, \eqref{eq:semid2} can be discretized in time with a \emph{symplectic Euler} method, for example:
\begin{eqnarray*}
  u_j^{m+1} &=& u_j^m\ +\ \Delta t\, v_j^{m+1}, \\
  v_j^{m+1} &=& v_j^m\ +\ \Delta t\, \biggl[\frac{w_{j+1}^m\ -\ w_j^m}{\Delta x}\ -\ \sin u_j^m\biggr],
\end{eqnarray*}
where $\Delta t > 0$ and $u_j^m := u(x_j, t_m)$, $t_m := m\Delta t$, $m = 1,2,\ldots$ After eliminating $v_j^n$ and $w_j^n$ from these equations we obtain the classical leap-frog scheme as a fully discrete anlog of the \acs{sg} equation
\begin{equation*}
  \frac{u_{j}^{m+1}\ -\ 2u_{j}^{m}\ +\ u_{j}^{m-1}}{\Delta t^2}\ -\ \frac{u_{j+1}^{m}\ -\ 2u_{j}^{m}\ +\ u_{j-1}^{m}}{\Delta x^2}\ +\ \sin u_j^m\ =\ 0.
\end{equation*}
After simple algebraic manipulations we can obtain the following discrete dynamical system for interior nodes of the lattice:
\begin{equation*}
  u_{j}^{m+1}\ =\ 2u_{j}^{m}\ -\ u_{j}^{m-1}\ +\ \Bigl(\frac{\Delta t}{\Delta x}\Bigr)^2\bigl[u_{j+1}^{m}\ -\ 2u_{j}^{m}\ +\ u_{j-1}^{m}\bigr]\ -\ \Delta t^2\,\sin u_j^m, \quad m\ =\ 1, 2, \ldots.
\end{equation*}
The treatment of nodes at junctions will be discussed below.

In general, one cannot expect to have the fully discrete energy conservation law similar to \eqref{eq:ed}, in contrast to the semi-discrete schemes. The reason is that a symplectic scheme aims to preserve the symplectic form and it does not guarantee anything about the Hamiltonian. However, the backward error analysis explains why, in general, the symplectic discretizations of PDEs show satisfactory energy conservation properties \cite{Moore2003}.

\section{Formulation on graphs}\label{sec:graphs}

In this Section we will describe the assembling procedure of 1D lattices $\ell_i := \bigl\{\x_j\in[\a_i = \x_1,\, \b_i = \x_{n_i}]\, |\, j=1\ldots,n_i\bigr\}$ into a network whose mathematical description is usually given on the language of the graph theory. In the sequel we will denote by $\a_i$, $\b_i$ the starting and terminal points of the lattice $\ell_i$ respectively.

Consider a simple oriented network-shaped weighted graph $G := (V,E)$. The finite set of vertices $V$ is basically the union of lattice initial and terminal nodes:
\begin{equation*}
  V\ :=\ \{\v_j\}_{j=1}^{m}\ \equiv\ \bigcup\limits_{i=1}^{|E|}\ \{\a_i, \b_i\}.
\end{equation*}
The finite set of edges $E \subseteq \Bigl\{(\a, \b) \in V^2\ |\ \a \neq \b \Bigr\}$. Every edge $e_i := (\ell_i,\, \omega_i)$ consists of 1D lattice segments whose orientation is naturally determined by the enumeration of discrete lattice points (or equivalently the choice of the first and last points $\a_i$ and $\b_i$). The \emph{length} of the edge $\ell_i$ can be prescribed through its weight $\omega_i$:
\begin{equation*}
  \omega_i\ :=\ |\ell_i|\Delta x_i\ \equiv \ n_i, \qquad i\ =\ 1,\ldots, |E|.
\end{equation*}
where $|\ell_i|$ is the number of points in the lattice and $\Delta x_i$ is the spacing between two consecutive points.

The internal organization of the network $G$ is traditionally given by graph-theoretical data structures such as the incidence and adjacency matrices \cite{Gould2012}. In the present study we will privelege the incidence matrix representation. By definition, the incidence matrix $\In = \bigl(b_{ij}\bigr)_{n\times m}\in\Mat_{n\times m}(\Z)$, $n := |V|$, $m := |E|$ has the following elements
\begin{equation*}
  b_{ij}\ =\ \begin{dcases}
  +1, & \mbox{ if the edge } e_j \mbox{ enters the vertex } v_i, \\
  -1, & \mbox{ if the edge } e_j \mbox{ leaves the vertex } v_i, \\
  0, & \mbox { otherwise.}
  \end{dcases}.
\end{equation*}
For the sake of illustration let us consider the graph $G_0 = (V,E)$ represented on Figure~\ref{fig:merc}. It is composed of four vertices $V = \{v_1, v_2, v_3, v_4\}$ and six edges $E = \{e_1,\ldots,e_6\}$. It is straightforward to check that the incidence matrix $\In_0$ of the graph $G_0$ is
\begin{equation*}
  \In_0\ =\ 
  \begin{pmatrix}
    -1 &  0 &  0 &  0 &  1 & -1 \\
     1 & -1 & -1 &  0 &  0 &  0 \\
     0 &  1 &  0 & -1 &  0 &  1 \\
     0 &  0 &  1 &  1 & -1 &  0
  \end{pmatrix}.
\end{equation*}
However, the most important information for us is the correspondance between the vertices $v_i$ with starting/terminal points of the lattices which compose the edges $e_j$. This correspondance is given as a list:
\begin{eqnarray*}
  v_1\ =\ \{\a_1, \a_6, \b_5\}, & v_3\ =\ \{\a_4, \b_2, \b_6\}, \\
  v_2\ =\ \{\a_2, \a_3, \b_1\}, & v_4\ =\ \{\a_5, \b_3, \b_4\}.
\end{eqnarray*}

\begin{figure}
  \centering
  % Generated with LaTeXDraw 2.0.8
  % Fri Aug 01 14:38:51 CEST 2014
  % \usepackage[usenames,dvipsnames]{pstricks}
  % \usepackage{epsfig}
  % \usepackage{pst-grad} % For gradients
  % \usepackage{pst-plot} % For axes
  \scalebox{0.85} % Change this value to rescale the drawing.
  {
  \begin{pspicture}(0,-5.088125)(9.321875,5.088125)
  \definecolor{color1}{rgb}{0.3215686274509804,0.3215686274509804,0.3215686274509804}
  \definecolor{color162b}{rgb}{0.21176470588235294,0.21176470588235294,0.21176470588235294}
  \psline[linewidth=0.051999997cm,linecolor=color1,arrowsize=0.193cm 2.0,arrowlength=1.4,arrowinset=0.4]{>-}(4.48,-2.4503126)(4.48,-0.2703125)
  \psline[linewidth=0.06cm,linecolor=color1,arrowsize=0.193cm 2.0,arrowlength=1.4,arrowinset=0.4]{->}(4.48,-0.3103125)(6.366702,1.383194)
  \psline[linewidth=0.06cm,linecolor=color1,arrowsize=0.193cm 2.0,arrowlength=1.4,arrowinset=0.4]{->}(4.48,-0.2903125)(2.6952674,1.3510033)
  \psline[linewidth=0.06cm,linecolor=color1](6.2889,1.3204715)(8.019998,2.9094405)
  \psline[linewidth=0.06cm,linecolor=color1](2.78,1.2896875)(1.0597715,2.7896874)
  \psline[linewidth=0.06cm,linecolor=color1](4.48,-2.2703125)(4.5,-4.5103126)
  \pscircle[linewidth=0.06,linecolor=color1,dimen=outer](4.55,0.0){4.55}
  \pscircle[linewidth=0.06,linecolor=color162b,dimen=outer,fillstyle=solid,fillcolor=color162b](4.49,-0.2603125){0.13}
  \pscircle[linewidth=0.06,linecolor=color162b,dimen=outer,fillstyle=solid,fillcolor=color162b](8.03,2.8996875){0.13}
  \pscircle[linewidth=0.06,linecolor=color162b,dimen=outer,fillstyle=solid,fillcolor=color162b](4.51,-4.4803123){0.13}
  \pscircle[linewidth=0.06,linecolor=color162b,dimen=outer,fillstyle=solid,fillcolor=color162b](1.05,2.7996874){0.13}
  \psline[linewidth=0.06cm,linecolor=color1](4.28,4.5296874)(4.82,4.6096873)
  \psline[linewidth=0.06cm,linecolor=color1](4.28,4.5096874)(4.82,4.4296875)
  \psline[linewidth=0.06cm,linecolor=color1](1.1,-2.9103124)(0.98,-2.5903125)
  \psline[linewidth=0.06cm,linecolor=color1](1.12,-2.9103124)(0.82,-2.7303126)
  \psline[linewidth=0.06cm,linecolor=color1](8.34,-2.4503126)(8.08,-2.6703124)
  \psline[linewidth=0.06cm,linecolor=color1](8.34,-2.4503126)(8.24,-2.7903125)
  \usefont{T1}{ptm}{m}{n}
  \rput(4.5514064,-4.8653126){$v_1$}
  \usefont{T1}{ptm}{m}{n}
  \rput(5.0314064,-0.4653125){$v_2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(8.511406,3.1946876){$v_3$}
  \usefont{T1}{ptm}{m}{n}
  \rput(0.83140624,3.0746875){$v_4$}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.981406,-2.6853125){$e_1$}
  \usefont{T1}{ptm}{m}{n}
  \rput(8.841406,-2.8453126){$e_6$}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.541406,4.8946877){$e_4$}
  \usefont{T1}{ptm}{m}{n}
  \rput(6.621406,0.9346875){$e_2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(2.4814062,1.1546875){$e_3$}
  \usefont{T1}{ptm}{m}{n}
  \rput(0.52140623,-2.8453126){$e_5$}
  \end{pspicture} 
  }
  \caption{\small\em A sample graph used in our study for the sake of illustration.}
  \label{fig:merc}
\end{figure}

\subsection{Conditions on junctions: the continuous case}

To fix the notations, consider a $Y$--junction composed of three semi-infinite strings $\S_{1,2,3}$ embedded into the Euclidean space $\R^2$:
\begin{equation*}
  \Y\ :=\ \Bigl\{\x\in\R^2: \exists i\in\{1,2,3\} \mbox{ such that } \x\in\S_i\Bigr\}
\end{equation*}
The intersection of all three strings is located at the unique point $C\in\R^2$ defined as (see Figure~\ref{fig:sketch1} for the illustration):
\begin{equation*}
  C\ :=\ \S_1\cap\S_2\cap\S_3.
\end{equation*}

A simple topological argument can be applied to show that the set $\Y$ is not homeomorphic to any Euclidean space. Indeed, let us remove virtually the junction point $C$ from this set $\Y$. It will be decomposed in three disjoint components. Clearly, any Euclidean space $\R^n$, $n \geq 1$ does not have such a point with the same property.

The first natural condition to be satisfied at the junction point $C$ is the continuity of the solution
\begin{equation}\label{eq:cont}
  \lim\limits_{\x\to C,\; \x\in\S_1} u(x,t)\ =\ \lim\limits_{\x\to C,\; \x\in\S_2} u(x,t)\ =\ \lim\limits_{\x\to C,\; \x\in\S_3} u(x,t).
\end{equation}
Moreover, this condition has to be completed by the ``\emph{charge}'' conservation property adopted also in some previous studies of Klein--Gordon \cite{Mehmeti2003, Bibikov2009} and \acs{sg} \cite{Caputo2014} equations studies:
\begin{equation}\label{eq:kirch}
  \left.\partial^1 u \right|_{\x\to C,\; \x\in\S_1}\ +\ \left.\partial^2 u \right|_{\x\to C,\; \x\in\S_2}\ +\ \left.\partial^3 u \right|_{\x\to C,\; \x\in\S_3}\ =\ 0,
\end{equation}
where $\partial^i$ denotes the first spatial derivative along the branch $\S_i$, $i=1,2,3$. Condition \eqref{eq:kirch} is a continuous analogue of the celebrated Kirchhoff's circuit law.

\subsubsection{Phase space structure}%\label{}

It can be beneficial for our understanding and analysis of the \acs{sg} dynamics on the graphs to look at the structure of the phase space. For this purpose let us consider for simplicity the $Y$-junction represented in the center of Figure~\ref{fig:merc}. The augmented phase space for the unsteady \acs{sg} equation restricted to a branch (say $e_1$) will consist of $\P\times \In$, where $\P := (\ddot{u}, \dot{u}, u, u_x, u_{xx})$ is the usual configuration space and $\In$ is the spatial extent of a branch ($\In = [v_1, v_2]$ in this particular case). We have to augment the configuration space $\P$ by $\In$ since the joints are realized on the boundaries of the interval $\In$. This portion of the phase space is schematically represented on Figure~\ref{fig:phase}. In Geometry these structures are called the \emph{foliations} \cite{Haefliger1970}. We stress out that it is only the local structure of the global phase space corresponding to the graph shown on Figure~\ref{fig:merc} which is represented. In order to represent the global structure, the borders of the phase spaces of individual branches (leaves on Figure~\ref{fig:phase}) have to be glued together in accordance with the scheme prescribed by the graph $G$.

\begin{figure}
  \centering
  \includegraphics[width=0.75\textwidth]{figs/PhaseSpace.eps}
  \caption{\small\em A schematic representation of the phase spaces foliation at a junction point.}
  \label{fig:phase}
\end{figure}

\subsubsection{Back-to-Manifold approach}\label{sec:man}

Condition \eqref{eq:kirch} can be demonstrated by converting (``\emph{inflating}'') the $Y$-junction domain into a manifold $\Y_\delta$ of small thickness $\delta > 0$. Thus, $\Y_\delta$ becomes a tubular neighbourhood of our network. The sketch of the proof is shown on Figure~\ref{fig:sketch2}. On the boundary $\partial\Y_{\delta}$ we impose the homogeneous Neuman condition \cite{Gulevich2006, Caputo2014}:
\begin{equation}\label{eq:neuman}
  \left.\partial_n u\,\right|_{\x\in\partial\Y_\delta}\ =\ 0, \qquad
  \partial_n u\ :=\ \grad u\cdot \n,
\end{equation}

The main idea of demonstrating \eqref{eq:kirch} consists in integrating the 2D version of the \acs{sg} equation \ref{eq:sg} over the domain $\Omega_\delta\cap\Y_\delta$\footnote{The domain $\Omega_\delta := \D_{2\delta}(C)$ is a disc centered at the junction point $C$ with radius $2\delta$ depicted on Figure~\ref{fig:sketch2}.} and using the Stokes theorem in the divergence (Laplacian) term:
\begin{equation*}
  \iint\limits_{\Omega_\delta\cap\Y_{\delta}}\Bigl[u_{tt}\ +\ \grad^2 u\ +\ \sin u\Bigr]\,\ud\x\ =\ \underbrace{\iint\limits_{\Omega_\delta\cap\Y_{\delta}}\Bigl[u_{tt}\ +\ \sin u\Bigr]\,\ud\x}_{(i)} + \underbrace{\int\limits_{\partial(\Omega_\delta\cap\Y_{\delta})}\partial_n u\,\ud s}_{(ii)}\ =\ 0.
\end{equation*}
By assuming the solution $u(\x,t)$ to be smooth and bounded in $\Omega_\delta$, the first integral $(i)$ scales with $\O(\delta^2)$. After taking into account the boundary condition \eqref{eq:neuman}, the second integral $(ii)$ becomes just the sum of three line integrals over $\ell_i$ in the interior of the domain $\Y_\delta$ (represented with red color on Figure~\ref{fig:sketch2}):
\begin{equation*}
  \sum_{i=1}^3\underbrace{\int_{\ell_i}\partial_n u\,\ud s}_{\O(\delta)}\ +\ \O(\delta^2)\ =\ 0\quad \biggr|\times \delta^{-1}
\end{equation*}
Let us observe that under the same assumptions on the solution $u(x,t)$ the first integral scales as $\O(\delta)$. Thus, dividing the last identity by $\delta$ and taking the limit $\delta\to +0$ leads the desired result \eqref{eq:kirch}.

\begin{figure}
  \centering
  % Generated with LaTeXDraw 2.0.8
  % Mon Jul 28 15:12:52 CEST 2014
  % \usepackage[usenames,dvipsnames]{pstricks}
  % \usepackage{epsfig}
  % \usepackage{pst-grad} % For gradients
  % \usepackage{pst-plot} % For axes
  \scalebox{0.9} % Change this value to rescale the drawing.
  {
  \begin{pspicture}(0,-5.12)(8.962812,5.12)
  \definecolor{color282}{rgb}{0.3764705882352941,0.3764705882352941,0.3764705882352941}
  \psline[linewidth=0.1cm,linecolor=color282](4.3809376,0.48)(1.4009376,4.1)
  \psline[linewidth=0.1cm,linecolor=color282](4.3809376,0.5)(7.7809377,4.1)
  \psline[linewidth=0.1cm,linecolor=color282](4.3809376,0.48)(4.3809376,-4.06)
  \psline[linewidth=0.04cm,linecolor=color282,linestyle=dashed,dash=0.16cm 0.16cm](1.3609375,4.16)(0.5809375,5.1)
  \psline[linewidth=0.04cm,linecolor=color282,linestyle=dashed,dash=0.16cm 0.16cm](7.7809377,4.1)(8.580937,4.94)
  \psline[linewidth=0.04cm,linecolor=color282,linestyle=dashed,dash=0.16cm 0.16cm](4.3809376,-4.06)(4.3809376,-5.1)
  \pscircle[linewidth=0.146,linecolor=color282,dimen=outer](4.4009376,0.52){0.18}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.862344,0.185){$C$}
  \usefont{T1}{ptm}{m}{n}
  \rput(7.6723437,3.345){$\mathcal{S}_2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(1.2123437,3.585){$\mathcal{S}_3$}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.972344,-3.195){$\mathcal{S}_1$}
  \usefont{T1}{ptm}{m}{n}
  \rput(1.7623438,-2.695){$\mathbb{R}^2$}
  \pscircle[linewidth=0.146,linecolor=color282,dimen=outer](5.0109377,1.13){0.17}
  \pscircle[linewidth=0.146,linecolor=color282,dimen=outer](3.8109374,1.17){0.17}
  \pscircle[linewidth=0.146,linecolor=color282,dimen=outer](4.3909373,-0.21){0.17}
  \usefont{T1}{ptm}{m}{n}
  \rput(5.4823437,1.125){$\boldsymbol{x}_2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(3.2823439,1.185){$\boldsymbol{x}_3$}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.762344,-0.415){$\boldsymbol{x}_1$}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.362344,0.96){$\boldsymbol{x}_0$}
  \end{pspicture} 
  }
  \caption{\small\em Sketch of a $Y$--junction composed of three strings $\S_{1,2,3}$ with the center located at the point $C$.}
  \label{fig:sketch1}
\end{figure}

\begin{figure}
  \centering
  % Generated with LaTeXDraw 2.0.8
  % Sat Jul 05 06:45:01 JST 2014
  % \usepackage[usenames,dvipsnames]{pstricks}
  % \usepackage{epsfig}
  % \usepackage{pst-grad} % For gradients
  % \usepackage{pst-plot} % For axes
  \scalebox{0.85} % Change this value to rescale the drawing.
  {
  \begin{pspicture}(0,-5.426683)(10.822812,5.426269)
  \definecolor{color183}{rgb}{0.6745098039215687,0.0,0.0}
  \definecolor{color1867b}{rgb}{0.6745098039215687,0.6745098039215687,0.6745098039215687}
  \definecolor{color1868}{rgb}{0.6352941176470588,0.6352941176470588,0.6352941176470588}
  \definecolor{color1882}{rgb}{0.00392156862745098,0.00392156862745098,0.00392156862745098}
  \definecolor{color1918}{rgb}{0.6745098039215687,0.00784313725490196,0.00784313725490196}
  \psframe[linewidth=0.002,linecolor=color1867b,dimen=outer,fillstyle=solid,fillcolor=color1867b](6.2209377,0.19331692)(4.4209375,-5.426683)
  \rput{-45.0}(0.23230292,6.167463){\psframe[linewidth=0.002,linecolor=color1868,linestyle=dotted,dotsep=0.16cm,dimen=outer,fillstyle=solid,fillcolor=color1867b](8.4609375,5.613317)(6.6609373,-0.006683079)}
  \rput{45.0}(2.8846302,-1.3574793){\psframe[linewidth=0.002,linecolor=color1868,linestyle=dotted,dotsep=0.16cm,dimen=outer,fillstyle=solid,fillcolor=color1867b](3.9809375,5.613317)(2.1809375,-0.006683079)}
  \psline[linewidth=0.002,linecolor=color1867b,fillstyle=solid,fillcolor=color1867b](6.2209377,0.17331693)(4.4209375,0.17331693)(5.3209376,1.153317)(6.2209377,0.19331692)
  \usefont{T1}{ptm}{m}{n}
  \rput(9.532344,4.0783167){$\mathcal{S}_2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(4.7723436,-5.061683){$\mathcal{S}_1$}
  \usefont{T1}{ptm}{m}{n}
  \rput(1.2123437,4.118317){$\mathcal{S}_3$}
  \psline[linewidth=0.02cm,linecolor=color1882,linestyle=dashed,dash=0.16cm 0.16cm,arrowsize=0.05291667cm 2.0,arrowlength=1.4,arrowinset=0.4]{<->}(7.0009375,3.493317)(8.240937,2.2533169)
  \usefont{T1}{ptm}{m}{n}
  \rput(7.8923435,3.078317){$\delta$}
  \pscircle[linewidth=0.02,linecolor=color1882,linestyle=dashed,dash=0.16cm 0.16cm,dimen=outer](5.2809377,0.8933169){2.0}
  \pscircle[linewidth=0.022,linecolor=color1882,dimen=outer,fillstyle=solid,fillcolor=color1882](5.3409376,0.75331694){0.08}
  \usefont{T1}{ptm}{m}{n}
  \rput(5.9423437,1.0183169){$C$}
  \usefont{T1}{ptm}{m}{n}
  \rput(2.2223437,-4.041683){$\mathbb{R}^2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(7.4023438,-0.14168307){$\Omega_\delta$}
  \psline[linewidth=0.02cm,linecolor=color1882,linestyle=dashed,dash=0.16cm 0.16cm,arrowsize=0.05291667cm 2.0,arrowlength=1.4,arrowinset=0.4]{<->}(2.3609376,2.2733169)(3.5609374,3.5733168)
  \usefont{T1}{ptm}{m}{n}
  \rput(2.7123437,3.1383169){$\delta$}
  \psline[linewidth=0.02cm,linecolor=color1882,linestyle=dashed,dash=0.16cm 0.16cm,arrowsize=0.05291667cm 2.0,arrowlength=1.4,arrowinset=0.4]{<->}(4.4409375,-1.8666831)(6.2409377,-1.846683)
  \usefont{T1}{ptm}{m}{n}
  \rput(5.3323436,-2.141683){$\delta$}
  \psline[linewidth=0.02cm,linecolor=color1882,linestyle=dashed,dash=0.16cm 0.16cm,arrowsize=0.05291667cm 2.0,arrowlength=1.4,arrowinset=0.4]{->}(5.2609377,0.8333169)(3.7009375,2.113317)
  \usefont{T1}{ptm}{m}{n}
  \rput(4.762344,1.678317){$2\delta$}
  \usefont{T1}{ptm}{m}{n}
  \rput(5.242344,-3.421683){$\mathcal{Y}_\delta$}
  \rput{-45.0}(0.9004505,4.0028687){\psarc[linewidth=0.06,linecolor=color183](5.282115,0.9144943){1.9798989}{53.020363}{108.954254}}
  \rput{45.0}(2.214508,-3.3712766){\psarc[linewidth=0.06,linecolor=color1918](5.176745,0.98750955){1.8596908}{67.62315}{126.45429}}
  \rput{180.0}(10.521875,1.7466339){\psarc[linewidth=0.06,linecolor=color183](5.2609377,0.87331694){1.98}{65.897766}{119.43849}}
  \usefont{T1}{ptm}{m}{n}
  \rput(5.2923436,-1.4216831){\color{color183}$\ell_1$}
  \usefont{T1}{ptm}{m}{n}
  \rput(7.132344,2.298317){\color{color183}$\ell_2$}
  \usefont{T1}{ptm}{m}{n}
  \rput(3.4523437,2.3783169){\color{color183}$\ell_3$}
  \end{pspicture} 
  }
  \caption{\small\em The $Y$-junction from the previous Figure~\ref{fig:sketch1} converted into a manifold $\Y_\delta$ by extending all the branches to some finite thickness $\delta > 0$.}
  \label{fig:sketch2}
\end{figure}

\subsubsection{Conservation laws approach}\label{sec:law}

A weaker version of condition \eqref{eq:kirch} can be derived from variational considerations by following \cite{Bibikov2009}. Consider a complex-valued field $u(x,t): \R\times\R^+ \mapsto \C$ whose behaviour is described by Lagrangian density \eqref{eq:Lc}. The energy-momentum tensor $\Tt = (\T^{\alpha\beta})$ of the field $u(x,t)$ satisfies the following conservation laws
\begin{equation}\label{eq:cons}
  \partial_\alpha\T^{\alpha\beta}\ =\ 0,
\end{equation}
where $\partial_0 := \partial_t$, $\partial_1 := \partial_x$. The components of tensor $\Tt$ are given in \cite{Ryder1996}
\begin{equation}\label{eq:tensor}
  \T^{\alpha\beta}\ =\ \frac{\delta \L^{c}_{\mathrm{sG}}}{\delta (\partial_\alpha u)}\;\partial^{\,\beta} u\ +\ \frac{\delta \L^{c}_{\mathrm{sG}}}{\delta (\partial_\alpha u^*)}\;\partial^{\,\beta} u^*\ -\ g^{\alpha\beta}\L^{c}_{\mathrm{sG}}, \qquad \alpha,\; \beta\ =\ 0,1,
\end{equation}
where $g^{\alpha\beta} := \diag\{1, -1\}$ is the Minkowski tensor. The contravariant derivative operator $\partial^{\,\alpha}$ is related to $\partial_\beta$ by $\partial^{\,\alpha} := g^{\alpha\beta}\partial_\beta$. It is straightforward to compute the components of $\Tt$ from \eqref{eq:tensor}:
\begin{equation}\label{eq:comp}
  \T^{00}\ =\ \half\partial_0 u\,\partial_0 u^*\ +\ \half\partial_1 u\,\partial_1 u^*\ +\ 1\ -\ \cos u, \quad \T^{01}\ =\ -\bigl(\partial_0u\,\partial_1u^*\ +\ \partial_0u^*\partial_1u\bigr).
\end{equation}
One can easily identify $\T^{00}$ with the energy density and $\T^{01}$ with the energy flux (see also equation \eqref{eq:cons}). Postulating the energy conservation on the $Y$-junction one obtains the following condition (see \cite{Bibikov2009} for more details):
\begin{equation*}
\left.\T^{10} \right|_{\x\to C,\; \x\in\S_1}\ +\ \left.\T^{10} \right|_{\x\to C,\; \x\in\S_2}\ +\ \left.\T^{10} \right|_{\x\to C,\; \x\in\S_3}\ =\ 0,
\end{equation*}
which can be expanded according to \eqref{eq:comp} to give finally the following analogue of the Kirchhoff condition:
\begin{equation}\label{eq:kirch2}
  \partial_t u^*\bigl(\partial^1u\ +\ \partial^2u\ +\ \partial^3u\bigr)\ +\ \partial_t u\bigl(\partial^1u^*\ +\ \partial^2u^*\ +\ \partial^3u^*\bigr)\ =\ 0 \mbox{ at } \x\ =\ C.
\end{equation}
The last condition is weaker than \eqref{eq:kirch} in the sense that all solutions to \eqref{eq:kirch} satisfy \eqref{eq:kirch2}. However, the inverse generally is not true.

In the sequel of the paper we will require that continuity \eqref{eq:cont} and ``\emph{charge}'' conservation \eqref{eq:kirch} conditions are satisfied at junction points.

\begin{remark}
One can notice that this approach can be generalized straightforwardly to any general nonlinearity in equation \eqref{eq:sg}, not only $f(u) = \sin u$:
\begin{equation*}
  u_{tt}\ -\ u_{xx}\ +\ f(u) \ =\ 0.
\end{equation*}
The reason is that the components of tensor $\Tt$ depend only on the differential terms in \eqref{eq:sg}.
\end{remark}

\begin{remark}
The compatibility conditions \eqref{eq:cont}, \eqref{eq:kirch} at the graph vertices can be straightforwardly generalized to the situations where any (finite) number of strings meet at one junction point.
\end{remark}

\begin{remark}
It goes without saying that the initial condition $u(\x,t=0) = u_0(\x)$ on a graph $G$ has to be compatible with (\ie satisfy) conditions \eqref{eq:cont}, \eqref{eq:kirch} as well.
\end{remark}

\subsection{Conditions on junctions: the discrete case}

It is straightforward to obtain the discrete version of compatibility conditions \eqref{eq:cont}, \eqref{eq:kirch} by following the approach proposed in \cite{Caputo2014}. By using the continuity condition \eqref{eq:cont}, we can employ for simplicity the forward finite differences written on adjacent nodes (see Figure~\ref{fig:sketch1} for the illustration). Let us denote the values of the solution at neighbouring points $\x_i$, $i=0,\ldots,3$ as $u_i := u(\x_i)$. Then, the discrete compatibility condition reads:
\begin{equation*}
  \frac{u_1 - u_0}{\Delta x}\ +\ \frac{u_2 - u_0}{\Delta x}\ +\ \frac{u_3 - u_0}{\Delta x}\ =\ 0 \quad \Rightarrow \quad u_0\ =\ \third\bigl(u_1\ +\ u_2\ +\ u_3\bigr).
\end{equation*}
Obviously, as in the continuous case, the last discrete condition can be generalized to any finite number of adjacent strings.

\section{Numerical results}\label{sec:num}

In order to represent graphically the solution, one has to specify also an embedding of the graph $G$ on $\R^2$ (for planar graphs and $\R^3$ in the general case), \ie a family of regular maps $g_i: e_i\in E \mapsto \R^2$ which satisfy the natural compatibility conditions at the vertices. For the graph $G_0$ we chose a natural embedding shown on Figure~\ref{fig:init}. The initial condition consists of three kinks \eqref{eq:kink0} initially placed on the edges $e_1$, $e_5$, $e_6$ propagating vertically upwards and connecting $0$ to $2\pi$. Solution values on other edges are chosen in order to satisfy the continuity condition \eqref{eq:cont}. The values of all numerical parameters are given in Table~\ref{tab:params}. The total energy of this system is equal to
\begin{equation*}
  \E(0)\ =\ 3\times 8\gamma\ \approx\ 76.86151382644181,
\end{equation*}
where $\gamma$ is defined in \eqref{eq:kink}. For the sake of comparison, the energy at the end of the simulation was equal to $\E(T) \approx 76.97$, which shows good conservative properties of the scheme (the relative error is less than $1.5\%$).

\begin{figure}
  \centering
  \includegraphics[width=0.65\textwidth]{figs/InitCond.eps}
  \caption{\small\em Embedding of the graph $G_0$ in the plane $\R^2$ along with the initial condition (represented with the color of the edges $E$) for three kinks propagating upwards in the network. Graph vertices $V$ are represented with bold blue points.}
  \label{fig:init}
\end{figure}

\begin{table}
  \centering
  \begin{tabular}{||>{\columncolor[gray]{0.85}}c||>{\columncolor[gray]{0.85}}c||}
  \hline\hline
  \textit{Parameter} & \textit{Value} \\
  \hline\hline
    Kink speed, $c_0$ & $0.95;\ 0.5$ \\
    Time step, $\Delta t$ & $0.01$ \\
    Final simulation time, $T$ & $25.0$ \\
    Number of time steps, $N_t$ & $4000$ \\
    Number of points, $N$ & $500$ \\
    Spatial discretization step, $\Delta x$ & $0.02$ \\
  \hline\hline
  \end{tabular}
  \caption{\small\em Parameters used in numerical simulations presented in this manuscript.}
  \label{tab:params}
\end{table}

\begin{figure}
  \centering
  \subfigure[$t\ =\ 6.4$]{% a
  \includegraphics[width=0.48\textwidth]{figs/SG_T_6_40.eps}}
  \subfigure[$t\ =\ 9.45$]{% b
  \includegraphics[width=0.48\textwidth]{figs/SG_T_9_45.eps}}
  \subfigure[$t\ =\ 16.3$]{% c
  \includegraphics[width=0.48\textwidth]{figs/SG_T_16_30.eps}}
  \subfigure[$t\ =\ 19.2$]{% d
  \includegraphics[width=0.48\textwidth]{figs/SG_T_19_20.eps}}
  \subfigure[$t\ =\ 26.1$]{% e
  \includegraphics[width=0.48\textwidth]{figs/SG_T_26_10.eps}}
  \subfigure[$t\ =\ 33.0$]{% f
  \includegraphics[width=0.48\textwidth]{figs/SG_T_33_00.eps}}
  \caption{\small\em Evolution of three \acs{sg} kinks (with $c = 0.95$) on the graph $G_0$.}
  \label{fig:dyn}
\end{figure}

The evolution of this initial condition on the interval $[0, T]$ under the \acs{sg} dynamics is shown on Figure~\ref{fig:dyn}. When the kinks arrive at vertices $v_2$, $v_3$ and $v_4$ (see Figure~\ref{fig:dyn}(\textit{a})) they split in six kinks (see Figure~\ref{fig:dyn}(\textit{b})), which collide again right in the middle of the edges $e_2$, $e_3$ and $e_4$. We observe a topological change at the moment of the collision (see Figure~\ref{fig:dyn}(\textit{c}, \textit{d})), since all the kinks switch from the ground state $0\ \leadsto\ 2\pi$ to $0\ \leadsto\ -2\pi$ as a result of the mutual reflection. Then, the newly generated kinks propagate vertically downwards along the graph edges $e_1$, $e_5$ and $e_6$ (see Figure~\ref{fig:dyn}(\textit{e})). Finally, at the end of the simulation the three kinks collide again in the vicinity of the vertice $v_1$. At the moment of collision there is another topological change from $0\ \leadsto\ -2\pi$ to $-2\pi\ \leadsto\ -4\pi$. The corresponding video illustration can be watched at this URL:

  \url{http://youtu.be/rwZ4d_T7nTs}

Let us perform another simulation, where we take the same set-up as described above, but the kinks are initialized to have the speed $c = 0.5$. One such kink is depicted on Figure~\ref{fig:kink}. The dynamics of this initial condition is shown on Figure~\ref{fig:dyn2}. In agreement with previous investigations \cite{Gulevich2008, Caputo2014} this kink does not possess enough energy to go through a junction. Consequently, the dynamics is confined only to the subgraph whose edges were initialized with kinks. Here again we observe a similar phenomenon to the previous case. When the kinks collide at the vertice $v_1$ (around $t = 30$), there is a topological change and all three solutions shift from $0\ \leadsto\ 2\pi$ to $2\pi\ \leadsto\ 4\pi$. The corresponding video illustration can be watched at this URL:

  \url{http://youtu.be/CtglMe0IcBk}

\begin{figure}
  \centering
  \subfigure[$t\ =\ 2.0$]{% a
  \includegraphics[width=0.48\textwidth]{figs/SG2_T_2_00.eps}}
  \subfigure[$t\ =\ 9.0$]{% b
  \includegraphics[width=0.48\textwidth]{figs/SG2_T_9_00.eps}}
  \subfigure[$t\ =\ 16.0$]{% c
  \includegraphics[width=0.48\textwidth]{figs/SG2_T_16_00.eps}}
  \subfigure[$t\ =\ 24.0$]{% d
  \includegraphics[width=0.48\textwidth]{figs/SG2_T_24_00.eps}}
  \subfigure[$t\ =\ 34.0$]{% e
  \includegraphics[width=0.48\textwidth]{figs/SG2_T_34_00.eps}}
  \subfigure[$t\ =\ 40.0$]{% f
  \includegraphics[width=0.48\textwidth]{figs/SG2_T_40_00.eps}}
  \caption{\small\em Evolution of three \acs{sg} kinks (with $c = 0.5$) on the graph $G_0$.}
  \label{fig:dyn2}
\end{figure}

\section{Conclusions and perspectives}\label{sec:concl}

In the present study we considered a discrete formulation of the celebrated \acf{sg} equation on domains which are not manifolds. More precisely, the 1D lattices can be assembled into arbitrary graphs (networks) and the coupled dynamics can be efficiently simulated using a simple symplectic numerical scheme. The edges of the graph deserve a special treatment based on local conservation laws. The performance of this formulation was illustrated on a practical example of a connected graph involving four cycles. A sequence of topological changes was detected during the collisions of elementary \acs{sg} kinks at the graph vertices. We considered two different situations when the kinks initially were super- and sub-critcal. In the latter case the system dynamics is restricted only to a sub-graph since the kinks do not possess enough energy to go through the junctions.

In this study we considered only simple graphs $G = (V,E)$. However, we believe that the generalization to multi\footnote{A multi-graph may contain multiple edges connecting two vertices.}- and pseudo\footnote{A pseudo-graph may contain loops connecting an edge with itself. Their effect on kink dynamics on a network is to be investigated.}-graphs should be straightforward.

Using two complementary approaches we derived the coupling conditions in junction points for the \acs{sg} equation. One of the methods relies heavily on the variational structure of these equations (see Section~\ref{sec:law}). Consequently, its application is limited to models which derive from a Lagrangian functional. The other method (back-to-manifold approach) is more general, since it consists in ``\emph{inflating}'' the channels, integrating the PDE on the resulting manifold and taking the limit $w\to 0$, where $w$ is the characteristic width of the ``\emph{tubes}''. The main findings of our study can be summarized as follows:
\begin{itemize}
  \item Two complementary methodologies to derive compatibility conditions for the \acs{sg} equation in junction points were proposed. They lead to equivalent results in our context
  \item A variational symplectic scheme for the \acs{sg} equation on graphs was described
  \item The scheme performance was illustrated on two test-cases of a system of kinks propagating in a closed network
 \end{itemize}
The application of the techniques developed in the present study to the classical fluids in the framework of \acs{nswe} equations is in progress.

\section*{Acknowledgments}
\addcontentsline{toc}{section}{Acknowledgments}

The authors acknowledge the Centre de Ressources Informatiques de Haute Normandie where most of the calculations were done. Moreover, we would like to thank Professor Michel \textsc{Raibaut} (LAMA UMR \#5127, Universit\'e Savoie Mont Blanc) for stimulating discussions on the geometrical matters.

%%% Bibliography
\addcontentsline{toc}{section}{References}
\bibliographystyle{abbrv}
\bibliography{biblio}

\end{document}

SG_Trees.bbl
\begin{thebibliography}{10}

\bibitem{Arnold1997}
V.~I. Arnold.
\newblock {\em {Mathematical Methods of Classical Mechanics}}.
\newblock Springer, New York, 2nd edition, 1997.

\bibitem{Bibikov2009}
P.~N. Bibikov and L.~V. Prokhorov.
\newblock {Mechanics not on a manifold}.
\newblock {\em J. Phys. A: Math. Gen}, 42(4):045302, Jan. 2009.

\bibitem{Bona2008}
J.~L. Bona and R.~Cascaval.
\newblock {Nonlinear dispersive waves on trees}.
\newblock {\em Canadian Applied Mathematics Quarterly}, 16(1):1--18, 2008.

\bibitem{Bressan2014}
A.~Bressan, S.~Cani\'{c}, M.~Garavello, M.~Herty, and B.~Piccoli.
\newblock {Flows on networks: recent results and perspectives}.
\newblock {\em EMS Surveys in Mathematical Sciences}, 1(1):47--111, 2014.

\bibitem{Caputo2014}
J.-G. Caputo and D.~Dutykh.
\newblock {Nonlinear waves in networks: model reduction for sine-Gordon}.
\newblock {\em Phys. Rev. E}, 90:022912, 2014.

\bibitem{Dashen1974}
R.~F. Dashen, B.~Hasslacher, and A.~Neveu.
\newblock {Nonperturbative methods and extended-hadron models in field theory.
  I. Semiclassical functional methods}.
\newblock {\em Phys. Rev. D}, 10(12):4114--4129, 1974.

\bibitem{Faddeev1987}
L.~D. Faddeev and L.~Takhtajan.
\newblock {\em {Hamiltonian Methods in the Theory of Solitons}}.
\newblock Springer, Berlin Heidelberg New York, 1987.

\bibitem{Fock1926}
V.~Fock.
\newblock {On the invariant form of the wave equation and the equations of
  motion for a charged point mass}.
\newblock {\em Z. Phys.}, 39(2-3):226--232, 1926.

\bibitem{Gould2012}
R.~Gould.
\newblock {\em {Graph Theory}}.
\newblock Dover Publications Inc., dover edition, 2012.

\bibitem{Gulevich2006}
D.~Gulevich and F.~Kusmartsev.
\newblock {Flux Cloning in Josephson Transmission Lines}.
\newblock {\em Phys. Rev. Lett.}, 97(1):017004, July 2006.

\bibitem{Gulevich2008}
D.~Gulevich, F.~Kusmartsev, S.~Savel'ev, V.~Yampol'skii, and F.~Nori.
\newblock {Shape Waves in 2D Josephson Junctions: Exact Solutions and Time
  Dilation}.
\newblock {\em Phys. Rev. Lett}, 101(12):127002, Sept. 2008.

\bibitem{Haefliger1970}
A.~Haefliger.
\newblock {Feuilletages sur les vari\'{e}t\'{e}s ouvertes}.
\newblock {\em Topology}, 9(2):183--194, May 1970.

\bibitem{Leimkuhler2004}
B.~Leimkuhler and S.~Reich.
\newblock {\em {Simulating Hamiltonian Dynamics}}, volume~14 of {\em Cambridge
  Monographs on Applied and Computational Mathematics}.
\newblock Cambridge University Press, Cambridge, 2004.

\bibitem{Mehmeti2003}
F.~A. Mehmeti and V.~R\'{e}gnier.
\newblock {Splitting of energy of dispersive waves in a star-shaped network}.
\newblock {\em ZAMM}, 83(2):105--118, Feb. 2003.

\bibitem{Moore2003}
B.~Moore and S.~Reich.
\newblock {Backward error analysis for multi-symplectic integration methods}.
\newblock {\em Numerische Mathematik}, 95(4):625--652, 2003.

\bibitem{Ryder1996}
L.~H. Ryder.
\newblock {\em {Quantum Field Theory}}.
\newblock Cambridge University Press, Cambridge, 2nd edition, 1996.

\bibitem{Scott2003}
A.~Scott.
\newblock {\em {Nonlinear Science: Emergence and Dynamics of Coherent
  Structures}}.
\newblock Oxford University Press, 2nd edition, 2003.

\bibitem{Takhtadzhyan1974}
L.~A. Takhtadzhyan and L.~D. Faddeev.
\newblock {Essentially nonlinear one-dimensional model of classical field
  theory}.
\newblock {\em Theor. Math. Phys.}, 21(2):1046--1057, 1974.

\bibitem{Vassilevskii2011}
Y.~Vassilevskii, S.~Simakov, V.~Salamatova, Y.~Ivanov, and T.~Dobroserdova.
\newblock {Numerical issues of modelling blood flow in networks of vessels with
  pathologies}.
\newblock {\em Russ. J. Numer. Anal. Math. Modelling}, 26(6):605--622, 2011.

\bibitem{Zuazua2013}
E.~Zuazua.
\newblock {Control and stabilization of waves on 1-d networks}.
\newblock {\em Lecture Notes in Mathematics}, 2062:463--493, 2013.

\end{thebibliography}
