text:
# 1 Introduction

## 1 Introduction

Analysis of nonlinear integrable systems is a subject of immense importance both from the physical and mathematical points of view. The inverse problem formulated on the basis of Lax operators can solve a wide class of problems. On the other hand, it is also important that problems be formulated with pre-assigned boundary conditions and solutions of the corresponding Cauchy problem be obtained. However at times the imposition of finite boundary conditions may lead to a loss of integrability. An important developement in this regard was the discovery by Sklyanin , who showed how non-periodic boundary conditions could be imposed on an integrable model without destroying its integrability. In this communication we have analysed a discrete nonlinear system, (also known as the DST model), initially studied by Ragnisco et al and also recently by Sklyanin et al from the view point of its integrability in presence of open boundary conditions. Our approach is different from that of Sklyanin in the sense that we have changed the form of the Lax operator at the end points where the boundary conditions are imposed. We have then shown that our approach can also accomodate the conditions laid down by Sklyanin for the system to be integrable under such circumstances. The present formalism was first used by Zhou  for a fermionic spin system. Both the classical and quantum hamiltonians have been derived explicitly. The quantum $`R`$ matrix has also been deduced so as to formulate the quantum inverse problem for the model under finite boundary conditions. Diagonalisation of the hamiltonian is carried out by using the method of separation of variables, as the boundary matrices are non diagonal in character. Of late there has been a great deal of interest in the study of canonical Bácklund transformations within the framework of classical $`r`$ matrix theory, specially under periodic boundary conditions . We have discussed this issue in case of quasi -periodic and finite boundary conditions. Finally in the case of the less stringent quasi-periodic boundary conditions we have shown how a connection with the eigenvalue of Baxter’s $`Q`$ operator may be established  by quantising the Bäcklund transformation.
The model under analysis is described by the following equations of motion.

$$\dot{q}_n=q_{n+1}q_n^2r_n$$

$$\dot{r}_n=r_{n1}+q_nr_n^2$$
$`(1.1)`$
where $`\dot{q}_n`$ stands for the time derivative of $`q_n`$. It was originally proposed by Ragnisco et al  in their analysis of Lie-Bäcklund symmetries of discrete systems. The Lax pair associated with (1.1) may be written as

$$\mathrm{\Psi }_{n+1}=L_n\mathrm{\Psi }_n\mathrm{\Psi }_{nt}=M_n\mathrm{\Psi }_n$$
$`(1.2)`$
where

$$L_n(\lambda )=\left(\begin{array}{cccc}\lambda +q_nr_n& q_n& & \\ r_n& 1& & \end{array}\right)M_n(\lambda )=\left(\begin{array}{cccc}\frac{\lambda }{2}& q_n& & \\ r_{n1}& \frac{\lambda }{2}& & \end{array}\right)$$
$`(1.3)`$
Note that consistency of (1.2) yields the equations of motion only for periodic boundary conditions. If however certain nontrivial boundary conditions are to be introduced then one has to adopt a different strategy as will be explained in the sequel. The system given by (1.1) possesses a hamiltonian structure with the following symplectic form .

$$\left(\begin{array}{ccc}q_n& & \\ r_n& & \end{array}\right)_t=\left(\begin{array}{cccc}0& 1& & \\ 1& 0& & \end{array}\right)\left(\begin{array}{cccc}\frac{\delta H}{\delta q_n}& & & \\ \frac{\delta H}{\delta r_n}& & & \end{array}\right)H=\frac{1}{2}\underset{n}{}(q_{n+1}r_n+q_nr_{n1}q_n^2r_n^2)$$
$`(1.4)`$
For the purpose of introducing the boundary conditions we modify the Lax operator in (1.2) in the following manner,

$$\mathrm{\Psi }_{j+1}=L_j(\lambda )\mathrm{\Psi }_j,j=1,2,\mathrm{}.N$$

$$\frac{d\mathrm{\Psi }_j}{dt}=M_j(\lambda )\mathrm{\Psi }_jj=2,3,\mathrm{}.N1$$

$$\frac{d\mathrm{\Psi }_{N+1}}{dt}=W_{N+1}(\lambda )\mathrm{\Psi }_{N+1}$$

$$\frac{d\mathrm{\Psi }_1}{dt}=W_1(\lambda )\mathrm{\Psi }_1$$
$`(1.5)`$
where $`W_{N+1},W_1`$ are two new $`2\times 2`$ matrices depending on the spectral parameter $`\lambda `$ and on the dynamical variables. The usual consistency condition of (1.2) viz

$$\frac{dL_j(\lambda )}{dt}=M_{j+1}L_j(\lambda )L_j(\lambda )M_j(\lambda )$$
$`(1.6)`$
is now replaced by the following set of equations:

$$\frac{dL_j(\lambda )}{dt}=M_{j+1}(\lambda )L_j(\lambda )L_j(\lambda )M_j(\lambda )(j=1,2\mathrm{}.N1)$$

$$\frac{dL_N(\lambda )}{dt}=W_{N+1}L_N(\lambda )L_N(\lambda )M_N(\lambda )$$

$$\frac{dL_1(\lambda )}{dt}=M_2(\lambda )L_1(\lambda )L_1(\lambda )W_1(\lambda )$$
$`(1.7)`$
That is the consistency conditions at the two ends are different. In the present case we have formally.

$$W_{N+1}=\left(\begin{array}{cccc}\frac{\lambda }{2}& q_{N+1}& & \\ r_N& \frac{\lambda }{2}& & \end{array}\right)W_1=\left(\begin{array}{cccc}\frac{\lambda }{2}& q_1& & \\ r_0& \frac{\lambda }{2}& & \end{array}\right)$$
$`(1.8)`$
so that upon imposing the boundary conditions $`r_0=\theta _{},q_{N+1}=\theta _+`$ the equations for the two ends of the discrete chain turn out to be as follows.

$$\dot{q}_1=q_2q_1^2r_1\dot{r}_1=(\theta _{}q_1r_1^2)$$

$$\dot{q}_N=(\theta _+q_N^2r_N)\dot{r}_N=(r_{N1}q_Nr_N^2)$$
$`(1.9)`$
We shall refer to the above system as an open chain.

## 2  Boundary conditions and classical $`r`$ matrix

To ascertain the hamiltonian associated with the open chain system it is always useful to compute the so called classical r-matrix through the poisson bracket given by the symplectic form in (1.4) . It is straight forward to show that.

$$\{L_n(\lambda )_,L_m(\mu )\}=[r(\lambda ,\mu ),L_n(\lambda )L_n(\mu )]\delta _{nm}$$
$`(2.1)`$
where $`r(\lambda ,\mu )=\frac{𝒫}{\lambda \mu },𝒫`$ standing for the permutation operator. The mondromy matrix is then defined by

$$T_N(\lambda )=\underset{n=1}{\overset{N}{}}L_n(\lambda )$$
$`(2.2)`$
From (2.1) and (2.2) it follows that

$$\{T_N(\lambda )_,T_N(\mu )\}=[r(\lambda ,\mu ),T_N(\lambda )T_N(\mu )]$$
$`(2.3)`$
We shall refer to this poisson algebra as the classical inverse scattering method (CISM-I) algebra. It follows that

$$\frac{dT_N}{dt}(\lambda )=\underset{n=1}{\overset{N}{}}L_N\mathrm{}L_{n+1}\dot{L}_nL_{n1}\mathrm{}.L_1(\lambda )=W_{N+1}T_N(\lambda )T_N(\lambda )W_1$$
$`(2.4)`$
from which it is easy to conclude that $`J(\lambda )=trT_N(\lambda )`$ is a constant of motion in the periodic case i.e when $`q_{m+N}=q_m`$ and $`r_{m+N}=r_m`$.

### 2.1 Quasiperiodic case

Instead of the periodic boundary conditions that are customarily assumed, one can impose the so called quasiperiodic boundary conditions on the system, without destroying its integrability. This requires the existence of a matrix $`C(\lambda ,\xi )`$ such that

$$\overline{J}(\lambda )=tr[C(\lambda ,\xi )T_N(\lambda )]$$
$`(2.5)`$
is time independent so that in the present situation one has

$$\frac{d\overline{J}(\lambda )}{dt}=tr[C(\lambda ,\xi )(M_{N+1}T_NT_NM_1)]=tr[C(\lambda ,\xi )M_{N+1}T_N]tr[M_1C(\lambda ,\xi )T_N]$$
$`(2.6)`$
from which it follows that

$$C(\lambda ,\xi )M_{N+1}=M_1C(\lambda ,\xi )$$
$`(2.7)`$
As a specific example one can take

$$C(\xi )=\left(\begin{array}{cccc}\xi ^{\frac{1}{2}}& 0& & \\ 0& \xi ^{\frac{1}{2}}& & \end{array}\right)$$
$`(2.8)`$
To analyse the structure of $`\overline{J}(\lambda )`$ in detail, we note the following asymptotics of the monodromy matrix in our case.

$$T_N^{11}(\lambda )\lambda ^N+\lambda ^{N1}S+p_2\lambda ^{N2}+\mathrm{}.$$

$$T_N^{12}(\lambda )\lambda ^{N1}q_1+O(\lambda ^{N2})$$

$$T_N^{21}(\lambda )\lambda ^{N1}r_N+O(\lambda ^{N2})$$

$$T_N^{22}\lambda ^{N2}p_2^{}+O(\lambda ^{N2})$$
$`(2.9)`$
whence we get

$$\overline{J}(\lambda )=tr(C(\xi )T_N(\lambda ))=\xi ^{\frac{1}{2}}\lambda ^N+(\xi ^{\frac{1}{2}}S)\lambda ^{N1}+(\xi ^{\frac{1}{2}}p_2)\lambda ^{N2}+(\xi ^{\frac{1}{2}}p_2^{})\lambda ^{N2}+\mathrm{}.$$
$`(2.10)`$
where

$$S=\underset{i=1}{\overset{N}{}}s_i,p_2=\underset{i=1}{\overset{N1}{}}(q_{i+1}r_i)+\underset{i<j}{}s_is_j,p_2^{}=r_Nq_1$$
$`(2.11)`$
One can verify, using the equations of motion, that each coefficient of $`\lambda ^j`$ is a constant. The following nonlinear combination of these constants then gives the Hamiltonian in the quasiperiodic case.

$$H=\frac{1}{2}\xi ^{\frac{1}{2}}(S^22p_2)+\xi ^{\frac{1}{2}}p_2^{}=\xi ^{\frac{1}{2}}[\underset{i=1}{\overset{N1}{}}q_{i+1}r_i\frac{1}{2}\underset{i}{}q_i^2r_i^2+\xi r_Nq_1]$$
$`(2.12)`$
### 2.2 Classical open chain

As mentioned earlier, the open chain in which independent boundary conditions are imposed at the two ends, represents a more general case compared to the quasi-periodic boundary conditions. The equations of motion are given by ( 1.1) and ( 1.9). Now by following the arguments given by Sklyanin in  one can show that

$$\tau (\lambda )=tr[K_+(\lambda )U(\lambda )],U(\lambda )=T_N(\lambda )K_(\lambda )T_N^1(\lambda )$$
$`(2.13)`$
represents the generator of the infinite conserved quantities with the matrices $`K_\pm (\lambda )`$ incorporating the effects of the nontrivial boundary conditions. In our case the classical $`r`$ matrix being symmetric, both $`K_\pm `$ satisfy the same equation namely

