\section{Introduction}
The interest in magneto-elastic materials finds its motivation in the growing
variety of new materials among which magneto-rheological elastomers or
magneto-sensitive polymeric composites \cite{Hossain-et-al-2015a, Hossain-et-al-2015b}
may be mentioned.  A whole Special Issue devoted to {\it Magnetoelastic Materials}
is going to be published soon \cite{SI2020} in the Journal {\it Materials}.
Many applications of magneto-elastic materials, covering a wide area of interest  from
technological to biomedical devices, see e.g. \cite{Ren},
can be listed. In particular, also
two dimensional problems are subject of applicative investigations \cite{Hadda}.
The model we consider is a two dimensional simplified one, however, we believe that, it might
open the way to further applications, possibly, via perturbative methods \cite{Bernard}.

We study the functional energy of a magneto-elastic material, that is
a material which is capable of deformation and magnetisation. The
magnetisation is a phenomenon that does not appear at a macroscopic
level, it is characterised by the {\it magnetisation vector} whose magnitude
is independent of the position while its direction which can vary from one point to another.\\
In this context,  the magnetisation vector ${\bf m}$ is a map from $\Omega$
(a bounded open set of $\mathbb{R}^2$) to $S^2$ (the unit sphere of
$\mathbb{R}^3$). In particular, here we assume $\Omega$ is the unit
disk of $\mathbb{R}^2$. The magnetisation distribution is well
described by a free energy functional which we assume composed of
three terms, namely the {\it exchange} energy ${E}_{\rm ex}$, the
{\it elastic} energy $E_{\rm el}$ and the {\it elastic-magnetic}
energy $E_{\rm em}$. In Section~\ref{model} we detail the three
energetic terms and, after some simplifications, derive the proposed
functional for describing some phenomena. Assuming the hypothesis of
radially symmetric maps, i.e.
$$ {\bf m} = ( \cos \theta \sin h(r), \sin \theta \sin h(r), \cos h(r)
),$$ we get to the analysis of a one-dimensional energy functional
that can be expressed in terms of the only scalar function $h$. The
effect of the elastic deformation reveals through a positive
parameter $\mu$ which characterizes the connection between the
magnetic and elastic processes. In Section~\ref{minproblem}
 the minimisation of the energy functional, namely
$$E(h)= \pi \int_0^1 \left[h_r^2+\left(\frac{\sin h}{r}\right)^2 -
\frac{\mu}{2} (\sin2h )^2 \right] rdr,$$ is the aim of  our paper.
In particular, we prove that there exists a critical value $\mu^0$
such that for $\mu \le \mu^0$ the functional energy is not negative
and there is only a global minimiser that is the trivial solution
$h\equiv 0$; for $\mu
> \mu^0$ other nontrivial minimisers appear, moreover the energy
takes negative values. The local
bifurcation analysis is carried out. More precisely we prove that at
the point $\mu^0$, two branches of  minimisers, with small norm,
bifurcate from the trivial stable solution. This local analysis does
not exclude the existence of other solutions of the minimisation
problem even for $\mu = 0$ (see also the results by Brezis and Coron
in \cite{BC} concerning the solutions of harmonic maps from the unit
disk in $\mathbb{R}^2$ to the sphere $S^2$).

