\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2007(2007), No. 133, pp. 1--20.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2007 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2007/133\hfil General solution] {General solution of overdamped Josephson junction equation in the case of phase-lock} \author[S. I. Tertychniy\hfil EJDE-2007/133\hfilneg] {Sergey I. Tertychniy} \address{Sergey I. Tertychniy \newline VNIIFTRI, Mendeleevo, Moscow Region, 141570, Russia} \email{bpt97@mail.ru} \thanks{Submitted May 25, 2007. Published October 12, 2007.} \subjclass[2000]{33E30, 34A05, 34A25, 34B30, 34B60, 34M05, \hfill\break\indent 34M35, 70K40} \keywords{Overdamped Josephson junction; nonlinear first order ODE; \hfill\break\indent linear second order ODE; double confluent Heun equation; phase lock; general solution} \begin{abstract} The first order nonlinear ordinary differential equation $\dot \varphi(t) + \sin\varphi(t)=B+A\cos\omega t$, which is commonly used as a simple model of an overdamped Josephson junction in superconductors is investigated. Its general solution is obtained in the case known as phase-lock where all but one solution converge to a common `essentially periodic' attractor. The general solution is represented in explicit form in terms of the Floquet solution of a double confluent Heun equation. In turn, the latter solution is represented through the Laurent series which defines an analytic function on the Riemann sphere with punctured poles. The series coefficients are given in terms of infinite products of $2\times2$ matrices with a single zero element. The closed form of the phase-lock condition is obtained and represented as the condition for existence of a real root of a transcendental function. The efficient phase-lock criterion is conjectured, and its plausibility is confirmed in numerical tests. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{conjecture}[theorem]{Conjecture} \section{Introduction} The nonlinear first order ordinary differential equation \begin{equation} \dot \varphi(t) + \sin\varphi(t)=q(t), \label{eq1} \end{equation} is commonly used in applied physics as the simple mathematical model describing the basic electric properties of Josephson junction (JJ) in superconductors \cite{R1,R2}. Here the right-hand side function $q(t)$, assumed to be known, specifies the external impact representing the appropriately normalized \textit{bias current} (or simply \textit{bias}) supplied by an external current source. The unknown real valued function $\varphi(t)$ called {\it the phase\/} describes the macroscopic quantum state of JJ. In particular, it is connected with the instantaneous voltage $V$ applied across JJ in accord with the equation $V=(\hbar/2e)\mathrm{d}\varphi/\mathrm{d}\tau$, where $\hbar$ is the Plank constant, $e$ is the electron charge, $\tau$ is the (dimensional) current time. The dimensionless variable $t$ entering \eqref{eq1} is defined as $t=\omega_c\tau$, where $\omega_c$ is a constant parameter depending on the junction properties and named JJ characteristic frequency. See \cite{R5} for more details of JJ physics. Equation \eqref{eq1} arises as the limiting case of the second order ODE utilized in more general Resistively Shunted Junction (RSJ) model \cite{R3,R4}. The reduction is legitimate if the role of the junction capacitance proves negligible. In practice, if JJ can be described by \eqref{eq1} it is named {\it overdamped}. Summarizing the aforementioned relationships, we shall name \eqref{eq1}, for brevity, {\em overdamped Josephson junction equation} (\eqref{eq1}--\eqref{eq2}). Under concordant conditions, the theoretical modelling applying \eqref{eq1}--\eqref{eq2} is in excellent agreement with experiments. It is also worth noting that nowadays electronic devices based on the Josephson effect play the important role in various branches of measurement technology. In particular, JJ arrays serve the heart element of the modern DC voltage standards \cite{R6}. The development of JJ-based synthesizers of AC voltage waveforms is currently in progress \cite{R7,R8,R9}. These and other successful applications stimulate the growing interest to the theoretical study and the modelling of JJ properties including investigation of capabilities of RSJ model and its limiting cases and the predictions they lead to. To be more specific, the case most important from viewpoint of applications and simultaneously distinguished by the wealthiness of the underlaid mathematics is definitely the one of the bias function representing harmonic oscillations. Without loss of generality, it can be represented in the form \begin{equation}\label{eq2} q(t)=B+A\cos\omega t, \end{equation} where $A,B,\omega>0$ are some real constant parameters. Hereinafter, the abbreviation \eqref{eq1}--\eqref{eq2} introduced above will refer to the couple (\ref{eq1}), (\ref{eq2}). In spite of apparent simplicity of \eqref{eq1}--\eqref{eq2}, few facts of its specific analytic theory had been available until recently. Some preliminary results concerning the problem of derivation of analytic solutions of \eqref{eq1}--\eqref{eq2} in the general setting had been obtained in \cite{R11a}. The approach put forward therein is elaborated in the present work. The focus is made on the case of manifestation by \eqref{eq1} of the {\em phase-lock\/} property which is one of its most important features from viewpoint of applications. The phase-lock property is formalized as follows: in the case of phase-lock any solution $\varphi=\varphi(t)$ of \eqref{eq1} either yields a periodic exponent $e^{\mathrm{i}\varphi}$, exactly two such distinguished solutions existing, or, as the time parameter grows, the exponent $e^{\mathrm{i}\varphi}$ exponentially converges to the similar exponent for the one (common for all $\varphi$'s) of the periodic functions just noted (another one plays the role of repeller). (Here and in what follows we shall not distinguish phase functions which differ by a constant equal to $2\pi$ times an integer.) The corresponding period coincides with one of $q(t)$; i.e., in the case (\ref{eq2}), $2\pi\omega^{-1}$. This behavior is stable with respect to weak parameter perturbations; i.e., the subset of parameter values leading to phase-lock is open. It is worth noting for completeness that in the opposite (no phase-lock) case no {\it stable\/} periodicity in the behavior of $e^{\mathrm{i}\varphi}$ is observed. There is also a third, intermediate, type of the phase behavior, where the asymptotic attracting and repelling solutions are, in a sense, merged \cite{R11}. It is realized on the lower-dimensional subset of the space of parameter values. In the present work, the complete analytic solution of \eqref{eq1}--\eqref{eq2} is obtained under assumption of the parameter choice ensuring phase-lock. The closed form of the phase-lock criterion in the form of the constraint imposed on the problem parameters is conjectured. \section{Overdamped Josephson junction equation against reduced double confluent Heun equation}\label{s1} The analytic theory of \eqref{eq1}--\eqref{eq2} can be based on its reduction to the following system of two {\em linear\/} first order ODEs, \cite{R11a}, \begin{equation}\label{eq3} \begin{gathered} 4\mathrm{i}\omega z^2 \;x'(z) = 2z x(z) +\left[2B z+A\left(z^2+1\right)\right] y(z),\\ -4\mathrm{i}\omega z^2\; y'(z) =\left[2B z+A\left(z^2+1\right)\right] x(z)+ 2z y(z), \end{gathered} \end{equation} where $z$ is the free complex variable. Indeed, on the universal covering $\Omega_1\simeq\mathbb{R}\ni t $ of the unit circle $S^1$ in $\mathbb{C}$; i.e., for $z=\exp(\mathrm{i}\omega t)$, any non-trivial solution of \eqref{eq3} determines a solution of \eqref{eq1}--\eqref{eq2} in accordance with the equation \begin{equation}\label{eq4} \exp(\mathrm{i}\varphi)=\frac{\Re x-\mathrm{i} \Re y}{\Re x+\mathrm{i} \Re y} \end{equation} ($\Re$ means the real part) supplemented with the continuity requirement. [The fulfillment of \eqref{eq1}--\eqref{eq2} follows from a straightforward computation taking into account the equality $\overline z=z^{-1}$ holding true on $\Omega_1$ which is utilized for the demonstration that, for real $A,B,\omega$, the functions $\Re x,\Re y$ also verify \eqref{eq3} on $\Omega_1$.] Conversely, any real valued solution $\varphi(t)$ of \eqref{eq1}--\eqref{eq2} induces through \eqref{eq4} some `initial data' $x(0),y(0)$ (with arbitrary norm $(x^2(0)+y^2(0))^{1/2}>0$) for \eqref{eq3}. Having solved the latter on $\Omega_1$ for $x,y$, one obtains, applying \eqref{eq4}, another phase function obeying \eqref{eq1}--\eqref{eq2}. It however must coincide, due to the identical initial value assumed at $t=0$, with the original $\varphi(t)$. Finally, since these $x,y$ defined on $\Omega_1$ simultaneously obey the linear ODEs (\ref{eq3}) with meromorphic coefficients which have the only singular points $z=0$ and $z^{-1}=0$, they can be extended, integrating (\ref{eq3}), say, along radial directions, to analytic functions defined on the whole $\Omega=\Omega_1\times \mathbb{R}_+$, the universal covering of the Riemann sphere with the punctured poles $z=0$ and $z^{-1}=0$. It is worth noting that on $\Omega_1$ the real valued functions $\tilde{x} =\tilde{x}(t)=\Re x(e^{\mathrm{i}\omega t}), \tilde{y}=\tilde{y}(t)=\Re y(e^{\mathrm{i}\omega t})$ satisfy the equations \begin{equation}\label{eq5} 2\mathrm{d} \tilde{x}/\mathrm{d} t = \tilde{x}+q \tilde{y},\quad -2\mathrm{d} \tilde{y}/\mathrm{d} t = q \tilde{x}+ \tilde{y}, \end{equation} leading, together with \eqref{eq4}, to the equation $({\mathrm{d}/\mathrm{d} t})\log(\tilde{x}^2+\tilde{y}^2)=\cos\varphi$ which implies \begin{equation}\label{eq5a} \mathrm{const}_1\,e^{-t}\le |\tilde{x}+\mathrm{i}\tilde{y}|^2\le \mathrm{const}_2\,e^t \quad \mbox{for }t>0. \end{equation} Note that if $x,y\not\equiv 0$ then $\mathrm{const}_1,\mathrm{const}_2$ may be assumed to be strictly positive. We shall refer to these bounds later on. The key observation enabling one to radically simplify the problem of description of the space of solutions of \eqref{eq3} is as follows \cite{R11a}. Let us introduce the analytic function $v(z)$ which satisfies the equation \begin{equation} \label{eq6} \Big[ z^2\frac{\mathrm{d}^2}{ \mathrm{d}^2 z} +\left( \mu (z^2+1) -n z\right) \frac{\mathrm{d}}{\mathrm{d} z} +(2\omega)^{-2} \Big] v =0, \end{equation} where \begin{equation} \label{eq7} n= -\big(\frac{B}{\omega}+1\big),\quad \mu=\frac{A}{2\omega}, \end{equation} are the constants replacing original $A,B$ which will be used below whenever it proves convenient. Then a straightforward calculation shows that the functions $x,y$ determined by the equations \begin{gather} \label{eq8} v = \hphantom{(2\omega z)^{-1}} \mathrm{i} z^{{n+1\over2}} \exp \Big( \frac{1}{2}\mu \left(-z+z^{-1}\right) \Big) (x-\mathrm{i} y), \\ \label{eq9} v' = (2\omega z)^{-1} z^{{n+1\over2}} \exp\Big(\frac{1}{2}\mu \left(-z+z^{-1}\right)\Big)(x+\mathrm{i} y) \end{gather} verify \eqref{eq3}. Conversely, defining the function $v(z)$ through the solution $x,y$ of \eqref{eq3} in accordance with \eqref{eq8}, a straightforward computation proves satisfaction of \eqref{eq9} and, then, \eqref{eq6} follows. Thus, \eqref{eq3} are equivalent to (\ref{eq6}) and \eqref{eq8}, \eqref{eq9} represent the corresponding one-to-one transformation. Equation \eqref{eq6} coincides, after appropriate identification of the constant parameters, with \cite[Eq.\ 1.4.40]{R14}. It represents therefore a particular instance of the double confluent Heun equation which can be shown to be non-degenerated in all cases where \eqref{eq1}--\eqref{eq2} is non-trivial. It is also worth reproducing here the {\it canonical form\/} of the ``generic'' double confluent Heun equation as it is given in \cite[Eq. 4.5.1]{R14}. It reads \begin{equation}\label{eq10} z^2\frac{\mathrm{d}^2y}{\mathrm{d} z^2} +(-z^2+c z+t)\frac{\mathrm{d} y}{\mathrm{d} z}+(-a z+\lambda)y=0 \end{equation} where $a,c,t,\lambda$ are some constants. To adjust it to our case, a single term has to be eliminated setting $a=0$. For brevity, we shall name this subclass of double confluent Heun equation's {\it reduced.} Besides, some obvious rescaling of the free variable $z$ is to be carried out. After these, the three residuary constant parameters exactly correspond to our constant parameters $n,\mu,\omega$ (we shall not need and so omit the reproducing of the concrete form of this transformation). The general analytic theory of double confluent Heun equation is given in the chapter 8 of the treatise \cite{R13}. Its solutions are there represented, up to nonzero factors given in explicit form, through the Laurent series whose coefficients are assumed to be computable through the `endless' chain of 3-term linear homogeneous equations (`recurrence relations'). In the present work, we derive the solution of reduced double confluent Heun equation in a cognate but more explicit form. As a technical limitation, we also stipulate in the present work for the additional condition to be imposed on the free constant problem parameters claiming of them the ensuring of the phase-lock property. On the base of practice of numerical computations, it can be conjectured that such parameter values fill up a non-empty open subset ({\em phase-lock area}) in the whole parameter space (see also the Conjecture \ref{conjA} below). The case where the parameters belong to its complement is left beyond the scope of the present work. \section{Formal solution of reduced double confluent Heun equation} Let us introduce yet another unknown function $E(z)$ replacing $v(z)$ by means of the transformation \begin{equation} \label{eq11} v(z)= z^{{n+\epsilon\over2}-\mathrm{i}{\kappa}} e^{-{\mu}z} E(z), \end{equation} where the discrete {\em parity parameter $\epsilon$\/} may assume one of the two values, either $\epsilon=0$ or $\epsilon=1$ (i.e.\ obeys the equation $\epsilon^2=\epsilon$), and $\kappa$ is some real positive constant which will be determined latter. For $v$ obeying (\ref{eq6}), $E(z)$ verifies the equation \begin{equation} \label{eq12} \begin{aligned} 0&= z^3 E'' +z\left[ \left(\epsilon- 2\mathrm{i}{\kappa}\right) z - \mu \left(z^2-1\right) \right] E' \\ &\quad +\Big[ \mu \big( {n-\epsilon\over2}+\mathrm{i} {\kappa}\big) z^2 +\Big((1-\epsilon)\big({1\over 4} + \mathrm{i} {\kappa}\big) -\kappa^2 -\big({n+1\over2}\big)^2 +\lambda\Big) z \\ &\quad +\mu\big({n+\epsilon\over 2}-\mathrm{i} {\kappa}\big) \Big]E, \end{aligned} \end{equation} where $\lambda=(2\omega)^{-2}-\mu^2$. Conversely, (\ref{eq11}),(\ref{eq12}) imply the fulfillment of \eqref{eq6}. At first glance, \eqref{eq12} seems `much worse' than the original double confluent Heun equation representation. Nevertheless, it is this equation which we shall attempt to solve searching for its solution in the form of Laurent series \begin{equation} \label{eq13} E=\sum^\infty_{k=-\infty} a_k z^{k} \end{equation} `centered' in the points $z=0$ and $z^{-1}=0$ (which are the only singular points for \eqref{eq12}) with unknown $z$-independent coefficients $a_k$. Then, carrying out straightforward substitution, one gets a sequence of 3-term recurrence relations binding triplets of neighboring $a$'s %series coefficients which can be written down %represented either as \begin{equation} \begin{aligned} 0&=- \mu \left(k-1-{n-\epsilon\over2}- \mathrm{i} {\kappa}\right) a_{k-1}\\ &\quad +(Z_{k}+\lambda) a_k + \mu \left(k+1+{n+\epsilon\over2} - \mathrm{i} {\kappa}\right) a_{k+1}, \end{aligned} \label{eq14} \end{equation} where \begin{equation} Z_k= \left(k+{\epsilon-1\over2} -\mathrm{i} {\kappa}\right)^2-\left({n+1\over2}\right)^2, \label{eq15} \end{equation} or as \begin{equation} \begin{aligned} 0&= - \mu \left(k-1-{n+\epsilon\over2}+\mathrm{i} {\kappa}\right)a_{-(k-1)} \\ &\quad +(\tilde {Z_{k}}+\lambda) a_{-k} + \mu \left(k+1+{n-\epsilon\over2} +\mathrm{i} {\kappa}\right) a_{-(k+1)}, \end{aligned} \label{eq16} \end{equation} where \begin{equation} \tilde Z_k=\left(k+{1-\epsilon\over2}+ \mathrm{i} {\kappa}\right)^2-\left({n+1\over2}\right)^2 \label{eq17} \end{equation} and $k=0,\pm1,\pm2\dots$. The sets of \eqref{eq14}) and \eqref{eq16} are exactly equivalent and each of them separately covers the whole set of equations the coefficients $a_k$ have to obey. However, in the approach utilized in the present work, we shall consider them both in conjunction, employing (\ref{eq14}) for coefficients $a_k$ with indices $k\ge -1$ and (\ref{eq16}) for $a_k$ with indices $k\le 1$. Thus, \eqref{eq14} and (\ref{eq16}) will be considered separately but on the common index variation `half-interval' $k\ge 0$ (remaining legitimate, in principle, for arbitrary integer $k$). Obviously, these two equation sets cover the complete set of conditions imposed to the coefficients $a_k$ and are `almost disjoined' coupling % intersecting in their `boundary' $k=\pm0$-members alone. Let us further consider for $k\ge 0$ the following {\it formal infinite products\/} \begin{equation} \label{eq18} R_{k}=\prod_{j=k}^\infty M_j,\quad \tilde R_{k}=\prod_{j=k}^\infty \tilde M_j \end{equation} of the $2\times2$ matrices \begin{gather} \label{eq19} M_j= \begin{pmatrix} 1+\lambda/ Z_{j} & \mu^2/ Z_{j} \\ 1 & 0 \end{pmatrix}, \\ \label{eq20} \tilde M_j= \begin{pmatrix} 1+\lambda/\tilde Z_{j} & \mu^2/\tilde Z_{j} \\ {\tilde Z_{j-1}}/{\tilde Z_{j}} & 0 \end{pmatrix}. \end{gather} It is assumed throughout that the matrices $M_j,\tilde M_j$ with {\it larger\/} indices $j$ are situated in products {\it to the right\/} with respect to ones labelled with {\it lesser\/} index values. Notice that in the case $\kappa=0$ and integer $n$, zero may appear in denominators of $M$-factors, making the above definitions meaningless. This apparent fault admits a simple resolution (see \eqref{eq34} and the discussion following it). For a while, we temporary leave out consideration of such specific parameter choices. It is also worth noting that the above definitions of $R_{k},\,\tilde R_{k}$ may be understood as a concise form of representation of the `descending' recurrence relations \begin{gather}\label{eq21} R_k = M_{k} R_{k+1},\\ \label{eq21a} \tilde R_k = \tilde M_{k} \tilde R_{k+1},\quad k=0,1,\dots \end{gather} among the neighboring $R$'s. These are the only dependencies which will be actually used below in derivations involving $R_{k},\tilde R_k$. The formulas (\ref{eq18}) are `formal' since neither the issue of the convergence of such sequences nor how one should understand the `initial values' $R_\infty, \tilde R_\infty$ necessary for the actual determination of the `finite index value' $R$-matrices by means of \eqref{eq21}--\eqref{eq21a} are here addressed. Now a straightforward calculation applying \eqref{eq21}--\eqref{eq21a} allows one to show that the following formulas \begin{gather} \label{eq22} \tilde a_k = \mu^k {\Gamma\left(1+{n+\epsilon\over2}- \mathrm{i} {\kappa}\right) \over \Gamma\left(k+1+{n+\epsilon\over2}- \mathrm{i} {\kappa}\right) } (0,1)\cdot R_k \cdot \begin{pmatrix}1\\ 0\end{pmatrix} \\ \label{eq23} \tilde a_{-k}= {\mu^k \over {\tilde Z_{k-1}}} {\Gamma\left(1+{n-\epsilon\over2}+ \mathrm{i} {\kappa}\right) \over \Gamma\left(k+1+{n-\epsilon\over2}+ \mathrm{i} {\kappa}\right) } (0,1)\cdot \tilde R_k \cdot \begin{pmatrix}1\\0 \end{pmatrix} \end{gather} (where $\Gamma$ is the Euler gamma-function) yield {\it the formal solutions\/} to \eqref{eq14} and \eqref{eq16}, respectively. \section{Validation of the formal solution} In this section we show that the formal solution of \eqref{eq12} presented in the form of expansion (\ref{eq13}) with coefficients given by \eqref{eq22},\eqref{eq23} represents a well defined analytic function of $z$. This means, first of all, that the infinite matrix products $R_k, \tilde R_k$ it involves converge. Moreover, the convergence takes place for any constant parameter values. The key auxiliary result which may be utilized for the proof of this assertion is as follows. \begin{theorem} \label{thm1} Let us consider the sequences of complex numbers $\alpha_j,\beta_j,\gamma_j,\delta_j$ satisfying the following `ascending' matrix recurrence relation \begin{equation} \label{eq24} \begin{pmatrix} \alpha_j &\beta_j\\ \gamma_j & \delta_j \end{pmatrix} = \begin{pmatrix} \alpha_{j-1} &\beta_{j-1}\\ \gamma_{j-1} & \delta_{j-1} \end{pmatrix} \begin{pmatrix} 1+\lambda/ Z_{j} & \mu^2/ Z_{j} \\ \sigma_j & 0 \end{pmatrix}, \end{equation} where \begin{equation} \label{eq25} \mbox{either $\sigma_j= 1$ or $\sigma_j= Z_{j-1}/ Z_{j}$.} \end{equation} Then all they converge as $j\to \infty$. Moreover, $\lim\beta_j=0=\lim\delta_j$ whereas for $\alpha,\gamma$-sequences there exist positive quantities $N_\alpha$, $N_\gamma$ and positive integer $j_0$, all depending at most on $n,\mu,\lambda,\kappa,\epsilon$, such that for all $j> j_0$ \begin{equation} \label{eq26} \begin{gathered} |\alpha_{j}-\lim_{j'\to \infty} \alpha_{j'}|< {N_\alpha \max_{j'>j_0}(|\alpha_{j'}|) j^{-1}},\\ |\gamma_{j}-lim_{j'\to \infty} \gamma_{j'}| < {N_\gamma \max_{j'>j_0}(|\gamma_{j'}|) j^{-1}}, \end{gathered} \end{equation} where the maxima are finite. \end{theorem} The outline of the proof can be found in Appendix. \smallskip \noindent{\it Remark}: Formally, we need not include in the theorem stipulation the requirement that either $\kappa\not =0$ or $n$ is non-integer (which would a priori evade possibility arising of contribution with vanishing $Z_{*}$ in the denominator) because with fixed constant parameters and sufficiently large $j_0$ no zero $Z_{*}$ may appear. Let us return to \eqref{eq18} and consider the four double-indexed sets of complex numbers $\{\alpha,\beta,\gamma,\delta\}^{(j_0)}_j$ defined as follows: \begin{equation} \label{eq26a} \begin{pmatrix} \alpha_{j}^{(j_0)} &\beta_{j}^{(j_0)}\\ \gamma_{j}^{(j_0)} & \delta_{j}^{(j_0)} \end{pmatrix} =R_{j_0}^{(j)} =\prod_{k=j_0}^{j} M_k. \end{equation} It is straightforward to verify that the sequences obtained by the picking out the elements with common value of the upper index $j_0$ obey the recurrence relations (\ref{eq24}) for the first choice in (\ref{eq25}). Hence it follows from the theorem that all they converge. We denote the corresponding limits as $\alpha^{(j_0)}$ etc. We have therefore the consistent definition for the infinite matrix products (\ref{eq18}), \begin{equation} \label{eq26b} R_{j_0}= \begin{pmatrix} \alpha^{(j_0)} &\beta^{(j_0)}\\ \gamma^{(j_0)} & \delta^{(j_0)} \end{pmatrix}. \end{equation} Let us further note that, increasing $j_0$, the `starting' sequence elements ($\alpha_{j_0}^{(j_0)},\alpha_{j_0+1}^{(j_0)}$ for $\alpha$ etc) tend to the corresponding elements of the idempotent matrix \begin{equation} \label{eq26c} M_{\infty}= \begin{pmatrix}1&0\\1&0 \end{pmatrix},\quad M_{\infty}^2=M_{\infty}, \end{equation} the discrepancy decreasing like $O(j_{0}^{-2})$. On the other hand, in accordance with (\ref{eq26}) the elements of $R_{j_0}^{(j)}$ differ from their $j$-limits constituting $R_{j_0}$ by $O(j_{0}^{-1})$-order quantities. This means that $R_{j_0}$ tends to $M_{\infty}$ as $j_0$ goes to infinity with the difference going to zero as $O(j_{0}^{-1})$. In other words, we have the following result. \begin{corollary} \label{cor0} 1. $R_{j_0}-M_{\infty}=O(j_{0}^{-1})$. Furthermore, in view of the convergence, 2. The modules of the elements of all the matrices $R_{j_0}$ are bounded from above in total --- provided `no-zeroes-in~denominators' condition for the parameter choice is respected, of course. \end{corollary} In according with the above relations, decomposing $R_{j}$ into the product of the two factors, $R_{j}=R_j^{j_0}\cdot R_{j_0}$, and approximating $R_{j_0} $ by $M_\infty$, one obtains a simple but important algorithm of computation of the products (\ref{eq18}) \begin{equation} \label{eq18a} R_{j}\approx\prod_{k=j}^{j_0} M_k\cdot M_\infty.\; \end{equation} The approximation is the better, the larger $j_0>j$ is selected. In the limit, one gets \begin{equation} \label{eq28a} R_j= \begin{pmatrix} \alpha^{(j)} &0\\ \gamma^{(j)} &0 \end{pmatrix}. \end{equation} This interpretation resolves (in a quite obvious way, though) the aforementioned uncertainty in %problem of specification of the `initial value' for the sequence $R_j$ treated through the `descending' recurrence relation `$[R_\infty]\dots\to R_j\to R_{j-1}\to \dots$' implied by (\ref{eq21}). After introducing a consistent representation of the matrices (\ref{eq18}), it is straightforward to do the same for the matrices $\tilde R_j$ (\ref{eq19}). The above speculation applies to them with minor modifications as well. The only distinction is the making use of the second choice in (\ref{eq25}) and the operating with complex conjugated quantities (equivalent in our case to the replacing $\kappa$ by $-\kappa$) throughout. We shall mark the elements of $\tilde R_j$ obtained in this way with tildes over the corresponding $\alpha$'s and $\gamma$'s. Now, with well defined $R_j$ and $\tilde R_j$ (\ref{eq18}) in hand, one is able to calculate the coefficients $a_k,a_{-k}$ for $k=0,1,2\dots$ in accordance with \eqref{eq22} \eqref{eq23}. The triple matrix products reduce to separate elements of $R$'s (or $\tilde R$'s) denoted above as $\gamma^{(k)}$ (for (\ref{eq22})), and $\tilde\gamma^{(k)}$(for \eqref{eq23}, respectively) which are the functions of the parameters $n,\mu,\lambda,\kappa,\epsilon$. Therefore, the sequences \begin{gather} \label{eq31} \tilde a_k = \mu^k \gamma^{(k)} { \Gamma\left(1+{n+\epsilon\over2}- \mathrm{i} {\kappa}\right) \over \Gamma\left(k+1+{n+\epsilon\over2}- \mathrm{i} {\kappa}\right) } \\ \label{eq32} \tilde a_{-k} ={\mu^k \tilde \gamma^{(k)} \over {\tilde Z_{k-1}}} { \Gamma\left(1+{n-\epsilon\over2}+ \mathrm{i} {\kappa}\right) \over \Gamma\left(k+1+{n-\epsilon\over2}+ \mathrm{i} {\kappa}\right) },\quad k=0,1,2,\dots \end{gather} are well defined and solve \eqref{eq14}, \eqref{eq16}, respectively, for $k=1,2,\dots$ The important feature of the expressions (\ref{eq31}), (\ref{eq32}) which is used below is their asymptotic behavior for large values of the index $k$ which is easy to derive in explicit form. Specifically, in accordance with inequalities (\ref{eq26}), the set of modules of $\gamma$- and $\tilde\gamma$-factors involved in \eqref{eq31}--\eqref{eq32} is bounded in total from above and each of these sequences converge to a finite limit. These imply in particular the validity of the estimates \begin{equation} \label{eq31a} |\tilde a_k| \propto {\zeta_1^k \over k!},\quad |\tilde a_{-k}| \propto {\zeta_2^k \over k^2 k!}, \end{equation} asymptotically, in the leading order, for some $k$-independent $\zeta_1,\zeta_2$. \section{Matching condition and phase-lock criterion} By construction, the $\tilde a$-coefficients defined by \eqref{eq22} and \eqref{eq23} obey the linear homogeneous equations (\ref{eq14}) and (\ref{eq16}), respectively, which are `the same', essentially, differing only in the associated intervals of the variation of the index, consisting of the positive integers for the former and negative ones for the latter. However, these two sequences cannot be joined, automatically, since they are `differently normalized'. This means, in particular, that their `edge elements' indexed with zeroes, generally speaking, differ. We will denote them as $\tilde a_0$ and $\tilde a_{-0}$, respectively, distinguishing here, in notations, the index `$0$' from the index `$-0$'. Now, referring to \eqref{eq31},(\ref{eq32}), one notes that in view of the factors of $\Gamma$-functions present in denominators and leading to asymptotic behaviors (\ref{eq31a}), the following power series in $z$ \begin{equation} \label{eq33} E_+(z)=1+\sum_{k=1}^\infty {\tilde a_k\over\tilde a_0}z^k,\quad E_-(z)=1+\sum_{k=1}^\infty {\tilde a_{-k}\over\tilde a_{-0}}z^{-k}, \end{equation} admit absolutely converging majorants. (We assumed above $\tilde a_{\pm0}\not=0$. Otherwise, i.e.\ if $\tilde a_{\pm0}=0$, $\tilde a_{\pm1}$ may not vanish and the series with the coefficients $\tilde a_{\pm k}/\tilde a_{\pm1}$ can be utilized instead.) Indeed, the Maclaurin series for the exponent function can play this role. Therefore, the series $E_+(z)$ and $E_-(z^{-1})$ define some {\it entire functions} of $z$. As a consequence, the expression \begin{equation} \label{eq34} E(z)= {{4\over\pi^2}\sin\left({\pi\over2}(n+\epsilon- 2\mathrm{i} {\kappa})\right) \sin\left({\pi\over2}(n+\epsilon+ 2\mathrm{i} {\kappa})\right) \over (n+\epsilon- 2\mathrm{i} {\kappa}) (n+2-\epsilon+2\mathrm{i} {\kappa}) } (E_+(z)+E_-(z)-1) \end{equation} represents a single-valued function of $z$ analytic everywhere on the Riemann sphere except the poles $z=0,z^{-1}=0$. They are the essential singular points for $E$. It has to be noted that the additional $z$-independent fractional factor in (\ref{eq34}) may be regarded as a specific common `normalization' of the (\ref{eq33})-type series which may be, in principle, arbitrary. However, its given form is, essentially, unique being fixed (up to an insignificant nonzero numerical factor) in view of the following reasons. The two sine-factors in the numerator regarded as holomorphic functions of $n+\epsilon\pm 2\mathrm{i} \kappa$ are introduced for the cancelling out zeroes in denominators arising due to the poles in the factors $Z_{j}^{-1}$ and $\tilde Z_{j}^{-1}$ involved in the products (\ref{eq21}) and regarded as the functions of the same parameters. The set of these (vicious, essentially) singularities constitute a homogeneous grid which is just covered by the grid of roots of the sine-factors in (\ref{eq34}) --- with the two exceptions. These two `superfluous' sine-factor roots are, in turn, `neutralized' by the two linear factors in the denominator in (\ref{eq34}) which are therefore also uniquely determined. As the result, in vicinity of any zero in denominators in coefficients of the power series defining $E(z)$ considered as the function of $n+\epsilon\pm 2\mathrm{i} \kappa$ (a root of some $Z_{*}$ or $\tilde Z_{*}$), the resulting expression takes the form of the ratio $\sin x/x\;(x\simeq 0)$ and is not now associated with any irregular behavior. Thus, as a matter of fact, the fractional factor involved in (\ref{eq34}) is distinguished (up to a numerical factor) by the claims (i) to cancel out the poles in the original expressions of the $\tilde a$-coefficients (\ref{eq31}), (\ref{eq32}) considered as the functions of $n+\epsilon\pm 2\mathrm{i} \kappa$ and (ii) to introduce neither more zeroes nor more poles as a result of such a `renormalization'. Now, when plugging the function (\ref{eq34}) in \eqref{eq12} in order to verify its fulfillment, we may provisionally drop out, sparing the space, $z$-independent fractional factor (restoring it afterwards). It is important to emphasize that the expressions (\ref{eq31}), (\ref{eq32}) verify, by construction, all the 3-term recurrence relations (\ref{eq14}),(\ref{eq16}) which bind the $a$-coefficients with indices {\em of a common sign}, either non-negative or non-positive. The only equation which does not fall into these categories, and, accordingly, has not been automatically fulfilled, is the `central' one binding the coefficients $a_{-1},a_0=1=a_{-0},a_1$; i.e., the equation \begin{equation} \label{eq35} \mu \left(1+{n-\epsilon\over2}+\mathrm{i} {\kappa}\right) a_{-1} +(Z_{0}+\lambda) a_0 + \mu \left(1+{n+\epsilon\over2}-\mathrm{i} {\kappa}\right) a_{1}=0. \end{equation} With normalization adopted in (\ref{eq33}), one has $a_0=1, a_1=\tilde a_1/\tilde a_0, a_{-1}=\tilde a_{-1}/\tilde a_{-0}$. Further, in accordance with (\ref{eq31}),(\ref{eq32}): \begin{gather*} % \label{eq36} \tilde a_0 = \gamma^{(0)},\quad \tilde a_1= {\mu\over 1+{n+\epsilon\over2}-\mathrm{i} {\kappa}}\gamma^{(1)} \\ \tilde a_{-0} = {\tilde\gamma^{(0)}\over\tilde Z_{-1}},\quad \tilde a_{-1}= {\mu\over 1+{n-\epsilon\over2}+\mathrm{i} {\kappa}}{\tilde\gamma^{(1)}\over\tilde Z_0} \end{gather*} Besides, one has $\gamma^{(0)}=\alpha^{(1)} $, $\tilde\gamma^{(0)}=\tilde\alpha^{(1)}\tilde Z_{-1}/\tilde Z_{0} $. Combining these dependencies, the following representation of (\ref{eq35}) arises \begin{equation}\label{eq37} 0 = \mu^2{\gamma^{(1)} \over\alpha^{(1)} } +(Z_{0}+\lambda) + \mu^2 {\tilde\gamma^{(1)}\over \tilde\alpha^{(1)}}. \end{equation} Accordingly, it is convenient to introduce the following function of the parameters $\kappa$ and $n,\lambda,\mu,\epsilon$ \begin{equation} \label{eq38} \begin{aligned} \Xi&= {{4\over\pi^2} \sin\left({\pi\over2}(n+\epsilon- 2\mathrm{i} {\kappa})\right) \sin\left({\pi\over2}(n+\epsilon+ 2\mathrm{i} {\kappa})\right) \over (n+\epsilon- 2\mathrm{i} {\kappa}) (n+2-\epsilon+2\mathrm{i} {\kappa}) }\\ &\quad \times \left( \mu^2{\gamma^{(1)}\tilde\alpha^{(1)}} +(Z_{0}+\lambda)\tilde\alpha^{(1)}\alpha^{(1)} + \mu^2 {\tilde\gamma^{(1)}\alpha^{(1)}} \right) \end{aligned} \end{equation} where $Z_0=\left({(\epsilon-1)/2} -\mathrm{i} {\kappa}\right)^2-\left({(n+1)/2}\right)^2$ (see (\ref{eq15})) and $\alpha$'s, $\gamma$'s are defined as the elements of the convergent matrix products as follows \begin{equation} \label{eq39} \begin{pmatrix} \alpha^{(1)} &0\\ \gamma^{(1)} &0 \end{pmatrix} = \prod_{j=1}^\infty M_j, \quad \begin{pmatrix} \tilde \alpha^{(1)} &0\\ \tilde\gamma^{(1)} &0 \end{pmatrix} = \prod_{j=1}^\infty \tilde M_j \end{equation} (see \eqref{eq18}-(\ref{eq20})). The fractional multiplier in the first line of \eqref{eq38} coincides with the one entering \eqref{eq34} and plays the identical role: It eliminates the vicious singularities arising for specific values of the parameters $n,\kappa$. We shall name $\Xi=\Xi(\kappa,n,\mu,\lambda,\epsilon)\equiv\Xi(\kappa;\dots)$ the {\em discriminant function\/} for brevity. The following statement holds true. \begin{proposition} \label{prop1} Restricting $\kappa$ to real values, the vanishing of the discriminant function is necessary and sufficient for the single valued analytic function \eqref{eq34} to verify \eqref{eq12} everywhere on the Riemann sphere except its poles $z=0,z^{-1}=0$. \end{proposition} Indeed, the vanishing of $\Xi$ implies the fulfillment of \eqref{eq35} (where $a$'s are expressed through $\tilde a$'s), the last equation among those binding coefficients of the expansion (\ref{eq13}) which has not been fulfilled as the result of the very coefficients definition. Now all the 3-term recurrence relations for $a$-coefficients, which \eqref{eq12} is equivalent to, are satisfied and the analytic function (\ref{eq34}) verifies \eqref{eq12} everywhere except of its own singular points $z=0$ and $z^{-1}=0$. The equation \begin{equation} \label{eq40} \Xi(\kappa;\dots) =0 \end{equation} referred to in the above proposition can be named {\it the matching condition\/} since it enforces the sequences of the coefficients $a_k$, $a_{-k}$, separately obeying the corresponding `halves' of the equation chain (\ref{eq14}) (equivalently, (\ref{eq16})) to be `matched' in their `edge' elements $a_{\pm0}=1,a_{\pm1}$. It is worth emphasizing that, up to this point, the parameter $\kappa$ (absent in \eqref{eq6}) has not been restricted in any way (apart of the claim to be real). Now and in what follows we regard the condition (\ref{eq40}) as $\kappa$ definition eliminating this odd `degree of freedom'. Hereinafter, it is a well defined function of the other parameters. It seems interesting enough that the addition of unspecified constant $\kappa$ to the transformation (\ref{eq11}) and its subsequent `fine tuning' by means of the claim of fulfillment of \eqref{eq40}) is necessary for the representation of solution of double confluent Heun equation (\ref{eq6}) in terms of convergent Laurent series. More precisely, it is clear that \eqref{eq15} can be solved {\em for any\/} (including trivial zero) choice of $\kappa$, choosing loosely $a_{j_0},a_{j_0+1}$ for arbitrary $j_0$ and then calculating, term by term, all the other coefficients $a_j$, advancing in parallel in both directions of $j$-index variation `off $j_0$ towards $\pm\infty$'. Then (\ref{eq13}) immediately yields a ($\kappa$-dependent!) formal solution to \eqref{eq12} and hence, through transformation (\ref{eq11}), to \eqref{eq6}. However, {\em it can be only formal and will necessarily diverge\/} for any $z$ unless the matching condition (\ref{eq40}) is fulfilled -- just in view of the uniqueness of solution with the analytic properties presupposed. On the other hand, considering separately the `halves' of the set of \eqref{eq15} and resolving them `in index variation directions' opposite to the ones assumed above (in a sense, `off $\pm\infty$ towards $\pm0$'), we obtain the {\em always converging\/} series (\ref{eq33}). However, as we have seen, we again have no solution, in this case even formal, unless the matching condition fixing $\kappa$ is fulfilled. Obviously, the uniqueness property implies that the introduction of the `branching' power function factor, as in (\ref{eq11}), is the only way to obtain a solution to \eqref{eq6} admitting representation in terms of convergent power series. Now, tracking back the relationship connecting \eqref{eq12} with the primary \eqref{eq1} and invoking the general theory of the latter applicable in the case of arbitrary continuous periodic $q(t)$ \cite{R11}, one can infer the following statements which however, in the present context, are only of the status of {\em conjectures\/} in view of the lack of their proof `from the first principles'; i.e., on the base of the properties of the discriminant function $\Xi$ following from its explicit definition. \begin{conjecture} \label{conjA} 1. There exists an open non-empty subset of the space of the problem parameters $n,\mu,\omega,\epsilon$ where \eqref{eq40} admits a real valued positive solution. 2. If real $\kappa$ solves \eqref{eq40}, $-\kappa$ is the solution as well. No more real roots exist. 3. Real roots of \eqref{eq40} obey the condition $|\kappa|\le (2\omega)^{-1}$. \end{conjecture} {\it Remark}: The last statement which is very important from viewpoint of applications, is nothing else but the form of the limitation on the rate of growth or decreasing of the functions $\tilde{x},\tilde y=x,y|_{z=e^{\mathrm{i}\omega t}}$ of the real variable $t$ implied by the inequalities (\ref{eq5a}) and the equation $$ x-\mathrm{i} y=-iz^{{\epsilon-1\over2}-\mathrm{i} \kappa} e^{-{\mu\over2}\left(z+{1\over z}\right)} E, $$ following from definitions. Notice that the latter clarifies the role of the {\em parity\/} parameter $\epsilon\; (=0\;\mbox{or}\; 1)$ which determines the multiplicity of the inverse to the map $S^1\ni z\to (x+\mathrm{i} y)/|x+\mathrm{i} y|\in S^1$ induced by the solution (\ref{eq34}). If $\epsilon=0$, the revolution along the circle $|z|=1$ leads to the reversing of the direction of the vector with components $(x,y)$ whereas in the case $\epsilon=1$ its direction is preserved. The above inverse map is double-valued in the former case and one-to-one in the latter one. More properties of the discriminant function can be inferred from the numerical experiments although, as opposed to the assertions of the Conjecture \ref{conjA}, they have, to date, no analytic arguments in their favor yet, even indirect. Nevertheless, the first item below is important enough from viewpoint of applications (seeming also plausible enough) to be {\em explicitly formulated\/} here. \begin{conjecture} \label{conjB} 1. {\bf Phase-lock criterion.} Equation \eqref{eq40} admits a real non-zero solution if and only if \begin{equation} \Xi(0;\dots)>0. \end{equation} This means in particular that $\Xi(0;\dots)$ is real. Moreover, the numerical study makes evidence that 2. $\Xi(\kappa;\dots)$ is real for real $\kappa$ (assuming the other parameters to be also real, of course). \end{conjecture} \section{Floquet solutions of double confluent Heun equation and involutive solution maps} Let us assume now that there exists a real positive solution $\kappa$ of \eqref{eq40}. With this $\kappa$, the function $E(z)$ defined by \eqref{eq34} verifies \eqref{eq12}. Let us consider the function $E_\#(z)$ defined through $E(z)$ as follows \begin{equation} \label{eq41} E_\#(z)= z^{2\mathrm{i}{\kappa}-\epsilon} \Big[ E'\big({1\over z}\big) +\Big(\big({n+\epsilon\over2}-\mathrm{i}{\kappa}\big)z-\mu\Big) E\big({1\over z}\big)\Big]. \end{equation} Then a straightforward computation shows that it verifies \eqref{eq12}, provided $E(z)$ does. It is worthwhile to note that the repetition of the transformation (\ref{eq41}) yields no more solutions to \eqref{eq12}. As a matter of fact, one has $_{\#}\circ_{\#}=(2\omega)^{-2} Id$. Thus $(2\omega)^{-1}{}_\#$ is the involutive map on the space of its solutions. Next, the functions $E(z)$ and $E_\#(z)$ are linearly independent for nonzero $\kappa$. Indeed, utilizing (\ref{eq12}), one finds % \begin{equation} \label{eq42} \begin{aligned} E_\# '(z) &= z^{ 2\mathrm{i}{\kappa}-\epsilon-1} \Big[ \Big(-{n+\epsilon\over2}+\mathrm{i}{\kappa} +\mu z\Big) E'({1\over z}) \\ &\quad +\Big( \mu\big({n+\epsilon\over2}-\mathrm{i}{\kappa}\big)(z^2+1) + \big( -({n+\epsilon\over2}-\mathrm{i}{\kappa})^2 +\lambda \big)z \Big) E({1\over z}) \Big]. \end{aligned} \end{equation} This reduction allows one to calculate the determinant of the linear transformation binding the pairs of functions $E_\#,E_\# '(z)$ and $E,E'(1/z)$ which proves equal to $(2\omega)^{-2}z^{2(-\epsilon + 2\mathrm{i}\kappa)}$ and is therefore nonzero. Hence $E_\#(z)$ is not zero, identically, (and may vanish at isolated points, at most, as well as $E(z)$). Finally, $E(z)$ is periodic on the unit circle centered at zero whereas $E_\#(z)$, for real $\kappa\not=0$, is not. Hence they are linear independent. The functions $E(z)$ and $E_\#(z)$ constitute therefore the fundamental system for \eqref{eq12} and any its solution %of the latter can be expanded in this basis with constant expansion coefficients. Consequently, for $\kappa\not=0$, the domain of general solution to \eqref{eq12} is $\Omega$, the universal covering of Riemann sphere with punctured poles. The analytic properties of $E,E_\#$ identify these solutions as (the unique pair of) the {\it Floquet solutions\/} of the reduced double confluent Heun equation under consideration. See \cite[section 2.4]{R13}. The function $E(z)$ obeys the important functional equation which can be derived as follows. A straightforward calculation shows that the right-hand side expression of \eqref{eq41} with the `branched' factor $z^{2\mathrm{i}{\kappa}}$ removed (i.e., $z^{-2\mathrm{i}{\kappa}} E_\#(z)$) satisfies the ODE which coincides with (\ref{eq12}) {\em up to the opposite sign of the parameter $\kappa$}. This means that, for real $\kappa$, $n,\mu,\omega$, the analytic function \begin{equation} \label{eq43} \hat E(z)=z^{-\epsilon} \Big[ \bar E'\big({1\over z}\big) +\Big(\big({n+\epsilon\over2}+\mathrm{i}{\kappa} \big)z-\mu\Big)\bar E\big({1\over z}\big) \Big], \end{equation} where $\bar E(z)=\overline{E(\overline z)}$, is the solution of \eqref{eq12} itself. As opposed to $E_\#$, this `yet another' solution has {\it the same\/} analytic properties as $E(z)$ and hence must coincide with it up to some numerical factor $C_C$; i.e., \begin{equation} \label{eq44} C_C E(z)= z^{-\epsilon} \Big[ \bar E'\big({1\over z}\big) +\Big(\Big({n+\epsilon\over2}+\mathrm{i}{\kappa} \big)z-\mu\Big)\bar E\big({1\over z}\big) \Big]. \end{equation} ($C_C$ may not vanish since otherwise $E_\#$ would also be identical zero.) This is the generalization of the similar property of the so called `Heun polynomials' established in Ref.~\cite{R11a}. The complex valued constant $C_C$ actually reduces to a single real constant. To show that, let us notice at first that if $E(z)$ verifies \eqref{eq12}, it follows from the latter and (\ref{eq44}), \begin{equation} \label{eq45} \begin{aligned} C_C E'(z)&= z^{-\epsilon-1} \Big[ \big( -{n+\epsilon\over2}-\mathrm{i}{\kappa} +\mu z \big)\bar E'(\frac{1}{z}) \\ &\quad + \Big( \mu {\epsilon+n\over2}(z^2+1) -\mathrm{i}\mu{\kappa}(z^2-1)\\ &\quad +\big( -{(n+\epsilon)^2\over4} -{\kappa^2}+\lambda \big)z \Big) \bar E(\frac{1}{z}) \Big]. \end{aligned} \end{equation} Evaluating \eqref{eq44}--\eqref{eq45} together with their complex conjugated versions at the point $z=z^{-1}=1$, one obtains four linear homogeneous equations binding the quantities $E(1),E'(1),\bar E(1),\bar E'(1)$ which may not vanish simultaneously. The corresponding consistency condition reads $ |C_C|^2=(2\omega)^{-2}$ implying \begin{equation} \label{eq46} C_C=(2\omega)^{-1}e^{\mathrm{i} C_c}, \end{equation} where $C_c$ is some {\em real\/} constant (actually, the function of the parameters $n,\mu,\omega,\epsilon$). It encodes all the monodromy data for \eqref{eq12}, essentially. It is straightforward to show that the transformation (\ref{eq44}) is also involutive. It manifests the specific symmetry in the behaviors of the function $E(z)$ in vicinities of the essentially singular points $z=0$ and $z^{-1}=0$. Remarkably, the possessing of this symmetry suffices itself for the fulfillment of \eqref{eq12}. Indeed, differentiating (\ref{eq44}) and taking into account (\ref{eq46}), one arrives at \eqref{eq12}. In a sense, (\ref{eq44}) together with stipulation for the analyticity of $E(z)$ can be considered as the equivalent to \eqref{eq12}. Additionally, (\ref{eq44}) implies anti-linear (involving complex conjugation) dependencies among the `distant' Laurent series coefficients $a_{-k}$ and $a_{k},a_{k+1}$ (whereas (\ref{eq14}) (or (\ref{eq16})) binds `nearby' $a_k,a_{k\pm1}$). In particular, it suffices to find all $a_k$ for $k>0$ and then $a_{-k}$ can be computed from the latter by means of a simple (anti-)linear transformation. \section{Essentially periodic and general solutions of overdamped Josephson junction equations} The connection between the functions $E(z),E_\#(z),\hat E(z)$ pointed out above is important for the consistent projecting the results concerning solutions of \eqref{eq12} to the level of original \eqref{eq1}--\eqref{eq2}. This procedure applying equations \eqref{eq4}, (\ref{eq8}), (\ref{eq9}), (\ref{eq11}) is straightforward and leads to the following conclusions. At first, the representation of the two special (and the most important) solutions to \eqref{eq1}--\eqref{eq2} for which the exponents $\exp(\mathrm{i}\varphi)$ are {\it periodic\/} (for brevity, we shall call such phase functions {\em essentially periodic}) follows. It reads \begin{gather} \label{eq47} \exp(-\mathrm{i}\varphi) = 2\mathrm{i}\omega \Big(z {E'(z)\over E(z)} +{n+\epsilon\over2}-\mathrm{i}{\kappa}-\mu z\Big),\\ \label{eq48} \exp(\mathrm{i}\varphi)=-2\mathrm{i}\omega \Big(z^{-1}{E'(z^{-1})\over E(z^{-1})}+{n+\epsilon\over2}-\mathrm{i}{\kappa} -\mu z^{-1}\Big), \end{gather} where \begin{equation} z=\exp(\mathrm{i}\omega t). \label{eq48a} \end{equation} For $\kappa>0$, the first of these formulas determines the asymptotic limit (the attractor) of a generic solution to \eqref{eq1}--\eqref{eq2} whereas the second solution is unstable (the repeller). It is important to emphasize that the functions $\varphi(t)$ defined by \eqref{eq47}, \eqref{eq48}, and \eqref{eq48a} are {\em real\/} and \eqref{eq44}, \eqref{eq46} are the crucial ones utilized in the calculation establishing this fact. At second, it is straightforward to carry out the `nonlinear superposition' of solutions (\ref{eq47}), (\ref{eq48}) operating with linear problem -related counterparts. The result is represented by the formula \begin{equation} \label{eq49} \begin{aligned} &\exp(\mathrm{i}\varphi) \\ &= -{\mathrm{i}\over2} \Big\{\cos\psi\cdot E(z) +\sin\psi\cdot z^{-\epsilon+2\mathrm{i}{\kappa}} \Big[ E'(z^{-1}) +\Big( \big({n+\epsilon\over2}-\mathrm{i}{\kappa} \big)z -\mu \Big)E(z^{-1}) \Big] \Big\} \\ &\times \Big\{\omega\cos\psi\cdot \Big[ z E'(z)+\big({n+\epsilon\over2}-\mathrm{i}{\kappa}-\mu z\big)E(z)\Big] + {1\over4\omega} \sin\psi\cdot z^{-\epsilon+1+2\mathrm{i}{\kappa}} E(z^{-1})\Big\}^{-1}, \end{aligned} \end{equation} where $\psi$ is an arbitrary real constant. More exactly, the set of all functions $\varphi$ %solutions described by (\ref{eq49}) is parameterized by points of $S^1$. %the unit circle. As opposed to (\ref{eq47}) and (\ref{eq48}) defined on the Riemann sphere with punctured poles, the function (\ref{eq45}) is defined on its universal covering, $\Omega$. Continuous (and then real analytic) function $\varphi(t)$ determined by this equation on $\Omega_1\in\Omega$, where $z$ is understood as $e^{\mathrm{i}\omega t}$, is just the general solution of \eqref{eq1}--\eqref{eq2} in the case of phase-lock. In particular, \eqref{eq47}, \eqref{eq48} represent the special instances of \eqref{eq49} arising for $\psi=0$ and $\psi=\pi/2$, respectively. As a consequence, asymptotic properties of the general solution mentioned above immediately follow. Indeed, as $t$ increases, the exponent (\ref{eq49}) is converging to (\ref{eq47}) and is taking off (\ref{eq48}) (unless it coincides with the latter). The two solutions described by (\ref{eq47}),(\ref{eq48}) are the only ones which are not affected by the translations %shifts $t\to {}t+2\pi\omega^{-1}$ (in the sense %that the exponents (\ref{eq47}),(\ref{eq48}) are kept unchanged) and preserve their functional form in asymptotics. Finally, at third, considering $\varphi$ defined by (\ref{eq47}) as analytic function of $z$ and taking in account \eqref{eq12}, one obtains \begin{equation} \label{eq50} \begin{aligned} {\mathrm{d}\varphi\over\mathrm{d} z} &= -\mathrm{i} z^{-2} \Big\{ z^3 \Big({E'(z)\over E(z)}\Big)^2 +z\big((1-z^2)\mu +z\big((\epsilon-1) -2\mathrm{i}{\kappa}\big)\big) {E'(z)\over E(z)} \\ &\quad -\Big(z+z\big({n-\epsilon\over2}+\mathrm{i}{\kappa}\big)-\mu\Big) \Big(\big({n+\epsilon\over2}-\mathrm{i}{\kappa}\big)-\mu z\Big) +{z\over4\omega^2}\Big\} \\ &\quad \times \Big\{ z {E'(z)\over E(z)} +\big({n+\epsilon\over2}-\mathrm{i}{\kappa}\big) -\mu z \Big\}^{-1} \end{aligned} \end{equation} On the unit circle, this $\varphi$ verifies \eqref{eq1}--\eqref{eq2}. It is therefore smooth (even real analytic) therein. Hence (\ref{eq50}) is continuous on $S^1$. Finally, since $\exp\mathrm{i}\varphi|_{z=\exp\mathrm{i}\omega t }$ is periodic, the following proposition holds. \begin{proposition} \label{prop2} The quantity \begin{equation} \label{eq51} k=(2\pi)^{-1} \oint {\mathrm{d}\varphi\over\mathrm{d} z} \mathrm{d} z, \end{equation} where ${\mathrm{d}\varphi/\mathrm d\, z}$ denotes the right-hand side in \eqref{eq50} and the integration is carried out over the circle $|z|=1$, is well defined and amounts to an integer. \end{proposition} This integer is the degree of the map ${S}^1\Rightarrow {S}^1$ induced by the function (\ref{eq47}). In physical applications, it is called {\it the phase-lock order\/} and is considered as the integer-valued function of the parameters which is obviously continuous and, thus, locally constant. Phase-lock order is involved in the formula representing the property of being `essentially periodic' for the phase function defined by (\ref{eq47}) (and asymptotically for a generic phase function) which reads $$ \varphi(t+2\pi\omega^{-1})=\varphi(t)+2\pi k \quad \forall t. $$ In a phase-lock state of JJ, the uniformly distributed discrete levels of averaged voltage equal to $k\cdot (\hbar\omega/2e)$ for some $k=0,\pm 1,\pm 2\dots$ are observed. \begin{conjecture} \label{conjC} Any integer map degree \eqref{eq51} is realized on some non-empty open subset of the space of the problem parameters $n,\mu,\omega,\epsilon$ (or, equivalently, $A,B,\omega,\epsilon$). \end{conjecture} This assertion is closely cognate to the item 1 of the above Conjecture \ref{conjA}. \section{Summary} It the present work, the general solution of the overdamped Josephson junction equation (\ref{eq1}) was obtained for the (co)sinusoidal right-hand side function (\ref{eq2}) in the case of one of three possible asymptotic behaviors known as the phase-lock mode. The solution is represented in explicit form in terms of the Floquet solution of the particular instance (arising in case of the vanishing of one of the four free constant parameters) of the double confluent Heun equation. The Floquet solution of double confluent Heun equation is represented in terms of the Laurent series whose coefficients are determined by explicit formulas involving convergent infinite products of $2\times 2$ matrices with a single zero element tending to the idempotent matrix (\ref{eq26c}). The derivation % method presupposes the existence of a real solution of the transcendental equation (\ref{eq40}) which is equivalent to the claim of realization of the phase-lock mode for the given parameter values. The plausible criterion of its existence (i.e., {\em the phase-lock criterion}) is conjectured. The derivation of general solution of \eqref{eq1}--\eqref{eq2} consists of a number of steps. In a concise form they can be described as follows. \begin{itemize} \item The investigation of the basic properties of \eqref{eq1} for arbitrary periodic (sufficiently regular) $q(t)$ allows one to establish the division of the space of the problem parameters into the two open subsets of which one corresponds to the phase-lock property of the \eqref{eq1}--\eqref{eq2} solutions whereas another corresponds to their pseudo-chaotic behavior revealing no %apparent stable periodicity. In the (lower-dimensional) complement to these areas the intermediate behavior is observed. The corresponding results are discussed in sufficient details in Ref.~\cite{R11}. \item The next important point is the intimate connection (first mentioned by Buchstaber; see, e.g., Ref.~\cite{R12}) between \eqref{eq1} and a simple {\em linear} system of the two first order ODEs (\ref{eq5}). For the (co)sinusoidal right-hand side function (\ref{eq2}), the latter takes the form (\ref{eq3}). \item At the next step, the transformation (\ref{eq8}) was found which converts the linear system (\ref{eq3}) to a particular instance of the double confluent Heun equation{} (\ref{eq10}). \end{itemize} Generally speaking, it could be solved by means of the expansion of the (modified) unknown function in Laurent series \cite{R13},\cite{R11a} centered at the singular points $z=0,z^{-1}=0$ but preliminarily the additional simple but important transformation has to be carried out: \begin{itemize} \item the `branched' power factor involving unspecified constant ($\kappa$-dependent contribution in (\ref{eq11})) is introduced. \item Additionally, the discrete `parity' parameter $\epsilon$, assuming either the value 0 or the value 1 (also absent in the equation to be solved), is introduced in the aforementioned power factor. This modification proves necessary for the subsequent exhaustive `indexing' of the solution space. \item After that, the standard technique of the power expansion leads to the `endless' sequence of the 3-term constraints (\ref{eq14}) (or, which is the same, (\ref{eq16})) imposed on the unknown series coefficients. \item The next step is the devision of the set of power series coefficients into the two subsets. The non-negative-index-value coefficients and non-positive-index-value ones are treated separately, solving the separate subsets of the equations (\ref{eq14}) and (\ref{eq16}), respectively, for $k\ge 1$. Analyzing them, the application of the continued fraction technique leads, after some transformations, to the `explicit' formulas (\ref{eq22}),\eqref{eq23} for the series coefficients involving infinite products of $2\times 2$ matrices converging for large index values to the idempotent matrix (\ref{eq26c}). This convergence is sufficiently fast to imply the convergence of the matrix products and, accordingly, the finiteness of the series coefficients. Moreover, the associated estimates make evident the existence of the absolutely converging majorants for the resulting Laurent series. Therefore, they actually determine an analytic function serving the Floquet solution of double confluent Heun equation. It proves representable as the sum of the two entire functions of the arguments $z$ and $z^{-1}$, respectively. \item The procedure producing Laurent series coefficients noted above proves however suffering from the improper introduction of a kind of vicious singularities arising as zeroes in denominator which appear for some special parameter values. They are eliminated my means of multiplication of the `raw' coefficient expressions by some $z$-independent (but parameter dependent) factors given in explicit form. \end{itemize} One more {\em remark\/} is here in order. The case where the above `regularizing' factor actually enters the play exactly corresponds to the situation where the problem solution can be built in terms of polynomials (it was discovered in Ref.~\cite{R11a}). In brief, the `vanishing' factors `neutralizing' diverging terms with vanishing denominators, eliminate, in pass, the terms which do not diverge. Ultimately, only finite number of `originally diverging' terms survive which lead to the reduction of the transcendental function describing generic solution to a polynomial. (The details of the realization of this reduction have not been properly elaborated yet.) \begin{itemize} \item Now the `solution candidate' for \eqref{eq12} can be represented as the analytic function (\ref{eq34}) which is well-defined for any parameter values. However, at the price of automatic convergence of the power series it has been built upon, generally speaking, it does not verify \eqref{eq12}. The equation is fulfilled if and only if the fulfillment of \eqref{eq40}, which is the transcendental equation for the still unspecified parameter $\kappa$, is stipulated. \end{itemize} At this stage, having solved \eqref{eq40}, a single solution (the Floquet solution) of double confluent Heun equation can be regarded as been explicitly constructed. \begin{itemize} \item The invariance of the space of solutions of double confluent Heun equation under consideration with respect to the transformation (\ref{eq41}) allows one to immediately obtain the fundamental system of its solutions in terms of the single Floquet solution noted above. \item The existence of automorphism represented by (\ref{eq44}) expresses the important intrinsic property of the Floquet solution of double confluent Heun equation. It is used for the derivation of the explicit representation of the exponent $\exp(\mathrm{i}\varphi)$ specifying the {\em real valued\/} phase function $\varphi$. It is obtained by means of the restriction of the analytic function (see (\ref{eq49})) from the universal covering of the Riemann sphere with punctured poles to the embedded universal covering of $S^1$. This trick yields the general solution of \eqref{eq1}, \eqref{eq2} {\em in the case of phase-lock}. \item Finally, employing the analytic properties and the periodicity (on $S^1$) of $\exp(\mathrm{i}\varphi)(z)$, the formula (\ref{eq51}) involving Floquet solution of double{} confluent Heun equation follows which gives the degree of the map ${S}^1\Rightarrow {S}^1$ it induces (the winding number) also known in application fields as the phase-lock order. \end{itemize} It is worth noting in conclusion that all the major constructions derived above admit a straightforward algorithmic implementation which have been used for the numeric verification of the relevant relationships. \section{Appendix: Outline of the proof of Theorem \ref{thm1}} Equation (\ref{eq24}) implies $ \beta_j=\mu^2\alpha_{j-1} Z_{j}^{-1}$, $\delta_j=\mu^2\gamma_{j-1} Z_{j}^{-1}$ and hence the asserted properties of the sequences $\beta_j,\delta_j$ follow from the existence of finite limits for the sequences $\alpha_j,\gamma_j$. As to the sequences $\alpha_j$ and $\gamma_j$, their elements have to obey the identical decoupled 3-term recurrence relations which, for $\alpha$'s, read \begin{equation} \label{eq52} \alpha_j= (1+\lambda Z_{j+(\epsilon-1)/2}^{-1})\alpha_{j-1} +\mu^2 Z_{j-\tilde\epsilon+(\epsilon-1)/2}^{-1}\alpha_{j-2}, \end{equation} where $\tilde\epsilon=1$ for the upper choice of $\sigma$ in (\ref{eq25}) and $\tilde\epsilon=0$ for the lower $\sigma$ choice therein. It suffices to consider the $\alpha$-sequence case. Evidently, for every integer $j_0>0$ and $l>0$ any solution of equations (\ref{eq52}) can be represented in terms of the decomposition \begin{equation} \label{eq53} \alpha_{j_0+l}=(1+p_{j_0,l})\alpha_{j_0-1}+q_{j_0,l}\alpha_{j_0-2}, \end{equation} for some coefficients $p_{j,l},q_{j,l}$ independent on the `starting' terms $\alpha_{j_0-1},\alpha_{j_0-2}$. Applying (\ref{eq52}), it is straightforward to derive the reduction \begin{equation} \label{eq:53a} q_{j_0,l+1}= (1+p_{j_0+1,l} )\mu^2 Z_{j_0-\tilde\epsilon}^{-1}, \end{equation} whereas for $p_{*,*}$ one gets the following recurrence relation: \begin{equation} \label{eq54} p_{j_0,l+1}= p_{j_0+1,l}+ \lambda(1+p_{j_0+1,l}) Z_{j}^{-1} +\mu^2(1+p_{j_0+2,l-1})Z_{j+1-\tilde\epsilon}^{-1}. \end{equation} Equation (\ref{eq54}) is equivalent to (\ref{eq52}), essentially, but it possesses the advantage of being endowed with the standard `initial conditions' \begin{equation} \label{eq55} p_{j_0,-1}= 0, p_{j_0,-2}=-1 \end{equation} independent of the specific `starting' values $\alpha_{j_0-1},\alpha_{j_0-2}$. Furthermore, one gets \begin{equation} \label{eq56} q_{j_0,-1}= 0, q_{j_0,-2}=1. \end{equation} It proves convenient to carry out one more rearrangement of unknowns introducing the differences \begin{equation} \label{eq57} \Delta p_{j_0,l}=p_{j_0,l+1}-p_{j_0,l}. \end{equation} which obey the own `initial conditions' \begin{gather} \label{eq58} \Delta p_{j_0,-2}= 1, \\ \label{eq58a} \Delta p_{j_0,-1}= \lambda Z^{-1}_{j_0}, \end{gather} and analogous recurrence relations \begin{equation} \label{eq59} \Delta p_{j_0,l+1}= \Delta p_{j_0+1,l}+ \lambda \Delta p_{j_0+1,l} Z_{j_0}^{-1} +\mu^2 \Delta p_{j_0+2,l-1}Z_{j_0+1-\tilde\epsilon}^{-1}. \end{equation} Now, summing up the subset of all the equations (\ref{eq59}) distinguished by the common sum of indices at the left and taking into account (\ref{eq58a}), all but two `free' $\Delta p$-terms cancel out and one obtains \begin{equation} \label{eq60} \begin{aligned} \Delta p_{j_0,l+1} &= \lambda Z^{-1}_{j_0+2+l} +\lambda \sum_{m=0}^{l+1} \Delta p_{j_0+1+m,l-m} Z_{j_0+m}^{-1} \\ &\quad +\mu^2 \sum_{m=0}^{l+1} \Delta p_{j_0+2+m,l-1-m}Z_{j_0+m+1-\tilde\epsilon}^{-1}. \end{aligned} \end{equation} In the sums, the second index of $\Delta p_{*,*}$ is everywhere less than the same index at the left that allows to apply the method of mathematical induction. For the `starting' values $-1,0$ of the second index one has \begin{gather*} Z^{}_{j_0}\Delta p_{j_0,-1} = \lambda, \\ Z_{j_0+1} \Delta p_{j_0,0} = \lambda(1+\lambda Z^{-1}_{j_0}) +\mu^2 Z_{j_0+1} Z^{-1}_{j_0+2-\tilde\epsilon}. \end{gather*} Therefore, for $l=-1,0$ there exist the finite limits $\lim_{j_0\to \infty}|Z^{}_{j_0+l+1}\Delta p_{j_0,l}|$. As a consequence, for these $l$'s, \begin{equation} \label{eq61} |\Delta p_{j_0,l}|< \tilde N |Z_{j_0+l+1}|^{-1} \end{equation} for appropriate constant $\tilde N$ which is convenient to choose $>1$. Let us consider this fact as the starting point of mathematical induction and assume that for some integer $l_0\ge 0$ and any integer $l$ from the interval $[-1,l_0]$ (\ref{eq61}) holds true. We may apply it for the estimating from above of the quantity $|\Delta p_{j_0,l_0+1}|$. This can be realized making use of the `decomposition' (\ref{eq60}) and the following elementary inequalities \begin{gather}\label{eq62} \sum_{m=j_0}^{L+j_0+1} |Z_{m}|^{-1}<{1+|n+1|^{-1} \over j_0-|n+1|/2+(\epsilon-1)/2 }, \\ \sum_{m=j_0+1-\tilde\epsilon}^{L+j_0+2-\tilde\epsilon} |Z_{m}|^{-1}< {1+|n+1|^{-1}\over j_0-|n+1|/2+(\epsilon-1)/2 }, \end{gather} where $L>0$ (and $n\not=-1$) . These imply the inequalities \begin{align*} | \Delta p_{j_0,l_0+1}| &\le |\lambda| |Z_{j_0+l_0+2}|^{-1} \\ &\quad +\tilde N |Z_{j_0+l_0+1}|^{-1} \Big( |\lambda| \sum_{m=0}^{l_0+1} |Z_{j_0+m}|^{-1} + |\mu^2| \sum_{m=0}^{l_0+1} |Z_{j_0+m+1-\tilde\epsilon}|^{-1} \Big) \\ &< |Z_{j_0+l_0+2}|^{-1} \Big(1+\tilde N{|Z_{j_0+l_0+2}|\over |Z_{j_0+l_0+1}|} {(|\lambda|+|\mu^2|)(1+|n+1|^{-1})\over (j_0-|n+1|/2+(\epsilon-1)/2)} \Big). \end{align*} Since we assumed $\tilde N>1$, there exists the bound of lower index values such that for any $j_0$ exceeding it the factor in brackets is less than $\tilde N $, and then the above inequalities imply $|\Delta p_{j_0,l_0+1}|<\tilde N |Z_{j_0+l_0+2}|^{-1}$. The induction step is accomplished. (\ref{eq61}) is therefore established for sufficiently large $j_0$ and arbitrary $l\ge0$. Increasing $\tilde N$ if necessary, (\ref{eq61}) proves valid for arbitrary $j_0$. In view of this property, one sees that the sum $\sum_{l=0}^\infty\Delta p_{j_0,l}$ has the majorant const$\times\sum_l|Z_{j_0+l+1}|^{-1}$ and thus converges itself. The sequence of its partial sums $\sum_{m=0}^l\Delta p_{j_0,m}=p_{j_0,l+1}-p_{j_0,-1}=p_{j_0,l+1}$ also converges as $l\to \infty$. Moreover, in view of (\ref{eq55}), (\ref{eq57}), (\ref{eq62}), (\ref{eq61}) one has the $l$-uniform bound \begin{equation}\label{eq62a} |p_{j_0,l+1}| <\tilde N\sum_{m=0}^l |Z_{j_0+m+2}|^{-1} <{\tilde N(1+|n+1|^{-1}) \over j_0+1+(\epsilon-1)/2}. \end{equation} It follows from \eqref{eq53}, \eqref{eq:53a} \begin{equation} \label{eq53b} \alpha_j= \alpha_{j_0+l}=(1+p_{j_0,l})\alpha_{j_0-1} + (1+p_{j_0+1,l-1} )\mu^2 Z_{j_0-\tilde\epsilon}^{-1}\alpha_{j_0-2} \end{equation} and the convergence of $\alpha$-sequence follows from the convergence of $p_{*,l}$. % as $l\to \infty$. Moreover, one has \begin{equation} \label{eq53c} \lim\alpha_j-\alpha_{j_0-1}= \lim_l p_{j_0,l}\cdot \alpha_{j_0-1} +(1+\lim_l p_{j_0+1,l} )\mu^2 Z_{j_0-\tilde\epsilon}^{-1}\cdot\alpha_{j_0-2}. \end{equation} The factors in front of the first and the second terms to the right scales as $j_0^{-1}$ and $j_0^{-2}$, respectively. We may therefore write down the inequality \begin{equation*} |\lim\alpha_j-\alpha_{j_0-1}|=N\max(|\alpha_{j_0-1}|, |\alpha_{j_0-2}|)j_0^{-1}, \end{equation*} where $N$ may depend on the parameters $n,\lambda,\mu,\kappa,\epsilon$ but not on the specific specimen of $\alpha$-sequence. This obviously implies the inequality (\ref{eq26}). It has also to be noted in conclusion that the case $n=-1$ formally falling off the above speculation does not actually correspond to an exceptional situation. Although inequalities (\ref{eq62}) formally fail, similar ones differing from \eqref{eq62} in the values of `constant' ($j_0$-independent) terms alone can be derived. The further reasoning holds true and leads to the same conclusions. \begin{thebibliography}{00} \bibitem{R1} B. D. Josephson, \emph{Possible new effects in superconducting tunelling}, {Phys.\ Lett.} {1}(1962) 251-253. \bibitem{R2} B. D. Josephson, Coupled superconductors, { Rev.\ Mod.\ Phys.} {36} (1964) 216-220. \bibitem{R3} D. E. McCumber, \emph{Effect of ac Impedance on dc Voltage-Current Characteristic of Superconductor Weak-Link Junctions}, {J.\ Appl.\ Phys.} {39}(1968) 3113-3118. \bibitem {R4} W. C. Stewart, \emph{Current-voltage characteristics of Josephson junctions}, {Appl.\ Phys\ Lett.} {12},(1968) 277-280. \bibitem{R5} A.\ Barone and G.\ Paterno, \emph{Physics and Applications of the Josephson Effect}, Willey, N.Y., 1982. \bibitem{R6} A. H. Hamilton, \emph{Josephson voltage standards}, {Rev.\ of Sci.\ Instr.}, {71} (2000) 3611-3623. \bibitem{R7} S. P. Benz and C. A. Hamilton, \emph{Application of the Josephson effect to Voltage Metrology}, {Proc IEEE} {92} (2004) 1617-1629. \bibitem{R8} J. Niemeyer, \emph{Josephson Voltage Metrology}, in: {Abstract Booklet EUCAS 2005}, Vienna (2005) 83-84. \bibitem {R9} A. Kemooinen, J. Nissil\"a, K. Ojasalo, J. Hassel, A. Manninen, P. Helist\"o and H. Sepp\"a, \emph{AC Voltage Standard based on an Externally-Shunted SIS Josephson Junction Array}, in: Abstract Booklet EUCAS 2005, Vienna (2005) 336. \bibitem{R10} R. L. Kautz; \emph{Noise, chaos, and Josephson voltage standard}, Rep. Prog. Phys., {59} (1996) 935-992. \bibitem{R11} S. I. Tertychniy, \emph{Long-term behavior of solutions to the equation $\dot \phi + \sin\phi=f$ with periodic $f$ and the modeling of dynamics of overdamped Josephson junctions}. Preprint {Arxiv:math-ph/0512058} (2005). \bibitem{R11a} S. I. Tertychniy, \emph{The modeling of a Josephson junction and Heun polynomials}. Preprint {Arxiv:math-ph/0601064} (2006). \bibitem{R12} V. M. Buchstaber, O. V. Karpov and S. I. Tertychniy, \emph{Quantum Josephson D/A converter driven by trains of short $2\pi$-pulses}, in: {Conference Digest CPEM 2002}, Ottawa, (2002) 502-503 (2002). \bibitem {R13} D. Schmidt, G. Wolf, \emph{Double confluent Heun equation, in: Heun's differential equations}, Ronveaux (Ed.) Oxford Univ. Press, Oxford, N.Y., (1995), Part C. \bibitem{R14} S. Yu. Slavyanov and W. Lay, \emph{Special functions: A Unified Theory Based on Singularities}, Nevskiy dialect, SPb, (2002), in Russian; English edition: Oxford Univ. Press, Oxford, N.Y., (2000). \end{thebibliography} \end{document}