$$[r(\lambda ,\mu ),K(\lambda )K(\mu )]+K^1(\lambda )r(\lambda +\mu )K^2(\mu )K^2(\mu )r(\lambda +\mu )K^1(\mu )=0$$
$`(2.14)`$
where $`K^1(\lambda )=K(\lambda )I`$ and $`K^2(\lambda )=IK(\lambda )`$. The following solutions of (2.14) are found to be admissible .

$$K_{}(\lambda )=\left(\begin{array}{cccc}\theta _{}& \lambda & & \\ 0& \theta _{}& & \end{array}\right)K_+(\lambda )=\left(\begin{array}{cccc}\theta _+& 0& & \\ \lambda & \theta _+& & \end{array}\right)$$
$`(2.15)`$
By explicit multiplication we find using (2.9) and ( 2.15) in ( 2.13)

$$\tau (\lambda )=\lambda ^{N+2}+2\lambda ^N(a_2\frac{1}{2}S^2+q_1\theta _{}+r_N\theta _+)+O(\lambda ^{N2})$$
$`(2.16)`$
where $`a_2=_{i<j=1}^Ns_is_j+_{i=1}^{N1}q_{i=1}r_i`$ and $`S=_{i=1}^Ns_i`$ so that the hamiltonian is given by the expression.

$$H=\underset{i=1}{\overset{N1}{}}q_{i+1}r_i\frac{1}{2}\underset{i=1}{\overset{N}{}}q_i^2r_i^2+q_1\theta _{}+r_N\theta _+$$
$`(2.17)`$
To ascertain if $`\dot{\tau }(\lambda )=0`$ we observe that

$$\frac{d\tau }{dt}(\lambda )=tr[K_+(\lambda )\dot{T}_N(\lambda )K_{}(\lambda )T_N^1(\lambda )+K_+(\lambda )T_N(\lambda )K_{}(\lambda )(T_N^1(\lambda )\dot{T}_N(\lambda )T_N^1(\lambda ))]$$

$$\frac{d\tau }{dt}(\lambda )=tr[K_+(\lambda )W_N(\lambda )U(\lambda )]tr[K_+(\lambda )T_N(\lambda )W_1(\lambda )K_{}(\lambda )T_N^1(\lambda )]$$

$$tr[K_+(\lambda )U(\lambda )W_N(\lambda )]+tr[K_+(\lambda )T_N(\lambda )K_{}(\lambda )W_1(\lambda )T_N^1(\lambda )]$$
$`(2.18)`$
hence $`\dot{\tau }(\lambda )=0`$ if and only if

$$K_+(\lambda )W_{N+1}(\lambda )=W_{N+1}(\lambda )K_+(\lambda )$$
$`(2.19)`$
$$W_1(\lambda )K_{}(\lambda )=K_{}(\lambda )W_1(\lambda )$$
$`(2.20)`$
which are found to be identical to the conditions deduced by Sklyanin .

## 3 Bäcklund transformation

In this section we derive a Bäcklund transformation (BT) for the model under consideration, using the recently developed approach of Sklyanin . This relies on the algebra of the local lax operator as given by (2.1). For the sake of clarity we shall begin by retracing the methods used by Sklyanin in  and then consider the case where the nonlinear variables do not obey periodic boundary conditions. Sklyanin’s formalism requires us to look for another representation of the local algebra given in (2.1). Let $`g_i(\lambda \sigma )`$ be such a representation satisfying

$$\{g_i(\lambda \sigma )_,g_i(\mu \sigma )\}=[r(\lambda \mu ),g_i(\lambda \sigma )g_i(\mu \sigma )]$$
$`(3.1)`$
where $`\lambda `$ and $`\mu `$ are spectral parameters and $`\sigma `$ is a complex parameter of the transformation. Moreover $`g_i(\lambda \sigma )`$ is assumed to depend on a canonical set of variables $`(t_i,T_i)`$. In (3.1) $`r(\lambda \mu )=\frac{𝒫}{\lambda \mu }`$. To construct the BT

$$B_\sigma :T_N(x,X;\lambda )T_N(y,Y;\lambda )$$
$`(3.2)`$
we shall first construct a local transformation

$$b_\sigma ^{(i)}:L_i(x_i,X_i;\lambda )L_i(y_i,Y_i;\lambda )$$
$`(3.3)`$
defined by the following

$$g_i(\lambda \sigma ;s_i,S_i)L_i(\lambda ;x_i,X_i)=L_i(\lambda ;y_i,Y_i)g_i(\lambda \sigma ;t_i,T_i)$$
$`(3.4)`$
here $`L_i(\lambda ;x_i,X_i)`$ is the same as in (1.3), we have simply relabelled the variables $`q_i,r_i`$ by $`x_i,X_i`$ for notational convenience. In (3.4) we see that in addition to the canonical variables $`(x_i,X_i)_{i=1}^N`$ and $`(y_i,Y_i)_{i=1}^N`$ defining the BT we have two additional sets of canonical variables $`(s_i,S_i)_{i=1}^N`$ and $`(t_i,T_i)_{i=1}^N`$. It is assumed that all the four sets of canonical variables are mutually commuting quantities. The transformation $`b_\sigma ^{(i)}`$ is therefore a transformation in an extended phase space. Now by imposing suitable constraints it is possible to eliminate the auxillary variables $`(s_i,Si)`$ and $`(t_i,T_i)`$ from (3.4) thereby obtaining a local transformation $`b_\sigma ^{(i)}`$. We get from (3.4) the following set of independent equations, taking $`g_i(\lambda \sigma ;s_i,S_i)`$ to be of the form.

$$g_i(\lambda \sigma );s_i,S_i)=(\begin{array}{cccc}1& s_i& & \\ S_i& \lambda \sigma s_iS_i& & \end{array})$$
$`(3.5)`$
$$(x_i+s_i)X_i=y_i(Y_iT_i),S_iX_i(x_i+s_i)=Y_iT_i+\sigma X_iy_i$$

$$t_i=y_i,X_i=S_i$$
$`(3.6)`$
Imposition of the constraints $`s_i=t_{i+1}`$ and $`S_i=T_{i+1}`$ in these equations leads to

$$X_i=\frac{1}{y_i}\frac{\sigma }{x_iy_{i+1}}Y_i=X_{i1}+\frac{x_iy_{i+1}}{y_i}X_i$$
$`(3.7)`$
which defines the local transformation $`b_\sigma ^{(i)}:(x_i,X_i)(y_i,Y_i)`$. It should be pointed out that this local BT may be obtained from a generating function $`G_\sigma `$ where

$$G_\sigma =\underset{i=1}{\overset{N}{}}[\frac{x_iy_{i+1}}{y_i}+\sigma \mathrm{log}(\frac{x_iy_{i+1}}{y_i})]$$
$`(3.8)`$
$$X_i=\frac{G_\sigma }{x_i},Y_i=\frac{G_\sigma }{y_i}$$
$`(3.9)`$
Moreover as a consequence of the constraints $`s_i=t_{i+1},S_i=T_{i+1}`$ and (3.6) we see that

$$g_i(\lambda \sigma ;s_i,S_i)=g_{i+1}(\lambda \sigma ;t_i,T_i)$$
$`(3.10)`$
and hence (3.4) becomes

$$g_{i+1}(\lambda \sigma ;y_{i+1},X_i)L_i(\lambda ;x_i,X_i)=L_i(\lambda ;x_iX_i)g_i(\lambda \sigma ;y_i,X_{i1})$$
$`(3.11)`$
Now it is well known that the generating function for the integrals of motion are given by the transfer matrix $`J(\lambda )`$ defined as the trace of the mondromy matrix, i.e

$$J(\lambda )=trT_N(x,X;\lambda )$$
$`(3.12)`$
where $`T_N(\lambda )`$ obeys the CISM-I algebra.

$$\{T_N(\lambda )_,T_N(\mu )\}=[r(\lambda \mu ),T_N(\lambda )T_N(\mu )]$$
$`(3.13)`$
Under the BT defined by (3.11) the transfer matrix assumes the form.

$$J(\lambda )=tr[g_{N+1}^1(\lambda \sigma ;y_{N+1},X_N)T_N(y,Y;\lambda )g_1(\lambda \sigma ;y_1,X_0)]$$
$`(3.14)`$
Clearly invariance of the integrals of motion and thus of the hamiltonian generated by $`J(\lambda )`$ requires that the variables in $`g_i(\lambda \sigma )`$ be periodic i.e

$$X_{N+m}=X_m,y_{N+m}=y_m$$
$`(3.15)`$
In such cases the BT will be canonical in nature since the hamiltonian is invariant.
In the quasi periodic case the boundary conditions on the dynamical variables are determined by (2.7), so that upon using the form of $`C(\xi )`$ as given in (2.8) we find that the variables in $`L_i(\lambda )`$ obey

$$x_{N+1}=\xi X_N,X_0=\xi X_N$$
$`(3.16)`$
We refer to these as the quasi periodic boundary conditions. Now in view of (2.5) which defines the generator of the integrals of motion $`\overline{J}(\lambda )`$ we see that as a consequence of the BT given by (3.12), it is modified to the following.

$$\overline{J}(\lambda )=tr[C(\lambda ,\xi )g_{N+1}^1(\lambda \sigma ;y_{N+1},X_N)T_N(y,Y;\lambda )g_1(\lambda \sigma ;y_1,X_0)]$$
$`(3.17)`$
Clearly invariance of $`\overline{J}(\lambda )`$ requires that

$$g_1(\lambda \sigma ;y_1,X_0)C(\lambda ,\xi )g_{N+1}^1(\lambda \sigma ;y_{N+1},X_N)=C(\lambda ,\xi )$$
$`(3.18)`$
which gives.

$$y_{N+1}=\xi y_1,X_0=\xi X_N$$
$`(3.19)`$
Structurally these equations are consistent with (3.16) because the Lax operator $`L_i(\lambda ;y_i,Y_i)`$ in (3.11) is obtained simply by replacing the variables $`(x_i,X_i)_{i=1}^N`$ with $`(y_i,Y_i)_{i=1}^N`$ . Therefore the BT generated by (3.11) through the matrix $`g_I(\lambda \sigma )`$ is a canonical transformation even when the nonlinear dynamical variables obey quasiperiodic boundary conditions. It is now natural to enquire if a BT can be constructed in case of open chains. Infact the most general type of boundary conditions for Liouville integrable lattices was formulated some years back by Sklyanin in terms of the representations of a new algebra known as the reflection equation algebra (REA). The classical limit of the REA is given by the following Poisson algebra, which shall be referred to as the CISM-II algebra.

$$\{T^{(1)}(\lambda ),T^{(2)}(\mu )\}=[r(\lambda \mu ),T^{(1)}(\lambda )T^{(2)}(\mu )]+T^{(1)}(\lambda )r(\lambda +\mu )T^{(2)}(\mu )T^{(2)}(\mu )r(\lambda +\mu )T^{(1)}(\lambda )$$
$`(3.20)`$
Now if $`K_\pm (\lambda )`$ be representations of the above algebra then it may be shown that $`U(\lambda )=T_N(\lambda )K_{}(\lambda )T_N^1(\lambda )`$ also satisfies the CISM-II algebra, where $`T_N(\lambda )`$ is defined by (2.2). Consequently it follows that