For the modelling of magneto-elastic interactions see also \cite{BPGV}, \cite{Bw},
\cite{CPV}, \cite{csvv}, \cite{csvv1}, \cite{He}, \cite{VV}. Magneto-viscoelastic problems are studied in \cite{GVS2010},  \cite{GVS2012} and  \cite{MGVS2017}.
Moreover we recall that the phenomenon of bifurcation of minimising
harmonic maps has been studied by Bethuel, Brezis, Coleman, H{\rm
\'{e}}lein (see \cite{BBCH}) in a different physical context.

\medskip

\section{The model} \label{model}
\setcounter{equation}{0}

 We start with the general three-dimensional theory.
We assume $\Omega \subset \mathbb{R}^3$ is the volume of the
magneto-elastic material and $\partial \Omega$ its boundary. Let
$x_i,\, i=1,2,3$ be the position of a point ${\bf x}$ of $\Omega$ and
denote by
$$u_i=u_i({\bf x}),\qquad i=1,2,3$$
the components of the displacement vector ${\bf u}$ and by
$$\varepsilon_{kl}({\bf u})=\frac{1}{2}(u_{k,l}+u_{l,k}),\qquad k,l=1,2,3$$
the deformation tensor where, as a common praxis, $u_{k,l}$ stands
for ${\frac{\partial u_k}{\partial x_l}}$. Moreover we denote by
$$m_{j}=m_{j}({\bf x}),\qquad j=1,2,3$$
the components of the magnetisation vector ${\bf m}$ that we assume of
unit modulus, i.e. $|{\bf m}|=1$. \\
In the sequel, where not specified,
the Latin indices vary in the set \{1,2,3\} and the summation over repeated indices is assumed.
We first define the exchange energy which arises from exchange
neighbourhood interactions as
\begin{equation} \label{Eex}
E_{\rm ex}({\bf m}) =\frac{1}{2} \int_{\Omega} a_{ij} m_{k,i}
m_{k,j}d\Omega
\end{equation}
where $a_{ijkl} = a_1\delta_{ijkl} + a_2\delta_{ij}\delta_{kl}$ with $a_1,a_2\ge0$ and  $\delta_{ijkl}=\delta_{ik}\delta_{jl}$ is the fourth-order identity tensor.
This integral represents the interface energy between magnetised domains with different orientations. For most magnetic materials $div\, {\bf m}=\delta_{ij}{m}_{i,j}=0$, so hereafter we assume $a_1=a>0$ and $a_2=0$ (see \cite{LL}).
The magneto-elastic energy is due to the coupling between the magnetic moments and the elastic lattice. For cubic crystals it is assumed to be
\begin{equation} \label{Eem}
E_{\rm em}({\bf m},{\bf u}) =\frac{1}{2} \int_{\Omega} \lambda_{i j k l}
m_{i} m_{j} \varepsilon_{kl}({\bf u})d\Omega
\end{equation}
where  $\big.\mathbb{L}= \{{\lambda}_{klmn}\}$ denotes the \textit{ {magneto-elasticity tensor}} whose entries $ \lambda_1, \lambda_2, \lambda_3\ge0$, and $\lambda_{i j k l}=\lambda_{1} \delta_{i j k l}+ \lambda_{2}
\delta_{i j }\delta_{kl}+ \lambda_{3}(\delta_{i k}\delta_{j
l}+\delta_{i l}\delta_{j k})$ with $\delta_{i j k l}=1$ if
$i=j=k=l$ and $\delta_{i j k l}=0$ otherwise. Moreover we
introduce the elastic energy
\begin{equation} \label{Eel}
E_{\rm el}({\bf u})= \frac12 \int_\Omega \sigma_{ijkl}
 \varepsilon_{ij}({\bf u})\varepsilon_{kl}({\bf u})\,d\Omega
 \end{equation}
where $\Big.\mathbb{E}= \{ \epsilon_{lm}\}$ indicates the  \textit{ {strain tensor}}  $\sigma_{ijkl}$
satisfying the following symmetry property
$$\sigma_{ijkl}=\sigma_{klij}=\sigma_{jilk}$$
and moreover the inequality
$$\sigma_{ijkl}\varepsilon_{ij}\varepsilon_{kl}\ge\beta\varepsilon_{ij}\varepsilon_{ij}$$
holds for some $\beta>0$. In the isotropic case
$$\sigma_{ijkl} = \tau_1\delta_{ijkl} + \tau_2\delta_{ij}\delta_{kl},\qquad \tau_1, \tau_2\ge0.$$
The resulting energy functional ${E}$ is given by
\begin{equation}
E({\bf m},{\bf u})= E(m,u) = E_{ex} (m) + E_{em} (m,u) + E_{ve} (u),
\end{equation}
which {{after some manipulations \cite{BPGV, GVS2012}}}, under the assumption the material is isotropic, reads
\begin{multline}
E({\bf m},{\bf u})=
 \frac{1}{2}\int_{\Omega} a|\nabla {{\bf m}}|^2 d\Omega +
\frac{1}{2}\int_{\Omega} \left[\tau_1
|\nabla {{\bf u}}|^2+\tau_2 (div\, {{\bf u}})^2\right] d\Omega +\\
+\frac{1}{2}\int_{\Omega} \left[\lambda_{1}\delta_{k l i j
u_{j,i}
m_{k}m_{l} +\lambda_{2} |{{\bf m}}|^2 div \,{{\bf u} }+2\lambda_3(\nabla
u_i\cdot {{\bf m}}) m_{i} \right] d\Omega~.
\end{multline}
\subsection{A simplified 2D model}

To get the proposed model we make some approximations. First of all
we assume
 $\Omega \subset \mathbb{R}^2$
 and neglect the components in plane of the
 displacement vector ${\bf u}$, i.e.  we assume  ${{\bf u}}=(0,0,w)$, which implies $div\, u=0$ since $w$ depends only on the plane coordinates. Let $\lambda_3= \lambda$ be a positive constant,

setting\footnote{No need to prescribe  nor $ \lambda_2$ nor $\tau_2\ge0$ since they both appear only as  factors of $div\, u$;  also $\lambda_1$ can be left arbitrary; indeed,   $\delta_{klij}u_{j,i}=0$ since $u_{j,i}\neq0$ only if $j=3$ and $i=1,2$ but $\delta_{klij}=0$ when $ j\neq i$.}
$\tau_1= 1$  and $a=1$,
the functional
${E}$ reduces to
\begin{equation} \label{E}
\displaystyle{{E}({\bf m},w)= \frac{1}{2}\int_{\Omega} \left( |\nabla {{\bf m}}|^2
+ 2\lambda m_3 (m_{\alpha} w_{,\alpha})+|\nabla w|^2\right)d\Omega}
\end{equation}
where the Greek indices vary in the set \{1,2\}.

Setting $\Omega \equiv D = \{ (x,y) \in \mathbb{R}^2 : \, x^2+y^2 <
1\}$ and assuming radial symmetry, further to $w=w(r)$, we can express the components of the  vector ${\bf m}$
in terms of $r$, that is, of the form
$$ {\bf m} = \left( \frac{x}{r} \sin h(r), \frac{y}{r} \sin h(r), \cos h(r)
\right), \qquad r=\sqrt{x^2+y^2},$$
where $h: (0,1)\subset {{\rm I}\!{\rm R}} \to {{\rm I}\!{\rm R}}$ is an unknown regular function.
Using the fact that $\partial_x r = \frac{x}{r}$ and $\partial_y r = \frac{y}{r}$ we deduce by the chain rule, where  $h_r:=\displaystyle{{ d h} \over {dr}}$  denotes the derivatives of the $h$ with respect to  the variable $r$, it follows
\begin{equation*}
\partial_x {\bf m} = \left(\frac{\sin h}{r}+ \frac{x^2}{r}\left(\frac{\sin h}{r}\right)_{\!r}, ~\frac{xy}{r}\left(\frac{\sin h}{r}\right)_{\!r} , ~\frac{x}{r}\left(\cos h\right)_{\!r}\right)
\end{equation*}
\begin{equation*}
\partial_y {\bf m} = \left({\frac{xy}{r}\left(\frac{\sin h}{r}\right)_{\!r} ,~ \frac{\sin h}{r}+ \frac{y^2}{r}\left(\frac{\sin h}{r}\right)_{\!r} , ~\frac{y}{r}\left(\cos h\right)_{r} }\right).
\end{equation*}
Thus we get
\begin{equation*}\displaystyle
\begin{aligned}
|\nabla {\bf m} |^2&= \left[\frac{\sin h}{r}+ \frac{x^2}{r}\left(\frac{\sin h}{r}\right)_{\!r} \right]^2 +\left[{\frac{\sin h}{r}+ \frac{y^2}{r}\left(\frac{\sin h}{r}\right)_{\!r} }\right]^2 +
2\left[{\frac{xy}{r}\left(\frac{\sin h}{r}\right)_{\!r}}\right]^2 +[\left(\cos h\right)_{\!r} ]^2\\
&=2\left(\frac{\sin h}{r}\right)^2 +
\frac{x^4+2x^2y^2 +y^4}{r^2}\left[{\left(\frac{\sin h}{r}\right)_{\!r}}\right]^2 +2 \frac{x^2+y^2}{r^2}\sin h \left(\frac{\sin h}{r}\right)_{\!r} + h_{r}^2 (\sin h)^2 \\
&= 2\left(\frac{\sin h}{r}\right)^2 + {r^2}\left[\left(\frac{\sin h}{r}\right)_{\!r}\right]^2 +2 \sin h\left(\frac{\sin h}{r}\right)_{\!r} + h_{r}^2(\sin h)^2 \\
&= \left(\frac{\sin h}{r}\right)^2 + \left[\frac{\sin h}{r}+{r}\left(\frac{\sin h}{r}\right)_{\!r}\right]^2
 +h_{r}^2 (\sin h)^2 \\
&= \left(\frac{\sin h}{r}\right)^2 + \left[\frac{\sin h}{r}+{r}\left(\frac{r h_{\!r} \cos h -
\sin h}{r^2}\right)\right]^2 + h_{r}^2 (\sin h)^2 \\
&= \left(\frac{\sin h}{r}\right)^2  +h_r^2.
\end{aligned}
\end{equation*}`

So the energy (\ref{E}), when we recall the assumed radial symmetry implies also $w=w(r)$, adopting the notation $w_r:=\displaystyle{{ d w} \over {dr}}$, becomes
$$
E(h,w)= \pi \int_0^1 \left[h_r^2 +\left(\frac{\sin h}{r}\right)^2
+\lambda \sin2h \, w_r + w_r^2\right] rdr$$
 and from that we deduce the governing
equations
\begin{equation} \label{eqs}
\left \{ \begin {array}{l}\displaystyle{ h_{rr} +\frac{h_r}{r}
-\frac{\sin2h}{2r^2} -\lambda \cos2h \,w_r =0}\\
\\
\displaystyle{ w_{rr} + \frac{w_r}{r} + \frac{\lambda}{2} \left[
(\sin 2h)_r + \frac{\sin 2h}{r}\right]=0}.
\end{array}\right.
\end{equation}
We prescribe the following boundary conditions
\begin{equation} \label{bcw}
w_r(0)=0,\,\,w(1)=0,
\end{equation}
where the first condition is motivated by the symmetry assumptions, while the second one corresponds
to prescribe  the boundary of $\Omega$ is fixed, and
\begin{equation} \label{bch}
 h_r(1)=0.
 \end{equation}
 Solving the second equation of (\ref{eqs}) which can be written
 \begin{equation*}
(r w_r)_{r} + \frac{\lambda}{2}(r \sin 2h)_{r} =0 \Leftrightarrow w_r = -\frac{\lambda}{2}\sin 2h
\end{equation*}
where the double implication is guaranteed when we set $h(0)=0$. Then, letting
$\mu= \lambda^2/2$ we get  the equation
\begin{equation} \label{eqh}
h_{rr} +\frac{h_r}{r} -\frac{\sin 2h}{2r^2} +\mu \sin 2h \cos2h  =0
\end{equation}
and the energy E becomes
\begin{equation} \label{Eh}
E(h)= \pi \int_0^1 \left[h_r^2+\left(\frac{\sin h}{r}\right)^2 -
\frac{\mu}{2} (\sin2h )^2 \right] rdr.
\end{equation}

The variational analysis of the functional $E(h)$ is the objective
of the following section.

\section{The minimisation problem} \label{minproblem}

\begin{lemma} \label{lemma1}
 Let us define
\begin{equation}
V=\{v ~\vert ~ v_r, {v\over r}\in L^2(0,1;r dr)\}.
\end{equation}
$V$ is a Hilbert space equipped with the norm
\begin{equation}
\displaystyle{\vert\!\vert v \vert\!\vert^2=
\int_0^1{ ({ v_r}^2+ {v^2\over r^2}) r dr~.}}
\end{equation}
\end{lemma}

\noindent\begin{proof}

Let $v_n$ be a Cauchy sequence in $V$, $\{(v_n)_r\}, \displaystyle{\left\{{v_n\over r}\right\}}$ are Cauchy sequences in
 $$L^2(r dr)=L^2(0,1;r dr)$$ and there exist $h, g$ such that
\begin{equation}\label{m1}
\{(v_n)_r\}, \displaystyle{\left\{{v_n\over r}\right\}} \to g,h~~\text{in}~~L^2(r dr)
\end{equation}
Set $\tilde h=h r$. Since
\begin{equation}\label{1m}
\int_0^1 \left({v_n\over r}-h\right)^2  r dr~\to 0
\end{equation}
one has
\begin{equation}
{v_n\to rh ~~\text{in}~~L^2 \big(0,1;{dr\over r}\big)~.}
\end{equation}
but also in ${\cal D}^{\prime}(0,1)$ so that
\begin{equation}
{(v_n)_r\to \tilde h_r ~~\text{in}~~{\cal D}^{\prime}~.}
\end{equation}
We deduce from \eqref{1m} that $\tilde h_r=g$ and thus $ \tilde h\in V$ and since
\begin{equation}
(v_n)_r, {v_n\over r} \to \tilde h_r, {\widetilde h\over r}~~\text{in}~~L^2(r dr)
\end{equation}
one has $v_n\to \tilde h \in V$. This completes the proof of  Lemma \ref{lemma1}.

\end{proof}\hfill\hfill $\Box$

\begin{lemma} \label{lemma2}

$$ V \subset \{v\in C([0,1]) ~\vert ~v(0)=0 \} $$
\end{lemma}

\noindent\begin{proof}

For $x,y\in (0,1]$ one has

\begin{equation}
\vert x v(x)- y v(y) \vert = \left\vert \int_x^y {(r v)_r dr }\right\vert =  \left\vert \int_x^y { r (v_r +{v\over r}) dr }\right\vert
\\ \le \int_x^y { r (|v_r| +{{|v|}\over r}) dr ~.}
\end{equation}
Using the Cauchy-Young  inequality $a \le {1\over 2} a^2 + {1\over 2}$ one gets

\begin{equation}
\vert x v(x)- y v(y) \vert \le \int_x^y \left\{ { {r\over 2} \left({v_r}^2 +{{v^2}\over {r^2}}\right) }
+ r \right\}dr \to 0 ~~
\text{when}~~ y \to x.
\end{equation}
It follows that $v$ is continuous at any point where $r\neq0$ on $[0,1]$.
Now, one has also
\begin{equation}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
\displaystyle v(x)^2 -  v(y)^2 =\!\!  \int_x^y \!\!  { {d\over {dr}} {v(r)}^2 } dr = \!\!  \int_x^y \!\!  { 2 {v_r \, v} }~ dr =
\!\!
 \int_x^y \!\!  {2 \sqrt{r} ~v_r {v\over {\sqrt{r}}} }dr \\ \\
 \displaystyle \le \int_x^y\!\!  \left\{ { {r} ({v_r}^2 +{{v^2}\over {r^2}})} \right\} ~dr \le \varepsilon
\end{array}
\end{equation}
for $x,y$ small enough (we used again  the Cauchy-Young  inequality).
\smallskip
Thus, when $x\to0$, $v(x)^2$ is a Cauchy sequence and there exist  $l\ge 0$ such that
\begin{equation}
\displaystyle{\lim_{x\to 0} {v(x)^2}=l}.
\end{equation}
If $l>0$ one has for $\varepsilon$ small enough
\begin{equation}
\Vert v \Vert^2\ge \int_0^1{ {v^2\over r} dr }\ge
\int_{\varepsilon^2}^\varepsilon{ \left({l\over 2}\right)^2{{dr}\over r} }
={l^2\over 4}(\ln \varepsilon- 2 \ln \varepsilon) = - {l^2\over 4} \ln \varepsilon
\end{equation}
and a contradiction when $\varepsilon\to 0$. Thus, $l=0$ and this completes the proof of the Lemma \ref{lemma2}.
\end{proof} \hfill $\Box$

\bigskip
\noindent{\bf Remark}
Since $V\subset H^1( \varepsilon,1)$, it follows that $V \subset C^{1/2}( \varepsilon,1)$ for every $\varepsilon$.

One sets

\begin{equation}
E(h)= \pi \int_0^1 \left\{ h_r^2+  ({{\sin h}\over{r}})^2 - \frac{\mu}{2} (\sin 2h)^2 \right\} rdr.
\end{equation}

\bigskip\noindent
One would like to show that $E(h)$ possesses a minimiser on $V$ for any $\mu$.

\noindent\begin{lemma} \label{lemma3}

The energy $E(h)$ is bounded from below on $V$ and one can find a minimising sequence $v_n$ such that
\begin{equation}\label{2m}
0\le v_n\le {\pi \over 2}~.
\end{equation}
\end{lemma}

\noindent\begin{proof}

One has clearly for every $h\in V$
\begin{equation}
E(h)\ge -\pi {|\mu|\over 2}\int_0^1 r dr =  -\pi {|\mu|\over 4}~.
\end{equation}
Thus $$ I= \inf_{h\in V} E(h) $$
exists. Let us denote by $v_n$ a sequence such that
$$ E(v_n) \to I~. $$

If $v_n\in V$, then also   $\vert v_n \vert\in V$ and one has $$ E(v_n)=  E(|v_n|)$$
so, without loss of generality, we assume $v_n\ge 0$.

\bigskip
\begin{figure}[H]
\centering
\epsfig{file=FIG-1-ver8.png, scale=0.31}
\caption{graphical representation.}
\label{fig4}
\end{figure}

\noindent
Then on $v_n> {\pi \over 2}$, we replace $v_n$ by $-v_n +\pi$ (cfr. Fig. 1).
It is clear that
\begin{equation}
\tilde v_n=v_n X_{\{v_n\le {\pi \over 2}\}} + (-v_n +\pi)X_{\{v_n> {\pi \over 2}\}}
\end{equation}
satisfies $\tilde v_n\in V$ and $$E(\tilde v_n)=E(v_n).$$
This completes the proof of the Lemma \ref{lemma3}.
\end{proof}

\hfill\hfill $\Box$

\begin{remark}\label{rem-p4} It could be that $-v_n +\pi$ achieves negative values, but clearly, after a finite number of operations like the one we just did we get a $v_n$ satisfying \eqref{2m}.
\end{remark}

\begin{lemma} \label{Lemma 4} There exists a minimiser $\tilde h$ of $E$ in $V$ satisfying
\begin{equation}
0 \le \tilde h\le  {\pi \over 2}~.
\end{equation}
\end{lemma}

\noindent\begin{proof}
We consider the sequence $\{v_n\}$ constructed in Lemma \ref{lemma3}.
We claim that $\{v_n\}$  is bounded in $V$ independently of $n$. Indeed, one has, since
for some constant $\lambda>0$ one has $\displaystyle{\left({{\sin x}\over{x}}\right)^2\ge \lambda}$,
$\forall x\in[0,{\pi \over 2}]$,
\begin{equation}  \label{m3}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
\displaystyle  \int_0^1 {r \left\{ {({v_n})_r}^2 +{{{v_n}^2}\over {r^2}} \right\} dr \le
\int_0^1 {r \left\{  {({v_n})_r}^2  +{1\over \lambda}\left({{\sin v_n}\over {r}}\right)^2 \right\} dr }}\\
\displaystyle \le \left(1 \vee {1\over \lambda}\right) \int_0^1 {r \left\{  {({v_n})_r}^2  +
\left({{\sin v_n}\over {r}}\right)^2 \right\} dr \le C}
\end{array}
\end{equation}
where $C$ is a constant independent of $n$ and $\vee$ denotes the maximum of two numbers.
Recall that since ${v_n}$ is a minimising sequence one has, for $n$ large enough,
$$ E({v_n}) \le E(0) =0$$
i.e. see the definition of $E$
\begin{equation}
\pi \int_0^1 {\left\{ {({v_n})_r}^2 + \left({{\sin v_n}\over {r}}\right)^2 \right\} r dr}\le
\pi {|\mu|\over 2}\int_0^1{{\sin^2 (2 v_n)} r dr} \le \pi {|\mu|\over 4}~.
\end{equation}
Since  $\displaystyle{\{({v_n})_r}\}, \{{{v_n}\over r}\}$ are bounded in $L^2(r dr)$ one finds a subsequence, still labelled by $n$, such that
$$  {{v_n}\over r} {{\rightharpoonup h}} ~~, ~~ {({v_n})_r} \rightharpoonup  g ~ \text{in} ~  L^2(r dr)~.$$
Set $\tilde h =h r$. The first weak convergence above reads
\begin{equation*}
\displaystyle \int_0^1 {{v_n}\over r} \Psi r dr \to  \int_0^1 h \Psi r dr~~~~, ~~~\forall  \Psi\in L^2 (rdr).
\end{equation*}
In particular, taking $\Psi \in {\cal D}(0,1)$ one see that
$$ {v_n} \to \tilde h =h r  ~ \text{in} ~ {\cal D}^{\prime}(0,1)$$
and thus, by the continuity of the derivative in ${\cal D}^{\prime}$
$$({v_n})_r \to \tilde h_r  = g ~~\text{in}~~  {\cal D}^{\prime}(0,1)~.$$
Thus, we have $ \tilde h\in V$. For any $k\ge 2$ one has also, thank to \eqref{m3}, that
$ {v_n}$ is bounded in $\displaystyle H^{1}\left({1\over k},1\right)$. Thus, by induction, one can find a subsequence
$\{n_{k}\}$ extracted from $\{n_{k-1}\}$ such that
$$ {v_{n_{k}} }\to \tilde h   ~ \text{in} ~ L^{2}\left({1\over k},1\right)~ \text{and a. e..}$$
Then clearly
$$ v_{n_{k}} \to \tilde h  ~~ ~  \text{a.e. on} ~ (0,1).$$
By the dominated Lebesgue theorem one has then that
$$ r \sin 2 {v_{n_{k}} } \to  r \sin 2  \tilde h ~~ ~ \text{in} ~ L^{2}(0,1) $$
$$  \sin {v_{n_{k}} } \to   \sin    \tilde h ~~ ~ \text{a.e.~ on} ~ (0,1)~. $$
Then, since $ x\mapsto x^2$ is convex by the Fatou lemma one has

\begin{equation}\label{weak-eps}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
{\displaystyle{ I={\underline{\lim}} ~E({v_{n_k} }) =
\pi\, {\underline{\lim}} \int_0^1 {\left\{ {({v_{n_k} })_r}^2 + \left({{\sin {v_{n_k} }}\over {r}}\right)^2 \right\} r dr}-
\pi \int_0^1{{\mu\over 2}({\sin (2 {v_{n_k} })})^2 r dr}\ge}}\\ \\
\displaystyle{\ge \pi ~ {\underline{\lim}} \int_0^1 {{({v_{n_n} })_r}^2 r dr}+  \pi ~{\underline{\lim}} \int_0^1 {\left({{\sin {v_{n_k} }}\over {r}}\right)^2 r dr}-
{\pi\mu\over 2} \int_0^1{({\sin (2 \tilde h)})^2 r dr}\ge}\\ \\
{\displaystyle{\ge \int_0^1 {({\tilde h_r})^2 r dr} +  \pi  \int_0^1 {{\underline{\lim}}{\left({{\sin {v_n}}\over {r}}\right)^2} r dr} -
{\pi\mu\over 2}\int_0^1{{(\sin (2 \tilde h))^2} r }dr} =}\\ \\
\displaystyle{\hphantom{1\over1}=E({\tilde h}) = I}~.
\end{array}
\end{equation}
This shows that ${\tilde h}$ is the minimiser that we are looking for.
\end{proof} \hfill $\Box$

 \begin{lemma}\label{Euler-eq}{
The Euler equation of the minimising problem is given by
\begin{equation} \label{m1b}
\left \{ \begin {array}{l}\displaystyle{ -h_{rr}
-\frac{h_r}{r}
+\frac{\sin 2h}{r^2} = \mu \sin 2h \cos 2h } ~~~\text{in}~ (0,1)\\
\\
h(0) =h_r(1)=0
\end{array}\right.
\end{equation}}
\end{lemma}

\noindent\begin{proof}
If $h$ is a minimiser of $E$ on $V$ one has
$$ {d \over{d \lambda}} E(h+\lambda v) \vert_0=0~~~,~~ \forall v\in V $$
Since
\begin{equation}
\displaystyle{ E(h+\lambda v)=
\pi   \int_0^1 {\left\{ (h+\lambda v)_r^2 +{ {{\sin(h+\lambda v)^2}}\over{{r}^2}}
-{\mu\over 2} {\sin (2 (h+\lambda v))}^2
\right\} }r dr}~.
\end{equation}
One gets $\forall v$
\begin{equation}\label{m2b}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
\displaystyle{
 \int_0^1 {\left\{ 2 h_r v_r + 2 {{ \sin h\cos h }\over {{r}^2}}v
-2{\mu}\, {\sin (2 h)\cos (2 h)}\, v
\right\} r dr}=0 }\\ \\
\displaystyle{ \Longleftrightarrow
 \int_0^1 {\left\{  h_r v_r +  {{\sin  (2 h)}\over {2{r}^2}} v
-{\mu}\, {\sin (2 h)\cos (2 h)}\, v
\right\} r dr}=0~~~, ~~~\forall v\in V}~.
\end{array}
\end{equation}
Thus, in the distributional sense
\begin{equation}\label{m3b}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
\displaystyle{- ( r h_r )_r +  {{\sin  (2 h)}\over {2{r}}}
-{\mu}\,r\, {\sin (2 h)\cos (2 h)} =0} \\ \\
\displaystyle{ \Longrightarrow  - r h_{rr} -h_r +  {{\sin  (2 h)}\over {2{r}}}
-{\mu}\,r\, {\sin (2 h)\cos (2 h)} =0~.}
\end{array}
\end{equation}
Dividing by $r$ we get the first equation of \eqref{m1b}. Integrating by parts in \eqref{m2b}
and using \eqref{m3b} we get
\begin{equation*}
\displaystyle{  \int_0^1 { ( r h_r v )_r - ( r h_r)_r v+  {{\sin  2 h v}\over {2{r}}}
-{\mu} {\sin (2 h)\cos 2 h r v}} =0~~~, ~~\forall v \in V }
\end{equation*}
i.e.
\begin{equation*}
\displaystyle{  \int_0^1 { ( r h_r v )_r } =0~~~, ~~\forall v \in V }
\end{equation*}
which gives
$$ h_r(1)=0~. $$
(in a weak sense) $h(0)=0$ follows from $h \in V$. This completes the proof of the Lemma.

 \hfill $\Box$  \end{proof}

 \begin{lemma}\label{p2b}
 If $h\neq 0$ is a nonnegative minimiser of $E$ on $V$ then $h>0$ on $(0,1)$.
\end{lemma}
\begin{proof}
Indeed, if $h$ vanishes at $r_0\in(0,1)$ then, since $h$ is smooth and $r_0$ is a minimum
for $h$, one would have
$$ h(r_0)=  h_r(r_0)=0$$
then from the theory of o.d.e's (see [1]), $h\equiv 0.$
\hfill $\Box$  \end{proof}

 \begin{lemma}\label{4.8}
 If $h$ is a positive minimiser of $E$  then $0<h\le {\pi\over 2}$.
\end{lemma}
\begin{proof}
If not then $h$ constructed as in the figure before (Fig.1) is a minimiser but it has a jump in the derivative unless this one is $0$. But then $h= {\pi\over 2}$ is solution of the o.d.e. on $h >{\pi\over 2}$  and a contradiction follows. Note that the solution of the elliptic equation (\ref{m1}) is smooth on $(0,1)$.
\end{proof}

 \begin{lemma}\label{4.9b}
A minimiser cannot vanish  on $(0,1)$ unless it vanishes identically.
\end{lemma}
\begin{proof}
If  $h$ is  a minimiser,  $\vert h \vert $ is also a minimiser. But, then, $\vert h \vert $ would
have a jump discontinuity in its derivative unless when it vanishes so does $h_r$.
This implies (theory of o.d.e's),  $h=0$.
\end{proof}

\rightline \hfill $\Box$

\begin{lemma}\label{4.9}
If $h \in V$ then   $\sin (kh) \in V ~~, ~~ \forall k\in{{\rm I}\!{\rm R}}$.
\end{lemma}
\begin{proof}
One has
\begin{equation}
\sin (kh)_r= k h_r\cos (kh)~~,~~\vert \sin (kh)\vert \le \vert kh \vert.
\end{equation}
Therefore one has
\begin{equation}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
\displaystyle
\Vert \sin (kh) \Vert^2 =  \int_0^1 \left\{  \sin (kh)_r^2 + \left( { \sin (kh) \over r}\right)^2 \right\} r dr\\ \\
\displaystyle{ \le  \int_0^1 \left\{ k^2 \cos (kh)^2  h_r^2+ k^2 {h^2 \over {r^2}} \right\} r dr \le  k^2 \Vert h\Vert ^2 }
\end{array}
\end{equation}
\end{proof}
\bigskip

\noindent It easy to check that $h\equiv 0$ solves (\ref{bch}),
(\ref{eqh}) and hence it is a stationary point of the functional
(\ref{Eh}). \\
Let $\gamma_0$ be the first eigenvalue of the problem
\begin{equation} \label{eigen}
\left \{ \begin {array}{l}\displaystyle{ -\phi_{rr}
-\frac{\phi_r}{r}
+\frac{\phi}{r^2} = \gamma \phi }\\
\\
\phi(0) =0, \quad \phi_r(1)=0 ~.
\end{array}\right.
\end{equation}
\begin{lemma}\label{gamma-0}
$$\gamma_0>1~~.$$
\end{lemma}
\begin{proof}
Suppose not,  i.e. $\gamma_0\le1$.
Let $\phi$ be the corresponding positive (or nonnegative) eigenfunction.
One has
\begin{equation}
 \begin {array}{l} \displaystyle{ -\phi_{rr}
-\frac{\phi_r}{r}= \phi (\gamma_0- {1\over {r^2}}) \le 0
~~~\text{since}~r\in (0,1)}\\ \\
\displaystyle{ (r \phi_r)_r \ge 0 \Longrightarrow r \phi_r \nearrow ~~
\Longrightarrow ~~r \phi_r \le 0 ~~\text{since}~~   \phi_r (1)=0~.}
\end{array}
\end{equation}
Thus, the maximum of $\phi$ is achieved at $0$ but, since $\phi(0)=0$,
we get a contradiction i.e. $\phi\equiv0$.
\hfill\hfill $\Box$
\end{proof}

\medskip\noindent
We have the following bifurcation lemma.

\begin{lemma}\label{biforcation}
If $\mu \le \gamma_0/2 $ we have $E(h) \ge 0$ and the global minimum
is attained  only for $h\equiv 0$. For $\mu > \gamma_0/2 $ the
global minimum is negative.
\end{lemma}
\begin{proof} The first equation of \eqref{eigen}  can also be written after a multiplication by $r$ as
\begin{equation*}
\displaystyle{ -(r\phi_r)_{r}+\frac{\phi}{r} = \gamma \phi r}.
\end{equation*}
Multiplying by $\phi$ and integrating over $(0,1)$ we derive by definition of $\gamma_0$ that
\begin{equation}\label{pointca}
\int_0^1 \left(\phi_r^2 + \frac{\phi^2}{r^2} \right)rdr \geq \gamma_0 \int_0^1  \phi^2 rdr ~~\forall \phi ~\text{ with }~\phi(0) =0, \quad \phi_r(1)=0.
\end{equation}

We divide the proof in two parts:

\bigskip
\noindent ({\it i} ) \quad $\mu \le \gamma_0/2$ \\

In this case we have (using \eqref{pointca} with $ \phi = \sin h$)
\begin{equation*}
\begin{aligned}
E(h) &=\pi \int_0^1
\left[(\cos h)^2 h_r^2+\left(\frac{\sin h}{r}\right)^2 - 2\mu (\sin h )^2(\cos h)^2 + (1-(\cos h)^2) h_r^2 \right] rdr \\
&\geq \int_0^1
\left[\gamma_0 (\sin h )^2 - 2\mu (\sin h )^2(\cos h)^2 + (1-(\cos h)^2) h_r^2 \right] rdr \\
& = \int_0^1
\left[ (\gamma_0  - 2\mu ))(\sin h )^2 + (1-(\cos h)^2) (2\mu (\sin h )^2 + h_r^2) \right] rdr \\
&\geq 0 =E(0)
\end{aligned}
\end{equation*}
the equality taking place only for $h=0$.

\bigskip
\noindent({\it ii }) \quad $\mu > \gamma_0/2$ \\

Let us denote by $\phi_0$ the first positive normalised eigenfuntion to \eqref{eigen}.

One has for $\epsilon >0$
\begin{equation*}
\begin{aligned}
E(\epsilon \phi_0) &\leq \pi \int_0^1
\left[(\epsilon \phi_0)_r^2+\frac{(\epsilon \phi_0)^2}{r^2} - \frac{\mu}{2} (\sin (2\epsilon \phi_0))^2\right] rdr \\
&= \pi \int_0^1
\left[\gamma_0(\epsilon \phi_0)^2 - \frac{\mu}{2} (\sin (2\epsilon \phi_0))^2\right] rdr. \\
\end{aligned}
\end{equation*}
Using with $x= 2\epsilon \phi_0$ the formula
\begin{equation*}
\sin x = x - \int_0^1 (1 - \cos (tx)) xdt
\end{equation*}
$E(\epsilon \phi_0) $ can be written as
\begin{equation*}
\begin{aligned}
E(\epsilon \phi_0)
&= \pi \int_0^1
\left[(\epsilon \phi_0)^2 \{\gamma_0- 2\mu (1-\int_0^1 (1 - \cos (2t\epsilon \phi_0))dt )^2\}\right] rdr \\
&< 0 =E(0)
\end{aligned}
\end{equation*}
for $\epsilon$ small since
\begin{equation*}
\int_0^1 (1 - \cos (2t\epsilon \phi_0))dt \to 0
\end{equation*}
when $\epsilon \to 0$.
 \hfill $\Box$  \end{proof}

\bigskip

\noindent
{\bf Alternative proof of $(ii)$
}

\medskip\noindent
Suppose $h\neq 0$ is a minimiser of $E$ one has
\begin{equation}
E(h)< E(0)
\end{equation}
i.e.
\begin{equation}
\begin{array}{cl@{\hspace{0.5ex}}c@{\hspace{1.0ex}}l}
\displaystyle \int_0^1
\left\{  h_r^2+\left(\frac{\sin h}{r}\right)^2 \right\} rdr < {\mu \over 2} \int_0^1
 \left({\sin 2 h}\right)^2  rdr = 2\mu \int_0^1 {\sin  h}^2 {\cos  h}^2 rdr  \\
\\
\displaystyle \Longrightarrow~~  \int_0^1
\left\{  {\cos  h}^2 h_r^2+\left(\lambda\frac{ \sin h}{r}\right)^2 \right\}rdr <\gamma_0
 \int_0^1 {\sin  h}^2  rdr \\ \\
\displaystyle \Longrightarrow~~ \gamma_0 >
{ \int_0^1\left\{   {\sin  h}_r^2 +\left(\frac{\sin h}{r}\right)^2 \right\}  rdr \over
\int_0^1 {\sin  h}^2  rdr }
\end{array}
\end{equation}
and a contradiction since ${\sin  h}\in V$ with the definition of $\gamma_0$.

\medskip\noindent
Consider the problem
\begin{equation} \label{one}
\left \{ \begin {array}{l}\displaystyle{ -h_{rr}
-\frac{h_r}{r}
+\frac{\sin 2h}{r^2} = \mu \sin 2h \cos 2h } ~~~\text{on}~ (0,1)\\
\\
h(0) =h_r(1)=0
\end{array}\right.
\end{equation}

\noindent
\begin{lemma}
If $\mu\le \gamma_0/2$ the only solution of (\ref{one})
such that $ \displaystyle{h\in \left[ -{\pi\over 2}, {\pi\over 2}\right]}$ is $h\equiv0$.
\end{lemma}

\noindent
\begin{proof}
Recall that for any $\phi \in V$ one has by definition of  $\gamma_0$
\begin{equation}\label{due}
\displaystyle{\gamma_0 \int_0^1 \phi^2 r dr \le
 \int_0^1 \left(\phi_r^2 + {\phi^2\over{r^2}}\right) r dr } ~~.
\end{equation}
Let us write the
equation  (\ref{one}) as
\begin{equation}
({r h_r})_{r} +\frac{\sin 2h}{2 r} = \mu \sin 2h \cos 2h \, r~.
\end{equation}
Multiply both sides by $\sin 2h$ and integrate on $(0,1)$. It comes
\begin{equation}
\displaystyle{\int_0^1 r \left\{ h_r  ({\sin 2h})_{r} +\frac{(\sin 2h)^2}{2 r^2} \right\} dr=
\mu \int_0^1 (\sin 2h)^2 \cos 2h ~ r dr} ~.
\end{equation}
One has
\begin{equation}
\displaystyle{({\sin 2h})_{r} =2 \cos 2h ~ h_r \Longleftrightarrow
h_r= {({\sin 2h})_{r} \over {2 \cos 2h}}}~.
\end{equation}
Thus, the equation above becomes
\begin{equation}
\displaystyle{\int_0^1 r \left\{   ({\sin 2h})^2_{r} {1\over {\cos 2h}} + \frac{(\sin 2h)^2}{ r^2} \right\} dr=
2 \mu \int_0^1 (\sin 2h)^2 \cos 2h r dr} ~.
\end{equation}
Suppose that $h$ is such that  $ \displaystyle{h\in \left[ -{\pi\over 2}, {\pi\over 2}\right]}$
 then since $-1\le \cos 2h\le 1$ one gets
 \begin{equation}
\displaystyle{\int_0^1 r \left\{   ({\sin 2h})^2_{r}  + \frac{(\sin 2h)^2}{ r^2} \right\} dr<
\gamma_0 \int_0^1 (\sin 2h)^2  r dr} ~
\end{equation}
i.e $\sin 2h \in V$ and satisfies an inequality contradicting (\ref{due}),
except if $h\equiv 0$.
\end{proof} \hfill\hfill $\Box$

\bigskip
\noindent Each minimiser of $E(h)$  solves the problem (\ref{bch}),
(\ref{eqh}). For the solutions of this problem  we can give the
following existence result around the bifurcation point.

\begin{lemma} \label{existence}
There exist two positive numbers $\rho_0$ and $\delta_0$ such that,
the problem (\ref{bch}), (\ref{eqh}) does not have non-zero
solutions for $\mu \in (\gamma_0/2 - \delta_0, \gamma_0/2]$ and
$\|h\|_0 \le \rho_0$. The problem has exactly two solutions $h_1$
and $h_2 = - h_1$ in the sphere $\|h\|_0 \le \rho_0$ for $\mu \in
(\gamma_0/2, \gamma_0/2 + \delta_0)$.
\end{lemma}
\begin{proof}
The proof follows from \cite[Theorem 6.12]{K}. Indeed the equation
(\ref{eqh})  can be written in the form
\begin{equation}\label{eqK}
 2\mu h =  L(h, r) + C(h, r, \mu) + D(h, r, \mu)
\end{equation}
 where $L$ is the linear operator
 $$L(h,r)=-h_{rr} -\frac{h_r}{r} +\frac{h}{r^2}$$
 and $C$, $D$  are given by
 $$C(h, r, \mu)=-\frac{2}{3}\frac{h^3}{r^2}+\frac{16}{3} \mu {h^3},$$
$$D(h, r, \mu)=-(\frac{2h}{2r^2} -\frac{\sin2h}{2r^2}) + \frac{2}{3}\frac{h^3}{r^2}
+\frac{\mu }{2}( 4h - \sin 4h) - \frac{16}{3} \mu {h^3}.
$$

\noindent It is easy to check that
\begin{equation} \label{omo}
C(th,r,\mu) = t^3 C(h,r,\mu),\quad (-\infty < t < \infty)
\end{equation}
and
\begin{equation}\label{D}
\|D(h,r,\mu)\|_0  = o (\| h\|^3).\end{equation} Moreover we
have
\begin{equation}\label{CF}
\displaystyle{((C(\phi^0,r,\mu),\phi^0))_0 =
\int_0^1 \left[-\frac{2}{3}\frac{(\phi^0)^4}{r}+\frac{16}{3} \mu
{(\phi^0)^4} r\right] dr}
> 0,\quad {\rm for}\,\, \mu\ge \frac{\gamma_0}{8}.
\end{equation}
Indeed from $((L(\phi^0,r)-\gamma_0 \phi^0, (\phi^0)^3))_0 = 0$ it
follows that
\begin{equation}
\displaystyle{\int_0^1 \left[{-\frac{d}{dr}{{(\phi^0_r r)}}(\phi^0)^3}+ {{(\phi^0)^4}\over{r}}
-\gamma {(\phi^0)^4} r \right] dr =0}
\end{equation}
The latter, on by parts integration
\begin{equation*}
\displaystyle{-{(\phi^0_r r)}(\phi^0)^3 \big\vert_{0}^{1}
+ \int_0^1 3 {(\phi^0_r)^2}  {(\phi^0)^2 } r dr  +
 \int_0^1\left[{{(\phi^0)^4}\over{r}} -\gamma_0 {(\phi^0)^4} r \right] dr =0}
\end{equation*}
that is
$$ \int_0^1 \frac{(\phi^0)^4}{r} - \gamma_0 (\phi^0)^4r \, dr \le 0,$$
and the inequality (\ref{CF}) can be easily derived.

\noindent The statements (\ref{eqK})- (\ref{CF}), together to the
local Lipschitz condition on the operators C and D, assure (see
\cite{K} ) the existence of exactly two branch of non-zero solutions
bifurcating from the point $\gamma_0/2$. Finally we remark that the
existence of two opposite branch follows from the odd functions in
(\ref{eqK}).

 \hfill $\Box$  \end{proof}

\begin{remark}\label{stability}
In order to establish the stability of the solutions to (\ref{bch}),
(\ref{eqh}) around the point $\mu_0=\gamma_0/2$, we perform  a
qualitative analysis of   the bifurcation equation to the lowest
order (see equation (\ref{loweq}) below). From (\ref{eqK}) setting
$$G(h,r,\mu)= -2\mu h + L(h, r) + C(h, r, \mu) + D(h, r, \mu)= 0 $$
and
$$ 2 \mu = \gamma_0 + \delta, \qquad |\delta| << 1,$$
assuming that each element $h \in \mathcal{H}^1(0,1)$ has the unique
representation
$$ h= \beta \phi^0 + P h,\qquad (P h, \phi^0)_0 = 0,\quad \beta \in
\mathbb{R},$$ we have
$$(G(h, r,\mu), \phi^0)_0= - \delta \beta  +
(C(\beta \phi^0, r,\mu),\phi^0)_0 + \dots$$ Moreover from
(\ref{omo}), (\ref{CF}) we can get to the simple l.o. bifurcation
equation, namely
\begin{equation} \label{loweq} - \delta \beta +
\beta^3 \bar{C}=0,\qquad \bar{C} =(C(\phi^0,r,\mu),\phi^0)_0 \ge 0.
\end{equation}
It is easy to check that:\\
for $\delta \le 0$ there is the only solution $\beta=0$ and this
solution is stable (indeed in this case: $-\delta + 3 \beta^2
\bar{C} \ge 0)$;\\
for $\delta > 0$ the trivial solution is no more stable but other
two stable solutions appear, i.e. $\beta = \pm
\sqrt{\delta/\bar{C}}$.
\end{remark}
\subsection*{Acknowledgements}
Under the financial support of G.N.F.M.-I.N.d.A.M.,  I.N.F.N. and Universit\`a di Roma
 \textsc{La Sapienza}, Rome, Italy.
M. Chipot acknowledges  Dip. S.B.A.I., Universit\`a di Roma  \textsc{La Sapienza}, for the kind hospitality.

\bibliographystyle{plain}

\nocite{*}

\hfill


\section{Introduction}

In monolayer honeycomb materials, including  graphene \cite{gra}, Xenes  \cite{Molle17} (silicene \cite{Aufray10,Liu11}, germanene \cite{cha,Liu11}, stanene \cite{sta}),
transition metal dichalcogenides \cite{tmdc}, or bismuthene \cite{Reis17}, the Fermi level appears in two non-equivalent valleys of the Brillouin zone.
In graphene ribbons \cite{nrbr} with zigzag edges there is a strict correspondence between
the valley state and the direction of the current flow \cite{zz1,zz2,zz3}.
The valley protects the chiral \cite{chiral} electron transport from backscattering by long
range potential disorder, which led to the
valley electronics \cite{rycerz,ryre} or application of the valley degree of freedom to information processing \cite{ryre,vq1,vq2}.

In bilayer \cite{morpugo} and staggered monolayer graphene \cite{niu}  a topological confinement inside the sample, far from the edges, was  found \cite{morpugo,niu} along a zero line of inversion-symmetry-breaking potential.
For bilayer graphene \cite{bgr} this potential is introduced by an electric field perpendicular to the layers \cite{morpugo}. The bias opens the energy gap in the band structure \cite{bgr,wry}.
 For an inhomogeneous electric field that is inverted at the zero line
in space,  topologically protected chiral currents have been found \cite{morpugo}
 confined in space to a region of the electric field flip. The flip line provides a one-dimensional confinement or quantum-wire-type channels \cite{macdo,peet,muktu,muktu2}.
The one-dimensional confinement of currents is also found in the quantum Hall conditions
at the n-p junctions \cite{np} induced by electric fields. In contrast to the currents localized at the n-p interface \cite{rick,tay,snejk}, the confinement in zero-line channels does not require external magnetic field \cite{morpugo}. { Note, that formation of an energy gap was also found in epitaxial monolayer graphene
due to the sublattice symmetry breaking by the substrate \cite{sic}. However, the gap of this
origin cannot be easily inverted for the topological confinement.}

A perpendicular electric field opens the energy gap for the buckled Xenes monolayers \cite{Molle17,ni,Drummond12}, and an inversion of the field induces topological confinement similarly as
in bilayer graphene \cite{Ezawa12a}. In contrast to bilayer graphene,
(i) in Xenes the chiral energy level
that passes from the valence to the conduction band
is a linear function of $k_y$.
In Xenes a (ii) single topological state per current direction is
present instead of two as in bilayer graphene.
As a consequence of (i) and (ii) the wave packet formed at
zero line in a Xenes should be stable against excitations and should travel
with a constant shape due to independence of the velocity on the wave vector.

In this paper we study the dynamics of the chiral wave packets along the zero lines
of the electric field in silicene \cite{Molle17,rivi,chow,Ezawa,s1,s2,s3}, which is the most advanced Xenes material, with
a successful application for the field effect transistor \cite{Tao15}.
The chiral wave packets move with the Fermi velocity and the motion of the topological packets is free from
Zitterbewegung \cite{ci2,citer,frolo,naturd} which is characteristic to the solution of the Dirac equation
and was found also for Weyl fermions in graphene \cite{citer,rome,frolo}.
Moreover, the topological electron packets move with a constant shape, similarly to solitons,
that in other conditions require interaction
with the environment --  see the electron solitons self-focused
with interaction to the metal gates in heterostructures \cite{bed} or the
 Trojan wave packets formed by carefully prepared electromagnetic field \cite{tro1,tro2}.

We show that the chiral electron packets can be transferred from one valley to the other by backscattering from the edge of the crystal and we find that the packet appears in the opposite valley in a restored shape.
Potential profiles that act as beam splitters and quantum rings \cite{qr}
 are proposed.
An interference of the split parts of the wave packet can be controlled by external
magnetic field threading the quantum ring. The ring stores the wave packet for a time that is
a periodic function with the period the  flux quantum. The Aharonov-Bohm \cite{ab} conductance oscillations are also demonstrated  for the system with an intervalley scattering present.
The topological rings are a new form of quantum rings for Dirac electrons, with respect
to previously considered systems with structural confinement \cite{structure0,structure,structure2,structure3,structure4},
mass confinement \cite{massconf,massconf2,massconf3} or the
confinement at circular n-p junctions
in the quantum Hall conditions \cite{amk}.

{
In the long wave packet limit, i.e. when the width of the packet is comparable with the
diameter of the ring, the interference within the ring is translated into an imbalance of the
scattering density in the left and right arms of the ring that becomes a periodic function
with doubled Aharonov-Bohm period. The result is also found in the stationary transport
for silicene and bilayer graphene. }

This paper is organized as follows.
In Section II we provide an analytic solution to the Dirac equation for the
topologically protected currents at the inversion
of the electric field (II.A). In section II.B the solution of the continuum approximation (II.A)
is translated to atomistic tight-binding description which naturally accounts
for the intervalley scattering, and the time stepping procedure is explained.
In the Section  III we  test the stability of the wave packet
against the intervalley transition (III.A) and  describe the Aharonov-Bohm \cite{ab} oscillations
of the electron storage time by the ring (III.B). Section III.C describes the conductance oscillations for intervalley scattering within the ring.
{ The long wave packet limit is disussed in III.D, and the limit
 is confronted with the standard  calculation based on the stationary electron scattering in Section III.E. Section III.F shows the results of the Landauer approach for
quantum rings defined by the zero lines in bilayer graphene.
 Section IV contains the summary.}

\begin{figure}
\begin{tabular}{llll}
(a) &\includegraphics[width=0.5 \columnwidth]{b-crop.pdf} &
(b) &\includegraphics[width=0.3\columnwidth]{a-crop.pdf} \\
(c) &\includegraphics[height=0.6\columnwidth]{wynikff_dopub-crop.pdf} &
(d) &\includegraphics[height=0.6\columnwidth]{reldypxl4-crop.pdf
\end{tabular}
\caption{(a) Schematic side view of a silicene monolayer embedded in a dielectric sandwiched
between top and bottom gates with the A sublattice on top (red dots in the inset)  closer to a positively $V_G>0$ biased top gate (red color) gated inducing a negative potential energy  for $x<0$. (b) Top view of the system. The $K'$ ($K$) valley  current flows leaving the negative (positive) potential energy  on the A sublattice on the left hand side. The flake has a zigzag (armchair) termination at the edges with constant $y$ ($x$).
(c)  The potential profile on the A and B sublattices (right axis) and the wave function on the
A sublattice as calculated in the continuum approximation [Eq. (6)]. The blue solid (red dashed) lines
show the wave functions with the positive (negative) product of the valley index $\eta$
and the spin $\sigma_z$ eigenvalue for $\lambda=4$ nm and $V_g=0.2$ eV. (d) The dispersion relation
calculated numerically.
 \label{sze}}
\end{figure}

\begin{figure}
\includegraphics[width=.6\columnwidth]{odzzodaa.png}
\put(-50,225){(a)} \put(-50,100){(b)}
\caption{Time evolution of the absolute value of the electron packet ($|\phi|$) backscattered by the zigzag (a) and armchair (b)
edge of the flake. The square flake has a zigzag edge  at $y=300$ nm (a) and an armchair edge at $x=300$ nm (b). The plots show the absolute value of the wave function
at  the inversion of the energy band that is introduced at $x=0$ (a) and at $y=0$ (b).
For the zigzag edge the wave packet is temporarily localized at the edge.
The vertical fringes observed
at the backscattering  result from the superposition of the incoming and outgoing waves from the opposite valleys  $K$ and $K'$.
}
 \label{aring}
\end{figure}

\section{Theory}

\subsection{Chiral currents confined by the energy gap inversion}
The chiral currents at the band inversion in silicene similar to the ones in bilayer graphene  \cite{morpugo} were found in Ref. \cite{Ezawa12a}.
For buckled monolayers as well as for bilayer systems the energy gap can be tailored in space
using multiple split gates with inverted polarization. The idea of the local manipulation of the energy gap by dual split gates  was proposed for topological confinement \cite{morpugo}
and pseudospin electronics in bilayer systems \cite{re2a,re2b,re2c}.

Here, we consider a buckled silicene monolayer sandwiched in between top and bottom gates  [Fig. \ref{sze}(a,b)].
  The gates are split, so that the electric field changes orientation along the $y$ axis (i.e. for $x=0$).
We model the potential at the A sublattice using an arctangent function
\begin{equation}
 V_A(x)=\frac{2V_g \arctan(x/\lambda)}{\pi}. \label{pote}
\end{equation}
We assume that the silicene is embedded symmetrically between the gates, so that
on the B sublattice we have  $V_B(x)=-V_A(x)$  [Fig.\ref{sze}(c)].
The potential bias between the sublattices opens an energy gap in the band structure \cite{ni,Drummond12}.
For potential of Eq. (\ref{pote}) the energy gap is inverted at $x=0$ by the flip of the electric field orientation

For the wave function components defined on sublattices
$\psi=\left( \psi_A , \psi_B\right)^T$  the low-energy approximation to the atomistic tight-binding Hamiltonian  reads \cite{Ezawa12a}
\begin{equation} H_\eta = \hbar v_F \left(k_x \tau_x -\eta k_y \tau_y \right)+V({\bf r})\tau_z
 -\eta \tau_z \sigma_z  t_{SO},\label{hamu} \end{equation}
where $\eta$ is the valley index ($\eta=1$ for the $K$ valley and -1  for the $K'$ valley),
 $\tau_x$, $\tau_y$ and $\tau_z$ are the Pauli
matrices in the sublattice space, ${\bf k}=-i\nabla$, the Fermi velocity $v_F$ is determined
by  the  nearest distance between Si atoms
$d=2.25$ \AA\;  and the tight-binding hopping parameter $t=1.6$ eV \cite{Liu11,Ezawa}
 $v_F=3dt/2\hbar$.  In Eq. (\ref{hamu}) $t_{SO}=3.9$ meV is the intrinsic spin-orbit coupling constant \cite{Liu11,Ezawa}.
The intrinsic spin-orbit coupling is diagonal in the basis of eigenstates of the $z$ component of the spin and $\sigma_z=\pm 1$ is treated as a quantum number.

With potential given by Eq. (1) the Hamiltonian commutes with the $\frac{\partial}{\partial y}$ operator. The common eigenfunctions of the energy and $y$ momentum component can be put in form
$\psi^{k_y}(x,y)=\exp(ik_yy)\left(\begin{matrix} \psi_A^{k_y}(x) & \psi_B^{k_y}(x)\end{matrix}\right)^T$.
The Hamiltonian eigenfunctions
fulfill the
 system of equations
\begin{eqnarray} \frac{(V_A(x)-E-\eta\sigma_z t_{SO})\Psi_A(x)}{i\hbar v_F}&=&(\Psi_B'-\eta k_y \Psi_B), \\
\frac{(-V_A(x)-E+\eta\sigma_z t_{SO})\Psi_B(x)}{i\hbar v_F}&=&(\Psi_A'+\eta k_y \Psi_A),
\end{eqnarray}
where the prime stands for $x$ derivative.
 We plug in Eqs. (3,4) a relation $\Psi_B=i\Psi_A$, i.e.  a guess based on a numerical solution that allows us to derive an analytical solution for states localized near the zero-line  area.
A sum of the resulting equations relates the potential and the wave function  \begin{equation} V_A(x)-\eta \sigma_z t_{SO}=-\hbar v_F \frac{\partial}{\partial x}\ln \Psi_A(x).\end{equation}
For the specific form of potential given by Eq. (1) the (unnormalized) wave function is found by a standard integration technique
\begin{equation}\Psi_A=(\lambda^2+x^2)^{\frac{\lambda V_g}{\pi\hbar v_F}}e^{(-\frac{2V_g x}{\pi \hbar v_F}\arctan\frac{x}{\lambda}+\frac{\eta\sigma_z t_{SO} x}{\hbar v_F})}.\end{equation}
The term $x\arctan(x/\lambda)$ in the exponent for low $x$ introduces a Gaussian-like dependence
to the wave function and keeps it localized near the band inversion area [Fig. \ref{sze}(c)].
The average value of $\langle x ^2\rangle^{1/2}$ for $\lambda= 1$, 2, 4, 12, and 20 nm is
2.45, 2.87, 3.55, 5.45 and 6.83 nm, respectively.
In presence of the spin-orbit coupling the wave function is not ideally symmetric
with respect to the center of the gap inversion line: the term including
the spin-orbit coupling shifts the wave function at left (right) of the inversion line for
negative (positive) product of valley index $\eta$ and the $\sigma_z$ eigenvalue.
For the applied parameters these shifts are not very strong -- see the blue and red lines in Fig. \ref{sze}(c).
 In the calculations below we set $\sigma_z=1$.

\begin{figure}
\begin{tabular}{ll}
(a) &\includegraphics[height=0.4\columnwidth]{od_zigzaga-crop.pdf}  \\
(b) &\includegraphics[height=0.4\columnwidth]{kakaprim-crop.pdf
\end{tabular}
\caption{(a) The average position in the $y$ (black line, left axis) direction for the wave packet
scattered by the armchair edge of the flake. The spread of the wave packet $\sqrt{\langle (y-\langle y \rangle )^2\rangle}$ multiplied by 4 (blue line, left axis), and the average position at the $x$ axis
(right axis, red line). The variation of $\langle x\rangle $ are due to the spin-orbit coupling and the valley transitions that follow the backscattering by the edge [see. Eq.(6)]. The electron spin is set at $\sigma_z=+1$.
The parameters of the system are the same as in Fig. 2(a), only the packet is started with the initial
average position at $y_0=100$ nm.
(b) The Fourier transform of the packet calculated for the $K'$ (red) and $K$ valley (black)
} \label{wyna}
\end{figure}

The energy of the states localized at the gap inversion can be calculated by adding Eqs. (3) and (4) still with the relation $\Psi_B=i\Psi_A$, which gives  \begin{equation} E=-\eta \hbar k_y v_F. \end{equation}
The entire dispersion relation calculated numerically with a finite difference approach \cite{sweep} is given in Fig. 1(d).
The linear band energy is independent of the $V_g$ or $\lambda$ which only affect
the transverse  wave function localization at the energy gap inversion line.
 Above the energy gap a continuous spectrum is found with a parabolic threshold as a function of $k_y$.
Near the zero energy only the localized reflectionless currents flow, and for $V_g=200$ meV,  the gap is wide enough to make the currents stable at room temperature. For the Fermi wave vector $k_F\simeq 0.25$ nm$^{-1}$ ($E_F\simeq 100$  meV) the continuum appears still 100 meV above the linear band.

The sign of the electron velocity within the linear band $v(k_y)=\frac{1}{\hbar}dE/dk_y=-\eta v_F$ depends
on the valley index.
 Hence, the transport at the inversion line is chiral, i.e. the electron states of the valley $K'$ ($K$) go up (down) along the electric field flip line [Fig. \ref{sze}(b)].
Generally, in the states localized along the field flip the current in  the $K'$ ($K$) valley flows with the negative (positive) potential at the A sublattice at the left-hand side.

\subsection{Chiral solitons in the atomistic tight-binding}
\subsubsection{Atomistic Hamiltonian}

Since the electron velocity in the linear chiral band is independent of $k_y$ the wave packet localized  at the flip of the electric field should move with an unchanged shape along the $x=0$ line.   The description of the electron wave packets stabilized by the valley degree of freedom calls for an approach that takes into account the intervalley scattering. A natural choice is the atomistic tight-binding approach.
The positions of the ions of the A sublattice
${\bf r}_{m}^A=m_1 {\bf a}_1+m_2 {\bf a}_2$  are generated
with the crystal lattice vectors
${\bf a}_1=a \left(\frac{1}{2},\frac{\sqrt{3}}{2},0\right)$
and ${\bf a}_2=a \left(1,0,0\right)$, where $a=3.89$ \AA\; is the silicene lattice constant, and $m_1$, $m_2$ are integers.
The B sublattice ions are generated by ${\bf r}_{m}^B={\bf r}_{m}^A+(0,d,\delta)$,
with the vertical shift of the sublattices $\delta=0.46$ \AA.
The coordinates of the center of the valleys in the reciprocal space
are ${\bf K}_\eta=\left(\frac{4\pi\eta}{3a},0\right)$ \cite{rivi}. The valleys for $\eta=1$ (-1) are
referred to as $K$ ($K'$).

We use
 the Hamiltonian \cite{Liu,Ezawa,chow}
\begin{eqnarray}
H_{TB}&=&-t\sum_{\langle m,l\rangle } p_{ml} c_{m}^\dagger c_{l}
+it_{SO} \sigma_z\sum_{\langle \langle m,l\rangle \rangle  } p_{ml} \nu_{ml} c^\dagger_{m} c_{l} \nonumber
\\ && +\sum_{m} V({\bf r}_m) c^\dagger_{m}c_{m}+\frac{g\mu_B B}{2}\sigma_z,  \label{hb0}
\end{eqnarray}
where $\langle m,l\rangle $ stands for the nearest neighbor ions,
$\langle\langle m,l\rangle\rangle $ for the  next-nearest-neighbor ions.
For the potential $V({\bf r}_m)$ we take $V_A({\bf r}_m)$  or $V_B({\bf r}_m)$.
The sign  $\nu_{ml}=\pm 1$ is  plus (minus)  for   the next nearest neighbor hopping path
 via the common neighbor ion that turns
 counterclockwise (clockwise). In Eq. (8) $p_{ml}$ is the Peierls phase
 $p_{ml}=e^{i\frac{e}{\hbar}\int_{\vec{r_m}}^{\vec{r_l}}\vec {A}\cdot \vec {dl}}$,
where $\vec{A}=(0,B x,0)$ is the vector potential, and $B$ is the value of the magnetic
field that is oriented perpendicular to the silicene plane. The last term in Eq. (8) is the spin Zeeman term with the Land\'e factor $g=2$, and Bohr magneton $\mu_B$.

\subsubsection{Initial condition and the time-stepping}

In the calculations to follow for the initial condition
we use  the solution of the continuum Hamiltonian [Eq. (2)] and localize
the packet along the band inversion using an envelope of form $\frac{1}{1+\frac{(y-y_0)^4}{D^4}}$,
with $D=80$ nm and $y_0$ that sets the center of the packet. We set the valley momentum with a plane wave and the $K$ or $K'$ coordinates.
Accordingly, for the atoms of the A sublattice we set as the initial condition
\begin{equation}
 \psi({\bf r}^A_m,t=0)=\exp\left(i ({\bf K_\eta}+{\bf k}_y) \cdot {\bf r}_m^A\right)\chi(y_m^A), \end{equation}
where ${\bf k}_y=(0,k_y,0)$ is the wave vector of the packet calculated with respect to the valley center, and $\Psi_A$ is given by Eq. (6).
For the atoms on the $B$ sublattice we take
\begin{equation}
 \psi({\bf r}^B_m,t=0)=i\eta \exp\left(i ({\bf K_\eta}+{\bf k}_y) \cdot {\bf r}_m^B\right)\chi(y_m^B), \end{equation}
where $\chi(y)=\frac{\Psi_A(y)}{1+\frac{(y-y_0)^4}{D^4}}$.
We set $k_y=0$ unless stated otherwise.

We solve the Schr\"odinger equation on the atomic lattice
$i\hbar \frac{\partial \psi} {\partial dt}=H\psi$, using the time step of $dt=10$ atomic units
or $dt=2.418\times 10^{-4}$ ps.
The wave function at the first step is calculated with the
implicit Crank-Nicolson scheme $\psi(dt)=\psi(0)+\frac{dt}{2i\hbar}H_{TB}\left(\psi(0)+\psi(dt)\right).$ The subsequent time steps are calculated with the
explicit Askar-Cakmack scheme $\psi(t+dt)=\psi(t-dt)+\frac{2dt}{i\hbar}H_{TB}\psi(t)$.
In presence of the external magnetic field the eigenstates of (3,4) need to be calculated numerically \cite{sweep}.
However, for the discussed range of the magnetic field ($B <1$ T) and the applied narrow flip area ($\lambda=4$ nm) no significant difference between the numerical eigenstates and formula of Eq. (6) used for the initial condition were found
in the wave packet evolution.

\section{Results}
\subsection{Wave packet motion}
We first consider a square flake $(x,y)\in[-300$ nm$,300$ nm$]\times [-300$ nm$,300$ nm$]$
with a zigzag edge at $y=\pm 300$ nm and and armchair edge at $x=\pm 300$ nm.
The packet is set in the $K'$ valley  ($\eta=-1$)
to make it move upwards (to increasing $y$ values -- see Fig. 1(b)).
Figure 2(a) shows the cross section of the packet along the $x=0$ line.
The packet indeed moves up in a stable form until it reaches the zigzag
edge of the flake.

The same result -- as long as the packets does not reach the edge -- is obtained for the solution of the time dynamics with the continuum Hamiltonian [by Eq.(2)].
The absence of Zitterbewegung \cite{ci2,citer,frolo,naturd} for the wave packet that follows the Dirac equation calls for a comment.
The velocity operator is obtained as ${\hat{v}_y}=\frac{1}{\hbar}\frac{\partial{H_\eta}}{\partial k_y}=- \eta \tau_y v_f $ \cite{citer}. For the Dirac equation this operator does not commute with the Hamiltonian  $[H_\eta,v_y]\neq 0$ which is usually invoked in the interpretation of the trembling motion \cite{citer}. However,
 the wave function localized at the
zero-line that follows the specific form used in Section II.A $\Psi=(\Psi_A, i\Psi_A)^T$ happens to be an eigenfunction of the $\hat{v}_y$ operator, with the eigenvalue $v_y=-\eta v_f$. By, the Ehrenfest theorem for we have $\frac{d}{dt}\langle \Psi |{\hat{v}_y}|\Psi\rangle=\frac{1}{i\hbar}\langle \Psi| {\hat{v}_y}H-H {\hat{v}_y}| \Psi\rangle=
\frac{1}{i\hbar}\left(\langle H \hat {v_y}\Psi|\Psi\rangle-\langle \Psi|H \hat{v_y}\Psi\rangle\right)=\frac{1}{i\hbar}\left(-\eta v_f \langle H\Psi|\Psi\rangle+\eta v_f \langle \Psi|H \Psi\rangle\right)=0$.

\subsection{Backscattering by the flake edge}
The zigzag edge of the crystal supports edge localized states. The incident packet couples to these states and in Fig. 2(a) we find formation
of a high peak of the absolute value of the wave function (see the red region at the right edge of the plot for $t\simeq$ 0.4 ps). The packet is backscattered and moves to the left with a shape restored to its original form, only in the opposite direction.

The average position $\langle y\rangle$, $\langle x\rangle$, and the size of the packet along the junction $\Delta y\equiv \sqrt{\langle (y-\langle y\rangle)^2 \rangle}$
is plotted in Fig. 3(a). The system is the same as in Fig. 2(a) only the packet
is started at $y_0=100$ nm.
In Fig. \ref{wyna}(b) we additionally plot
the square of the absolute value of the Fourier transform of the wave packet calculated
for $K$ and $K'$ valleys,
\begin{equation}
\psi({\bf k},t)=\frac{1}{2\pi} \int \psi(x,y,t) \exp(-i{\bf k}\cdot {\bf r}) {\bf dr},
\end{equation}
for ${\bf k}={\bf K}_{-1}$ (i.e. the $K'$ valley) and ${\bf k}={\bf K}_1$ ($K$ valley).
We find that as the packet reaches the edge of the flake it is scattered to the other valley
that makes it travel in the opposite direction.
In Fig. 3(a) we can also see that the average $\langle x \rangle$
oscillates with a small magnitude which results from the valley flips
at the backscattering by the edge which
changes the sign of the exponent with $t_{SO}$ in the wave function ($\sigma_z$ is set to 1).

We find that the restoration of the soliton upon scattering is observed for both the armchair
and the zigzag termination of the flake.
For backscattering by the armchair edge  we flipped the potential by $-\pi/2$ angle. The packet which is still started in the $K'$ valley moves toward larger $x$ values
until it reaches the edge [Fig. 2(b)]. Here no peak of the wave function at the edge
is observed since the armchair termination does not support the edge localized states \cite{Wakabayashi2010}.

\begin{figure}
\includegraphics[width=.5\columnwidth]{rysy-1.png}
\caption{The potential at the A sublattice for a quantum ring of radius 100 nm.
The potential at the B sublattice is assumed opposite. The arrows indicate the
orientation of the $K'$ currents. The currents in the $K$ valley flow in the opposite direction.} \label{ring}
\end{figure}

\begin{figure*}
\begin{tabular}{lllllll}
\includegraphics[width=0.24\columnwidth]{oo-7.png}\put(-58,51){0.084ps} \put(-75,5){(a)} \put(-85,30){\rotatebox{90}{$B=0$}}&
\includegraphics[width=0.24\columnwidth]{oo-17.png}\put(-58,51){0.204ps}&
\includegraphics[width=0.24\columnwidth]{oo-27.png}\put(-58,51){0.326}&
\includegraphics[width=0.24\columnwidth]{oo-37.png}\put(-58,51){0.447}&
\includegraphics[width=0.24\columnwidth]{oo-47.png}\put(-58,51){0.568}&
\includegraphics[width=0.24\columnwidth]{o57.png}\put(-58,51){0.689}&
\includegraphics[width=0.24\columnwidth]{oo-67.png}\put(-58,51){0.81} \\
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp30.png}\put(-58,51){0.36ps}
\put(-75,5){(b)} \put(-85,10){\rotatebox{90}{$B=0.066$T}}&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp50.png}\put(-58,51){0.6ps}&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp60.png}\put(-58,51){0.72}&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp77.png}\put(-58,51){0.924}&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp88.png}\put(-58,51){1.056}&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp98.png}\put(-58,51){1.176}&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{pp108.png}\put(-58,51){1.296} \\

\end{tabular}
\caption{The snapshots of the time evolution of the wave packet for $B=0$ (a) and $B=0.066$ T (b).
The magnetic field used in (b) corresponds to half the flux quantum threading the ring of radius 100 nm.
The time from the start of the simulation is given in picoseconds in the frames.
} \label{snap}
\end{figure*}

\begin{figure}
\begin{tabular}{l}
\includegraphics[width=0.7\columnwidth]{wyniki_z_kolnierzem-pub3-crop.pdf}\put(1,100){(a)} \\
\includegraphics[width=0.7\columnwidth]{czasy_do_pub-crop.pdf}\put(1,100){(b)}\\
\includegraphics[width=0.7\columnwidth]{w_dopub-crop.pdf}\put(1,100){(c)}
\end{tabular}
\caption{(a) The part of the wave packet before the ring (black lines),
inside the ring (blue lines) and above the ring (red lines)
for zero magnetic field (solid lines), half of the flux quantum ($B=0.066$T)
and the magnetic field corresponding to the flux equal to $\Phi_0/8.5$ ($B=0.0155$ T).
(b) The time at which less than 10\% (blue line), 50\% (red line), and less than 75\% of the electron density stays within the ring as a function of the external magnetic field. Local extrema	 of all the three
times are found for integer and half the flux quanta which correspond to even and odd multiples
of $0.066$ T, respectively.
(c) The part of the wave packet inside the ring for the selected moments in time from the start of the simulation.
} \label{inzain}
\end{figure}

\begin{figure}
\includegraphics[width=0.9\columnwidth]{backscatter.pdf}
\caption{The backscattering probability for a sharp perturbation to the potential
given by Eq. (13). The red, green  and blue lines correspond to $k_y=0$, $0.1$/nm and $0.15$/nm, respectively. For each value of $k_y$ two lines are plotted: the lower and upper bound to the
backscattering probability. The spacing between these two lines is given by the probability density
localized within the ring at the end of the simulation (6.1 ps from the start).  } \label{R}
\end{figure}

\begin{figure*}
\begin{tabular}{llllll}
\includegraphics[width=0.24\columnwidth]{l13.png}\put(-58,51){0.145ps} &
\includegraphics[width=0.24\columnwidth]{l22.png}\put(-58,51){0.266ps}\put(-42,7){$\leftarrow$}\put(-28,7){$\rightarrow$}\put(-42,.5) {\small K'}&
\includegraphics[width=0.24\columnwidth]{l32.png}\put(-58,51){0.386}
\put(-52,12){$\nwarrow$}\put(-52,1) {\small K'}\put(-24,17){$\nearrow$}
&
\includegraphics[width=0.24\columnwidth]{l42.png}\put(-58,51){0.508}
\put(-56,26){$\uparrow$}\put(-58,16) {\small K'}\put(-21,45){$\nwarrow$}
\put(-24,17){$\nearrow$} \color{teal} \put(-14,13){$\swarrow$} \put(-10,4) {K} &
\includegraphics[width=0.24\columnwidth]{l52.png}\put(-58,51){0.628}
\put(-44,39){$\nearrow$}\put(-21,45){$\nwarrow$}
\put(-29,52){$\uparrow$}
\put(-29,42){$\downarrow$}
 \color{teal} \put(-14,13){$\swarrow$} \put(-10,4) {K}
\put(-29,3){$\downarrow$}
\put(-29,15){$\uparrow$}
&
\includegraphics[width=0.24\columnwidth]{l62.png}\put(-58,51){0.749}

\put(-29,52){$\uparrow$}
\put(-29,35){$\downarrow$}
 \color{teal} \put(-24,5){$\leftarrow$}
\put(-29,3){$\downarrow$}
\put(-29,22){$\uparrow$}

\\
\includegraphics[width=0.24\columnwidth]{l72.png}\put(-58,51){0.870}

 \put(-29,35){$\downarrow$}\put(-29,52){$\uparrow$}\color{teal} \put(-30,3){$\downarrow$}\put(-29,22){$\uparrow$}
&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{l82.png}\put(-58,51){0.991ps}

\put(-42,7){$\leftarrow$}\put(-27,7){$\rightarrow$}
\color{teal}
\put(-29,37){$\uparrow$}
\put(-20,45){$\searrow$}
\put(-45,38){$\swarrow$}
&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{l92.png}\put(-58,51){1.11ps}
\put(-56,26){$\uparrow$}\put(-6,26){$\uparrow$}
\color{teal}\put(-47,26){$\downarrow$}\put(-19,26){$\downarrow$}
&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{l102.png}\put(-58,51){1.23}
\put(-44,39){$\nearrow$}\put(-18,42){$\nwarrow$}\put(-29,52){$\uparrow$} \put(-29,42){$\downarrow$}\
\color{teal}\put(-52,12){$\searrow$} \put(-30,3){$\downarrow$}
&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{l112.png}\put(-58,51){1.35}
\put(-29,52){$\uparrow$} \put(-29,42){$\downarrow$}\put(-18,42){$\nwarrow$}
\color{teal}\put(-52,12){$\searrow$}\put(-30,3){$\downarrow$}
&
\includegraphics[width=0.24\columnwidth, trim= 00 0 00 0,clip]{l122.png}\put(-58,51){1.47}

\put(-29,35){$\downarrow$}\put(-29,52){$\uparrow$}
 \color{teal}
\put(-29,3){$\downarrow$}
\put(-29,22){$\uparrow$}
\end{tabular}
\caption{Snapshots of the electron density for $B=0$, $k_y=0.1$/nm (the green lines in Fig. 7(a)).
The blue (green) arrows indicate the motion of the parts of the packet moving
in the $K (K')$ valley. The time from the start of the simulation is given in picoseconds in the frames. } \label{wiring}
\end{figure*}
\subsection{Quantum ring}
A form of a quantum ring with the chiral waveguides for the electron flow can be defined with an  engineering of the electric field. For that purpose
one needs a local inversion of the electric field  which requires additional top and bottom gates in the system. Here we consider a circular area of radius $R$,
and set the model potential at the $A$ sublattice to
\begin{align}
V_A=  \frac{4V_g}{\pi^2}\arctan \left(\frac{r-R}{\lambda}\right)\arctan \left(\frac{x}{\lambda}\right),
\end{align}
where $R=100$ nm is the circle radius with the center at the origin and $r=\sqrt{x^2+y^2}$ is the distance from the origin. The potential at the $A$ sublattice is plotted in Fig. \ref{ring}. As above we take $V_B({\bf r})=-V_A({\bf r})$.
Fig. \ref{ring} shows also the direction of the flow for currents in the $K'$ valley -- with
the negative (positive) potential at left (right) side of the current flow.
When the line of $V_A=0$ meets the ring at $y=-100$ nm the $K'$ current can go to either
the left or the right arm of the ring. This potential profile introduces  a beam splitter
for the electron wave packet in this way.  The central bar of the ring is inaccessible for the $K'$ current going up.

In quantum rings \cite{qr} Aharonov-Bohm \cite{ab} oscillations of coherent conductance
are observed that in the Landauer-B\"uttiker \cite{labu} approach are explained as due to variation of the electron transfer probability
across the ring with phase shifts accumulated from the vector potential of the magnetic field.
In the present system the electron backscattering from the ring is prohibited by the topological protection of the  valley current, so the transfer probability of the soliton packet is 1 for any magnetic field.
However, we find that the time that the electron spends within the ring
changes due to the phase shifts introduced by the vector potential.

In this section and in the rest of the paper we neglect the intrinsic spin-orbit coupling ($t_{SO}=0$) that introduces a weak asymmetry in the electron injection to the ring due to the spin-valley dependent shift off the zero line [cf. Fig. 3(a) and Eq. (6)]. { The calculations for the quantum ring
require long leads to prevent return of the packet to the scattering area upon reflection
from the edge of the crystal. For long wave packets the entire computational box is taken long up to 6 $\mu$m. Systems these long are treated with the scaling method of Ref. \cite{scaling}
for which the crystal constant is scaled as $a_s=a s_f$ with the hopping parameter $t_s=t/s_f$.
In Hamiltonian (8) $t_s$ replaces $t$, and $a_s$ replaces $a$ while the Si ions are generated in the computational box. We use the scaling factor $s_f=3$ or $4$ in the calculations for silicene below. }

Figure 5(a) shows the snapshots of the simulation of the packet transfer across the ring
for $B=0$. For $t=0$ the packet is started 350 nm below the center of the ring of radius 100 nm.
The snapshots taken at 0.204 ps and 0.325 ps show that the packet is split into two parts
at the entrance to the ring.
In both the left and right arms of the ring the $K'$ current moves up leaving the
negative potential at the A sublattice at the left-hand side [see Fig. 4].
The split packets meet at the exit of the ring [$t=0.568$ ps and $t=0.689$ ps]
with the same phase and the packet of its original size is restored [$t=0.810$ ps].

Figure 6(a) shows the parts of the packet before ("b4", black lines) the ring,
within ("in", blue lines), as well as the transferred part ("passed", red lines) as a function
of time. The results for $B=0$ that correspond to Fig. 5(a) are plotted with the
solid lines. For $t=0.75$ ps the entire packet is transferred above the ring.

In Fig. 5(b) we plotted the snapshots of the simulation for $B=0.066$ T
which corresponds to half the flux quantum $B\pi R^2=\Phi_0/2=\frac{h}{2e}$ threading the ring with $R=100$nm.
The parts of the packet that meet at the exit ($t=0.6$ ps) acquire a relative $\pi$ phase
due to the  Aharonov-Bohm effect.
A node of the wave function is formed at the exit to the upper channel.
A nondestructive interference  is observed within the area  below the exit from the ring,
which directs the packet to the internal bar within the ring ($t=0.72$ps), i.e. to the only path where the $K'$ current can go for the exit to the upper channel blocked by the Aharonov-Bohm effect. The bar is transparent for the $K'$ current going down [Fig.4]. The packet is split again to the left and right arms of the ring at $0.924$ ps.
After the second round the parts of the packet meet in phase $2\pi$ at the exit and the packet
smoothly leaves the ring [$t=1.176$ ps, $t=1.296$ ps].

In Fig. 6(a) the results for half the flux quantum are plotted with the dotted line.
The packet is transferred to the exit with a delay but completely and in a single move.
For comparison in Fig. 6(a) the results for the magnetic field of $0.0155$ T
 which corresponds to $\Phi_0/8.5$ are plotted  with the dashed lines. Here a part of the packet passes to the
exit as fast as for $B=0$, but due to a phase difference
at the exit a part of the packet stays inside the ring and leaves it in portions at
subsequent rounds, which produces the steps in the dashed lines in Fig. 6(a).

In Fig. 6(b) we plotted the moments in time for which less than 75\%, less than 50\% and
less than 10\% stays within the ring as a function of the external magnetic field.
The plot is periodic with the Aharonov-Bohm period of $0.132$ T.
The ring is emptied the fastest for the integer and half quanta. The result of Fig. 6(a)
for $B=0.0155$ T corresponds to a local maximum of the time for which more than 10\% of
the packet stays within the ring.
Finally Fig. 6(c) shows the part of the packet contained within the ring for a fixed moment
in time as a function of the magnetic field. For $t=1$ps and $1.25$ps the ring-confined part is locally maximal for the magnetic field corresponding to odd multiples of half the flux quantum
$B=(2n+1)0.066$ T. For longer times these maxima are turned into minima due to
compensation of the phase difference after the second round of electron circulation [Fig. 5(b)] within the ring.

\subsection{Intervalley scattering and conductance oscillations}
For the potential profile plotted in Fig. 4 the transfer probability can fall below 1
only provided that a intervalley transition is present within the ring.
 The intervalley scattering is introduced by potential
variation that is short on the atomistic scale. For the modeling we introduced a point potential defect
of a Lorentzian  form
{
\begin{equation} V_d({\bf r})=\frac{V_g}{1+{({\bf r}-{\bf r}_d)^2}/{l^2}}, \end{equation}
where ${\bf r}_d=(R,0,0)$, and $V_g=7.2/s_f^2$ eV and $l=12a=4.7$ nm. In the scaling approach \cite{scaling} the smooth potential that changes slowly on the atomic scale, in particular the one given by Eq. (1) and (12) stays invariant with $s_f$. For the abrupt short-range defect potential we found that scaling of $V_g$ with $s_f$ is necessary to keep the same effectiveness of the intervalley scattering.
 The defect potential $V_d$ enters with the same sign to potential on both sublattices,
as $V_A+V_d$ on  sublattice A and $V_B+V_d=-V_A+V_d$ on sublattice B. }

In presence of the defect, the results start to change significantly with $k_y$.
The backscattering probability  $R$ as a function of the external field is given in Fig. 7
for $k_y=0$, $0.1$/nm and $0.15$ /nm.
For each value of $k_y$ two lines are plotted, which are the minimal and maximal bound for the
backscattering probability. The difference between the two is determined by the part of the electron
packet that stays within the ring at the end of the simulation (6.1 ps).
$R$ as a function of $B$ is approximately periodic with the flux quantum
 The behavior of the electron packet is displayed in Fig. 8 for $B=0$ and $k_y=0.1$/nm. The packet is incident in the $K'$ valley. For $t=0.508$ ps we spot the scattering center
near ${\bf r}_d$. A part of the packet passes across the defect moving still in $K'$ valley, and
a larger part is backscattered and move in the direction which is only allowed for the $K$ valley.
For $t=0.628$ ps the $K$ current reaches the entrance to the ring, a part of it is backscattered
to the input channel, and another goes up along the bar. The electron packets of opposite valley
meet within the bar for $t=0.749$ ps and 0.87 ps.
For $t=0.991$ ps we can see recycling of currents for both the valleys: the $K$ valley current cannot pass to the output channel and the $K'$ one to the input channel.
The opposite valley currents meet again, this time in the arms of the ring for $t=1.11$ ps. For $t=1.47$ ps the packet distribution
 is similar to $t=0.87$ ps, only a smaller portion of the electron packet is still
present within the ring.

\subsection{The limit of long wave packets}

The wave packet dynamics in the limit of a large size of the packet should approach the conditions of the stationary electron flow. In order to study the limit of a plane incident wave we
consider a Gaussian envelope of the packet, i.e. in the initial condition given by Eqs. (9) and (10) we set $\chi=\exp(-(y-y_0)^2/4\sigma ^2)/(2\pi\sigma^2)^{1/4}$. The Fourier transform of the packet
produces the probability density with the standard deviation of $\sigma_k=\frac{1}{2\sigma}$ in the wave vector space.

We solve  the time evolution and integrate the parts of the packet in the left and right arms
of the ring and in the center bar over time as the packet transfers the ring.
 In Fig. \ref{lmit} we plot the results
for increasing length of the packet in the initial condition.
The results for $\sigma=30$ nm and $\sigma=60$ nm are nearly identical. The packet passess equally
through the left and right arms of the ring. The integrals in Fig. \ref{lmit}(a) and  Fig. \ref{lmit}(b) are periodic functions of the magnetic field with the quantum of the  flux threading the ring of radius $R=100$ nm, i.e. $B=0.132$ T, as above. However, an asymmetry of the electron transfer
across the arms appears for $\sigma=120$ nm [Fig. 9(c)] and becomes very strong for $\sigma=240$ nm [Fig. 9(d)]. Moreover, as the result of this asymmetry, the period of the integrals doubles.
For explanation of this effect the solution of the stationary scattering problem of the next subsection is helpful.

\begin{figure}
\begin{tabular}{l}
\includegraphics[width=0.9\columnwidth]{_ry60-crop.pdf}\put(1,80){(a)}\put(-200,85){$\sigma=30$ nm} \\
\includegraphics[width=0.9\columnwidth]{_ry120-crop.pdf}\put(1,80){(b)}\put(-200,85){$\sigma=60$ nm}\\
\includegraphics[width=0.9\columnwidth]{_ry240-crop.pdf}\put(1,80){(c)}\put(-200,85){$\sigma=120$ nm} \\
\includegraphics[width=0.9\columnwidth]{_rys480-crop.pdf}\put(1,80){(d)}\put(-200,85){$\sigma=240$ nm}
\end{tabular}
\caption {Probability density in the left (red dashed line), and right (black line) arms
of the ring and in the central bar (blue line) integrated over time, for a Gaussian enveloppe of
the wave packet in the initial condition $\exp(-(y-y_0)^2/4\sigma ^2)$, with $\sigma$ given
in the figure. The wave vector applied in the initial condition is $k_y=0.01$/nm.
} \label{lmit} \color{black}
\end{figure}

\subsection{Stationary electron flow}

We solve the standard stationary quantum scattering problem for the Fermi level electron
using the wave function matching technique \cite{wfm} for the atomistic tight binding Hamiltonian. We set  $E_F=6.4$ meV for which the Fermi wave vector is displaced by $k_y=0.01$/nm from the $K'$ Dirac point
The integral over the scattering density in space is plotted in Fig. \ref{ciajm}(a).
The result corresponds very well with Fig. \ref{lmit}(d), only the features are more abrupt
in the stationary case, which is due to the presence of a finite range of $k_y$ in the wave packet dynamics.

Figures \ref{ciajm}(b,c,d) show the scattering density for 0.12, 0.13 and 0.14 T, respectively.
For 0.13 T [Fig.\ref{ciajm}(c)] the parts of the electron density passing through both the arms of the ring meet in phase at the exit to the ring, and the electron wave function does not enter the central bar. However, for slightly different magnetic field [Fig.\ref{ciajm}(b,d)] a phase difference appears, the
part of the wave function is injected to the central bar from above, and the interference
within the ring promotes right or left arm of the ring. In the extreme conditions of Figs. \ref{ciajm}(b,d) the
 electron circulates around a {\it half } of the entire ring, which is the origin of the period doubling
on the magnetic field scale, as the area for the magnetic field flux is halved.

In the time-dependent dynamics for short packets the injection to both the arms of the ring is nearly ideally symmetric. Figure \ref{ciajm}(e) gives the integrals for the stationary flow
for $B=0.075$ T as a function of $k_y$.  The scattering density switches very fast from the left
to the right arm of the ring.  The horizontal bars at the top of the Figure
show the segments of $k_y$ from $0.01 \frac{1}{\mathrm{nm}} \pm \sigma_k=0.01 \frac{1}{\mathrm{nm}} \pm \frac{1}{2\sigma}$, for $\sigma=30,60,120$ and 240 nm (from top to bottom). Only for large $\sigma$ the asymmetry survives averaging over $k_y$ range contained within the packet, hence
the symmetric transfer for shorter packets.

Figure \ref{ciajm}(f) shows the backscattering probability $R$
for the ring with the defect given by Eq. (13). The conductance $G=(1-R)\frac{2e^2}{h}$
has the period of $0.264$ T -- corresponding to the flux quantum
threading half of the ring
in consistence with the period of the scattering density integrals of Fig. \ref{ciajm}(a).

\begin{figure}
\begin{tabular}{l}
\includegraphics[width=0.9\columnwidth]{sf4_niezaburzony_calka.pdf}\put(1,80){(a)}\\
\includegraphics[width=0.3\columnwidth]{map_120mT-crop-1.png} \put(-50,80){(b) 0.12 T}\;
\includegraphics[width=0.3\columnwidth]{map_130mT-crop-1.png} \put(-50,80){(c)  0.13 T}\;
\includegraphics[width=0.3\columnwidth]{map_140mT-crop-1.png} \put(-50,80){(d) 0.14 T }\\
\includegraphics[width=0.8\columnwidth]{bu.pdf}\put(1,80){(e)}\\
\includegraphics[width=0.6\columnwidth]{sf4_zaburz_bezSO.pdf}\put(1,80){(f)}\\
\end{tabular}
\caption { Solutions of the stationary scattering problem for $E_F=6.4$ meV,
which corresponds to the wave vector $k_y=0.01/$nm as calculated with respect to the $K'$ valley.
(a) Scattering  density integrated over  the left (red dashed line), the right (black line) arm
of the ring, and in the central bar (blue line).
(b-d) The scattering density for $B=0.12$, $0.13$, $0.14$ T.
(e) Same as (a) only as a function of $k_y$ for $B=0.075$T. The horizontal bars show
the segments $[k_y-\sigma_k,k_y+\sigma_k]$ with the standard deviation of the Gaussian packet
$\sigma_k=\frac{1}{2\sigma}$ for $\sigma=30$, 60, 120, and 240 nm, from top to bottom.
(f) Backscattering probability obtained for the defect potential given by Eq. (13).
} \label{ciajm}
\end{figure}

\subsection{Quantum rings defined by zero lines in bilayer graphene}
Qualitatively similar results for the transfer across the quantum rings
defined by the zero lines of the symmetry breaking electric field are
found for the bilayer graphene.

For bilayer graphene we use the atomistic tight-binding Hamiltonian spanned by $p_z$ orbitals,
\begin{equation}
   H=\sum_{\langle i,j\rangle }\left(t_{ij} c_i^\dagger c_j+h.c.\right)+\sum_i V({\bf r}_i) c_i^\dagger c_i,
\label{eq:dh}
\end{equation}
where  $V({\bf r}_i)$ is the external potential
at the $i$-th site at position $\mathbf{r}_i$, and in the first term we sum over the nearest neighbors. We use the tight-binding parametrization of
Bernal stacked layers \cite{Partoens}, with
$t_{ij}=-3.12$ eV for the nearest neighbors within the same layer.
 For the interlayer coupling, we take $t_{ij}=-0.377$ eV for the A-B dimers, $t_{ij}=-0.29$ eV for skew interlayer hoppings \cite{Partoens}  between atoms of the same sublattice (A-A or B-B type), and $t_{ij}=0.12$ eV for skew interlayer hopping between atoms of different sublattices.

For simulation of the ring, we assume potential of the form given by Eq. (12) on the upper
sublattice and an opposite potential on the lower sublattice [see Fig. \ref{bila}(a)].
We set $V_g=0.2$ eV, $\lambda=4$ nm as above, but  the radius of the ring is taken equal to $R=50$ nm.
The magnetic field period corresponding to the flux quantum threading the circle
of this radius is $0.53$ T.

The dispersion relation for armchair nanoribbon is displayed in Fig. \ref{bila}(b).
For calculations we take the the Fermi energy is $E_F=0.1$ eV.
For bilayer graphene we have two energy bands instead of the single one moving up the ribbon towards the ring. The integrals of the scattering density are plotted in Fig. \ref{bila}(c)
and display the periodicity with doubled period of $2\times 0.53$ T, as in silicene.

In order to produce the backscattering we removed an atom of the upper graphene layer
from the center of the right arm. We selected an atom that does not form a vertical dimer with the lower layer.
The backscattering probability -- the sum of probabilities for each of the incident subbands --
is given in Fig. \ref{bila}(d) and display the periodicity corresponding to the flux
through half the ring, as found above for silicene.

\begin{figure}
\begin{tabular}{l}
\includegraphics[width=0.4\columnwidth]{scheme.pdf}\put(1,30){(a)}
\includegraphics[width=0.6\columnwidth]{rel_dysp_bilayer.pdf} \put(1,30){(b)}\;\\
\includegraphics[width=0.8\columnwidth]{dens_frac_pristine.pdf} \put(1,80){(c)}\;\\
\includegraphics[width=0.8\columnwidth]{R_vacancy.pdf}\put(1,80){(d)}\\
\end{tabular}
\caption {(a) Schematics of the potential on the upper and lower layers
of bilayer Bernal stacked graphene. Blue and red colors correspond to opposite sign of
the potential.
We consider an armchair nanoribbon of width $W=49.2$ nm, with 1604 atoms in the elementary cell (on both layers). The zero line forms the ring of radius $R=50$ nm. The computational box in the ring area covers an area of radius  74.5 nm.
(b) The dispersion relation for the amrchair nanoribbon feeding the current to the ring. The dots indicate the Fermi wave vectors
at the bands with positive velocity for $E_F=0.1$ eV. $a_{cc}$ is the nearest neighbor distance for graphene.
(c) Integral of the scattering density in the left and right arms of the ring and in the bar for $E_F=0.1$ eV.
(d) Backscattering probability for a vacancy: a carbon atom removed from the center of the right
arm from a position that does not form a vertical dimer with the bottom layer.
} \label{bila}
\end{figure}

\color{black}
\section{Summary}
We studied the dynamics of electron wave packets in buckled silicene in inhomogeneous vertical electric field
that breaks the symmetry between the sublattices using an atomistic tight-binding approach.
We have demonstrated that the line of the electric field flip in silicene supports a smooth
untrembling motion of unspreading wave packets that are topologically protected from backscattering.
We proposed a form of a quantum ring that uses branching of the zero line to split the wave packets
and to make them interfere again. The ring stores the packet for a finite time that can be controlled
with the external magnetic field. { For short wave packets the time spent by the electron in the left and right arms of the ring is a periodic function of the flux with the period
of the flux quantum threading the ring. We found that for long packets, close to the plane waves, the electron transport across the rings becomes asymmetrical with an imbalance of the electron transfer across the left and right halves of the ring. In consequence the magnetic period is doubled.
 This  result is reproduced by stationary scattering calculations.
We demonstrated that the same effect is found for rings defined in  bilayer graphene.
The point defects produce backscattering probability that has a period of the flux quantum threading
ring for short packets. For long packets and in the stationary transport the period of the backscattering probability is doubled.
}

\section*{Acknowledgments}
The results of Section III.F were obtained by B.R. (stationary transport silicene) that is supported by
by the Polish Government budget for science 2017--2021 within the Diamentowy Grant project (Grant No. 0045/DIA/2017/46).
The results of Section III.G were calculated by A.M.-K. (stationary transport bilayer graphene) that is supported by National Science Centre (NCN) according to decision DEC-2015/17/B/ST3/01161.
The rest of the results were provided by B.S. that is supported by the National Science Centre (NCN) according to decision DEC-2016/23/B/ST3/00821.
B.R. acknowledges the support of EU PhD Project POWER.03.02.00-
00-I004/16. The calculations were performed
on PL-Grid Infrastructure at ACK Cyfronet AGH.