\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 255, pp. 1--18.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/255\hfil Unification of $q$-difference equations]
{Unification of integrable $q$-difference equations}

\author[B. S\.il\.ind\.ir, D. Soyo\u{g}lu  \hfil EJDE-2015/255\hfilneg]
{Burcu S\.il\.ind\.ir,  Duygu Soyo\u{g}lu}

\address{Burcu S\.il\.ind\.ir \newline
Department of Mathematics,
Dokuz Eyl\"ul University, Tınaztepe Campus, 
35160, Buca, \.Izmir, Turkey}
\email{burcusilindir@gmail.com}

\address{Duygu Soyo\u{g}lu \newline
Department of Mathematics, Izmir University of Economics,
35330, Bal\c{c}ova, Izmir, Turkey}
\email{duygusoyoglu@gmail.com}


\thanks{Submitted September 3, 2015. Published October 2, 2015.}
\subjclass[2010]{37K10, 37K40, 39A13, 39A14}
\keywords{Integrability; $q$-soliton solutions; $q$-difference KdV equation;
\hfill\break\indent  $q$-difference-$q$-difference Toda equation; $q$-difference
sine-Gordon equation; \hfill\break\indent Hirota direct method}

\begin{abstract}
 This article presents a unifying framework for
 $q$-discrete equations. We introduce a generalized $q$-difference
 equation in Hirota bilinear form and develop the associated
 three-$q$-soliton solutions which are described in polynomials of
 power functions by utilizing Hirota direct method. Furthermore, we
 present that the generalized $q$-difference soliton equation
 reduces to $q$-analogues of  Toda, KdV and sine-Gordon equations
 equipped with their three-$q$-soliton solutions by appropriate
 transformations.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks


\section{Introduction}

The concept of integrability possesses a key position in the field
of theoretical and mathematical physics. In the landmark article
\cite{hirota1}, Hirota introduced a very essential method, the
so-called \emph{Hirota direct method} which allows not only to
construct multi-soliton solutions or some special type of
solutions, but also to investigate the integrability criteria of a
given nonlinear evolution equation \cite{grammaticos, gurses,
hierentina0, hierentina1, hierentina2}. Another important hallmark
of Hirota's method over other methods; such as inverse scattering
transform \cite{Gardner}, or B\"{a}cklund transformation
\cite{wahl}, is the fact that it is algebraic rather than
analytic. The intrinsic feature of the method is to convert a
nonlinear partial differential or difference equation to
\emph{Hirota bilinear form} which is expressed by means of a
polynomial in \emph{Hirota-$D$ derivative operator}. In the
literature, it is conjectured that all integrable nonlinear
evolution equations can be revealed in Hirota bilinear forms while
the converse is not true. As an aside, notice that the equations
in Hirota bilinear form equipped with three soliton solutions are
defined to be \emph{Hirota integrable} and they are widely
considered to be integrable. \cite{Ma,Ma0}. In the present paper,
we stick to Hirota integrability definition.

B\l{a}szak et al. \cite{bgss} accomplished that all
discrete systems that are generated by distinct vector fields are
not globally equivalent. Besides, it is concluded that
$q$-difference systems on $\mathbb{R}^-$, are not isomorphic to
lattice systems on $\mathbb{R}$. This in-equivalence beget a
deeper analysis on the $q$-discretization of Toda lattice system
of equations. In the literature, Hirota direct method was applied
to a vast variety of differential or difference type of equations.
In \cite{burcu5}, the method was also shown to be applicable to
$q$-difference equations such as $q$-difference-$q$-difference and
differential-$q$-difference Toda equations to obtain multi-soliton
solutions. We observed that the constructed solutions obey
classical soliton attitudes as well as they have power
counterparts for $q$-discrete variables. We defined such solutions
as \emph{q-soliton solutions}.


In this same vein, one can present $q$-analogues of several
soliton equations. Instead, it is of great interest to intensify
on a single soliton equation that gathers various $q$-discrete
type of equations under one roof. The aim of this paper is to
create a unifying framework for $q$-discrete equations and analyze
the applicability of Hirota direct method to develop their
multi-soliton solutions. This formalism is based on the
$q$-discretization of equations determined by $q$-forward jump
operator. The fundamental feature of the framework is to introduce
appropriate $q$-deformed Hirota bilinear forms in a way that they
recover continuous Hirota bilinear forms. Significantly,
$q$-deformed Hirota bilinear forms enjoy a key position as they
provide not only $q$-analogues of corresponding equations but also
their multi-soliton solutions. For this purpose, we present a
generalized $q$-difference soliton equation which comprises
$q$-analogues of various soliton equations such as Toda, KdV and
sine-Gordon equations. We develop its three soliton solutions  by
the use of Hirota direct method and we stress that the constructed
solutions appear to be $q$-solitons. Unlike the discrete
generalized Toda equation \cite{Hirota10} whose solutions are of
exponential type, this generalized $q$-difference equation admits
soliton solutions that are expressed in terms of a polynomial in
power functions. This is a consequence of non-commutativity
between $q$-forward shift and exponential transformation.

It is possible to present $q$-discretization of  a continuous
equation in several approaches.  They can be derived  by the frame
of $q$-derivative operator, or analogously as in the present
article, they can be constructed by the use of $q$-forward jump
operator. In \cite{burcu5}, a counter-example to Hirota's theorem
\cite{hirotatheorem} is revealed by examining
$q$-differential-$q$-difference version of Toda equation which is
expressed in terms of $q$-derivative operator. Although it can be
presented in Hirota bilinear form and satisfies the sufficient
conditions to admit at least two-soliton solutions, it possesses
only solitary wave like a solution determined by $q$-exponential
function. The nonexistence of further $q$-exponential type of
soliton solutions is a consequence of the lack of additive
property of $q$-exponentials \cite{sch} and lack of
time-independency on interaction terms. Even though Hirota direct
method is applicable to $q$-difference equations, it fails to
produce $q$-exponential type of multi-soliton solutions for
$q$-differential equations governed by $q$-derivative operator.
Accordingly, in Section 4, we introduce  $\Delta$-Hirota
$D$-operator to analyze $\Delta$-differential equations on
arbitrary time scales. In this more general case, classical Hirota
perturbation does not provide multi-soliton solutions for any
difference equation studied on discrete intervals with
non-constant \emph{graininess } (e.g. $q$-differential
equations) or on such time scales. To be more precise, here we
conjecture that other than the unifying framework proposed in the
present article, it is not possible to acquire another unifying
approach via classical Hirota perturbation for integrable
$\Delta$-differential equations on time scales with non-constant
graininess.

The current article is organized as follows: In Section 2, we
present the preliminary notions regarding Hirota $D$-operator,
$q$-forward jump operator and $q$-exponential identity which
allows to convert $q$-discrete equations into Hirota bilinear
forms. Section 3 is devoted to derive three-$q$-soliton solutions
of the generalized $q$-difference soliton  equation as existence
of three-soliton solutions is a benchmark for integrability. In
Section 4, we introduce the proper reductions on this generalized
equation yielding to the $q$-difference-$q$-difference Toda,
$q$-difference-$q$-difference KdV and $q$-difference sine-Gordon
equations. We intimately demonstrate the notion of continuous
limit arising in association between $q$-deformed and classical
Hirota bilinear forms. Furthermore, we present $q$-soliton
solutions of the considered equations explicitly, resulting from
the reductions on the findings of Section 3.

\section{Preliminaries}

 In this section,  to declare the source of