$$\tau (\lambda )=trK_+(\lambda )U(\lambda )$$
$`(3.21)`$
is a generating function of the integrals of motion. We have already seen that under the BT given by (3.11) the monodromy matrix undergoes the following transformation.

$$B_\sigma :T_N(x,X;\lambda )g_{N+1}^1(\lambda \sigma ;y_{N+1},X_N)T_N(y,Y;\lambda )g_1(\lambda \sigma ;y_1,X_0)$$
$`(3.22)`$
Next consider $`2\times 2`$ matrices $`V_\pm (\lambda )`$ and construct the following quantities:

$$Z(\lambda )T_N(x,X;\lambda )V_{}(\lambda )T_N^1(x,X;\lambda )\stackrel{~}{J}(\lambda )tr(V_+(\lambda )Z(\lambda ))$$
$`(3.23)`$
note that $`Z(\lambda )`$ per se might not satisfy the CISM -II algebra. But if we apply a BT as given by (3.11) then $`\stackrel{~}{J}(\lambda )`$ assumes the following form.

$$\stackrel{~}{J}(\lambda )=tr[\{g_{N+1}(\lambda \sigma )V_+(\lambda )g_{N+1}^1(\lambda \sigma )\}T_N(y,Y;\lambda )\times $$

$$\times \{g_1(\lambda \sigma )V_{}(\lambda )g_1^1(\lambda \sigma )\}T_N^1(y,Y;\lambda )]$$
$`(3.24)`$
It is evident that $`\stackrel{~}{J}(\lambda )`$ would be a generator of the integrals of motion provided that the quantities in curly brackets are representations of the CISM -II algebra. This infact opens up the possibility of deriving a hierarchy of suitable boundary conditions for integrable systems starting from appropriate matrices $`V_+`$. We will examine this possibility by deriving the usual scalar boundary matrices. First of all if the quantities in curly brackets were similarity transformations, then we could have taken the matrices $`V_\pm `$ to be simply the matrices $`K_\pm (\lambda )`$ as defined in ( 2.15). This is a consequence of the observation that the CISM-II algebra is invariant under a similarity transformation of its elements. However this possibility is spoilt by the fact that the arguments of $`g_{1(N+1)}`$ and $`g_{1(N+1)}^1`$ in ( 3.24) are different. Consider therefore the following forms of the matrices $`V_\pm (\lambda )`$

$$V_+(\lambda )=\frac{\lambda }{\lambda +\sigma }\left(\begin{array}{cccc}a+\frac{d}{\lambda }& b& & \\ 1& a+\frac{d}{\lambda }& & \end{array}\right)$$
$`(3.25)`$
$$V_{}(\lambda )=\frac{\lambda }{\lambda \sigma }\left(\begin{array}{cccc}\lambda A_1+A_0+\frac{\delta }{\lambda }& \beta \lambda ^2+B_0& & \\ C_0& \lambda A_1A_0+\frac{\delta }{\lambda }& & \end{array}\right)$$
$`(3.26)`$
so that

$$g_{N+1}(\lambda \sigma )V_+(\lambda )g_{N+1}^1(\lambda \sigma )=K_+(\lambda )=\left(\begin{array}{cccc}\theta _+& 0& & \\ \lambda & \theta _+& & \end{array}\right)$$
$`(3.27)`$
and

$$g_1(\lambda \sigma )V_{}(\lambda )g_1^1(\lambda \sigma )=K_+(\lambda )=\left(\begin{array}{cccc}\theta _{}& \lambda & & \\ 0& \theta _{}& & \end{array}\right)$$
$`(3.28)`$
Here

$$a=y_{N+1}\theta _+d=\sigma \theta _+,b=y_{N+1}(2\theta _+y_{N+1})$$
$`(3.29)`$
and $`A_1=X_0,A_0=(\sigma X_0+\theta _{}y_1X_0^2),\delta =\sigma \theta _{}`$

$$\beta =1,C_0=X_0^2,B_0=(y_1X_0\sigma )^2+2\theta _{}y_1$$
$`(3.30)`$
Disregarding the overall factor in $`V_+(\lambda )`$ we see that its structure is similar to the general scalar solution of the CISM -II algebra. While that of $`V_{}(\lambda )`$ is similar to that derived in  ( compare their eqns ( 2.3) -(2.5) with $`\alpha =\gamma =0,\beta =1)`$. The difference in our case is that all the entries in $`V_{}(\lambda )`$ are scalars. Thus while the former is by itself a solution of the CISM -II algebra, for $`V_{}(\lambda )`$ to do the same it is necessary that $`X_0=\theta _{}=0`$. The other pertinent issue is whether the matrices $`V_\pm (\lambda )`$ are compatible with the integrability of the model under consideration. the constraint of integrability requires these matrices to satisfy (2.20) which is a stringent condition. In case of $`V_+(\lambda )`$ the boundary condition $`y_{N+1}=\theta _+`$ ensures $`a=b=0`$ so that the integrability condition is not fulfiled unless $`Y_N=constant`$ which clearly is in contradiction to the fact that $`(y_i,Y_i)_{i=1}^N`$ are dynamical variables. Thus $`V_\pm (\lambda )`$ by themselves cannot be made to satisfy the requirements of integrabilty and CISM -II algebra simultaneously. However as a result of the BT the quantities $`Z(\lambda )`$ and $`\stackrel{~}{J}(\lambda )`$ are modified in a such a manner that the new boundary matrices $`K_\pm (\lambda )`$ appearing in lieu of $`V_\pm (\lambda )`$ are not only solutions of CISM -II algebra but are also compatible with the integrability of the model under consideration.

## 4 Quantum open chain

In this section we consider the quantum inverse problem for a open chain. Note that the Lax operator in (1.3) is essentially similar to the discrete self trapping (DST) model . The integrability and quantum inverse problem in absence of finite boundary conditions were studied in . From the hamiltonian structure of our system we have seen that $`(q_n,r_n)`$ are a conjugate pair of variables; so in the quantum version of the model we impose the following commutation relations.

$$[q_n,r_m]=\eta \delta _{nm}(\eta =i\mathrm{})[q_n,q_m]=[r_n,r_m]=0$$
$`(4.1)`$
As before the monodromy matrix is defined by $`T_N(\lambda )=_{n=1}^NL_n(\lambda )`$ and obeys the QISM -I algebra

$$R_{12}(\lambda \mu )T_N^1(\lambda )T_N^2(\mu )=T_N^2(\mu )T_N^1(\lambda )R_{12}(\lambda \mu )$$
$`(4.2)`$
with $`R_{12}(\lambda \mu )=+\frac{\eta }{\lambda \mu }𝒫`$. The matrices determining the boundary conditions $`K_\pm (\lambda )`$ obey the conditions .

$$R_{12}(\lambda \mu )K_{}^{(1)}(\lambda )R_{12}(\lambda +\mu )K_{}^{(2)}(\lambda )=K_{}^{(2)}(\mu )R_{12}(\lambda +\mu )K_{}^{(1)}R_{12}(\lambda \mu )$$
$`(4.3)`$
$$R_{12}(\lambda +\mu )K_+^{(1)t_1}(\lambda )R_{12}(\lambda \mu 2\eta )K_+^{(2)t_2}(\mu )$$

$$=K_+^{(2)t_2}(\mu )R_{12}(\lambda \mu 2\eta )K_+^{(1)t_1}(\lambda )R_{12}(\lambda +\mu )$$
$`(4.4)`$
The form of $`K_\pm (\lambda )`$ given in (2.15) is seen to satisfy these equations. The generators of the conserved densities are then built from.

$$U(\lambda )=\left(\begin{array}{cccc}A(\lambda )& B(\lambda )& & \\ C(\lambda )& D(\lambda )& & \end{array}\right)=T(\lambda )K_{}(\lambda \eta /2,\xi _{})\sigma _2T^t(\lambda )\sigma _2$$
$`(4.5)`$
which obeys the algebra

$$R_{12}(\lambda \mu )U^1(\lambda )R_{12}(\lambda +\mu \eta )U^2(\lambda )=U^2(\lambda )R_{12}(\lambda +\mu \eta )U^1(\lambda )R_{12}(\lambda \mu )$$
$`(4.6)`$
and is given by

$$\tau (\lambda )=tr[K_+(\lambda \eta /2,\xi _+)U(\lambda )]=\lambda ^{2N+2}+2\lambda ^{2N}\widehat{H}_q+\mathrm{}\mathrm{}\mathrm{}.$$
$`(4.7)`$
where the hamiltonian with quantum corrections is given by

$$\widehat{H}_q=\underset{i=1}{\overset{N1}{}}q_{i+1}r_i\frac{1}{2}\underset{i=1}{\overset{N}{}}(q_ir_i)^2\frac{\eta ^2}{8}+\xi _+r_N+\xi _{}q_1$$
$`(4.8)`$
On the other hand with $`U(\lambda )`$ given by (4.5), we can write the general form of $`\tau (\lambda )`$ as

$$\tau (\lambda )=\xi _+(A(\lambda )+D(\lambda ))+(\lambda +\eta /2)B(\lambda )$$
$`(4.9)`$
### 4.1 Separation of variables

From (4.6) we can find the operator algebra for the elements $`A(\lambda ),B(\lambda )`$ etc. It is however more convenient to define instead of $`D(\lambda )`$ the following operator.

$$D^{}(\lambda )=2\lambda D(\lambda )\eta A(\lambda )$$
$`(4.10)`$
The commutation relations relevant for our purposes are then obtainable from (4.3) and are as follows.

$$[B(\lambda ),B(\mu )]=0$$
$`(4.11)`$
$$A(\lambda )B(\mu )=\frac{(\lambda \mu \eta )(\lambda +\mu \eta )}{(\lambda \mu )(\lambda +\mu )}B(\mu )A(\lambda )+\frac{\eta (2\mu \eta )}{(\lambda \mu )2\mu }B(\lambda )A(\mu )\frac{\eta }{(\lambda +\mu )2\mu }B(\lambda )D^{}(\mu )$$
$`(4.12)`$
$$D^{}(\lambda )B(\mu )=\frac{(\lambda \mu +\eta )(\lambda +\mu \eta )}{(\lambda \mu )(\lambda +\mu )}B(\mu )D^{}(\lambda )+$$

$$+\frac{\eta (2\lambda +\eta )(2\mu \eta )}{(\lambda +\mu )2\mu }B(\lambda )A(\mu )\frac{\eta (2\lambda +\eta )}{(\lambda \mu )2\mu }B(\lambda )D^{}(\mu )$$
$`(4.13)`$
Furthermore from the definition of $`U(\lambda )`$ as given in (4.5) the operators $`A(\lambda ),B(\lambda )`$ etc are found to have the following forms:

$$A(\lambda )=()^Nr_N\{\lambda ^{2N}+(S+\eta /2)\lambda ^{2N1}+\mathrm{}\mathrm{}..\}$$
$`(4.14)`$
$$D(\lambda )=()^Nr_N\{\lambda ^{2N}+(S\eta /2)\lambda ^{2N1}+\mathrm{}\mathrm{}..\}$$
$`(4.15)`$
$$B(\lambda )=()^N(\lambda \eta /2)\{\lambda ^{2N}+\mathrm{}\mathrm{}\}()^N(\lambda \eta /2)\underset{\alpha =1}{\overset{N}{}}(\lambda \widehat{u}_\alpha )(\lambda +\widehat{u}_\alpha )$$
$`(4.16)`$
Where $`S=_{i=1}^Nq_ir_i`$. From ( 4.9) we note that the presence of the operator $`B(\lambda )`$ in the expression for the transfer matrix $`\tau (\lambda )`$ prevents us from setting up the usual algebraic Bethe ansatz. Consequently we have recast $`B(\lambda )`$ in the form given in (4.16) where $`\widehat{u}_\alpha `$ are the zeros of this operator. Henceforth $`\{\widehat{u}_\alpha \alpha =1,2\mathrm{}N\}`$ are themselves to be regarded as operators defining $`B(\lambda )`$ and it follows from their commutativity given by (4.11) that

$$[\widehat{u}_\alpha ,\widehat{u}_\beta ]=0$$
$`(4.17)`$
We next define operators $`\widehat{v}_\alpha ^\pm `$ by left substitution $`\lambda \widehat{u}_\alpha `$ into the operators $`A(\lambda )`$ and $`D^{}(\lambda )`$ as shown below :

$$\widehat{v}_\alpha ^+D^{}(\lambda \widehat{u}_\alpha ),\widehat{v}_\alpha ^{}A(\lambda \widehat{u}_\alpha )$$
$`(4.18)`$
It then follows from the algebra of the operators $`A(\lambda ),D^{}(\lambda )`$ that

$$[\widehat{v}_\alpha ^\pm ,\widehat{u}_\beta ]=\pm \eta \widehat{v}_\alpha ^\pm \delta _{\alpha \beta }[\widehat{v}_\alpha ^+,\widehat{v}_\beta ^+]=0[\widehat{v}_\alpha ^{},\widehat{v}_\beta ^{}]=0$$
$`(4.19)`$
Notice that another set of dynamical variables can be constructed from the zeros of $`B(\lambda )`$ at $`\lambda =\widehat{u}_\alpha ,\alpha =1,2\mathrm{}N`$ and obeying

$$[\widehat{w}_\alpha ^\pm ,\widehat{u}_\beta ]=\pm \eta \widehat{w}_\alpha ^\pm \delta _{\alpha \beta }[\widehat{w}_\alpha ^+,\widehat{w}_\beta ^+]=0[\widehat{w}_\alpha ^{},\widehat{w}_\beta ^{}]=0$$
$`(4.20)`$
where $`\widehat{w}_\alpha ^+=A(\lambda \widehat{u}_\alpha )`$ and $`\widehat{w}_\alpha ^{}=D^{}(\lambda \widehat{u}_\alpha )`$. Then by using Lagrange interpolation one can rewrite the operators $`A(\lambda )`$ and $`D^{}(\lambda )`$ in the following manner .

$$A(\lambda )=\underset{\alpha =1}{\overset{2N}{}}\underset{\beta =1,\beta \alpha }{\overset{2N}{}}\frac{\lambda \widehat{\lambda }_\beta }{\widehat{\lambda }_\alpha \widehat{\lambda }_\beta }\widehat{\mu }_\alpha ^++()^Nr_N\underset{\alpha =1}{\overset{2N}{}}(\lambda \widehat{\lambda }_\alpha )$$
$`(4.21)`$
$$D^{}(\lambda )=\underset{\alpha =1}{\overset{2N}{}}\underset{\beta =1,\beta \alpha }{\overset{2N}{}}\frac{\lambda \widehat{\lambda }_\beta }{\widehat{\lambda }_\alpha \widehat{\lambda }_\beta }\widehat{\mu }_\alpha ^{}+()^Nr_N\underset{\alpha =1}{\overset{2N}{}}(\lambda \widehat{\lambda }_\alpha )$$
$`(4.22)`$
$$\widehat{\lambda }_\alpha =\widehat{u}_\alpha ,1\alpha N$$

$$\widehat{\lambda }_\alpha =\widehat{u}_\alpha N+1\alpha 2N$$

$$\widehat{\mu }_\alpha ^\pm =\widehat{v}_\alpha ^{}1\alpha N$$

$$\widehat{\mu }_\alpha ^\pm =\widehat{w}_\alpha ^\pm N+1\alpha 2N$$

The transfer matrix then assumes the form:

$$2\widehat{u}_\alpha \tau (\widehat{u}_\alpha )=\xi _+(2\widehat{u}_\alpha +\eta )v_\alpha ^{}+\xi _+v_\alpha ^+$$
$`(4.23)`$
Now by extensively using the preperties of quantum determinants  it can be shown that from the commutation relations for the operators $`\widehat{u}_\alpha ,v_\alpha ^\pm ,w_\alpha ^\pm `$ one can define a Hilbert space where the eigenvalue problem reduces to

$$2u_\alpha \tau (u_\alpha )\mathrm{\Psi }(u_1,\mathrm{}..u_N)=\xi _+(2u_\alpha \eta )\mathrm{\Delta }^{}(u_\alpha )\mathrm{\Psi }(\mathrm{}..,u_\alpha \eta ,\mathrm{}.)+\xi _+\mathrm{\Delta }^+(u_\alpha )\mathrm{\Psi }(\mathrm{},u_\alpha +\eta ,\mathrm{}.)$$
$`(4.24)`$
where

$$\mathrm{\Delta }^{}(u)=\xi _{}+(u\frac{\eta }{2})\mathrm{\Delta }^+(u)=(2u\eta )(\xi _{}(u+\frac{\eta }{2}))$$
$`(4.25)`$
Setting

$$\mathrm{\Psi }(u_1,\mathrm{}u_N)=\underset{\alpha =1}{\overset{N}{}}\varphi (u_\alpha )$$
$`(4.26)`$
we have finally

$$2u_\alpha \tau (u_\alpha )\varphi (u_\alpha )=\xi _+(2u_\alpha +\eta )\mathrm{\Delta }^{}(u_\alpha )\varphi (u_\alpha \eta )+\xi _+\mathrm{\Delta }^+\varphi (u_\alpha +\eta )$$
$`(4.27)`$
## 5 Relation with Baxter’s Q-operator

In this section we consider a quantum mechanical version of the Bäcklund transformation and show how in the quasi periodic case it may be connected to Baxter’s Q-Operator. For this purpose we quantise the dynamical variables occurring in the Lax operator given in (1.3) which then asssumes the form.

$$L_i(\lambda )=\left(\begin{array}{cccc}\lambda \eta q_i_{q_i}& q_i& & \\ \eta _{q_i}& 1& & \end{array}\right)$$
$`(5.1)`$
where we have set $`r_i=\eta _{q_i}`$ and $`[q_i,r_j]=\eta \delta _{ij}`$. The quantum analog of the BT for integrable DST model which essentially has the same Lax operator as ours was studied in great detail by Kuznetsov, et al in ; others have also studied the same model from alternative points of view . The essential idea is to reproduce an operator satisfying the following conditions,

$$\widehat{J}(\lambda )\widehat{Q}(\lambda )=\mathrm{\Delta }_+(\lambda )\widehat{Q}(\lambda +\eta )+\mathrm{\Delta }_{}(\lambda )\widehat{Q}(\lambda \eta )$$
$`(5.2)`$
$$\widehat{Q}(\lambda )\widehat{Q}(\mu )=\widehat{Q}(\mu )\widehat{Q}(\lambda )$$
$`(5.3)`$
$$\widehat{J}(\mu )\widehat{Q}(\lambda )=\widehat{Q}(\lambda )\widehat{J}(\mu )$$
$`(5.4)`$
and having a common set of eigenfunctions

$$\widehat{Q}(\lambda )\mathrm{\Phi }=Q(\lambda )\mathrm{\Phi },\widehat{J}(\lambda )\mathrm{\Phi }=J(\lambda )\mathrm{\Phi }$$
$`(5.5)`$
with the eigenvalues satisfying a finite difference equation.

$$\tau (\lambda )Q(\lambda )=\mathrm{\Delta }_+(\lambda )Q(\lambda +\eta )+\delta _{}(\lambda )Q(\lambda \eta )$$
$`(5.6)`$
Incidentally the last equation arose in Baxter’s analysis of the XXX model where $`\tau (\lambda )`$ was the transfer matrix of the model with $`Q(\lambda )`$ being a suitable polynomial.
To prove (5.3) it is sufficient to demand the existence of an operator $`\widehat{R}_\sigma `$ such that

$$𝒢(\lambda \sigma ;t,_t)L(\lambda ;y,_y)\widehat{R}_\sigma =\widehat{R}_\sigma L(\lambda ;q,_q)𝒢(\lambda \sigma ;s,_s)$$
$`(5.7)`$
where we have dropped the indices $`{}_{}{}^{}i_{}^{}`$ on the variables for brevity. Note that here $`𝒢(\lambda \sigma )`$ is assumed to satisfy the following algebra.

$$R_{12}(\lambda \mu )𝒢^1(\lambda )𝒢^2(\mu )=)𝒢^2(\mu ))𝒢^1(\lambda )R_{12}(\lambda \mu )$$
$`(5.8)`$
and as before

$$R(\lambda \mu )=+\frac{\eta }{\lambda \mu }𝒫$$
$`(5.9)`$
regarding the explicit form of $`𝒢(\lambda \sigma )`$ we take it to be as follows.

$$𝒢(\lambda \sigma ;s,_s)=\left(\begin{array}{cccc}\lambda \sigma s_s& s& & \\ \eta _s& 1& & \end{array}\right)$$
$`(5.10)`$
On the basis of the above remarks proceeding as in  we finally arrive at a form of $`\widehat{Q}(\sigma )`$ as an integral operator with a kernel given by

$$𝒬(\sigma )=\underset{i=1}{\overset{N}{}}w_i(\sigma ;y_{i+1},y_i,q_i)$$
$`(5.11)`$
where

$$w_i(\sigma ;y_{i+1},y_i,q_i)=\frac{i}{2\pi }\mathrm{\Gamma }(\sigma /\eta +1)y_i^1(\frac{y_{i+1}q_i}{y_i})^{\frac{\sigma }{\eta }1}exp(\frac{y_{i+1}q_i}{\eta y_i})$$
$`(5.12)`$
To derive the last property of Baxter’s Q-operator we note that $`\widehat{J}(\lambda )=tr[C(\lambda ,\xi )T_N(\lambda )]`$ in the quasi periodic case. So we consider

$$[\widehat{Q}(\sigma )\widehat{J}(\sigma )\mathrm{\Phi }](\stackrel{}{y})=tr𝑑q^N𝒬(\sigma )T_N(\sigma )\mathrm{\Phi }(\stackrel{}{q})$$
$`(5.13)`$
which yields for the kernel of the joint operator $`\widehat{Q}(\sigma )\widehat{J}(\sigma )`$