Hirota bilinear forms of $q$-discrete equations, we intend to
review \emph{q-exponential identity}, stated in \cite{burcu5}.
For this purpose, we first present Hirota $D$-derivative
(operator)$D:S\times S\to S$
\begin{equation}
\begin{aligned}
&[D_x^{m_1}D_t^{m_2}\dots]\{f\cdot g\} \\
&=[(\partial _x-\partial
_{x'})^{m_1}(\partial _t-\partial
_{t'})^{m_2}\dots]f(x,t,\dots)\cdot g(x',t',\dots)|
_{x'=x,t'=t,\dots},
\end{aligned} \label{dop}
\end{equation}
where $S$ is a space of differentiable functions
$f:\mathbb{C}^n \to \mathbb{C}$, $x,t,\dots $ are independent variables and
$m_i\in\mathbb{Z}^{+}$, for  $i\geq 1$. Indeed, such
differential operator \eqref{dop} represents a novel calculus
obeying the following properties:

\begin{proposition}[\cite{hirota9}] \label{eq:1}
Let $f(x,t,\dots)$ and $g(x,t,\dots)$ be differentiable functions and
$P(D)$ be any polynomial in $D$, then
\begin{itemize}
\item[(i)] $P(D)\{f\cdot 1\} = P(\partial)f$,
$P(D)\{1\cdot f\} =P(-\partial)f$;
\item[(ii)]  $P(D)\{f\cdot g\} = P(-D)\{g\cdot f\}$,
\end{itemize}
hold, where $\partial$ denotes the ordinary differential
operator.
\end{proposition}

Notice that, Hirota $D$-derivative \eqref{dop} can also be
introduced by the frame of the exponential identity
\begin{equation}
\exp(hD_x)f(x)g(x) =f(x+h)g(x-h),\label{eq:4}
\end{equation}
 which is very beneficial in deriving Hirota bilinear forms of
differential-difference type of equations and B\"{a}cklund
transformations \cite{hirota9}. Here $f,g$ are smooth functions of
$x$ and $h$ is a parameter.

Suppose that we have continuous one-parameter group of
diffeomorphisms $\mathbb{R}\ni h\mapsto \sigma_h$, acting as
forward jump operators. There is one-to-one correspondence between
one-parameter group of transformations and their
\emph{infinitesimal generators}. Moreover, such diffeomorphisms
are determined by exponentiation of the infinitesimal generator as
they appear to be the solutions of the ordinary differential
equations \cite{olver}. Accordingly, for such $ \sigma_h$ we have
\cite{bgss}
\begin{equation}
\sigma_h(x)=e^{h\chi(x)\partial_x}x
\end{equation}
if and only if
\begin{equation}
e^{h\chi(x)\partial_x}f(x)=f(e^{h\chi(x)\partial_x}.x)=f(\sigma_h(x)).
\end{equation}
 Here $h$ is a positive deformation parameter, $f(x)$ is
a smooth function and the vector field $\chi(x)\partial_x$ is
infinitesimal generator.


It is possible to suggest the infinitesimal generators of the
form $\chi(x)\partial_x=x^{1-n}\partial_x$ on
$\mathbb{R}$.  The choice $n=1$ yields the \emph{forward jump
operator} of lattice type
\begin{equation}
\sigma_h(x)=e^{h\partial_x}x=x+h\quad\Leftrightarrow\quad
e^{h\partial_x}f(x)=f(x+h).
\end{equation}
In addition, the choice $n=0$ gives rise to the
\emph{$q$-forward jump operator}
\begin{equation}
\sigma_h(x)=e^{hx\partial_x}x=e^hx=qx\quad
(q\equiv e^h)\quad\Leftrightarrow\quad
e^{hx\partial_x}f(x)=f(qx).\label{shiftingrelation1}
\end{equation}

\begin{definition} \rm
 We define the $q$-forward jump operator $E_q$ acting on any smooth
function $f(x)$ as
\begin{equation}
E_q(f(x)):=e^{hx\partial_x}f(x)=f(qx),\label{q}
\end{equation}
where $x\in\mathbb{R}$ and $h$ is a deformation parameter. Similarly, $q$-backward jump operator is
introduced to act as
\begin{equation}
E_q^{-1}f(x):=e^{-hx\partial_x}f(x)=f(\frac{x}{q}).\label{qjump}
\end{equation}
\end{definition}

\begin{proposition} \label{limitoperator}
The continuous limit of $q$-forward jump
operator $E_q$ is
\begin{equation}
\lim_{q\to 1}E_q(x)=\lim_{h\to 0}\sigma_h(x)=x.
\end{equation}
\end{proposition}

\begin{proof}
The limit process can be presented by the expansion
into Taylor series with respect to $h$ near zero. That is
\begin{equation*}
\lim_{q\to 1}E_q(x)=\lim_{h\to
0}[x+hx\partial_x(x)+\frac{h^2}{2}(x\partial_x)^2(x)+O(h^3)]=x.
\end{equation*}
As well as
\begin{equation*}
\lim_{q\to 1}E_q(x)=\lim_{q\to 1}qx=x.
\end{equation*}
\end{proof}

We emphasize that all discrete systems generated by infinitesimal
generators $\chi(x)\partial_x$ are not equivalent.
To be more precise, if we consider $\chi(x)=x^{1-n}$(discrete case) where
$n\neq 0$ is odd, and $\chi'(x')=1$, it is possible to find a
local transformation
$x'= {\frac{1}{n}x^n}$
which is a bijection on $\mathbb{R}- \{0\}$. Thus, the discrete systems given
by $\chi(x)=x^{1-n}$ with odd $n$, turns out to be Toda lattice
type of equations. On the other hand, if we consider $\chi(x)=x$
($q$-difference case) and $\chi'(x')=1$, we have the
transformation $x=e^{x'}$  and it is not a bijection if
$x\in\mathbb{R}^-$.

Therefore, there does not exist an isomorphism between lattice
systems on $\mathbb{R}$ and $q$-difference systems on
$\mathbb{R}^-$ obtained by $q$-forward jump operators $E_q$ given
by \eqref{q}.

Inspired by this in-equivalence, throughout this work we present
$q$-difference equations that are determined by $q$-forward jump
operators $E_q$.  In order to construct Hirota bilinear forms of
$q$-difference type of equations, it is worthwhile to state the
$q$-analogue of exponential identity which is given in terms of
$q$-forward and $q$-backward jump operators.

\begin{theorem}[\cite{burcu5}]
 Let $f(x)$, $g(x)$ be continuously differentiable
functions, then the $q$-exponential identity
\begin{equation}
e^{hxD_x}f(x)g(x)=f(qx)g(\frac{x}{q})=E_qf(x)E_q^{-1}g(x), \quad
x\in \mathbb{R},\label{eq:5}
\end{equation}
holds where $h$ and $q$ are quantum parameters related
as $q=e^h$.
\end{theorem}

 For the proof of the above theorem we refer \cite{burcu5}.

\section{$q$-soliton solutions}

  In this section, we propose a generalized $q$-difference
soliton equation, namely a $q$-discrete analogue of Hirota-Miwa
equation
\begin{equation}
P(D_1,D_2,D_3)\{f\cdot f\}=\sum_{i=1}^3\lambda_i\cosh{(D_i)}\{f\cdot f\}=0,
\label{hirotabilinearform}
\end{equation}
where $\lambda_i$' s are arbitrary parameters, $D_i$' s are linear
combinations of the operators $xD_x$, $yD_y$, $tD_t$, i.e.,
\begin{equation}
D_i=a_itD_t+b_ixD_x+c_iyD_y,\quad a_i,b_i,c_i\in\mathbb{R},
\quad i=1,2,3.\label{operator}
\end{equation}
  Note that \eqref{hirotabilinearform} reproduces various
$q$-discretized soliton equations, by utilizing proper
identifications and reductions of parameters. The associated
transformations from \eqref{hirotabilinearform} into several
$q$-difference equations are intimately demonstrated in Section 4.

\begin{corollary} \label{cor3.1}
The $q$-exponential identity
in three variables $t,x,y\in \mathbb{R}$
\begin{equation}
\exp(a_itD_t+b_ixD_x+c_iyD_y)f(t,x,y)g(t,x,y)
 =f(q_it,p_ix,r_iy)g(\frac{t}{q_i},\frac{x}{p_i},\frac{y}{r_i})
\label{qexponentialidentity}
\end{equation}
holds, for any continuously differentiable functions $f$ and $g$, equipped with
the relations between quantum parameters $e^{a_i}=q_i$,
$e^{b_i}=p_i$ and $e^{c_i}=r_i$, for all $i=1,2,3$,
respectively.
\end{corollary}


To acquire multi-$q$-soliton solutions of the equation
\eqref{hirotabilinearform}, we utilize the so-called Hirota
perturbation. Upon substituting the finite perturbation expansions
of the dependent variable $f(t,x,y)$ around a formal perturbation
parameter $\varepsilon$
\begin{equation}
f(t,x,y)=1+\varepsilon f^{(1)}(t,x,y)+\varepsilon^2f^{(2)}(t,x,y)+\dots.\label{f}\end{equation}
  into the Hirota bilinear form
\eqref{hirotabilinearform}, we derive
\begin{equation}\label{perturbation}
\begin{aligned}
&P(D_1,D_2,D_3)\{f(t,x,y)\cdot f(t,x,y)\}\\
&=P(D_1,D_2,D_3)[\{1.1\}+\varepsilon
\{ 1\cdot f^{(1)}+f^{(1)}.1\}+\varepsilon^2
\{ 1\cdot f^{(2)}+f^{(2)}.1+f^{(1)}\cdot f^{(1)}\}
\\
&\quad +\varepsilon^3\{ 1\cdot f^{(3)}+f^{(3)}\cdot 1+f^{(1)}
\cdot f^{(2)}+f^{(2)}\cdot f^{(1)}\}\\
&\quad + \varepsilon^4\{ 1\cdot f^{(4)}+f^{(4)}\cdot1+f^{(1)}\cdot f^{(3)}
+f^{(3)}\cdot f^{(1)}+f^{(2)}\cdot f^{(2)}\}+\dots].
\end{aligned}
\end{equation}
The last step towards the method is to
analyze the conditions on the coefficients of $\varepsilon^i$,
for all $i\geq 0$ for multi-$q$-soliton solutions. From the
coefficient of the first term $\varepsilon^0$, we have
\begin{equation*}
P(D_1,D_2,D_3)\{1\cdot1\}=\lambda_1+\lambda_2+\lambda_3.
\end{equation*}

\begin{theorem}[\cite{hirotatheorem}]\label{constraint}
Any equation in Hirota bilinear form $P(D_t,D_x,D_y)f\cdot f=0$,
satisfying the sufficient conditions
\begin{gather}
P(0,0,0)=0,\\
P(D_t,D_x,D_y)=P(-D_t,-D_x,-D_y),
\end{gather}  
 admits at least two-soliton solutions.
\end{theorem}

To satisfy the conditions of Theorem
\eqref{constraint}, hereafter we need to have the constraint
\begin{equation}
P(0,0,0)=\lambda_1+\lambda_2+\lambda_3=0.\label{assumption}
\end{equation}
The coefficient of $\varepsilon^1$ implies
\begin{equation}
\begin{aligned}
&P(D_1,D_2,D_3)\{ 1\cdot f^{(1)}+f^{(1)}.1\} \\
&= 2P(\partial_1,\partial_2,\partial_3)f^{(1)} \\
&= 2[\sum_{i=1}^3\frac{\lambda_i}{2}(\exp{(a_it\partial_t+b_ix\partial_x
 +c_iy\partial_y)}
\exp{(-a_it\partial_t-b_ix\partial_x-c_iy\partial_y)})
]f^{(1)}\\
&=0.
\end{aligned} \label{eq:9}
\end{equation}
In the literature, soliton solutions of both differential
\cite{hirota1} or difference \cite{Hirota10} type of equations
tend to be of the exponential form.
 However, $q$-difference equations have an exclusive nature.
 We remark that the $q$-difference equation \eqref{eq:9} does not
admit exponential type of solutions and its solution needs to include power
counterparts for  the $q$-discrete variables, which is indeed a
consequence of change of variables. Therefore, as all variables
are $q$-discrete, the $q$-difference equation
 \eqref{eq:9} admits a starting solution of the power form
\begin{equation}
f^{(1)}(t,x,y)=\eta t^{\alpha}x^{\beta}y^{\gamma},\label{eq:10}
\end{equation}
where $\alpha,\eta,\beta$ are arbitrary
 constants. The solutions of the form \eqref{eq:10} provide
 \emph{$q$-soliton} solutions. The notion of $q$-solitons are introduced
in \cite{burcu5} as follows.

\begin{definition} \rm
A solution possessing usual soliton behaviors
 and having  power counterparts for $q$-discrete variables are called as
$q$-soliton  solution.
\end{definition}

Substituting such solution \eqref{eq:10} in
 \eqref{eq:9}, we obtain the so-called dispersion relation which determines
the relation among the  parameters as
\begin{equation}
P(v)=\sum_{i=1}^3\frac{\lambda_i}{2}(q_i^{\alpha}p_i^{\beta}r_i^{\gamma}
+q_i^{-\alpha}p_i^{-\beta}r_i^{-\gamma})=0,\label{dispersion}
\end{equation}
  where we denote the vector  $v=(\alpha,\beta,\gamma)$.
The coefficient of $\varepsilon^2$ resulted from
\eqref{perturbation} can be written as
\begin{equation}
P(D_1,D_2,D_3)\{f^{(1)}\cdot f^{(1)}\}=
-2P(\partial_1,\partial_2,\partial_3)f^{(2)}.\label{eq:13}
\end{equation}
  Substituting $f^{(1)}$ given in \eqref{eq:10} on the
left hand side of \eqref{eq:13}, yields as
\begin{equation*}
P(D_1,D_2,D_3)\{f^{(1)}\cdot f^{(1)}\}=(\lambda_1+\lambda_2
+\lambda_3)\eta^2t^{2\alpha}x^{2\beta}y^{2\gamma},
\end{equation*}
  which vanishes by the ansatz \eqref{assumption}.
Therefore, for all $j\geq 2$, we assume that $f^{(j)}=0$. As a
generalization, for $i$-$q$-soliton solution, we assume
$f^{(j)}=0$ for all $j\geq {i+1}$. Setting $\varepsilon=1$, we
express the solution describing one-$q$-soliton as
\begin{equation}
f(t,x,y)=1+\eta t^{\alpha}x^{\beta}y^{\gamma}.
\end{equation}
To obtain two-$q$-soliton solutions, we start with the following solution of
\eqref{eq:9}
\begin{equation*}
f^{(1)}=\sum_{i=1}^2\eta_it^{\alpha_i}x^{\beta_i}y^{\gamma_i},
\end{equation*}
where $\eta_i,\alpha_i,\beta_i$' s are constants for all
$i=1,2$. By the constraint \eqref{assumption}, the coefficient of
$\varepsilon^0$ vanishes and the coefficient of $\varepsilon^1$
implies the dispersion relation
\begin{equation}
P(v_j)=P(\alpha_j,\beta_j,\gamma_j)
=\sum_{i=1}^3\frac{\lambda_i}{2}(q_i^{\alpha_j}p_i^{\beta_j}r_i^{\gamma_j}
+q_i^{-\alpha_j}p_i^{-\beta_j}r_i^{-\gamma_j})=0,\quad\forall
j=1,2.\label{dispersion2}
\end{equation}

Subsequently, from the coefficient of $\varepsilon^2$ we derive
\begin{equation}
\begin{aligned}
-P(\partial)f^{(2)}
&=\eta_1\eta_2t^{\alpha_1+\alpha_2}
x^{\beta_1+\beta_2}y^{\gamma_1+\gamma_2}\sum_{i=1}^3
\frac{\lambda_i}{2}[q_i^{\alpha_1-\alpha_2}p_i^{\beta_1-\beta_2}
r_i^{\gamma_1-\gamma_2}\\
&\quad +q_i^{\alpha_2-\alpha_1}p_i^{\beta_2-\beta_1}
r_i^{\gamma_2-\gamma_1}],
\end{aligned}\label{ekare}
\end{equation}
which implies that $f^{(2)}$ is of the form
\begin{equation}
f^{(2)}=A(1,2)\eta_1\eta_2t^{\alpha_1+\alpha_2}x^{\beta_1+\beta_2}
y^{\gamma_1+\gamma_2}.\label{f2}
\end{equation}
  Substituting such $f^{(2)}$, given in \eqref{f2} into
\eqref{ekare}, we find the phase shift among two-$q$-solitons as
\begin{equation*}
A(1,2)=-\frac{P(v_1-v_2)}
 {P(v_1+v_2)},
\end{equation*}
  where the vector notation stands for
\begin{equation*}P(v_1\pm v_2)=\sum_{i=1}^3\frac{\lambda_i}{2}
[q_i^{\alpha_1\pm\alpha_2}p_i^{\beta_1\pm\beta_2}r_i^{\gamma_1\pm\gamma_2}
+q_i^{-(\alpha_1\pm\alpha_2)}p_i^{-(\beta_1\pm\beta_2)}
r_i^{-(\gamma_1\pm\gamma_2)}].
\end{equation*}
Because of the fact that all higher order terms of $f^{(i)}$,
$i\geq 3$ vanishes and the dispersion relation \eqref{dispersion2}
holds, the coefficient of $\varepsilon^j=0$ for all $j\geq 3$.
Therefore, two-$q$-solitons are given by
\begin{equation}
f(t,x,y)=1+\eta_1t^{\alpha_1}x^{\beta_1}y^{\gamma_1}
+\eta_2t^{\alpha_2}x^{\beta_2}y^{\gamma_2}
+A(1,2)\eta_1\eta_2t^{\alpha_1+\alpha_2}x^{\beta_1+\beta_2}
y^{\gamma_1+\gamma_2}.\label{2ss}
\end{equation}
For three-$q$-soliton solutions, we begin with
\begin{equation*}
f^{(1)}=\sum_{i=1}^3\eta_it^{\alpha_i}x^{\beta_i}y^{\gamma_i},
\end{equation*}
where $\alpha_i,\eta_i, \beta_i$ are constants for $i=1,2,3$. The
coefficient of $\varepsilon^1$ enables to have a similar
dispersion relation
\begin{equation}
P(v_j)=\sum_{i=1}^3\frac{\lambda_i}{2}(q_i^{\alpha_j}
p_i^{\beta_j}r_i^{\gamma_j}+q_i^{-\alpha_j}p_i^{-\beta_j}r_i^{-\gamma_j})=0,
\quad \forall j=1,2,3 \label{dispersion3}.
\end{equation}
Further,  $f^{(2)}$ follows
\[
f^{(2)}=\sum_{i<j}^{(3)}A(i,j)\eta_i\eta_jt^{\alpha_i+\alpha_j}
x^{\beta_i+\beta_j}y^{\eta_i+\eta_j},
\]
 from the coefficient of $\varepsilon^2$
\begin{align*}
-P(\partial)f^{(2)}&=\sum_{i<j}^{(3)}\eta_i\eta_jt^{\alpha_i+\alpha_j}
 x^{\beta_i+\beta_j}y^{\gamma_i+\gamma_j}
\Big[\sum_{k=1}^3\frac{\lambda_k}{2}(q_k^{\alpha_i-\alpha_j}
p_k^{\beta_i-\beta_j}r_k^{\gamma_i-\gamma_j}
\\
&\quad +q_k^{\alpha_j-\alpha_i}p_k^{\beta_j-\beta_i}r_k^{\gamma_j-\gamma_i})\Big]
,\end{align*}
 where $\sum_{i<j}^{(3)}$ denotes the
summation over all elements such that $i<j$ and $i,j=1,2,3$. Here
one can derive the associated interaction terms
\begin{equation}
A(j,k)=-\frac{P(v_j-v_k)}{P(v_j+v_k)},\quad j<k,\; j,k=1,2,3,
\end{equation}
  among three-$q$-soliton solutions. The coefficient of
$\varepsilon^3$ yields $f^{(3)}$ as
\[
f^{(3)}=A(1,2,3)\eta_1\eta_2\eta_3t^{\alpha_1+\alpha_2+\alpha_3}
x^{\beta_1+\beta_2+\beta_3}y^{\gamma_1+\gamma_2+\gamma_3},
\]
  where
\begin{equation}
\begin{aligned}
A(1,2,3)&=-\Big(A(1,2)P(v_3-v_1-v_2)+A(1,3)P(v_2-v_1-v_3)\\
&\quad +A(2,3)P(v_1-v_2-v_3)\Big) \big/ P(v_1+v_2+v_3).
\end{aligned}\label{b.2}
\end{equation}
  If the coefficient of $\varepsilon^4$ is under
consideration, having $f^{(4)}=0$, we encounter another expression
for $A(1,2,3)$
\begin{equation}
A(1,2,3)
=A(1,2)A(1,3)A(2,3).\label{b2.2}
\end{equation}
  Both expressions \eqref{b.2} and \eqref{b2.2} for
$A(1,2,3)$ are equivalent provided that the \emph{three-soliton
solution condition}
\begin{equation}
\sum_{\sigma_i=\pm
1}P(\sum_{i=1}^3\sigma_iv_i)\prod_{i<j}^{(3)}P(\sigma_iv_i-\sigma_jv_j)=0,
 \quad i,j=1,2,3,\label{hirotacondition}
\end{equation}
is satisfied (see \cite{hirota9}). To sum up, to guarantee the
existence of three-soliton solutions, we end up with the condition
\eqref{hirotacondition} arising as a constraint on $P$. Thus, we
present the three-$q$-soliton solutions
\begin{equation}\label{3.2ss}
\begin{split}
f(x,t)
&=1+\sum_{i=1}^3\eta_i t^{\alpha_i}x^{\beta_i}y^{\gamma_i}
+\sum_{i<j}^{(3)}A(i,j)\eta_i\eta_j t^{\alpha_i+\alpha_j}x^{\beta_i+\beta_j}
y^{\gamma_i+\gamma_j}\\
&\quad + A(1,2)A(1,3)A(2,3)\eta_1\eta_2\eta_3
t^{\alpha_1+\alpha_2+\alpha_3}x^{\beta_1+\beta_2+\beta_3}
y^{\gamma_1+\gamma_2+\gamma_3},
\end{split}
\end{equation}
which are expressed in terms of a polynomial in power
functions.

\section{Special cases}

In the literature, the concept of a dependent variable
transformation that converts a given nonlinear partial
differential or difference equation into the bilinear form is one
of the significant tools. The transformed new variables, that are
expressed as Wronski or Casorati type of determinants, are said to
be $\tau$-\emph{functions} and they solve equations in bilinear
forms \cite{sato,satsuma1}. In \cite{satsuma2}, two-dimensional
$q$-discrete Toda lattice equation is presented and the $\tau$
functions are determined by means of Wronski determinant. Two
different $q$-deformations of KdV hierarchies are introduced in
\cite{frenkel} and in \cite{klr} on which the notion of
integrability is discussed via bi-Hamiltonian structures. In
\cite{adler-horozov}, it is stated that any KdV $\tau$-function
provides a $q$-KdV $\tau$-function. However, to our knowledge no
research has been addressed to $q$-discretization of sine-Gordon
equation. In this section, we intend to present other
$q$-discretized versions of Toda, KdV equations and introduce
$q$-difference sine-Gordon equation resulted by proper reductions
of parameters on the generalized $q$-difference soliton equation
\eqref{hirotabilinearform}. Such reductions on the three
$q$-soliton solutions, presented in Section 3, provide three
$q$-soliton solutions of the considered equations in one move.

  Before embarking to the details we explain the source of
$q$-discretization as follows.

 There are several approaches to $q$-discretize continuous
equations. The $q$-deformed equations can be constructed by the
frame of $q$-forward jump operator $E_q$ given by \eqref{q} or
they can be obtained by means of $q$-derivative operator
$\partial_q$ \cite{jackson},
\begin{equation}
\partial_{q,x} f(x)=\frac{f(qx)-f(x)}{qx-x}.
\end{equation}
Here $f$ is a $q$-differentiable function. However, classical Hirota method
fails to produce $q$-exponential type of multi-soliton solutions
for $q$-differential equations with respect to $\partial_q$. In
\cite{burcu5}, $q$-analogue of Hirota $D$-operator acting on
$q$-differentiable functions $f,g$
\begin{equation}
D_{q,x}^{m}\{f.g\}:=(\partial_{q,x}-\partial_{q, x^{'}})^m
f(x)g(x^{'})| _{x'=x}, \quad
m\in\mathbb{Z}^{+},\label{qdop}
\end{equation}
 was introduced and developed where $q$-deformed Hirota bilinear
forms were expressed in terms of such $q$-Hirota $D$-operator.
$q$-differential-$q$-difference version of Toda equation was
analyzed. We note that even though it is in $q$-deformed Hirota
bilinear form
\begin{align*}
&P(D_{q,t},D_x)\{f(x,t)\cdot f(x,t)\}\\
&=[D_{q,t}^2-(\exp(hxD_x)+\exp(-hxD_x)-2)]\{f(x,t)\cdot f(x,t)\}=0,
\end{align*}
  and satisfies the sufficient conditions
\begin{align*}
P(D_{q,t},D_x)=P(-D_{q,t},-D_x),\quad P(0,0)=0,
\end{align*}
it admits only solitary wave like a solution determined
by $q$-exponential function as $f^{(1)}(x,t)=\eta x^\alpha
e_{q}^{\beta t}$. Indeed, classical Hirota perturbation \eqref{f}
does not provide further $q$-exponential type of solutions. The
restrictive condition  $$wz=qzw,\quad z,w\in q^{\mathbb{Z}},$$ to
satisfy additive property of $q$-exponentials
\[
e_q^ze_q^w=e_q^{z+w},\quad z,w\in q^{\mathbb{Z}},
\]
causes to detract one of the summands in two-$q$-soliton solutions
and the solution reduces to one-$q$-soliton solution. Rather than
the choice $f^{(2)}=A(1,2)x^{\alpha_1+\alpha_2}e_q^{\beta_1
t+\beta_2 t}$, one can analyze the function
$f^{(2)}=A(1,2)x^{\alpha_1+\alpha_2}e_q^{\beta_1 t}e_q^{\beta_2
t}$. But in this case such a choice yields to have time dependency
on the interaction constant $A(1,2)$. Note that, all other
admissible choices for the function $f^{(2)}$ lead similar
dead-ends.

To be more general, we can analyze Hirota approach on an arbitrary
time scale. For this purpose, we introduce the notion of
$\Delta$-Hirota $D$-operator.

\begin{definition} \rm
Let $\mathbb{T}$ be an arbitrary time scale.
Let $\mathbb{T}^{\kappa}$ denote Hilger's above truncated
set consisting of $\mathbb{T}$ except for a possible left
scattered maximal point. Let $f(x),g(x): \mathbb{T}\mapsto
\mathbb{R}$ be arbitrary functions and $x\in \mathbb{T}^{\kappa}$.
We introduce $\Delta$-Hirota $D$-operator as
\begin{equation}
\mathcal{D}_x^m\{f.g\}:=(\Delta_x-\Delta_{x^{'}})^mf(x).g(x^{'})| _{x'=x},
\quad m\in\mathbb{Z}^{+}.\label{deltahirota}
\end{equation}
\end{definition}

  Here delta derivative of $f$ is defined in
\cite{bohner}, to act as
\begin{equation}\label{del}
\Delta_x f(x) = \lim_{s\to x} \frac{f(\sigma (x))-f(s)}{\sigma
(x)-s}, \quad x\in \mathbb{T}^{\kappa},  s\in\mathbb{T},
\end{equation}
where $\sigma:\mathbb{T}\mapsto\mathbb{T}$ is forward jump
operator given by
\begin{equation}
\sigma(x)=\inf \{ y \in {\mathbb T}: y> x\}, \quad x\in\mathbb{T}.\label{forward}
\end{equation}

We note that when $\mathbb{T}=\mathbb{R}$, then $\sigma(x)=x$,
$\Delta_x = \frac{d}{dx}$ and $\Delta$-Hirota $D$-operator
\eqref{deltahirota} turns out to be usual Hirota $D$-operator
\eqref{dop}.  If $\mathbb{T}=\mathbb{K}_q:=q^{\mathbb{Z}}\bigcup
\{0\}$, $q\neq 1$ then $\sigma(x)=qx$, $\Delta_x=\partial_{q,x} $
and $\Delta$-Hirota $D$-operator \eqref{deltahirota} reduces to
$q$-analogue of Hirota $D$-operator \eqref{qdop}.

One can make use of the $\Delta$-Hirota $D$-operator
\eqref{deltahirota} to write $\Delta$-differential equations
\cite{metin} on arbitrary time scales   in Hirota bilinear form.
By utilizing classical Hirota pertubation \eqref{f} developed on
such $\Delta$-Hirota $D$-operator, one can  produce multi soliton
solutions only when the \emph{graininess function}
$\mu:\mathbb{T}\to [0,\infty)$
\begin{equation}
\mu(x):=\sigma(x)-x, \quad x\in \mathbb{T},
\end{equation}
is constant (e.g. difference equations on $h\mathbb{Z}$ or
$q$-difference equations). To be more precise, in
order to derive multi soliton solutions on arbitrary time scales
with arbitrary non-constant graininess $\mu$, one should construct
a deformed perturbation different than the classical Hirota
perturbation \eqref{f}. However, such an approach will be no more
classical Hirota direct method.

To sum up, we conjecture that other than the unifying framework
proposed in the present article, it is not possible to create
another unifying approach via classical Hirota perturbation
\eqref{f} for integrable equations on discrete intervals with
non-constant graininess $\mu$, or on arbitrary such time scales.

In this work, we present $q$-discretization of equations that are
expressed by $q$-forward jump operator $E_q$ as it recovers the
continuous case. Our framework is based on the inverse procedure.
The core of the procedure is to introduce appropriate $q$-deformed
Hirota bilinear forms (resulted from proper reductions on
\eqref{hirotabilinearform}) in a way that they reduce to
continuous Hirota bilinear forms as $q\to 1$ (equivalently
in the small limit of $h$). In addition, from such $q$-deformed
Hirota bilinear forms, we develop not only $q$-analogues of
corresponding equations but also their multi-$q$-soliton
solutions.

\subsection{The $q$-difference-$q$-difference Toda equation}

  It is proposed in \cite{burcu5} that
$q$-difference-$q$-difference Toda equation admits the Hirota
bilinear form
\begin{equation}
\begin{aligned}
&[h^{-1}(\exp(h \tau D_\tau)+\exp(-h\tau D_\tau)-2)
 -(\exp(\bar{h} yD_y) \\
 &+\exp(-\bar{h}yD_y)-2)]\{f(\tau,y)\cdot f(\tau,y)\}=0.
\end{aligned} \label{e}
\end{equation}
  In the present work, we determine \eqref{e} by the
reductions
\begin{equation}
D_1=h\tau D_{\tau}, \quad  D_2=\bar{h}yD_y,\quad
D_3=0,\quad \lambda_1=2h^{-1},\quad \lambda_2=-2, \quad
\lambda_3=2-2h^{-1}\label{red}
\end{equation}
on the generalized  $q$-difference equation \eqref{hirotabilinearform}.
Note that $\lambda_1+\lambda_2+\lambda_3=0$ to satisfy the
condition \eqref{assumption}. In order to construct the standard
form of the $q$-difference-$q$-difference Toda equation, we need
to introduce the following operator.
\begin{definition} \cite{burcu5}
 The central $q$-difference operator $\Lambda_{\tau}^2$ acting on a
function $f(\tau)$, is defined as
\begin{equation}
\Lambda_{\tau}^2f(\tau)=f(q\tau)+f(\frac{\tau}{q})-2f(\tau), \quad
q\neq 1,\; \tau\in\mathbb{R}.\label{centdifq}
\end{equation}
\end{definition}

  Using the inverse procedure, the standard form of the
$q$-difference-$q$-difference Toda equation is proposed written by
the language of \eqref{centdifq}
\begin{equation}
\Lambda_\tau^2\log(1+V(\tau,y))=\Lambda_y^2\log (1+h V(\tau,y)),
\label{qdif-qdiftodaeqn}
\end{equation}
  which follows from the transformation
\begin{equation}
V(\tau,y):=h^{-1}\big[\frac{f(q\tau,y)f(\frac{\tau}{q},y)}{f^2(\tau,y)}-1\big]
=\frac{f(\tau,\bar{q}y)f(\tau,\frac{y}{\bar{q}})}{f^2(\tau,y)}-1,
\label{transformation}
\end{equation}
on the bilinear form \eqref{e}. Here $e^h=q$, $e^{\bar{h}}=\bar{q}$.

\begin{proposition}
One-$q$-soliton solution of $q$-difference-$q$-difference Toda equation
\eqref{qdif-qdiftodaeqn} is
\begin{equation}
V(\tau,y)=\frac{\eta\tau^\alpha y^{\beta}[\bar{q}^\beta+\bar{q}^{-\beta}-2]}{(1+\eta\tau^{\alpha}y^{\beta})^2},
\label{21ss}
\end{equation}
provided that the dispersion relation
\begin{equation}
h^{-1}(q^\alpha+q^{-\alpha}-2)=\bar{q}^\beta
+\bar{q}^{-\beta}-2\label{dispersion21},
\end{equation}
is satisfied.
\end{proposition}

\begin{proof}
The identifications
\begin{equation}
t=\tau,\quad \gamma=0,\quad a_1=h,\quad b_2=\bar{h},\quad
  a_2=a_3=b_1=b_3=c_1=c_2=c_3=0\label{identification}
\end{equation}
resulting from the reductions \eqref{red} allows to obtain
$f^{(1)}=\eta\tau^{\alpha}y^{\beta}$. One-$q$-soliton solution
follows from straightforward calculation of using
$f=1+\eta\tau^{\alpha}y^{\beta}$  on the dependent variable
transformation \eqref{transformation}.
The identifications \eqref{identification} on the dispersion relation
\eqref{dispersion} lead to derive the dispersion relation
\eqref{dispersion21} which determines the relation among the
parameters in \eqref{21ss}.
\end{proof}

  If $\tau$, $y\in q^{\mathbb{Z}}$, namely $\tau=q^n$ and
$y=\bar{q}^m$, $n,m\in\mathbb{Z}$, then
$q$-difference-$q$-difference Toda equation
\eqref{qdif-qdiftodaeqn} can be rewritten
\begin{equation}
\begin{aligned}
&\frac{(1+V(q^{n+1},(\bar{q})^{m}))(1+V(q^{n-1},(\bar{q})^{m}))}
{(1+V(q^n,(\bar{q})^{m}))^2}\\
&= \frac{(1+hV(q^n,(\bar{q})^{m+1}))(1+hV(q^n,(\bar{q})^{m-1}))}
{(1+hV(q^n,(\bar{q})^{m}))^2},
\end{aligned} \label{newform2}
\end{equation}
whose one-$q$ soliton solution is explicitly deduced as
\begin{equation}
 V(\tau,y)=\frac{\eta q^{n\alpha} (\bar{q})^{m\beta}
[(\bar{q})^\beta+(\bar{q})^{-\beta}-2]}{(1+\eta
 q^{n\alpha}(\bar{q})^{m\beta})^2}.
 \label{21sss}
\end{equation}
Subsequently, it is possible to rewrite two-$q$-soliton
\eqref{2ss} and three-$q$-soliton solutions \eqref{3.2ss} using
the above reductions and identifications. Therefore, by the frame
of the reductions \eqref{red} on \eqref{hirotabilinearform} we
recover the results found in the Ref.\cite{burcu5} for
$q$-difference-$q$-difference Toda equation. Moreover, it is
observed that the plotted  waves associated to two-$q$-soliton and
three-$q$-soliton solutions of $q$-difference-$q$-difference Toda
equation obey the classical soliton behaviors. Additionally, since
solutions are presented by means of power functions, the lengths
of waves increase as space variable increases (We refer
\cite{burcu5} for the graphs). Therefore the solutions presented
in Section 3, appear to satisfy $q$-soliton conditions.

\begin{proposition}\label{prop}
Hirota bilinear form of the $q$-difference-$q$-difference Toda
equation \eqref{e} reduces to
Hirota bilinear form of the differential-$q$-difference Toda
equation
\begin{equation}
[{D_t}^2-(\exp(\bar{h}yD_y)+\exp(-\bar{h}yD_y)-2)]\{f(t,y)\cdot f(t,y)\}=0,
\label{d}
\end{equation}
in the small limit of $h$. 
\end{proposition}

  Indeed, if we replace $h$ by $h^2$ and let $\tau=\exp(ht)$, then the 
expression on \eqref{e}
\begin{equation}
4h^{-2}\sinh^2(\frac{h^2\tau D_{\tau}}{2})= h^{-2}(\exp(hD_t)+\exp(-hD_t)-2),
\end{equation}
tends to $D_t^2$ as $h\to 0$.   To be more
precise, Hirota bilinear form \eqref{e}, is a generalization of
Hirota bilinear form \eqref{d}, from which we establish the
standard form of the $q$-discretized Toda equation
\eqref{qdif-qdiftodaeqn}.

\subsection{The $q$-difference-$q$-difference KdV equation}
  We propose the $q$-difference-$q$-difference KdV
equation in Hirota bilinear form
\begin{equation}
\sinh(\frac{h^2\tau D_{\tau}+\bar{h}^2yD_y}{2})[h^{-1}\sinh(h^2 \tau
D_{\tau})+2\sinh(\bar{h}^2yD_y)]\{f(\tau,y)\cdot f(\tau,y)\}=0,
\label{qdifqdifkdv}
\end{equation}
that can be identified by the below reductions
\begin{equation}
D_1 = \frac{1}{2}(3h^2\tau D_{\tau}+\bar{h}^2yD_y) ,\quad
D_2 = \frac{1}{2}(h^2\tau D_{\tau}+3\bar{h}^2yD_y),
D_3 = \frac{1}{2}(h^2\tau D_{\tau}-\bar{h}^2yD_y),  \label{reductions2}
\end{equation}
\begin{equation} \lambda_1=1,\quad\lambda_2=2h,\quad
\lambda_3=-1-2h,\label{reductions2a}
\end{equation}
on \eqref{hirotabilinearform}.

\begin{proposition} \label{prop2}
Hirota bilinear form of $q$-difference-$q$-difference KdV equation
\eqref{qdifqdifkdv} reduces to Hirota bilinear form of the
continuous KdV equation \cite{hirota1}
\begin{equation}
[D_x(D_t+D_x^3)]\{f(t,x)\cdot f(t,x)\}=0,\label{kdv}
\end{equation}
as $h$, $\bar{h}$ tends to zero.
\end{proposition}

\begin{proof} Setting  $\tau=\exp(ht)$, Hirota
bilinear form of $q$-difference-$q$-difference KdV equation
\eqref{qdifqdifkdv} turns out to be
\begin{equation}
\sinh(\frac{h D_{t}+\bar{h}^2yD_y}{2})[h^{-1}\sinh(h
D_{t})+2\sinh(\bar{h}^2yD_y)]\{f\cdot f\}=0.\label{2qdifqdifkdv}
\end{equation}
  One can verify that as $h\to 0$, the equation
\eqref{2qdifqdifkdv} reduces to
\begin{equation}
\sinh(\frac{\bar{h}^2yD_y}{2})[D_t+2\sinh(\bar{h}^2yD_y)]\{f(t,y)\cdot f(t,y)\}=0.
\label{difqdifkdv}
\end{equation}
  By a similar fashion, we set $y=\exp(\bar{h}x)$ which
implies \eqref{difqdifkdv} as
\begin{equation}
\sinh(\frac{\bar{h}D_x}{2})[D_t+2\sinh(\bar{h}D_x)]\{f(t,x)\cdot f(t,x)\}=0.
\label{difqdifkdv1}
\end{equation}
As a final step, we fix  $D_t=\frac{\bar{h}^3D_t}{3}-2\bar{h}D_x$
in \eqref{difqdifkdv1} and we derive
\begin{equation}
\sinh(\frac{\bar{h}D_x}{2})[\frac{\bar{h}^3D_t}{3}-2\bar{h}D_x
+2\sinh(\bar{h}D_x)]\{f\cdot f\}=0.\label{2difqdifkdv}
\end{equation}
  Dividing \eqref{2difqdifkdv} with $\bar{h}^4$ and taking
the small limit of $\bar{h}$, we end up with the continuous KdV
equation in Hirota bilinear form \eqref{kdv}.
\end{proof}

To be more precise, $q$-deformed Hirota bilinear form
\eqref{qdifqdifkdv} is a generalization of Hirota bilinear form
\eqref{kdv} of usual KdV equation.

\subsection{The standard form of $q$-difference-$q$-difference KdV equation}
  In this section, we intend to introduce the standard
form of the $q$-difference-$q$-difference KdV equation. Before
embarking to the details we list some useful identities for
hyperbolic functions.
\smallskip



\noindent\textbf{Properties}: \cite{hirota9}
For continuously differentiable functions $f$, $g$, the following identities hold
\begin{gather}
\begin{aligned}
&\sinh\frac{1}{2}(D_2-D_3)\{\sinh\frac{1}{2}
(D_2+D_3)\cdot 2\sinh D_1f\cdot f\}\\
&.\{\cosh(\frac{1}{2}(D_2+D_3)-D_1)f\cdot f\} \\
&=\sinh D_1\{\cosh D_2 f\cdot f\}\cdot\{\cosh D_3 f\cdot f\},
\end{aligned} \label{a}
\\
\cosh D_1\{\cosh D_2f\cdot f\}\cdot\{\cosh D_2f\cdot f\}=\cosh D_2\{\cosh
D_1f\cdot f\}.\{\cosh D_1f\cdot f\},\label{b}
\\
\exp(\alpha \partial_1)(\frac{f}{g})=exp(\alpha
D_1)\{f\}.\{\frac{g}{\cosh(\alpha D_1) g.g}\},\label{c}
\end{gather}
where the operators $D_i$ are of the form
\eqref{operator} and $\alpha$ is an arbitrary constant.

\begin{definition} \rm
The $q$-difference operator $\delta_{\tau}$, operating on an arbitrary
function  $u(\tau)$, is defined as follows
\begin{equation}
\delta_{\tau}u(\tau):=u(q\tau)-u(\frac{\tau}{q}),\quad
\tau\in \mathbb{R},\quad q\neq 1 \label{Delta}.
\end{equation}
\end{definition}

We stress that the operators $\Lambda^2$, defined by \eqref{centdifq}
and  $\delta^2$ and are equivalent.

\begin{proposition} \label{prop4.7}
The standard form of the
$q$-difference-$q$-difference KdV equation is given as
\begin{equation}
\delta_{\tau}(\frac{1}{V(\tau,y)})
=-2h^{1/2}\delta_y V(\tau,y),\label{kdv2}
\end{equation}
provided that the dependent variable transformation
\begin{equation}
V(\tau,y):=-\frac{f(\tau,\bar{q}y)f(\tau,\frac{y}{\bar{q}})}{f(q\tau,
y)f(\frac{\tau}{q},y)},
\label{transformation2}
\end{equation}
  holds, where $e^h=q$ and $e^{\bar{h}}=\bar{q}$.
\end{proposition}

\begin{proof}
 We first adopt an equivalent version of Hirota bilinear form
\eqref{qdifqdifkdv} of $q$-difference-$q$-difference KdV equation
\begin{equation}
\begin{aligned}
&\sinh(\frac{h\tau D_{\tau}-\bar{h}yD_y}{2})\{\sinh(\frac{h\tau
D_{\tau}+\bar{h}yD_y}{2})[h^{-1/2}\sinh(h \tau
D_{\tau})\\
&+2\sinh(\bar{h}yD_y)]f\cdot f\} \cdot \{\cosh(\frac{h\tau
D_{\tau}-\bar{h}yD_y}{2})f\cdot f\}=0,
\end{aligned}\label{qdifqdifkdv*}
\end{equation}
where for convenience we interchange $h^2$ and $\bar{h}^2$ by $h$
and $\bar{h}$, respectively. Using the identity \eqref{a}, the
bilinear form \eqref{qdifqdifkdv*} can be rewritten in the form
\begin{equation}
\begin{aligned}
&h^{-1/2}\sinh(h\tau D_{\tau})\{\cosh(h\tau
D_{\tau})f\cdot f\}.\{\cosh(\bar{h}yD_y)f\cdot f\}\\
&= -2\sinh(\bar{h}yD_y)\{\cosh(h\tau
D_{\tau})f\cdot f\}.\{\cosh(\bar{h}yD_y)f\cdot f\}.
\end{aligned} \label{qdifqdifkdv**}
\end{equation}

  Upon dividing the form \eqref{qdifqdifkdv**} by
\eqref{b} equipped with $D_1=h\tau D_{\tau}$, $D_2=\bar{h}yD_y$,
we obtain
\begin{align*}
&\frac{h^{-1/2}\sinh(h\tau D_{\tau})(\cosh (h\tau
D_{\tau})f\cdot f).({\cosh(\bar{h}yD_y)f\cdot f})}{ \cosh(h\tau
D_{\tau})(\cosh(\bar{h}yD_y)f\cdot f)\cdot(\cosh(\bar{h}yD_y)f\cdot f)}\\
&=\frac{-2\sinh(\bar{h}yD_y)(\cosh
(h\tau D_{\tau})f\cdot f).({\cosh(\bar{h}yD_y)f\cdot f})}{\cosh(\bar{h}y
D_y)(\cosh(h\tau D_{\tau}f\cdot f))\cdot (\cosh(h\tau
D_{\tau})f\cdot f)},
\end{align*}
 which takes the form
\begin{equation}
 h^{-1/2}\sinh(h\tau \partial_{\tau})(\frac{\cosh(h\tau
D_{\tau})f\cdot f}{\cosh(\bar{h}yD_y)f\cdot f})
=-2\sinh(\bar{h}y\partial_y)(\frac{\cosh(\bar{h}yD_y)f\cdot f}{\cosh(h\tau
D_{\tau})f\cdot f}),\label{qdifqdifkdv****}
\end{equation}
Here we make use of the identity \eqref{c}. Armed with the $q$-exponential
identity \eqref{qexponentialidentity}, by the frame of $e^h=q$ and
$e^{\bar{h}}=\bar{q}$, then \eqref{qdifqdifkdv****} boils down to
\begin{equation}
\begin{aligned}
&[\frac{f(q^2\tau,y)f(\tau,y)}{f(q\tau,\bar{q}y)\cdot f(q\tau,\frac{y}{\bar{q}})}
-\frac{f(\tau,y)f(\frac{\tau}{q^2},y)}{f(\frac{\tau}{q},\bar{q}y)
\cdot f(\frac{\tau}{q},\frac{y}{\bar{q}})}] \\
&=2h^{1/2}[\frac{f(\tau,y)f(\tau,\frac{y}{\bar{q}^2})}{f(q\tau,\frac{y}{\bar{q}})
 \cdot f(\frac{\tau}{q},\frac{y}{\bar{q}})}
-\frac{f(\tau,\bar{q}^2y)f(\tau,y)}{f(q\tau,\bar{q}y)
 \cdot f(\frac{\tau}{q},\bar{q}y)}].
\end{aligned}\label{s}
\end{equation}
  Utilizing the dependent variable transformation
\eqref{transformation2}, we rewrite the equation \eqref{s} as
\eqref{kdv2} in the language of $q$-difference operator
\eqref{Delta}.
\end{proof}


\begin{proposition} \label{prop4.8}
The $q$-difference-$q$-difference KdV equation \eqref{kdv2} admits
one-$q$-soliton solution
\begin{equation}
V(\tau,y)=-\frac{[1+\eta\tau^\alpha (\bar{q}y)^{\beta}
+\eta\tau^\alpha (\bar{q})^{-\beta}y^{\beta}+\eta^2\tau^{2\alpha}y^{2\beta}]}
{(1+\eta (q\tau)^{\alpha}y^{\beta})(1+\eta
\tau^{\alpha}(q)^{-\alpha}y^{\beta})},\label{twoss}
\end{equation}
where the dispersion relation is
\begin{equation}
(\bar{q})^{\frac{\beta}{2}}[q^{\frac{3\alpha}{2}}-q^{-\frac{\alpha}{2}}]+
(\bar{q})^{\frac{-\beta}{2}}[q^{\frac{-3\alpha}{2}}-q^{\frac{\alpha}{2}}]+2h\{
q^{\frac{\alpha}{2}}[(\bar{q})^{\frac{3\beta}{2}}-(\bar{q})^{-\frac{\beta}{2}}]
+ q^{\frac{-\alpha}{2}}[(\bar{q})^{\frac{-3\beta}{2}}
 -(\bar{q})^{\frac{\beta}{2}}]\}=0.\label{dispersion31}
\end{equation}
\end{proposition}

\begin{proof} Using the identifications
\begin{gather}
t=\tau,\quad\gamma=0,\quad a_1=\frac{3h}{2},\quad
b_1=\frac{\bar{h}}{2},\label{identification2}\\
a_2=a_3=\frac{h}{2},\quad b_2=\frac{3\bar{h}}{2},
\quad b_3=-\frac{\bar{h}}{2},\quad
c_1=c_2=c_3=0,\label{identification2a}
\end{gather}
resulting from the reductions \eqref{reductions2},
\eqref{reductions2a} we obtain
$f^{(1)}=\eta\tau^{\alpha}y^{\beta}$. One can find one-$q$-soliton
solution \eqref{twoss} using $f=1+\eta\tau^{\alpha}y^{\beta}$ on
the transformation \eqref{transformation2}. Similarly, using
\eqref{identification2}, \eqref{identification2a} on the
dispersion relation \eqref{dispersion} yields the dispersion
relation \eqref{dispersion31}.
\end{proof}

  If $\tau$, $y\in q^{\mathbb{Z}}$, i.e., $\tau=q^n$ and
$y=\bar{q}^m$, $n,m\in\mathbb{Z}$, then one-$q$-soliton solution
\eqref{twoss} turns out to be
\begin{equation}
V=-\frac{[1+\eta
q^{n\alpha} (\bar{q})^{\beta(m+1)}+\eta q^{n\alpha}
(\bar{q})^{\beta(m-1)}+\eta^2\tau^{2n\alpha}(\bar{q})^{2m\beta}]}
{(1+\eta (q)^{\alpha(n+1)}(\bar{q})^{m\beta})(1+\eta
q^{\alpha(n-1)}(\bar{q})^{m\beta})}.
\end{equation}
  Furthermore, one can explicitly present  two and
three-$q$-soliton solutions using the reductions on \eqref{2ss}
and \eqref{3.2ss}, respectively.

\subsection{The $q$-difference sine-Gordon Equation}
  We propose Hirota bilinear form of $q$-difference
sine-Gordon equation
\begin{equation}
\begin{aligned}
&[2\sinh(\bar{h}^2yD_y)\sinh(h^2\tau
D_{\tau}) h\bar{h}\cosh(h^2 \tau
D_{\tau}+\bar{h}^2yD_y+kzDz) \\
&-h\bar{h}\cosh(h^2 \tau
D_{\tau}-\bar{h}^2yD_y)]\{f(\tau,y,z)\cdot f(\tau,y,z)\}=0,
\end{aligned} \label{qdifsinegordon1}
\end{equation}
which follows from the reductions below on the equation
\eqref{hirotabilinearform},
\begin{gather} \label{reductions3}
\begin{gathered}
D_1 = h^2\tau D_{\tau}+\bar{h}^2yD_y, \\
D_2 = h^2\tau D_{\tau}+\bar{h}^2yD_y+kzDz, \\
D_3 = h^2\tau D_{\tau}-\bar{h}^2yD_y,
\end{gathered} \\
\lambda_1=1,\quad\lambda_2=h\bar{h},\quad
\lambda_3=-1-h\bar{h}.\label{reductions3a}
\end{gather}
  To present the source of the limits which reveals Hirota
Bilinear form of  classical sine-Gordon equation
\cite{hirota-satsuma}, we need to rewrite \eqref{qdifsinegordon1}
in a proper decomposition. For this purpose, we adopt the
\emph{periodicity} definition on $q$-numbers.

\begin{definition}  \label{def4.10}\rm
A function $f(x)$ is said to be $q^{n}$-periodic if
\begin{equation}
f(q^nx)=f(x),  \quad q>1, \; n\in\mathbb{Z},\;  x\in \mathbb{K}_q.
\end{equation}
\end{definition}

\begin{proposition} \label{prop3}
Hirota bilinear form of $q$-difference sine-Gordon equation
\eqref{qdifsinegordon1} reduces to Hirota bilinear form of the
continuous sine-Gordon equation \cite{hirota-satsuma}
\begin{gather}
D_xD_t\{\bar{g}.\bar{f}\}=\bar{g}.\bar{f},\label{sinegordonhirotabilinearform1}\\
D_xD_t\{\bar{f}.\bar{f}-\bar{g}.\bar{g}\}=0,\label{sinegordonhirotabilinearform2}
\end{gather}
in the small limit of $h$ and $\bar{h}$.
\end{proposition}

\begin{proof}
  Assume that $f$ in \eqref{qdifsinegordon1} is
${\bar{q}}^2$-periodic function equipped with
\begin{equation*}
f(\bar{q}z)=f(\frac{z}{\bar{q}}),\end{equation*}
then the discussion includes the function $f$ and its $q$-shifted version,
 say $\tilde{f}$, i.e.
$e^{kz\partial_z}f(z)=f(\bar{q}z):=\tilde{f}(z)$,
provided that $e^k=\bar{q}$. Such periodicity condition on
\eqref{qdifsinegordon1} implies that
\begin{align*}
&[(h\bar{h})^{-1}\sinh(\bar{h}^2yD_y)\sinh(h^2\tau
D_{\tau})-\frac{1}{2}\cosh(h^2 \tau
D_{\tau}-\bar{h}^2yD_y)]\{f(\tau,y,z)\cdot f(\tau,y,z)\} \\
&=-\frac{1}{2}\cosh(h^2 \tau
D_{\tau}+\bar{h}^2yD_y)\{\tilde{f}(\tau,y,z).\tilde{f}(\tau,y,z)\},
\end{align*}
   provided that $h\bar{h}\neq 0$. Equivalently one can encounter
the bilinear form
\begin{equation}
\begin{aligned}
&[(h\bar{h})^{-1}\sinh(\bar{h}^2yD_y)\sinh(h^2\tau
D_{\tau})]\{f\cdot f\} \\
&=\frac{1}{2}\cosh(\bar{h}^2yD_y)\cosh(h^2 \tau
D_{\tau})]\{f\cdot f-\tilde{f}.\tilde{f}\},
\end{aligned}\label{qdifsinegordon2}
\end{equation}
whose continuous limit is intimately related to classical
sine-Gordon equation. Setting
\begin{equation}
f=\bar{f}+i\bar{g} \quad
\tilde{f}=\bar{f}-i\bar{g},\label{g}
\end{equation}
on \eqref{qdifsinegordon2}, we have
\begin{equation}
\begin{aligned}
&[(h\bar{h})^{-1}\sinh(\bar{h}^2yD_y)\sinh(h^2\tau
D_{\tau})]\{\bar{g}\cdot\bar{f}\}\\
&=\cosh(\bar{h}^2yD_y)\cosh(h^2 \tau
D_{\tau})]\{\bar{g}\cdot\bar{f}\}
\end{aligned}\label{qdifsinegordon3}
\end{equation}
and
\begin{equation}
\sinh(\bar{h}^2yD_y)\sinh(h^2\tau
D_{\tau})\{\bar{f}.\bar{f}-\bar{g}.\bar{g}\}=0.\label{qdifsinegordon4}
\end{equation}
Let $y=e^{\bar{h}x}$ and $\tau=e^{ht}$, then Hirota bilinear form
\eqref{qdifsinegordon3}, \eqref{qdifsinegordon4} turns out to be
Hirota bilinear form of sine-Gordon equation
\eqref{sinegordonhirotabilinearform1},
\eqref{sinegordonhirotabilinearform2} respectively, in the limit
$h,\bar{h}\to 0$.
\end{proof}

\subsection{Standard form of $q$-difference sine-Gordon equation}

To establish the ordinary form of
$q$-difference sine-Gordon equation, we make use of its decomposed
bilinear form \eqref{qdifsinegordon3}, \eqref{qdifsinegordon4} on
which for the sake of convenience we interchange $h^2$ and
$\bar{h}^2 $ by $h$ and $\bar{h}$, respectively as
\begin{gather}
[(h\bar{h})^{-1/2}\sinh(\bar{h}yD_y)\sinh(h\tau
D_{\tau})]\{\bar{g}.\bar{f}\}=\cosh(\bar{h}yD_y)\cosh(h \tau
D_{\tau})]\{\bar{g}.\bar{f}\},\label{qdifsinegordon5}
\\
\sinh(\bar{h}yD_y)\sinh(h\tau
D_{\tau})\{\bar{f}.\bar{f}-\bar{g}.\bar{g}\}=0.\label{qdifsinegordon6}
\end{gather}
  The restrictions \eqref{g} imply to assume
$\bar{f}:=\exp(\rho).\cos(\phi)$,
$\bar{g}:=\exp(\rho).\sin(\phi)$. Then we rewrite
\eqref{qdifsinegordon5} as
\begin{align*}
&(1-(h\bar{h})^{1/2})\exp(\rho(q\tau,py)+\rho(\frac{\tau}{q},\frac{y}{p})).
\sin(\phi(q\tau,py)+\phi(\frac{\tau}{q},\frac{y}{p})) \\
&=(1+(h\bar{h})^{1/2})\exp(\rho(q\tau,\frac{y}{p})+\rho(\frac{\tau}{q},py)).
\sin(\phi(q\tau,\frac{y}{p})+\phi(\frac{\tau}{q},py)),
\end{align*}
and \eqref{qdifsinegordon6} as
\begin{align*}
&\exp(\rho(q\tau,py)+\rho(\frac{\tau}{q},\frac{y}{p})).
\cos(\phi(q\tau,py)+\phi(\frac{\tau}{q},\frac{y}{p})) \\
&= \exp(\rho(q\tau,\frac{y}{p})+\rho(\frac{\tau}{q},py)).
\cos(\phi(q\tau,\frac{y}{p})+\phi(\frac{\tau}{q},py)),
\end{align*}
respectively by the frame of $e^h=q$ and
$e^{\bar{h}}=p$. Solving the resulting equation for $\phi$, we
derive
\begin{equation}
\begin{aligned}
&\sin(\phi(q\tau,py) + \phi(\frac{\tau}{q},\frac{y}{p})
-\phi(q\tau,\frac{y}{p})-\phi(\frac{\tau}{q},py)) \\
&=(h\bar{h})^{1/2}\sin(\phi(q\tau,py)+\phi(\frac{\tau}{q},\frac{y}{p})
+\phi(q\tau,\frac{y}{p})+\phi(\frac{\tau}{q},py)).
\end{aligned}\label{eqn}
\end{equation}

\begin{definition} \label{4.11} \rm
 The $q$-sum operator $\Gamma_{\tau}$, operating on any function
$u(\tau)$ is defined as
\begin{equation}
 \Gamma_{\tau}u(\tau):=u(q\tau)+u(\frac{\tau}{q}),\quad
\tau\in \mathbb{R}, \; q\neq 1.\label{gamma}
\end{equation}
\end{definition}

  We introduce the standard form of the $q$-difference
sine-Gordon equation as
\begin{equation}
\sin[\delta_y\delta_{\tau}\phi(\tau,y)]=(h\bar{h})^{1/2}
\sin[\Gamma_y\Gamma_{\tau}\phi(\tau,y)],\label{eqn1}
\end{equation}
by rewriting \eqref{eqn} in the frame of  $q$-sum operator \eqref{gamma}
and $q$-difference operator \eqref{Delta}.

\begin{proposition}
The $q$-difference sine-Gordon equation \eqref{eqn1} admits one-$q$-soliton
solution
\begin{equation}
\phi=4\tan^{-1}(\eta \tau^{\alpha}y^{\beta}z^{\gamma}),\label{thress}
\end{equation}
provided that the dispersion relation
\begin{equation}
(q^{\alpha}-q^{-\alpha})(p^{\beta}-p^{-\beta})+
h\bar{h}[q^{\alpha}(p^{\beta}(\bar{q})^{\gamma}-p^{-\beta})
+q^{-\alpha}(p^{-\beta}(\bar{q})^{-\gamma}-p^{\beta})]=0,\label{dispersion41}
\end{equation}
is satisfied.
\end{proposition}

\begin{proof} The identifications
\begin{gather}
t=\tau,\quad a_1=a_2=a_3=h,\quad b_1=b_2=\bar{h},\label{identification3}\\
c_2=k \quad b_3=-\bar{h},\quad c_1=c_3=0,\label{identification3a}
\end{gather}
resulting from the reductions \eqref{reductions3},
\eqref{reductions3a} leads the starting solution of
$f^{(1)}=\eta\tau^{\alpha}y^{\beta}z^{\gamma}$.
If we expand $\bar{f}$ and $\bar{g}$ in series of
$\varepsilon$,
\begin{gather*}
\bar{f}=1+\varepsilon^2f^{(2)}+\varepsilon^4f^{(4)}+\dots\\
\bar{g}=\varepsilon f^{(1)}+\varepsilon^3f^{(3)}+\dots
\end{gather*}
it is possible to write $\bar{g}/\bar{f}=f^{(1)}$, equipped with
$\varepsilon=1$. To be more intense, we derive one-$q$-soliton
solution
\begin{equation*}
\phi=4\arctan(\frac{\bar{g}}{\bar{f}})
=4\arctan(\eta\tau^{\alpha}y^{\beta}z^{\gamma}).
\end{equation*}
  Further, the dispersion relation \eqref{dispersion41}
results from the reductions \eqref{identification3},
\eqref{identification3a} on the dispersion relation
\eqref{dispersion}.
\end{proof}

In the same way, it is possible to rewrite two-$q$-soliton \eqref{2ss} and
three-$q$-soliton solutions \eqref{3.2ss} using the above
reductions and identifications.

\subsection*{Conclusion}
In the present article, we have introduced a general unifying
framework for integrable $q$-discrete equations and their
multi-soliton solutions. We presented a generalized $q$-difference
equation, which reproduces $q$-discretized soliton equations such
as Toda, KdV and sine-Gordon equations by proper transformations.
We showed that Hirota direct method  produces three-$q$-soliton
solutions of this generalized $q$-difference  equation. However,
the classical method fails to construct multi-soliton solutions
for $\Delta$-differential equations on arbitrary time scales with
non-constant graininess. We plan to devote our next investigation
to concerning Wronskian technique to construct soliton and other
special (\emph{e.g.} complexiton) solutions of bilinear equations
on arbitrary time scales. Such a work is in progress and will be
communicated in a separate paper.

\subsection*{Acknowledgements}
The authors are grateful to Prof. Maciej B\l{a}szak for useful
discussions. This work is partially  supported by the Scientific
and Technological Research Council of Turkey (T{\" U}B\.{I}TAK).

\begin{thebibliography}{00}

\bibitem{adler-horozov} M. Adler, E. Horozov, P. Moerbeke;
\emph{The solution to the $q$-KdV equation}, Phys. Lett. A. 242 (1998), 139-151.

\bibitem{bgss} M. B\l{a}szak, M. G\"urses, B.  Silindir, B. M. Szablikowski;
\emph{Integrable discrete systems on $\mathbb{R}$ and related
 dispersionless systems}, J. Math. Phys. 49 (2008), 1-17.

\bibitem{bohner} M. Bohner, A. Peterson;
\emph{Dynamic Equations on time scale: An introduction with applications},
 Birkhauser, Boston, 2001.

\bibitem{frenkel} E. Frenkel;
\emph{Deformations of the KdV hierarchy and related
soliton equations}, Int. Math. Res. Notices. 2 (1996), 55-76.

\bibitem{Gardner}  C. S. Gardner, J. M. Greene, M. D. Kruskal, R. M. Miura;
\emph{Method for solving the Korteweg-de Vries equation}, Phys.
Rev. Lett. 19 (1967), 1095-1097.

\bibitem{grammaticos} B. Grammaticos, A. Ramani, J. Hietarinta;
\emph{A search for integrable bilinear equations:
The Painlev\'{e} approach},  J. Math. Phys. 31 (1990), 2572-2578.

\bibitem{metin} M. G\"urses M, S.G. Guseinov, B. Silindir;
\emph{Integrable equations on time scales}, J. Math. Phys. 46 (2005), 1-22.

\bibitem{gurses} M.  G\"urses, A. Pekcan;
\emph{2+1 KdV(N) equations}, J. Math. Phys. 52 (2011), 1-9.

\bibitem{hierentina0} J. Hietarinta;
\emph{A search of bilinear equations passing Hirota's three-soliton condition:
I. KdV-type bilinear equations}, J. Math. Phys. 28 (1987), 1732-1742.

\bibitem{hierentina1} J. Hietarinta;
\emph{Searching for integrable PDE's by testing Hirota's
three-soliton condition}. In: Watt SM, editor. ISSAC' 91
Proceedings of the 1991 International Symposium on Symbolic and
Algebraic Computation, ACM, Bonn, New York, 1991, pp. 295-300.

\bibitem{hierentina2} J. Hietarinta, D. Zhang; \emph{Hirota's method and
the search for integrable partial difference equations. 1.
Equations on a $3\times 3$ stencil}, J. Diff. Eq. and Appl. 19
(2013), 1292-1316.

\bibitem{hirota1} R. Hirota;
\emph{Exact solution of the Korteweg-de Vries equation for
multiple collisions of solitons}, Phys. Rev. Lett. 27 (1971)
1192-1194.

\bibitem{hirota-satsuma} R. Hirota;
 \emph{Exact solution of the sine-Gordon equation for multiple collisions of
solitons}, J. Phys. Soc. Japan, 33 (1972), 1459-1463.

\bibitem{hirotatheorem} R. Hirota;
\emph{Direct methods in soliton theory}. In:
Bullough RK, Caudrey PJ, editors. Topics in Current
Physics Solitons,  Springer-Verlag, New York, 1980, pp. 157-176.

\bibitem{Hirota10} R. Hirota;
\emph{Discrete analogue of a Generalized Toda equation},
 J. Phys. Soc. Japan, 50 (1981), 3785-3791.

\bibitem{satsuma1} R. Hirota, Y. Ohta, J. Satsuma;
\emph{Wronskian structures of solutions for soliton equations},
Prog. Theor. Phys. Suppl. 94 (1988), 59-72.

\bibitem{hirota9} R. Hirota;
\emph{The Direct Method in Soliton Theory}, Cambridge University Press,
New York, 2004.

\bibitem{jackson} F. H. Jackson;
\emph{On q functions and certain difference operator},
Trans. R. Soc. Edinb. 46 (1908), 253-281.

\bibitem{satsuma2} K. Kajiwara, J. Satsuma;
\emph{q-difference version of two-dimensional Toda Lattice equation},
J. Phys. Soc. Japan. 60 (1991), 3986-3989.

\bibitem{klr} B. Khesin B, V. Lyubashenko, C. Roger C;
\emph{Extensions and contractions of the Lie algebra of $q$-pseudodifferential
symbols on the circle}, J. Funct. Anal. 143 (1997), 55-97.

\bibitem{Ma} W. X. Ma;
\emph{Integrability}. In: Scott A, editor. Encyclopedia of Nonlinear Science,
Taylor \& Francis,  New York, 2005, pp. 450-453.

\bibitem{Ma0} W.X. Ma, W. Strampp;
\emph{Bilinear forms and B\"{a}cklund transformations of the perturbation systems},
 Phys. Lett. A. 341 (2005), 441-449.

\bibitem{olver} P. J. Olver;
\emph{Applications of Lie Groups to Differential Equations},
Springer-Verlag, New York, 1993.

\bibitem{sato} M. Sato;
 \emph{Soliton equations as dynamical systems on infinite dimensional
Grassmann manifolds}, RIMS Kokyuroku. 439 (1981), 30-46.

\bibitem{sch} M. P. Sch\"{u}tzenberger;
\emph{Une interpr\'{e}tation de certaines solutions de l'\'{e}quation
fonctionnelle: F(x + y) = F(x)F(y)}, C. R.  Acad. Sci. Paris. 236 (1953),
 352-353.

\bibitem{burcu5} B. Silindir;
\emph{Soliton solutions of $q$-Toda lattice
by Hirota Direct Method},  Advances in Difference Equations.
Article Number: 121 (2012),  1-22.

\bibitem{wahl} H.D. Wahlquist,  F.B. Estabrook;
\emph{B\"{a}cklund transformation for solutions of the Korteweg-de Vries
equation}, Phys. Rev. Lett. 31 (1973), 1386-1390.

\end{thebibliography}

\end{document}