$$[\widehat{Q}(\sigma )\widehat{J}(\sigma )\mathrm{\Phi }](\stackrel{}{y}\stackrel{}{x})=tr\{C(\xi )T_N^{}(\sigma )\underset{i=1}{\overset{N}{}}w_i\}=tr\{C(\xi )\underset{i=1}{\overset{N}{}}L^{}(\sigma )w_i\}$$
$`(5.14)`$
where $`L^{}(\sigma )`$ represents the adjoint of the corresponding local Lax operator. The last equation may be recast into the following form

$$[\widehat{Q}(\sigma )\widehat{J}(\sigma )\mathrm{\Phi }](\stackrel{}{y}\stackrel{}{x})=𝒬(\sigma )tr[C(\xi )\stackrel{~}{T}_N(\sigma )]$$
$`(5.15)`$
with

$$\stackrel{~}{T}_N(\sigma )=\underset{i=1}{\overset{N}{}}\stackrel{~}{L}_i(\sigma )=\underset{i=1}{\overset{N}{}}\left(\begin{array}{cccc}\sigma +\eta +\eta q_i_{q_i}\mathrm{log}w_i& q_i& & \\ \eta _{q_i}\mathrm{log}w_i& 1& & \end{array}\right)$$
$`(5.16)`$
It is now possible to triangularise the matrix $`\stackrel{~}{L}_i(\sigma )`$ by the following transformation:

$$\stackrel{~}{L}_i(\sigma )S_{i+1}^1\stackrel{~}{L}_i(\sigma )S_i,S_i=\left(\begin{array}{cccc}1& y_{i+1}& & \\ 0& 1& & \end{array}\right)$$
$`(5.17)`$
This transformation leaves the quasiperiodic boundary matrix $`C(\xi )`$ invariant, while causing $`\stackrel{~}{L}_i(\sigma )`$ to become triangular as shown below.

$$\stackrel{~}{L}_i(\sigma )\left(\begin{array}{cccc}\sigma \frac{w_i(\frac{\sigma }{\eta }1)}{w_i(\frac{\sigma }{\eta })}& 0& & \\ & \frac{w_i(\frac{\sigma }{\eta }+1)}{w_i(\frac{\sigma }{\eta })}& & \end{array}\right)$$
$`(5.18)`$
It is then straightforward to show that

$$J(\sigma ,\xi )\underset{i=1}{\overset{N}{}}w_i(\frac{\sigma }{\eta })=\xi ^{\frac{1}{2}}\sigma ^N\underset{i=1}{\overset{N}{}}w_i(\sigma /\eta 1)+\xi ^{\frac{1}{2}}\underset{i=1}{\overset{N}{}}w_i(\sigma /\eta +1)$$
$`(5.19)`$
This relation ought to be compared with the eigenvalue deduced by the standard method of algebraic Bethe ansatz for the quasi periodic case which is of the form:

$$\mathrm{\Lambda }(\sigma )\underset{i=1}{\overset{m}{}}(\sigma \mu _i)=\xi ^{\frac{1}{2}}\sigma ^N\underset{i=1}{\overset{m}{}}(\sigma \mu _i\eta )+\xi ^{\frac{1}{2}}\underset{i=1}{\overset{m}{}}(\sigma \mu _i+\eta )$$
$`(5.20)`$
## 6 Discussion

We have analysed a discrete nonlinear integrable system in the presence of finite boundary conditions and have shown how the diagonalisation problem can be handled when the boundary matrices are upper (lower) triangular in nature. The construction of Bäcklund transformations under different types of boundary conditions was analysed and for the quasiperiodic case we have shown that a connectionbetween the quantised form of the BT and Baxter’s Q operator may be established.
One of us (AGC) wishes to acknowledge the financial support provided by UGC ( Govt. of India) PSW-026/99-00.

## 7 References

$`[1]`$ E.K. Sklyanin, (1988)J. Phys. A. Math. Gen. 21, 2375
$`[2]`$ I. Merola, O. Ragnisco and Tu G. Zhang (1994) Inverse Problems101315
$`[3]`$ H. Q. Zhou (1996)J. Phys. A math. Gen. 29 L607
$`[4]`$ H. Zhang, G.Z.Tu, W. Oevel and B. Fuchssteiner, (1991) J. Math. Phys 32 1908.
$`[5]`$ E.K. Sklyanin, (1987) Func. Anal. and Appl. 21(2), 86
$`[6]`$ V. B. Kuznetsov, and E. K. Sklyanin (1998) J. Phys. A Math. Gen31, 2241
$`[7]`$ A. N. W. Hone, V. B. Kuznetsov and O. Ragnisco , Solv-int /9904003
$`[8]`$ L.D. Faddeev and L.A. Takhtajan, (1979) Russian math. Sueveys 34 ( 5), 11; V. A. Korepin, N. M. Bogoliubov and A. G. Izergin Quantum Inverse Scattering Method and Correlation Functions, Cambridge University Press, 1993
$`[9]`$ E. K Sklyanin in Nankai lactures on Mathematical Physics ” Quantum groups and Quantum Integrable Systems” World Scientific Ed. Mo-Lin Ge pp 63 (1992)
$`[10]`$ V. B. Kuznetsov, M. F. Jorgensen and P. L. Christiansen 91995) J. Phys. A Math Gen. 28 4639
$`[11]`$ V. Pasquier, M. Gaudin (1992)J. Phys. A. Math. Gen 25, 5243; R. J. Baxter, Exactly Solved models in Statistical mechanics Academic Press, New York/London (1982)
$`[12]`$ V. Z. Enol’skii, M. salerno, N. A. Kostov and A. C. Scott (1991) Physica Scripta43, 229
$`[13]`$ V. Z. Enol’skii, V. B. Kuznetsov and M. Salerno (1993)Physica D 68, 138
$`[14]`$ K. Hikami, 91996)J. Phys. Soc. Jpn65(5), 1213
$`[15]`$ V. B. Kuznetsov, M. Salerno and E.K.Sklyanin 92000)J. Phys. A Math. Gen 33, 171
$`[16]`$ E.K. Sklyanin in ”Bäcklund transformation and Baxter’s Q operator” Lect. at the seminar: ”Integrable systems from Classical to Quantum” ( Univ. de Montreal, July 26- August 6, 1999) and S. E Derkachov solv-int/9902015


format: text

text:
# Nonparametric modeling and spatiotemporal dynamical systems

## 1 Introduction

One of the major goals of spatiotemporal data analysis consists in inferring a model from a given spatiotemporal data set. Conventional ways of modeling are purely theoretical arguing and a posteriori evaluation of possible models with suitable measures. In this article, statistical methods and applications to infer generalized additive models from data are presented. The attention is put on spatio-temporal systems, but there are certainly more possibilities of using the ideas in physics and engineering. Generalized additive models have been introduced to statistical analysis in the 1980’s . The development of solving algorithms and mathematical proofs has been developed contemporarily (see  and references therein).

During the last two decades, natural sciences profited a lot by the increasing computer power of modern machines. In nonlinear and statistical sciences, old techniques could be exploited to their full extent and new methods have been developed. The identification of nonlinearities is one of the important topics in data analysis of. e.g., chaotic or pattern forming systems . This in turn is feasible by statistical tools.

Careful application of algorithms and thorough interpretation of results leads to a direct way to obtain models from data. Hereafter, the basic ideas of the techniques to solve for a generalized additive model are explained and an application to spatio-temporal data is given. The intention is not to give a full overview of statistical methods, nor to report an experiment in detail. Rather, it is shown how a consistent procedure to infer models from data can be used as a feedback to theory which in turn can motivate experimental designs to yield new data and so on.

The paper is structured as follows: In Sec. 2 the algorithmic solutions for generalized additive models are pointed out. This is done by the example of ACE, the Alternating Conditional Expectation value algorithm, which uses the backfitting technique. This presentation is given some room; even though nothing new for statisticians, technical points need some explanation to be understood. Sec. 3 provides more information to be used in the context of spatio-temporal data analysis. In Sec. 4 problems arising from preprocessing of data are discussed along with the pattern formation example of vortex ripple dynamics. The paper concludes with a short discussion in Sec. 5.

## 2 Generalized additive models and backfitting

The standard tool of data analysis for physicists and engineers is the multivariate linear regression, where one fits the model

$$U_0=C+\underset{i=1}{\overset{M}{}}\alpha _iU_i+ϵ$$
(1)
to experimental data. The symbols $`ϵ`$ and $`U_i`$ (upper case letters), $`i=0,\mathrm{},M`$ denote a random process with given distribution, $`C`$ is a constant and the $`\alpha _i`$ are model parameters. One measures a realization $`u_i`$ (lower case letters), $`i=0,\mathrm{},M`$ of the random process and finds estimates $`\widehat{\alpha _i}`$ for the parameters, e.g., by least squares fitting . Throughout the text, estimations are denoted by a hat, only where ambiguity exists, otherwise additional symbols are omitted.

There are many ways of generalization. One often used way is the generalized linear regression model,

$$U_0=C+\underset{i=1}{\overset{M}{}}\alpha _if_i(U_i)+ϵ.$$
(2)
In this model, the data analyst chooses beforehand the functions $`f_i`$ due to some prior knowledge about the system or simply by guessing. Then, a linear regression procedure yields estimates for the parameters $`\alpha _i`$. So, the model is expressed in parametric form by the estimates $`\widehat{\alpha _i}`$. For details about linear regression, see, e.g., .

If no obvious choice for a function is at hand, it is desirable to fit the nonparametric model

$$U_0=C+\underset{i=1}{\overset{M}{}}f_i(U_i)+ϵ.$$
(3)
Given some data, estimates $`\widehat{f_i}`$ must result. With this kind of modeling no prior knowledge about functional dependences is needed as input; one can find very general functional forms. A disadvantage is the absence of even more general terms, like $`f_i(U_i,U_j)`$, $`ji`$. The reason for this lack is a practical one: as will become clear below, any implementation finding functions of $`d`$ variables has to compete with the curse of dimensionality. For the generalized additive model (3), estimates for the functions $`f_i`$ are found by the backfitting algorithm . This is an iterative procedure, working by the following rules:

1. Initialization: $`C=E(U_0)`$, $`f_i=f_i^0`$, $`i=1,\mathrm{},M`$.
2. Iteration: for $`i=1,\mathrm{},M`$ calculate

$$f_i=E\left(U_0C\underset{ki}{}f_k(U_k)|U_k\right)$$

until convergence,

where $`E()`$ denotes the expectation value and the $`f_i^0`$ are some appropriate initial settings. In applications, one has to obtain the functions from data as realizations of the generating process. Then, this operator has to be replaced by its estimator. In each iteration step we estimate the functions by $`\widehat{f_i}(U_i)=S_i(|U_i)`$, a smoothing operator which returns a function dependent on $`U_i`$.

The discussion of properties and choice of the smoother is a subtle issue and lies far beyond the scope of this paper, a detailed discussion is given in . The examples presented in this article have been obtained using the simple running mean smoother (moving average). This is not the best choice for many problems, e.g., points at the boundaries have to be considered with care or must be cut - throughout this article, there have been enough data to afford for the luxury to cut the boundaries, in spatio-temporal systems there are often many thousand of data points. Smoothing splines are good smoothers in many situations due to their differentiability properties . In any case, a careful check should be undertaken for a concrete set of data.

The iteration works by adjusting only for one function, subtracting all the others from the estimation of $`U_0`$, i.e. one smoothes the partial residual $`(U_0_kf_k)`$ against $`U_i`$. In spatiotemporal data analysis, on is faced with the task of finding equations of motion from data. These can be coupled map lattices, a set of ordinary coupled differential equations, partial differential equations or even more complicated models. In many of these formulae, the lhs. is linear (e.g. a time derivative) and Eq. (3) constitutes an appropriate additive, nonlinear model. If the measured data are nonlinear transforms through a measuring function, it can be desirable to model a nonlinear lhs, too:

$$f_0(U_0)=C+\underset{i=1}{\overset{M}{}}f_i(U_i)+ϵ.$$
(4)
The constant can be absorbed w.l.o.g. into the functions $`f_i`$. Now, the model is more symmetrical and indeed, in physics, often situations without a clear predictor (lhs) - response (rhs) structure are found. A typical example is pattern formation where one is interested in the nonlinear interaction of stationary states. A priori it is not clear which variables constitute the state space nor is it always possible to measure them. Data analysis provides then an estimation of the model as a nonlinear, possibly noninvertible, transformation of the measured data. As an example, take the relation $`U_0=G(U_1+U_2)`$ with $`G`$ some nonlinear invertible function. The inverted equation $`G^1(U_0)=U_1+U_2`$ cannot be found with (3), but with the model (5). An example for this case is given below.

An algorithm to solve (4) is ACE, (the Alternating Conditional Expectation value algorithm) , it works by minimizing the squared error

$$E\left[f_0(U_0)\underset{i=1}{\overset{M}{}}f_i(U_i)\right]^2.$$
(5)
The method shall be explained stepping from the one-dimensional estimation over two dimensional to the M-dimensional problem.

#### 1) $`U_0=f_1(U_1)+ϵ`$.

For this simple model, one obtains from Eq. (5)

$$E\left[U_0f_1(U_1)\right]^2=min,$$
(6)
where minimization is achieved by variation of $`f_1`$ in the space of measurable functions . After some elementary steps one finds the solution

$$f_1(U_1)=E(U_0|U_1).$$
(7)
For applications, again a smoother is used as estimator of the expectation value operator. The result of the smoothing on a set of numerically produced data is shown for the example $`U_0=U_1^2`$ in Fig. 1. The data has been generated by first drawing 10 000 equally distributed random numbers for $`u_1`$ in the interval $`(3,3)`$. These numbers have been squared to yield the values for $`u_0`$, after this Gaussian noise $`N(0,0.1^2)`$ has been added. The bandwidth of the running mean smoother has been set here and below to 200 points. In Fig. 1, the scatterplot of the points $`(u_1,u_0)`$ is shown (grey) together with the estimate for the function $`f_1`$, found by Eq. (7) (black, thick line). As indication for the estimation error the pointwise standard deviation has been calculated from the estimated function and the residuals (upper and lower black lines).

#### 2) $`f_0(U_0)=f_1(U_1)+ϵ`$.

In this case, the solution is found iteratively:

1. Initialization: $`f_0(U_0)=\left[U_0E(U_0)\right]/\sqrt{var(U_0)}`$.
2. Computation of $`f_1(U_1)=E[f_0(U_0)|U_1]`$
3. Computation of $`\stackrel{~}{f_0}(U_0)=E[f_1(U_1)|U_0]`$
4. Normalization: $`f_0(U_0)=\stackrel{~}{f_0}(U_0)/\sqrt{var(\stackrel{~}{f_0}(U_0))}`$.
5. Iteration of 2)-4) until convergence.

A smoother with nonzero bandwidth is contracting and thus the trivial solution is excluded by normalization (point 4) with the variance $`var(\stackrel{~}{f_0}(U_0))`$.

In order to show how ACE handles non-invertible relations, data ($`A\mathrm{cos}\theta `$, $`A\mathrm{sin}\theta `$) have been generated, with $`\theta =\phi +ϵ_\phi `$ and $`\phi `$ an equally distributed random number in $`[0,2\pi )`$ and $`A=1+ϵ_A`$. The noise $`ϵ_\phi ,ϵ_A`$ is Gaussian distributed, with $`N(0,0.05^2)`$. This corresponds to the noisy measurement of an amplitude and a phase with a nonlinear, noninvertible measurement transformation.

Without noise, one has $`U_0=\mathrm{cos}(\phi )`$, $`U_1=\mathrm{sin}(\phi )`$ and thus $`U_0^2=1U_1^2`$. With noise, the relation reads to first order $`U_0^22ϵ_AU_0=1U_1^2+2ϵ_AU_1`$. This means that multiplicative noise enters in predictor and response, touching the errors in variables problem which is relevant for many applications. The strong, multiplicative transformation of the noise is mirrored directly in the results.

In Fig. 2 a) the data are displayed in a scatterplot, in Fig. 2 b) and c), the resulting functions are shown on top of the respective residuals together with twice the pointwise standard deviation as an error indicator. Due to the nonlinear noise transformation the functions are distorted. The pointwise standard deviation is less significant as error indicator due to the very asymmetric local distribution of values, to be seen in the upper and lower curves of Fig. 2 b) and c). In this case asymmetric measures should be calculated. The above example underlines a conceptual problem when one is faced with errors in the measurement variables. The general treatment of this problem is not simple, for more details see  and references therein.

Even though the precise estimation of the functional dependence is not possible in the above case, one finds an approximation which would not be so bad for many experiments. Please note that linear tools would not be able to find this relation easily, rather one should undergo a trial-and-error procedure until a guess for the right form of the functions is found. Some problems with a non-Gaussian distribution, however, remain in linear regression, too.

#### 3)$`f_0(U_0)=_{i=1}^Mf_i(U_i)+ϵ`$

This equation defines a (noisy) hypersurface in $`M`$-dimensional space (in two dimensions, one finds a line). The functions $`f_i`$, fall on a hyperplane: the algorithm performs a transformation of the points $`\{U_1,\mathrm{},U_M\}`$ to the points $`\{f_1,\mathrm{},f_M\}`$ such that the best plane in the mean squared sense is found.

The modification in the algorithm has to be done in the steps 2. and 3. above. They are replaced by

1. Computation of $`f_i(U_i)=E[f_0(U_0)_{ki}f_k(U_k)|U_i]`$ by backfitting.
2. Computation of $`\stackrel{~}{f_0}(U_0)=E[_if_i(U_i)|U_0]`$

This shows that the ACE algorithm is based on backfitting. To demonstrate a higher-dimensional example, the relation $`U_0=\mathrm{exp}(U_1+U_2^2)`$ has been chosen. The variables have been generated using equally distributed values for $`U_1,U_2`$ in the interval $`(2,2)`$. From the resulting numbers, $`u_0`$ has been calculated, thus the smallest possible value of $`u_0`$ is $`\mathrm{exp}(2)0.14`$, the largest is $`\mathrm{exp}(6)400`$. After this procedure (in contrast to the previous example), Gaussian, $`N(0,0.2^2)`$, noise has been added to $`u_i`$, $`i=0,1,2`$. This avoids the problems with transformation of measurement noise, discussed above. But at this place, another interesting question for the modeler is raised instead: what happens if one tries to include some additionally measured variable in the model?

As an example, one can imagine a measurement of four variables, where it is not clear a priori which variables are needed in a minimal model. As a first guess, a generalized additive model $`f_0(U_0)=f_1(U_1)+f_2(U_2)+f_2(U_3)`$ will be tested. This situation shall be now related to the data generated above: an additional, equally distributed variable in the interval $`(0,100)`$, uncorrelated to $`U_0,U_1,U_2`$, is generated and fed as fourth variable into ACE.

The three-variable additive model is $`f_0(U_0)=f_1(U_1)+f_2(U_2)`$ and one expects ACE to find the estimates $`\widehat{f_0}(U_0)=\mathrm{ln}(U_0),\widehat{f_1}(U_1)=U_1,\widehat{f_2}(U_2)=U_2^2`$ due to the additive structure. For the four dimensional model, one expects an unchanged result for $`\widehat{f_0},\widehat{f_1},\widehat{f_2}`$ plus a zero function $`\widehat{f_3}(U_3)=0`$ in the additional variable.

In Fig. 3, the 3D scatterplot of the generated relevant data is shown by the points $`(\mathrm{ln}(U_0),U_1,U_2)`$, the $`U_0`$-axis is logarithmic. A two dimensional representation is not able to display the full relation, if one plots, e.g., $`U_0`$ vs. $`U_1`$ (Fig. 4), no functional dependence is cognizable. The transformations $`(\widehat{f_0},\widehat{f_1},\widehat{f_2},\widehat{f_3})`$ found by ACE are shown in Fig. 5 as functions of their arguments. The function $`f_0`$ is found to be logarithmic, $`f_0=0.5\mathrm{ln}(u_0)+c`$, $`f_1`$ is linear with a slope of $`0.5`$, and $`f_2=au_2^2+c`$ is quadratic also with a factor $`a=0.5`$. The factor $`0.5`$ reflects possible ambiguities in the result, as well as the addition of the constant, see Sec. 3.2.

After multiplying the results by $`2`$ and adding a constant to $`f_0`$ and $`f_2`$ one can speak of a good coincidence of expected and found result. The additional variable $`u_3`$ has no dependence on the others and yields $`f_3=O(10^2)0`$. Thus, the method is able to identify uncorrelated terms. A comparison of runs with and without the additional third term on the rhs showed that the functions are virtually unchanged, confirming the stability of the algorithm. Despite the seemingly convincing result it must be noted that results can be misleading if there exist correlations between an additional term and “true” variables. If the task is to find a minimal model, one has to consider this effect. The connection of correlation with ACE is discussed in Sec.3.3.

## 3 Some further topics

### 3.1 Smoothers and the bias-variance dilemma

As mentioned above, the choice of the smoother can be delicate. Typically, to each smoother belongs a smoothing parameter which characterizes its bandwidth (window width). For the running mean smoother, one obtains for problem (7)

$$\widehat{f}_1(U_{1,l})=1/N_l\underset{jB(U_{0,l})}{}U_{0,j},$$
(8)
where the hat denotes estimation and $`B(U_{0,l})`$ is a neighborhood of the point $`U_{0,l}`$. Increasing bandwidth enlarges the neighborhood and thus the number $`N_l`$ of data points in $`B`$. The respective expectation and variance are

$`E(\widehat{f}_1(U_{1,l}))`$ $`=`$ $`1/N_l{\displaystyle \underset{jB(U_{0,l})}{}}f(U_{0,j}),`$ (9)
$`var(\widehat{f}_1(U_{1,l}))`$ $`=`$ $`\sigma ^2/N_l,`$ (10)
where $`\sigma _l=\sigma `$ is assumed for simplicity. Thus, the variance shrinks with increasing bandwidth; the bias, however, increases because more terms different from $`f(U_{0,l})`$ are involved. If one wants to find the optimum choice of the bandwidth in order to have neither large bias nor large variance, one uses measures like the mean-squared error or the average predictive squared error; these are evaluated using cross-validation and yield a criterion for the bandwidth to be chosen. For more information it is referred to the literature .

The iteration procedure should run until convergence. In fact, convergence is not guaranteed for all classes of smoothers, asymptotic properties and convergence is discussed in a rigorous way in . Intuitively, it is clear that the iteration should stop when the errors on the functions fall below the noise $`ϵ`$ in the model (4). Then it is a question how to find properly the respective errors. There are different ways to estimate the error consistently, depending on the smoother and the error models used. Locally, one can choose twice the pointwise standard error as a characterization, global confidence criteria are harder to derive , an elegant and modern approach is given by the Bayesian formulation of the backfitting algorithm . In the case of errors in variables, care has to be taken to calculate the error correctly, a Bayesian approach is appealing in that case. In the examples, the most naive criterion, i.e., twice the pointwise standard error has been used.

### 3.2 Peculiarities of ACE

The ACE algorithm can be related to other algorithms, namely canonical variates and alternating least squares methods. Special properties of all three methods are discussed in , peculiarities of ACE are discussed in , here a few items shall be given without a claim for completeness.
i) ACE is symmetric in the predictor and the response, which is quite unusual for a regression tool, but suits well several setups in physical applications.
ii) Transformations of the variables are not reproduced by ACE. Take, e.g., the case $`U_0=f_1(U_1)+ϵ`$. If one applies a transformation to each of the sides, ACE does not necessarily find $`f_1`$ again. This reflects parts of the nature of the inverse problem, which is in general not uniquely solvable. Especially, adding a constant to either of the sides in the model is ambiguous in ACE.
iii) Different distributions in the model constituents can “distort” the resulting functions. E.g., it may happen that one term in a model with two variables is Gaussian distributed and another one is equally distributed in some interval. Then, the Gaussian variable is bent at the ends (cf. Fig. 5).

### 3.3 Maximal Correlation

The ACE algorithm can be regarded as a regression tool, but there exists another possible interpretation. Instead of solving the least squares error minimization problem

$$E[f_0(U_0)\underset{i=1}{\overset{M}{}}f_i(U_i)]^2=min$$
(11)
one can reformulate the above to

$$\mathrm{\Psi }=corr(f_0(U_0),\underset{i=1}{\overset{M}{}}f_i(U_i))=max,$$
(12)
where $`corr`$ denotes the correlation function. This reflects the principle of the maximal correlation . The functions that fulfill Eq. (12) are called optimal (in the sense of correlation) transformations, found by the ACE procedure. To have a global measure for the importance of a single term $`f_i(U_i)`$ , there are several possibilities, e.g., one can choose the overall variance, $`var(f_i)`$, but this does not reflect correlation with other terms and a noise term with large variance would be judged to be important. Here, the author suggests to use Eq. (12) in a symmetric way,

$$\mathrm{\Psi }_i=\left|\frac{E(f_i_{ji}f_j)}{\sqrt{E(f_i^2)E(_{ji}f_j)^2}}\right|,$$
(13)
to characterize the importance of a single term. The use of this measure is demonstrated for the data from the example used for Fig. 5. The values $`\mathrm{\Psi }_0=0.9779`$, $`\mathrm{\Psi }_1=0.9272`$, $`\mathrm{\Psi }_2=0.9690`$ and $`\mathrm{\Psi }_3=0.0426`$ are found, confirming again the unimportance of the additional uncorrelated variable, this time by the global measure $`\mathrm{\Psi }_3`$.

If however, the additional variable is correlated, e.g., by passively following another variable, a high correlation will follow, yielding a model with more components. In the case of spatiotemporal modeling one has to be specially careful, because one is interested in minimal (in the number of terms) models, but dynamical dependencies can result easily in correlations, like in the case of slaved variables.

### 3.4 Generalizations

One can generalize the procedure to include more complex, non-additive couplings, e.g. in the three dimensional case, one considers the model

$$U_0=F(U_1,U_2)+ϵ,$$
(14)
with $`F`$ some function, describing a two dimensional surface in a three dimensional space spanned by $`U_0,U_1,U_2`$. Like for Eq. (7), one finds

$$F(U_1,U_2)=E(U_0|U_1,U_2).$$
(15)
Going to higher dimensions requires in general more data due to the curse of dimensionality, additionally one has to decide again which estimator to use for correct results.

As an example, the relation $`U_0=U_1U_2`$ (Fig. 6a) is used. Ten thousand equally distributed data points in the range $`(1,1)`$ for $`U_1,U_2`$ have been generated, from these $`U_0`$ is obtained by multiplication. Gaussian, $`N(0,0.1^2`$, noise has been added afterwards to each of the three variables as noise realizations. The estimated function coincides well with the expectation (Fig. 6b). Note that in this case ACE yields as a result $`\mathrm{log}U_0=\mathrm{log}U_1+\mathrm{log}U_2`$ because of the additive structure of the model (5)

For this model, many applications with quasi-two-dimensional setting can be found, e.g., in geological and hydrodynamical systems (e.g. the two dimensional Euler equations).

## 4 Application to spatiotemporal data

Typical spatiotemporal models of interest for a physicist are coupled map (CML) systems , coupled ordinary differential equations (ODE’s) or partial differential equations (PDE’s). In all these models, derivatives or finite differences are involved. Here, the measured quantity shall consist of a space-time record, given by a movie, i.e. a series of consecutive pictures. From the time series in one point one builds finite differences in time, from the pictures one calculates differences in space. This procedure is numerically not trivial: addition (or subtraction, respectively) is numerically bad conditioned and can yield extinction of digits, measurement noise is then shifted to the first digits. Basically, one must use the apparatus known from numerical integration  to use correct sampling in space and time. In the case of a periodic setup, accuracy loss can be limited by using spectral methods .

As an example, let us take the Complex Ginzburg-Landau Equation. It describes a well pattern formating system above onset :

$$_tA=(1+ic_1)\mathrm{\Delta }A+(1ic_3)f_1(|A|^2)A.$$
(16)
The function $`f_1`$ is a complex nonlinear function. The task of data analysis -supposed that $`A`$ can be measured- is to determine the nonlinearity $`f_1`$. To apply the ideas of backfitting we put $`U_0(x,t)=_tA/A`$, $`U_1(x,t)=\mathrm{\Delta }A/A`$, $`U_2(x,t)=|A|^2`$. Real and imaginary part must be treated separately. Using the iteration procedure, one finally obtains estimates for the parameters $`c_1`$, $`c_3`$ and the full function $`f_1`$. In a slightly more complicated setup, this has been performed in . In fact, for this example the numerically critical step is rather the correct evaluation of derivatives than the application of backfitting.

If one wants to find a model, but there are no physical arguments which derivatives should occur or it is even unclear if the measured quantities are state variables or have undergone some measurement transformations, one can still try to find a generalized additive model. In that case, it can be worth to fit a discrete model (CML), that represents the system dynamics, avoiding a good part of the trouble with finite differences (one can choose mapping units big enough to avoid numerical problems).

A recently presented experiment considers underwater vortex ripples . The ripple formation is driven by a periodical motion of water over a sand bed with given amplitude and frequency. A picture of the ripples -viewed from the side- is given in Fig. 7 a), the resulting space time plot is shown in Fig. 7 b). In this paper, experimental details are not given, for those it is referred to the literature . But the treatment of the data and subsequent results are presented to demonstrate application of the nonparametric regression.

The model is formed by coarse graining of space and time. In time, one uses a mapping from period to period, in space, ripples are characterized by their length $`\lambda `$. The interaction between ripples is modeled by a nonlinear interaction function which must be concave by stability analysis. The model is obtained for ripples close to the final length; without discussing further the physical mechanisms, it reads:

$$\lambda _i(T+1)\lambda _i(T)=f(\lambda _{i1(T)})+f(\lambda _i(T))f(\lambda _{i+1}(T))+ϵ.$$
(17)
The subscript $`i`$ denotes the space index.

The process of data analysis involved the following steps:
i) calculating the ripple length, this has been done by fitting triangles to the raw profiles and calculating the distances from crest to crest (cf. Fig. 7. The accuracy amounted to two digits.
ii) Calculation of the rhs of Eq. (17), this has been done by first calculating the difference, the result has been very noisy due to the numerical extinction of digits. To compensate for this effect, the differences have been smoothed and small numbers have been discarded. Please note that no “unsuitable” points have been excluded (a common, but often wrong technique for outliers), rather it is a numerical requirement. In any case, the model has been developed from considerations of long ripples not too far from the final length and thus very small ripples are not expected to follow the model dynamics.
iii) To apply a backfitting algorithm, one has to identify the variables $`U_i`$ as they have been denoted in the previous sections. We put formally $`\mathrm{\Delta }_t\lambda =\lambda _i(T+1)\lambda _i(T)=U_0`$, $`\lambda _{i1}(T)=U_1`$, $`\lambda _i(T)=U_2`$, and $`\lambda _{i+1}(T)=U_3`$. The lhs shall be linear, so the model (3) has been used. For this example, the functions are denoted by $`f_1=f_l`$, $`f_2=f_c`$, $`f_3=f_r`$, for the left, center and right one, respectively. In the present context this allows easier interpretation in terms of the spatial structure. The data under consideration have been obtained by nine experimental runs, each between 1000 and 6000 data points.

The experimental runs have been stationary as far as experimental technique is concerned. Each data point $`(\mathrm{\Delta }_t\lambda _i(T),\lambda _{i1}(T),\lambda _i(T),\lambda _{i+1}(T))`$ can be considered as an independent realization of the underlying process. Thus, the data from all nine measurements have been analyzed together. The result is shown in Fig. 8. For each function, the importance criterion is calculated with $`\mathrm{\Psi }_0=0.732`$, $`\mathrm{\Psi }_1=0.560`$, $`\mathrm{\Psi }_2=0.8604`$ and $`\mathrm{\Psi }_3=0.6376`$. From these numbers one obtains basically two informations: 1) There is a large part of variance which cannot be modeled. This can have two reasons, the first is measurement noise, which has been reduced as much as possible, the second is that parts of the dynamical equations (higher order interactions or similar) are not included in the model. A priori it is not easy (or even impossible) to distinguish which is the reason, again the error in variables problem has to be treated with care. 2) The growth of a ripple depends mainly on the length of this ripple, but the neighbors cannot be neglected. This is seen by the lower correlations for the neighboring functions.

As a basic test, the result functions must show the expected concave shape. This is perfectly true in the range the model should hold. There are, however differences between the neighbor functions $`f_l`$, $`f_r`$ and the center one, $`f_c`$. Mainly one notes a very large standard deviation at certain regions of $`f_l,f_r`$ and a systematic deviation, a bend, at small ripple lengths. At this point the data analysis result gives feedback to the experiment and the theoretical modeling. It must be clarified if these differences have a physical background or if they are artificial effects due to the measurement procedure. This is ongoing work.

In Fig. 9 the left, center and right function are plotted together. Encouragingly, with respect to the established model, the functions fall close together except for small ripple lengths. The model has been developed in the region close to the final length, thus the result confirms the model. For small lengths, the result should be considered as input for further modeling. In the case of small ripple lengths, the deviations of the side functions to the center indicate that in the small ripple length region, a difference between center and sides exists. If one averages over left, center and right function, assuming that all three functions are equal and differences arise from measurement errors one obtains the thick black line in Fig. 9.

The model has been obtained as a great simplification of the highly complex granular-fluid system. Given additionally a certain amount of measurement and numerical errors from data processing, one can say that the model is well confirmed by direct, nonparametric data analysis. Further activity aims at a more general modeling in the sense of Eq. (14) to check if there are differences in the interaction functions of the sides and the center which is theoretically very improbable. The direct justification of the model (17) is a success not only for the modeler but as well for statistical data analysis in physical experiments.

## 5 Discussion and conclusion

Some of the facts about generalized additive models, backfitting algorithms and ACE have been presented. With a step by step explanation the basic ideas of algorithms and some technical details have been discussed. In the examples, problems with nonlinear transformation of data have been pointed out and the effect of additional, uncorrelated variables has been investigated. With data from an example for a complex pattern forming system, it has been shown how to apply the algorithms to spatio-temporal data and explained the main difficulties for applications. With data analysis a simple model has been validated, which has been developed by physical arguments and coarse graining of the highly complex system of underwater vortex ripples.

Generalized additive models occur often in physical and engineering applications. The idea to obtain models directly from data can be realized with the techniques developed in statistical science. The focus of the presented methods lies in the possibility to determine nonlinearities. Results from the nonparametric data analysis help the theoretician to improve analytical models as well as the applied physicist or engineer who just needs a working model to represent the system dynamics. These aspects can be closed to a circuit - measurement - data analysis - theoretical modeling - which can be iterated until a consistent model for the system under consideration is obtained. It can be expected that the use of nonparametric methods is helpful as well in situations where little is known about the system. There, working models can be inferred without further theoretical justification but as tools

While in this article additive models have been, there exist techniques to solve for either additive or multiplicative models, e.g. the marginal integration  or alternatively Breimans $`\mathrm{\Pi }`$-method  (without a claim for completeness). Implementations of the backfitting procedure or the ACE algorithm exist for the S-PLUS language or the public domain clone R-PLUS, or in the worldwide web .

Statistical approaches are promising tools for the analysis of complex systems, since nonlinearities are automatically fetched. Bayesian formulations promise to yield better error models and another, future application can be stochastic modeling, finding a representation of the residuals in form of a noise term with measured probability distribution.

## 6 Acknowledgements

I thank K. Andersen, J. Kurths, A. Pikovsky and A. Politi for support and discussion, a special thanks to H. Voss who led my attention to this kind of problems. F. Schmidt has provided computer wisdom, L. Barker-Exp. was a constant source of inspiration. I acknowledge support by the DFG (German research foundation), project number PI 200/7-01.


text:
# Inner density waves: the need for two-dimensional spectroscopy

## 1. Introduction

High resolution photometry shows the presence of non-axisymmetric structures in the inner parts of galaxies: spirals, bars and lopsided luminosity distributions. These are assumed to be small versions of the large-scale density waves observed in disk galaxies. although the involved formation and evolution processes may be significantly different (see e.g. Emsellem 2002 and references therein). We still do not know the respective roles of these modes in the evolution of the central kiloparsec, although possible links with the nuclear activity have been searched for. More importantly, the dynamical status of these photometric structures has not been properly assessed: we still don’t really know if these are true kinematic density waves, or transient density perturbations. Hence the need for detailed spectroscopic information.

I have discussed some features of inner density waves in a recent paper (Emsellem 2002). I wish here to further illustrate why classical long-slit spectroscopy is not well adapted to unveil the full complexity of the central regions of galaxies. These non axisymmetric structures are also strongly time dependent. This adds an extra dimension to the solution space to probe. We then (and always) need all the information we can gather to properly constrain the morphology and dynamics of these systems.

## 2. Where should we put the slit?

In this Section, I give two prototypical examples where the use of long-slit spectrography can lead to a restricted and biased interpretation of the observed structures.

### 2.1. A biased view

One of the best studied nearby nucleus belongs to the spiral galaxy M 31. With a distance of about 0.7 Mpc ($`3.4`$ pc/arcsec), it is the closest disk galaxy for which we have a rather privileged view of its central structures. HST photometry (Lauer et al. 1993) has revealed a double nucleus, the brightest peak (P1) being offset by about $`1.8`$ pc from the centre of the outer isophotes (P2). Ground-based stellar kinematics (Bacon et al. 1994) showed that the velocity field, at a spatial resolution of 0$`.^{\prime \prime }`$9, was nearly symmetric with respect to P2. The velocity dispersion however peaked on the side opposite to P1, with respect to P2. Subsequent ground-based long-slit data (Kormendy & Bender 1999, 0$`.^{\prime \prime }`$64 FWHM) confirmed the offset of the dispersion peak from the assumed velocity centre. Statler et al. (1999) then presented high resolution FOC long-slit data taken along the P1-P2 axis: although the signal-to-noise of these data was rather low, they were able to uncover a strong asymmetry in the stellar velocity profile, and a high central dispersion value.

We have recently published new high resolution data of the nucleus of M 31, obtained with the adaptive optics assisted integral field spectrograph OASIS (Bacon et al. 2001), as well as archival STIS kinematics. Taking advantage of the two-dimensional spatial coverage of the OASIS data, and helped by multi-band dithered HST/WFPC2 photometry, we have reconsidered the centering of the various published data sets. Using integral field spectroscopy, we can accurately define a reference centre a posteriori. Long-slit data of M 31 were alternatively taken along the P1-P2 axis, and along the kinematical major-axis of the nucleus, these axes differing by more than 10 degrees. Centres were also not necessarily well registered, and significant shifts between the different data sets were measured. These two point were important enough to confuse the interpretations, particularly during the modeling process (see Bacon et al. 2001 for details).

As emphasized in Fig. 1, it is thus not clear where a slit should be put. The fact that P1 does not lie on the nucleus kinematical major-axis must be taken into account when comparisons with models are performed. Only high resolution two-dimensional spectroscopy coupled with accurate photometry did allow a consistent picture to emerge.

### 2.2. The Sombrero: a double barred galaxy?

The Sombrero galaxy (M 104) is another famous and well studied nearby spiral. It is close to edge-on ($`i=84^{}`$), which tends to hide any departure from axisymmetry within its equatorial plane. M 104 has been (and is still) assumed to be an atypical (because of its huge spheroidal component) but axisymmetric early-type spiral. Long-slit spectroscopy has thus often been used to probe its major or minor-axis kinematics and line-strength gradients (e.g. Hes & Peletier 1993).

However, the presence of an outer prominent ring of HI, CO and dust, led Emsellem (1995) to suggest the existence of a tumbling potential, thus linking the ring to its Outer Lindblad Resonance. This was later confirmed by the discovery of an additional Ultra Harmonic ring at a radius of 43<sup>′′</sup>(Emsellem et al. 1996). But the real surprise came from the the TIGER integral field spectroscopic data (Fig. 2). The two-dimensional coverage allowed us to spatially deconvolve the full TIGER datacube (reaching an equivalent resolution of $`0.^{\prime \prime }5`$), something which would have been impossible using long-slit data.

The resulting emission line distribution shows the presence of a spiral-like structure, and the gas velocity map exhibits strong non-circular motions within the few central arcseconds (Emsellem & Ferruit 2000). We also detected the inner gaseous spiral in the narrow band WFPC2/HST (\[NII\]/H$`\alpha `$) imaging (see also Pogge et al. 2000), and a puzzling dust lane system using HST colour maps (Fig. 3). At a minimum inclination of 82, the deprojection stretching factor along the minor-axis is about 6. If we then assume that the detected gas lies within the equatorial plane of the galaxy, its distribution must therefore be significantly elongated. Hence we tentatively suggested the presence of an inner bar, the Sombrero galaxy being then a double barred spiral.

## 3. Time evolution

As mentioned in the introduction, density waves in the inner regions of galaxies imply a strongly time dependent morphology. This can for example be shown by numerical simulations with which we can take snapshots of observable quantities at different times.

### 3.1. Following a double bar system

Double bars are ideal cases to illustrate the issue of time dependency, as the two bars are generally thought not to share the same pattern speed. This implies a strongly evolving dynamics and morphology particularly at the interface region, as nicely emphasized by Maciejewski et al. (2002). A N-body$`+`$SPH simulation of a double barred galaxy (kindly provided by Daniel Friedli) is presented in Fig. 4. In this model, the inner bar rotates faster than its large-scale parent: this can easily be seen by tracing their position angle difference. The gas distribution at and within the Inner Lindblad Resonance ring exhibits a strongly evolving shape.

In real galaxies, these structures are projected on the sky. Strong departures from circular motions will have to be interpreted in the frame of a deprojected morphology. Only a full two-dimensional view can reconcile the observations with the models, and help constraining the true morphology and dynamics of the target.

### 3.2. More ingredients

The model of a double bar galaxy presented in the previous Section (Fig.4), did not include processes such as star formation or nuclear activity. When these ingredients are added, things get even more confused. Cases like NGC 1808 (Emsellem et al. 2001) are there to remind us that there are important couplings between various intervening processes which may lead to complex cycles (see e.g. Combes 2001).

I present here an example of such a complex case, namely the double barred galaxy NGC 470. It is a star forming inclined spiral, the photometry of which looks rather regular in the near-infrared (Friedli et al. 1995). We recently obtained OASIS integral field observations of the central region of NGC 470 (Fig. 5). The gaseous and stellar distribution in the optical are strongly asymmetric. This is partially due to the presence of dust on the NW part, but there seems to be a true residual lopsidedness in the central few arcseconds. Without these maps, it would have been difficult to determine where the kinematical centre lies.

## 4. Conclusions

Density waves seem to be ubiquitous in the central regions of galaxies, although we still do not have sufficient observational constraints to assess their existence. If they do exist, density waves certainly play an important role in redistributing the dissipative component and reshaping the inner parts of disk galaxies. Two-dimensional spectrographs are now becoming more common, and every large telescope will have its own. With their enhanced spatial resolution and the parallel increase in collecting aperture, these will be ideal instruments to probe the existence and role of inner density waves.

## References

Bacon, R., Emsellem, E., Combes, F., Copin, Y., Monnet, G., Martin, P., 2001, A&A 371, 409

Bacon, R., Emsellem, R., Monnet, G., Nieto, J.-L., A&A, 281, 691

Combes, F., 2001, INAOE, eds. D. Kunth & I. Aretxaga, astro-ph/0010570

Emsellem, E., 2002, in Disk of Galaxies, ASP Conference Series, Eds Athanassoula & Bosma, in press (astro-ph/0201176)

Emsellem, E., Bacon, R., Monnet, G., Poulain, P., 1996, A&A 312, 777

Emsellem, E., Greusard, D., Combes, F., Friedli, D., Leon, S., Pécontal, E., Wozniak, H., 2001, A&A 368, 52

Emsellem, E., Ferruit, P., 2000, A&A, 357, 111

Friedli, D., Wozniak, H., Rieke, M., Martinet, L., Bratschi, P., 1996, A&AS, 118, 461

Lauer, T. R., Faber, S. M., Groth, E. J., et al., 1993, AJ, 106, 1436

Pogge, R. W., Maoz, D., Ho, L. C., Eracleous, M., 2000, ApJ 532, 323

Statler, T. S., King, I. R., Crane, P., Jedrzejewski, R. I., 1999, AJ, 117, 894

Maciejewski, W., Teuben, P. J., Sparke, L. S., Stone, J. M., 2002, MNRAS, 329, 502

Hes, R., Peletier, R., F., 1993, A&A, 268, 539