\section{Introducci\'on}

En [3] desarrollamos la teor\'ia para poder resolver ecuaciones diferenciales de la forma
\begin{equation}
\phi^{\prime}=f_{1}(\phi)f_{2}(\phi)\cdots f_{n}(\phi)
\end{equation}
en t\'ermino de las soluciones de las ecuaciones m\'as simples $\phi^{\prime}=f_{i}(\phi)$.

Queremos hacer lo mismo con ecuaciones de diferencia de la forma

\begin{equation}
y_{n+1}=f_{1}(y_{n})f_{2}(y_{n})\cdots f_{n}(y_{n})
\end{equation}
esto es, expresar la soluci\'on $y_{n}$ en t\'ermino de las soluciones de las ecuaciones
$y_{n+1}=f_{i}(y_{n})$. Esto ser\'a muy \'util por ejemplo para resolver la ecuaci\'on de
diferencia log\'istica $y_{n+1}=\mu y_{n}(1-y_{n})$. Esta es la ecuaci\'on en diferencia no
lineal m\'as simple y hasta el momento no tiene una soluci\'on exacta.

Como el operador de derivaci\'on $\delta$ y el operador de diferencia forward $\vartriangle$
son ejemplos de operadores delta, entonces iremos un paso m\'a all\'a definiendo y resolviendo
sistemas din\'amicos uni-dimensionales tipo delta.

Este art\'iculo est\'a dividido de la siguiente forma. Primeramente establecemos conceptos
b\'asicos de sistemas din\'amicos unidimensionales sobre dominios de integridad y conceptos
y resultados de la teor\'ia de operadores delta y del c\'alculo umbral desarrollado por G.C. Rota.
Seguido definimos sistemas din\'amicos unidimensionales tipo delta sobre dominios de integridad y
establecemos algunos resultados b\'asicos sobre estos. Uno de los resultados importantes alcanzado
es que si fijamos la funci\'on $f$ del sistema y variamos los operadores delta obtenemos un grupo
de sistemas din\'amicos. Finalizamos este art\'iculo introduciendo algunos
sistemas din\'amicos tipo delta cuando delta es el operador de diferencia forward, el operador
de diferencia backward, el operador Abel y el operador Touchard.

En todo este artículo $R$ será un anillo de característica cero. Además, $\mathbb{N}$ será el monoide
de enteros no negativos y $\mathbb{N}_{1}$ será el semigrupo de enteros positivos.

\section{Preliminares}

Iniciamos este art\'iculo presentando la informaci\'on necesaria para construir la teor\'ia de
sistemas din\'amico tipo delta definidos sobre un dominio de integridad de caracter\'istica cero.
Primero son presentados los resultados b\'asicos de [3]. Necesitamos el siguiente isomorfismo de
anillos
\begin{equation}
\dhurw_{R}[[x]]\rightarrow\mathfrak{A}(\dhurw_{R}[[x]])\rightarrow\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])
\end{equation}
en donde $\dhurw_{R}[[x]]$ es el anillo de sucesiones de derivadas de funciones en $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$,
el anillo $\mathfrak{A}(\dhurw_{R}[[x]])$ es construido usando el operador aut\'onomo definido sobre
$\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Este operador es definido usando los polinomios de Bell, por lo tanto exhibe
su car\'acter combinatorial. Finalmente el anillo $\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])$ es el anillo
que contiene los semi flujos soluci\'on de las ecuaciones $\phi^{\prime}=f(\phi)$ para cada
$f\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Otros resultados que necesitaremos son obtenido del c\'alculo umbral de Rota,
ve\'ase [4], particularmente los operadores delta que son operadores lineales que generalizan a los operadores
de derivaci\'on.

\subsection{Ecuaciones diferenciales aut\'onomas de orden uno sobre dominios de integridad}

Denote $R$ un dominio de integridad de caracter\'istica cero y sea $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$ el anillo
de Hurwitz, ve\'ase [2], de las funciones generadoras exponenciales de las sucesiones definidas
sobre $R$ con suma y producto ordinarios de series de potencias. Denote $\delta$ la derivaci\'on
definida sobre $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$ y denote $\dhurw_{R}[[x]]$ el anillo expansi\'on de Hurwitz de
$\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$, esto es, el anillo de sucesiones de derivadas $(f(x),\delta f(x),\delta^{2}f(x),...)$ para cada $f(x)$ en $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$, donde la suma $+$ en $\dhurw_{R}[[x]]$ es definida
componente a componente y el producto $\ast$ es el producto Hurwitz de sucesiones dado por
\begin{equation}
\Delta f(x)\ast\Delta g(x)=\left(\sum_{k=0}^{n}\delta^{k}f(x)\delta^{n-k}g(x)\right)_{n\in\mathbb{N}}.
\end{equation}

Como $\Delta f+\Delta g=\Delta(f+g)$ y $\Delta f\ast\Delta g=\Delta(fg)$, entonces los anillos
$\dhurw_{R}[[x]]$ y $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$ son isomorfos.
El anillo $\dhurw_{R}[[x]]$ sirve de base para la definici\'on del operador aut\'onomo.
Sea $\mathop{\mbox{\textsl{H}}}\nolimits_{S}$ el anillo de sucesiones definidas sobre el anillo $S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. El
\textbf{operador aut\'onomo} $\mathfrak{A}$ es el mapa no lineal
$\mathfrak{A}:\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow\mathop{\mbox{\textsl{H}}}\nolimits_{S}$ definido por
\begin{equation}
\mathfrak{A}(\Delta f)=(A_{n}([\delta^{n-1}f]))_{n\in\mathbb{N}_{1}}
\end{equation}
con $[\delta^{n}f]=(f,\delta f,...,\delta^{n}f)$, en donde los $A_{n}$ son definidos
recursivamente por
\begin{eqnarray}
A_{1}([\delta^{0}f])&=&A_{1}([f])=f,\\
A_{n+1}([\delta^{n}f])&=&Y_{n}(A_{1}([\delta^{0}f]),A_{2}([\delta^{1}f])...,A_{n}([\delta^{n-1}f]);\delta^{1}f,...,\delta^{n}f),
\end{eqnarray}
$n\geq 1$.\\
Los polinomios $A_{n}$ en las indeterminadas $\delta^{0}f,\delta^{1}f,...,\delta^{n-1}f$
serán llamados \textbf{polinomios aut\'onomos}.
Denote $\mathfrak{A}(\dhurw_{R}[[x]])$ el conjunto de todas las sucesiones $\mathfrak{A}(\Delta f)$.
Definiendo la suma $\boxplus$ en $\mathfrak{A}(\dhurw_{R}[[x]])$ de la siguiente manera
\begin{equation}\label{def_bplus}
\mathfrak{A}(\Delta f)\boxplus\mathfrak{A}(\Delta g)=\mathfrak{A}(\Delta f)+\mathfrak{A}(\Delta g)+\left( H_{n}(f,g)\right)_{n\in\mathbb{N}_{1}}
\end{equation}
en donde $H_{1}(f,g)=0$ y
\begin{equation}\label{eqn_Hn}
H_{n+1}(f,g)=f\delta A_{n}([\delta^{(n-1)}g])+g\delta A_{n}([\delta^{(n-1)}f])+(f+g)\delta H_{n}(f,g)
\end{equation}
y el producto $\circledast$ en $\mathfrak{A}(\dhurw_{R}[[x]])$ de la siguiente manera
\begin{equation}\label{def_cdast}
\mathfrak{A}(\Delta f)\circledast\mathfrak{A}(\Delta g)=\left(\sum_{l=1}^{n}\mathcal{A}_{l}(f)\mathcal{A}_{n-l+1}(g)\right)_{n\in\mathbb{N}_{1}}
\end{equation}
en donde
\begin{equation}
\mathcal{A}_{l}(f)=\sum_{\vert p(n)\vert=l}\alpha_{\vert p(n)\vert} A_{j_{1}}([\delta^{j_{1}-1}f])A_{j_{2}}([\delta^{j_{2}-1}f])\cdots A_{j_{r}}([\delta^{j_{r}-1}f])
\end{equation}
y $\alpha_{\vert p(n)\vert}$ son n\'umeros apropiados, entonces $(\mathfrak{A}(\dhurw_{R}[[x]]),\boxplus,\circledast)$ es un anillo conmutativo con unidad. Este anillo ser\'a llamado \textbf{anillo
aut\'onomo} .

Como $\mathfrak{A}(\Delta f)\boxplus\mathfrak{A}(\Delta g)=\mathfrak{A}(\Delta f+\Delta g)$ y $\mathfrak{A}(\Delta f)\circledast\mathfrak{A}(\Delta g)=\mathfrak{A}(\Delta f\ast\Delta g)$, entonces los anillos $\dhurw_{R}[[x]]$ y
$\mathfrak{A}(\dhurw_{R}[[x]])$ son isomorfos.

Ahora mostraremos la relaci\'on entre anillo aut\'onomo y sistemas din\'amicos. Definimos un
\textbf{sistema dinámico} sobre el anillo $R$ como la terna $(R,\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Phi)$
donde $R$ es el \textbf{conjunto de tiempos}, $\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el \textbf{espacio de fases} y
$\Phi$ es el \textbf{flujo}
\begin{equation}\label{flujo_diferencial}
\Phi(t,x,\Delta f(x))=x+\sum_{n=1}^{\infty}A_{n}([\delta^{(n-1)}f(x)])\dfrac{t^{n}}{n!}.
\end{equation}
esto es, $\Phi$ es el mapa definido
por $\Phi:R\times\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ donde
$S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.

El flujo $\Phi$ satisface
\begin{enumerate}
\item $\Phi(0,x,\Delta f(x))=x$,
\item $\Phi(t,\Phi(s,x,\Delta f(x)),\Delta f(x))=\Phi(t+s,x\Delta f(x))$,
\item $f(x)\delta_{x}\Phi(t,x,\Delta f(x))=\delta_{t}\Phi(t,x,\Delta f(x))=f(\Phi)$,
\item Para todo $a\in R$ se cumple que $\Phi(t,x,a\Delta^{}f(x))=\Phi(at,x,\Delta^{}f(x))$.
\end{enumerate}

Las propiedades 1. y 2. le dan a $\Phi$ estructura de grupo abeliano. Este grupo act\'ua sobre
el espacio de fases $\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$. La propiedad 4. dota a $\Phi$ de estructura de $R$-m\'odulo.
Ponga $\Phi_{t,r}(x)\equiv\Phi(t,x,r\Delta f(x))$. Entonces
\begin{equation}
\Phi_{t,r}(x)=\Phi_{rt,1}(x)=\Phi_{r,t}(x)=\Phi_{1,rt}(x)
\end{equation}
Luego podemos definir
\begin{equation}
\Phi_{R,1}(x)=\{\Phi_{t,1}(x):t\in R\}.
\end{equation}

Defina el mapa $\star:R\times\Phi_{R,1}\rightarrow\Phi_{R,1}$ por
$s\star\Phi_{t,1}=\Phi_{st,1}$. Entonces $\Phi_{R,1}$ es un $R$-módulo.

Por otro lado, sea $\rho_{t}$ el mapa llevando $\mathfrak{A}(\dhurw_{R}[[x]])$ a
$\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])$ definido por
\begin{equation}
\mathfrak{A}(\Delta f(x))\mapsto\sum_{n=1}^{\infty}A_{n}([\delta^{n-1}f(x)])\frac{t^{n}}{n!}
\end{equation}

Definimos un \textbf{semi sistema dinámico} como la terna $(R,\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Psi)$
donde $R$ es el conjunto de tiempos, $\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el espacio de fases, con
$S=\mathfrak{A}(\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]])$ y $\Psi$ es el \textbf{semi flujo}
\begin{equation}
\Psi(t,x,\Delta f(x))=\Phi(t,x,\Delta f(x))-x=\rho_{t}\mathfrak{A}(\Delta f(x))
\end{equation}
esto es, $\Psi$ es el mapa
$\Psi:R\times\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ .

Luego es posible extender las operaciones del anillo $\mathfrak{A}(\dhurw_{R}[[x]])$ al conjunto
$\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])$. Tome $\Delta f(x)$ y $\Delta g(x)$ de $\dhurw_{R}[[x]]$. Definimos
la suma $\boxplus$ y el producto $\circledast$ en $\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])$ así
\begin{eqnarray}
\rho_{t}\mathfrak{A}(\Delta f)\boxplus\rho_{t}\mathfrak{A}(\Delta g)&=&\sum_{n=1}^{\infty}\left(A_{n}([\delta^{n-1}f(x)])+A_{n}([\delta^{n-1}f(x)])+H_{n}(f,g)\right)\\
\rho_{t}\mathfrak{A}(\Delta f)\circledast\rho_{t}\mathfrak{A}(\Delta g)&=&\sum_{n=1}^{\infty}\left(\sum_{l=1}^{n}\mathcal{A}_{l,n}(f)\mathcal{A}_{n-l+1,n}(g)\right)\frac{t^{n}}{n!}
\end{eqnarray}

Entonces $\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])$ con la suma $\boxplus$ y el producto $\circledast$ es
un anillo conmutativo con unidades $0$ y $\rho_{t}\mathfrak{A}(\textbf{e}_{\Delta})=t$. Claramente
los anillos $\mathfrak{A}(\dhurw_{R}[[x]])$ y $\rho_{t}\mathfrak{A}(\dhurw_{R}[[x]])$ son isomorfos. Luego podemos
factorizar semi flujos.
Sea $f=g_{1}g_{2}\cdots g_{m}$ un producto de funciones en $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Entonces
\begin{equation}
\rho_{t}\mathfrak{A}(\Delta f)=\rho_{t}\mathfrak{A}(\Delta g_{1})\circledast\cdots\circledast\rho_{t}\mathfrak{A}(\Delta g_{m})
\end{equation}

implica que
\begin{equation}
\Phi(t,x,\Delta f(x))=x+(\Phi(t,x,\Delta g_{1}(x))-x)\circledast\cdots\circledast(\Phi(t,x,\Delta g_{m}(x))-x)
\end{equation}
De este modo podemos expresar el flujo de un sistema din\'amico en t\'erminos de producto de
semi flujos. Lo cual es muy \'util para encontrar soluciones de ecuaciones diferenciales en
t\'ermino de soluciones de ecuaciones m\'as simples.

\subsection{C\'alculo finito de operadores}

Sea $\mathbb{F}$ un campo de caracter\'istica cero y $\mathbb{F}[t]$ el $\mathbb{F}$-\'algebra de polinomios con coeficientes
en $\mathbb{F}$ en la variable $t$.

Una \textbf{sucesi\'on polinomial} $(p_{n}(t))_{n\in\mathbb{N}}$ en $\mathbb{F}[t]$ es un sucesi\'on de polinomios
tal que $p_{n}(t)$ tiene exactamente grado $n$. Una sucesi\'on polinomial se dice de
\textbf{tipo binomial} si se cumple
\begin{equation}
p_{n}(t+s)=\sum_{k=0}^{n}\binom{n}{k}p_{k}(t)p_{n-k}(s)
\end{equation}
para todo $n\in\mathbb{N}$.

Si $a$ es un elemento de $\mathbb{F}$, entonces $E^{a}$ es el operador shift enviando un polinomio $p(t)$
a $p(t+a)$. Se puede mostrar que
\begin{equation}
E^{a}=e^{a\delta}=\sum_{k=0}^{\infty}\frac{a^{k}}{k!}\delta^{k}
\end{equation}

Un operador $Q$ es \textbf{shift invariante} si $QE^{a}=E^{a}Q$ para cualquier elemento $a$ en $\mathbb{F}$.
El operador $\delta$ es claramente un operador shift invariante. Denote $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[\delta]]$
el anillo de Hurwitz de operadores shift invariantes.

Un operador $Q:\mathbb{F}[t]\rightarrow\mathbb{F}[t]$ es un \textbf{operador delta} si $Q$ es shift invariante y
$Qt=c$, donde $c$ es una constante no cero. Una sucesi\'on polinomial $(p_{n}(t))_{n\in\mathbb{N}}$
es un \textbf{conjunto base} para el operador delta $Q$ si $p_{0}(t)=1$, $p_{n}(0)=0$ para todo
$n\geq1$ y $Qp_{n}(t)=np_{n-1}(t)$.

Algunas propiedades de operadores delta son:
\begin{enumerate}
\item Si $p(t)$ es un polinomio de grado $n$, entonces $Qp(t)$ es un polinomio de grado $n-1$.
\item Sean $(p_{n}(t))_{n\in\mathbb{N}}$ una sucesi\'on de polinomios y $Q:\mathbb{F}[t]\rightarrow\mathbb{F}[t]$ un operador
definido por $Qp_{n}(t)=np_{n-1}(t)$. Entonces son equivalente
\begin{enumerate}
\item La sucesi\'on $(p_{n}(t))_{n\in\mathbb{N}}$ es de tipo binomial
\item El operador $Q$ es un operador delta
\item El operador $Q$ tiene la expansión formal
\begin{equation}
Q=\sum_{k=0}^{\infty}a_{k}\frac{\delta^{k}}{k!}
\end{equation}
es decir, $Q\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[\delta]]$.
\end{enumerate}
\end{enumerate}

Lo anterior significa entre otras cosas que las sucesiones polinomiales de tipo binomial son
caracterizadas por operadores delta. Los siguientes son resultados importantes en la teor\'ia
de operadores delta:

\textbf{Primer teorema de expansi\'on}. Sea $T$ un operador shift invariante y $Q$ un operador delta
con conjunto base $(p_{n}(x))$. Entonces
\begin{equation}
T=\sum_{k=0}^{\infty}\frac{[Tp_{k}(x)]_{x=0}}{k!}Q^{k}
\end{equation}

\textbf{Teorema de isomorfismo}. Sea $Q$ un operador delta con conjunto base $(q_{n}(x))_{n\in\mathbb{N}}$.
Entonces el mapa $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[\delta]]\rightarrow\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[t]]$ dado por
\begin{equation}
T\mapsto\sum_{k=0}^{\infty}[Tq_{k}(x)]_{x=0}\frac{t^{k}}{k!}
\end{equation}
es un isomorfismo de anillos.

Sea $Q$ un operador delta con conjunto base $p_{n}(t)$ y sea $Q=q(\delta)$. Sea $q^{-1}(t)$ la
serie de potencias formal inversa. Entonces
\begin{equation}\label{eqn_exp_inver}
\sum_{n=0}^{\infty}\frac{q_{n}(t)}{n!}u^{n}=e^{tq^{-1}(u)}.
\end{equation}

Un operador $L:\mathbb{F}[t]\rightarrow\mathbb{F}[t]$ es un \textbf{operador umbral} si existen dos sucesiones
bases $(p_{n}(t))_{n\in\mathbb{N}}$ y $(q_{n}(t))_{n\in\mathbb{N}}$ tal que $Lp_{n}(t)=q_{n}(t)$, para todo $n\in\mathbb{N}$.
Existe un \'unico operador lineal $L$ actuando sobre $\mathbb{F}[t]$ tal que $L(t^{n})=q_{n}(t)$. Este
operador $L$ es conocido como la \textbf{representaci\'on umbral} de la sucesi\'on $q_{n}(t)$.
Si
\begin{equation}
p_{n}(t)=\sum_{k=0}^{n}a_{n,k}t^{k},
\end{equation}
entonces la \textbf{composici\'on umbral} es la sucesi\'on $(r_{n}(t))_{n\in\mathbb{N}}$ definida por
\begin{equation}
r_{n}(t)=\sum_{k=0}^{n}a_{n,k}q_{k}(t)
\end{equation}
es decir $r_{n}(t)=Lp_{n}(t)$, donde $L$ es la representaci\'on umbral de $(q_{n}(t))_{n\in\mathbb{N}}$.
Usaremos la notaci\'on $r_{n}(t)=p_{n}(\underline{q}(t))$ para $r_{n}(t)=Lp_{n}(t)$. Luego
definimos la composici\'on $(p\circ q)(t)$ como la sucesi\'on $(p_{n}(\underline{q}(t)))_{n\in\mathbb{N}}$.
Denote $\mathcal{U}$ el conjunto de sucesiones de base. Entonces $\mathcal{U}$ es un grupo de Lie con
la composici\'on umbral $\circ$ como operaci\'on en donde dos sucesiones $(p_{n}(t))_{n\in\mathbb{N}}$
y $(q_{n}(t))_{n\in\mathbb{N}}$ son inversas si $p_{n}(\underline{q}(t))=t^{n}$. Si $p_{n}(t)$ y $q_{n}(t)$
son sucesiones bases con operadores delta $P=f(\delta)$ y $Q=g(\delta)$, entonces
$p_{n}(\underline{q}(t))$ es un conjunto base con operador delta $f(g(\delta))$. Luego
$(p_{n}(t))_{n\in\mathbb{N}}$ y $(q_{n}(t))_{n\in\mathbb{N}}$ son inversas si y s\'olo si $f(g(x))=t$.

\section{Sistemas Dinámicos tipo Delta}

Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$. Ahora sea $L$ la
representaci\'on umbral de los polinomios base $q$, esto es, $L(t^{n})=q_{n}(t)$. El operador delta
se comporta como una derivada cuando actua sobre polinomios. Lo que queremos hacer es definir
una nueva clase de ecuaciones diferenciales en donde en vez de usar una derivada ordinaria usaremos
operadores delta y luego encontrar analog\'ias en los resultados encontrados en sistemas din\'amicos
unidimensionales definidos sobre dominios.

\begin{theorem}
Denote $\Phi_{Q}(t,x,\Delta f(x))=L[\Phi(t,x,\Delta f(x))]$. Entonces
\begin{equation}
Q\Phi_{Q}(t,x,\Delta f(x))=f(\Phi_{Q}(t,x,\Delta f(x))).
\end{equation}
\end{theorem}
\begin{proof}
Usamos la definici\'on de $\Phi(t,x,\Delta f(x))$ y aplicamos $QL$. Tenemos
\begin{eqnarray*}
QL[\Phi(t,x,\Delta f(x))]&=&QL\left(x+\sum_{k=1}^{\infty}A_{k}([\delta^{k-1}f(x)])\frac{t^{k}}{k!}\right)\\
&=&xQL(1)+\sum_{k=1}^{\infty}A_{k}([\delta^{k-1}f(x)])\frac{QL(t^{k})}{k!}\\
&=&\sum_{k=1}^{\infty}A_{k}([\delta^{k-1}f(x)])\frac{kL(t^{k-1})}{k!}\\
&=&\sum_{k=0}^{\infty}A_{k+1}([\delta^{k}f(x)])\frac{L(t^{k})}{k!}\\
&=&L\left(\sum_{k=0}^{\infty}A_{k+1}([\delta^{k}f(x)])\frac{t^{k}}{k!}\right)\\
&=&L\left(f(x)\delta_{x}\Phi(t,x,\Delta f(x))\right)
\end{eqnarray*}
Por las propiedades de flujo la \'ultima l\'inea arriba da $L[f(\Phi(t,x,\Delta f(x)))]$. Finalmente
poniendo $\Phi_{Q}(t,x,\Delta f(x))=L[\Phi(t,x,\Delta f(x))]$ obtenemos el resultado deseado.
\end{proof}

Luego de este resultado definimos ecuaciones diferenciales aut\'onomas de orden uno con
derivaci\'on delta de este modo

\begin{definition}
Sea $Q$ un operador umbral. Tome $f(x)$ en $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Una ecuaci\'on diferencial aut\'onoma
uni-dimensional \textbf{tipo delta} es una ecuaci\'on de la forma
\begin{equation}
Q\Phi=f(\Phi).
\end{equation}
\end{definition}

Esta definici\'on generaliza a las ecuaciones que ya conocemos. Cuando $Q=\delta$ obtenemos
las ecuaciones diferenciales ordinarias. Si $Q$ es el operador de diferencia $\vartriangle$, entonces
obtenemos las ecuaciones en diferencias de orden uno. Ahora bien esta definición permite obtener
nuevos tipos de ecuaciones según el operador delta $Q$ que usemos. Si $Q=\triangledown$, obtenemos
\textbf{ecuaciones de diferencia backward}. Si $Q=E^{a}\delta$, entonces obtenemos
\textbf{ecuaciones de tipo Abel}. Si $Q=\log(I+\delta)$, entonces obtenemos \textbf{ecuaciones
de tipo Touchard}. Las soluciones a este tipo de soluciones ser\'an estudiadas m\'as adelante.
Ahora queremos encontrar una soluci\'on que generalice a las soluciones encontradas en secciones
anteriores.

Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ con representaci\'on umbral
$L$. Tome $\mathfrak{A}(\Delta f(x))$ en
$\mathfrak{A}(\dhurw_{R}[[x]])$ y denote $\rho_{q}$ el mapa definido por
\begin{equation}
\rho_{q}(\mathfrak{A}(\Delta f(x)))=\sum_{n=1}^{\infty}A_{k}([\delta^{n-1}f(x)])\frac{q_{n}(t)}{n!}.
\end{equation}
Luego $\rho_{q}\mathfrak{A}(\Delta f(x))=L(\rho_{t}\mathfrak{A}(\Delta f(x)))$ y $L$ es la representaci\'on
umbral de $\rho_{q}\mathfrak{A}(\Delta f(x))$.

\begin{definition}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$. Tome $\Delta f(x)$ y
$\Delta g(x)$ de $\dhurw_{R}[[x]]$. Definimos la suma $\boxplus_{q}$ y el
producto $\circledast_{q}$ en $\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]])$ así
\begin{eqnarray}
\rho_{q}\mathfrak{A}(\Delta f)\boxplus_{q}\rho_{q}\mathfrak{A}(\Delta g)&=&\sum_{n=1}^{\infty}(A_{n}([\delta^{n-1}f(x)])+A_{n}([\delta^{n-1}f(x)])\nonumber\\
&&+H_{n}(f,g))\frac{q_{n}(t)}{n!}
\end{eqnarray}
donde
$$H_{n+1}(f,g)=f\delta A_{n}([\delta^{n-1}g])+g\delta A_{n}([\delta^{n-1}f])+(f+g)\delta H_{n}(f,g)$$
y
\begin{eqnarray}
\rho_{q}\mathfrak{A}(\Delta f)\circledast_{q}\rho_{q}\mathfrak{A}(\Delta g)&=&\sum_{n=1}^{\infty}\left(\sum_{l=1}^{n}\mathcal{A}_{l,n}(f)\mathcal{A}_{n-l+1,n}(g)\right)\frac{q_{n}(t)}{n!}
\end{eqnarray}
\end{definition}

Ahora dotamos de estructura de anillo al conjunto $\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]])$

\begin{theorem}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$.
El conjunto $\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]])$ con la suma $\boxplus_{q}$ y el producto
$\circledast_{q}$ es un anillo conmutativo con unidades $\rho_{q}\mathfrak{A}(0)=0$ y
$\rho_{q}\mathfrak{A}(\textbf{e}_{\Delta})=q_{1}(t)=t$.
\end{theorem}
\begin{proof}
De (\ref{def_bplus}) y (\ref{def_cdast}) se sigue que
\begin{eqnarray*}
\rho_{q}\mathfrak{A}(\Delta f)\boxplus_{q}\rho_{q}\mathfrak{A}(\Delta g)&=&\rho_{q}[\mathfrak{A}(\Delta f)\boxplus\mathfrak{A}(\Delta g)]\\
\rho_{q}\mathfrak{A}(\Delta f)\circledast_{q}\rho_{q}\mathfrak{A}(\Delta g)&=&\rho_{q}[\mathfrak{A}(\Delta f)\circledast\mathfrak{A}(\Delta g)]
\end{eqnarray*}
Luego la prueba sigue igual a la que se hizo en [?].
\end{proof}

As\'i el mapa $\rho_{q}$ es un isomorfismo de anillos.

Denote $L\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$ la imagen de $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$ por el operador umbral. Extenderemos la
definici\'on de semi sistemas din\'amicos a operadores delta

\begin{definition}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$.
Definimos un \textbf{semi sistema dinámico tipo delta} como la terna $(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Psi_{Q})$
donde $R$ es el conjunto de tiempos, $L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el espacio de fases, con
$S=\mathfrak{A}(\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]])$ y $\Psi_{Q}$ es el \textbf{semi flujo delta}
\begin{equation}
\Psi_{Q}(t,x,\Delta f(x))=\rho_{q}\mathfrak{A}(\Delta f(x))
\end{equation}
esto es, $\Psi_{Q}$ es el mapa
$\Psi_{Q}:R\times L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ .
\end{definition}

Esta definici\'on junto con que $\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]])$ es un anillo nos lleva a la siguiente
definici\'on

\begin{definition}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$.
Definimos el anillo de semi sistemas din\'amicos tipo delta como el conjunto
\begin{equation}
(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]]))=\{(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\Psi_{Q}):\Psi_{Q}=\rho_{q}\mathfrak{A}(\Delta f)\}
\end{equation}
donde
\begin{small}
\begin{eqnarray*}
(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\Delta f))\boxplus_{q}(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\Delta g))&=&(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\Delta(f+g)))
\end{eqnarray*}
y
\begin{eqnarray*}
(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\Delta f))\circledast_{q}(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\Delta g))&=&(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\Delta(f\ast g)))
\end{eqnarray*}
\end{small}
y en donde $(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],0)$ y $(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],t)$ son las unidades con respecto a
$\boxplus_{q}$ y $\circledast_{q}$ respectivamente.
\end{definition}

Como veremos m\'as adelante el anillo $(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[x]],\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]]))$ contiene
todas las soluciones a las ecuaciones tipo delta $Q\Phi=f(\Phi)$ para cada
función $f(x)$ en $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Gracias a este anillo será posible descomponer las soluciones
de una ecuación diferencial aut\'onoma tipo delta de orden uno en soluciones m\'as simples.\\

Un tipo importante de ecuaciones tipo delta se dan cuando $f(x)$ es un polinomio. Si $f(x)$
es separable, entonces podemos factorizar la soluci\'on de $Q\Phi=f(\Phi)$ usando la factorizaci\'on
de $\rho_{q}\mathfrak{A}(\Delta f(x))$ en el anillo $\rho_{q}\mathfrak{A}(\dhurw_{R}[[x]])$. Primero calcularemos
$\rho_{q}\mathfrak{A}(\Delta(ax+b))$

\begin{proposition}\label{prop_pol_lin}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ y sea $Q=p(\delta)$. Entonces
\begin{equation}
\rho_{q}\mathfrak{A}(\Delta(ax+b))=\frac{1}{a}(ax+b)\left(e^{tp^{-1}(a)}-1\right),
\end{equation}
$a\neq0$ y en donde $p^{-1}$ es la serie de potencia formal inversa de $p$.
\end{proposition}
\begin{proof}
Aplicamos los operadores $\Delta$, $\mathfrak{A}$ y $\rho_{q}$ a la funci\'on $f(x)=ax+b$. As\'i tenemos
\begin{eqnarray*}
\rho_{q}\mathfrak{A}(\Delta(ax+b))&=&\rho_{q}\mathfrak{A}(ax+b,a,0,...)\\
&=&\rho_{q}(ax+b,a(ax+b),a^{2}(ax+b),...)\\
&=&\sum_{n=1}^{\infty}a^{n-1}(ax+b)\frac{q_{n}(t)}{n!}\\
&=&\frac{1}{a}(ax+b)\sum_{n=1}^{\infty}a^{n}\frac{q_{n}(t)}{n!}\\
&=&\frac{1}{a}(ax+b)\left(e^{tp^{-1}(a)}-1\right)
\end{eqnarray*}
en donde hemos usado (\ref{eqn_exp_inver}).
\end{proof}

Ahora si $f(x)$ factoriza como $(a_{1}x+b_{1})(a_{2}x+b_{2})\cdots(a_{n}x+b_{n})$ entonces
tenemos

\begin{theorem}\label{theo_poly_prod}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ y sea $Q=p(\delta)$. Suponga
$f(x)=\prod_{k=1}^{n}(a_{k}x+b_{k})$ un polinomio separable en $\mathbb{F}[x]\subset\mathop{\mbox{\textsl{H}}}\nolimits_{\mathbb{F}}[[x]]$.
Entonces
\begin{equation}
\rho_{q}\mathfrak{A}(\Delta f(x))=\bcastq_{k=1}^{n}\frac{1}{a_{k}}(a_{k}x+b_{k})\left(e^{tp^{-1}(a_{k})}-1\right)
\end{equation}
en donde $p^{-1}$ es la serie de potencia formal inversa de $p$.
\end{theorem}
\begin{proof}
Por aplicaci\'on directa de la Proposici\'on \ref{prop_pol_lin}.
\end{proof}

Ahora mostraremos otra forma de calcular $Q\Phi=f(\Phi)$. Suponga que
$f(x)=\sum_{k=0}^{n}a_{k}x^{k}$ es un polinomio de grado $n$ en $R[x]\subset\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
Queremos calcular el polinomio $\rho_{q}\mathfrak{A}(\Delta f(x))$. Usaremos la siguiente notaci\'on
\begin{eqnarray}\label{eqn_pot_rho}
\rho_{q}\mathfrak{A}[\Delta x]_{\circledast_{q}}^{k}&=&\rho_{q}\mathfrak{A}[\Delta x]\circledast_{q}\cdots\circledast_{q}\rho_{q}\mathfrak{A}[\Delta x]
\end{eqnarray}

Primero calcularemos los monomios $\rho_{q}\mathfrak{A}[\Delta x]_{\circledast_{q}}^{k}$

\begin{proposition}\label{lemma_potencia}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ y sea $Q=p(\delta)$. Sea
$L$ la representaci\'on umbral de $q_{n}(t)$. Cuando $k=0$ tenemos que
$\rho_{q}\mathfrak{A}[\Delta x]_{\circledast_{q}}^{0}=at$ y cuando $k=1$, claramente
$\rho_{q}\mathfrak{A}[\Delta x]_{\circledast_{q}}^{1}=x(e^{tp^{-1}(a)}-1)$. Para $k\geq2$ se tiene
\begin{equation}
\rho_{q}\mathfrak{A}[\Delta x]_{\circledast_{q}}^{k}=L\left(\frac{x}{\sqrt[k-1]{1-a(k-1)x^{k-1}t}}-x\right)
\end{equation}
\end{proposition}
\begin{proof}
En [?] fue mostrado que
\begin{equation*}
\rho_{t}\mathfrak{A}[\Delta x]_{\circledast}^{k}=\frac{x}{\sqrt[k-1]{1-a(k-1)x^{k-1}t}}-x.
\end{equation*}
Como $L$ es la representaci\'on umbral de $q_{n}(t)$, entonces
$\rho_{q}\mathfrak{A}[\Delta x]_{\circledast}^{k}=L(\rho_{t}\mathfrak{A}[\Delta x]_{\circledast}^{k})$
y todo sigue de aqu\'i.
\end{proof}

Por aplicación directa de esta Proposici\'on tenemos el siguiente resultado

\begin{theorem}\label{theo_poly_sum}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ y sea $Q=p(\delta)$. Suponga
$f(x)=\sum_{k=1}^{n}a_{k}x^{k}$ un polinomio $R[x]\subset\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
Entonces
\begin{eqnarray}
&&\rho_{q}\mathfrak{A}(\Delta f(x))=\bboxq_{k=0}^{n}\rho_{q}\mathfrak{A}[a_{k}\Delta x]_{\circledast_{q}}^{k}\nonumber\\
&&=a_{0}t\boxplus_{q}x(e^{tp^{-1}(a_{1})}-1)\boxplus_{q}\bboxq_{k=2}^{n}L\left(\frac{x}{\sqrt[k-1]{1-a_{k}(k-1)x^{k-1}t}}-x\right)
\end{eqnarray}
\end{theorem}

Finalmente damos la definici\'on de sistemas din\'amicos tipo delta

\begin{definition}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ y denote $L$ la representaci\'on
umbral de $q_{n}(t)$. Definimos un \textbf{sistema dinámico tipo delta} sobre el anillo $R$ como
la terna $(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Phi_{Q})$ donde $R$ es el \textbf{conjunto de tiempos},
$L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el \textbf{espacio de fases} y $\Phi_{Q}$ es el \textbf{flujo delta}
\begin{equation}\label{flujo_delta}
\Phi_{Q}(t,x,\Delta f(x))=x+\sum_{n=1}^{\infty}A_{n}([\delta^{(n-1)}f(x)])\dfrac{q_{n}(t)}{n!}.
\end{equation}
esto es, $\Phi_{Q}$ es el mapa definido
por $\Phi_{Q}:R\times L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ donde
$S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
\end{definition}

De este modo

\begin{equation}
\Phi_{Q}(t,x,\Delta f(x))=x+\rho_{q}\mathfrak{A}(\Delta f(x))
\end{equation}
y junto con los Teoremas \ref{theo_poly_prod} y \ref{theo_poly_sum} encontramos la soluci\'on a la
ecuaci\'on $Q\Phi=f(\Phi)$ cuando $f(x)$ es un polinomio.

La representaci\'on umbral $L$ de $\Phi_{Q}(t,x,\Delta f(x))$ implica que

\begin{enumerate}
\item $\Phi_{Q}(0,x,\Delta f(x))=0$
\item $\Phi_{Q}(t+s,x,\Delta f(x))=\Phi_{Q}(t,\Phi_{Q}(s,x,\Delta f(x)),\Delta f(x))$
\item $Q\Phi_{Q}(t,x,\Delta f(x))=f(x)\delta_{x}\Phi_{Q}(t,x,\Delta f(x))=f(\Phi_{Q})$.
\end{enumerate}

Podemos reescribir (\ref{flujo_delta}) de la siguiente forma
\begin{equation}\label{eqn_delta_rep}
\Phi_{Q}(t,x,\Delta f(x))=\sum_{n=0}^{\infty}Q^{n}\Phi_{Q}(0,x,\Delta f(x))\frac{q_{n}(t)}{n!}
\end{equation}
Llamamos a (\ref{eqn_delta_rep}) una \textbf{delta representaci\'on formal} de
$\Phi_{Q}(t,x,\Delta f(x))$.

El \'ultimo resultado de esta secci\'on es mostrar que el conjunto de todos los flujos tipo delta
forman un grupo

\begin{theorem}
Fije una $f(x)$ en $\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Denote $\mathcal{F}$ el conjunto de flujos tipo delta $\Phi_{Q}$
soluci\'on de $Q\Phi_{Q}=f(\Phi_{Q})$ para todo $Q\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[\delta]]$. Tome operadores delta
$A,B\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[\delta]]$ con polinomios base $(a_{n}(t))_{n\in\mathbb{N}}$ y $(b_{n}(t))_{n\in\mathbb{N}}$,
respectivamente. Defina la composici\'on $\circledcirc$ de flujos deltas $\Phi_{A}$ y $\Phi_{B}$
por $\Phi_{A}\circledcirc\Phi_{B}=\Phi_{C}$, donde
\begin{equation}
\Phi_{C}=x+\sum_{n=1}^{\infty}A_{n}([\delta^{n-1}f(x)])\frac{c_{n}(t)}{n!}.
\end{equation}
y $c_{n}(t)=a_{n}(\underline{b}(t))$. Entonces $(\mathcal{F},\circledcirc)$ es un grupo de Lie
con unidad $\Phi$.
\end{theorem}
\begin{proof}
Como $(c_{n}(t))_{n\in\mathbb{N}}$ es una sucesi\'on de conjuntos base, entonces $\Phi_{C}\in\mathcal{F}$ y
de aquí el producto $\circledcirc$ es cerrado en $\mathcal{G}$. Si $\Phi_{A}\circledcirc\Phi_{B}=\Phi$, entonces $c_{n}(t)=t^{n}$ implica que $(a_{n}(t))$ y $(b_{n}(t))$ son inversos con la
composici\'on umbral. Luego $\Phi_{A}$ y $\Phi_{B}$ son inversos con respecto a $\circledcirc$.
La asociatividad de $\circledcirc$ sigue de la asociatividad de la composici\'on de umbral $\circ$.
De este modo $(\mathcal{F},\circledcirc)$ es un grupo. Ahora defina el mapa
$\xi:\mathcal{F}\rightarrow\mathcal{U}$ por $\xi(\Phi_{A})=(a_{n}(t))_{n\in\mathbb{N}}$. Entonces
$\xi(\Phi_{A}\circledcirc\Phi_{B})=(a_{n}(t))_{n\in\mathbb{N}}\circ(b_{n}(t))_{n\in\mathbb{N}}$. As\'i
$\xi$ es un isomorfismo de grupos y por tanto $\mathcal{F}$ es un grupo de Lie.
\end{proof}

Ahora hallaremos la representaci\'on en serie de potencias formal de $\Phi_{Q}(t,x,\Delta f(x))$.
Tendremos presente que $q_{n}(t)=\sum_{i=1}^{n}\beta_{n,i}t^{i}$.

\begin{theorem}\label{theo_matriz}
Sea $Q$ un operador delta con conjunto base $q=(q_{n}(t))_{n\in\mathbb{N}}$ con
$q_{n}(t)=\sum_{k=1}^{n}\beta_{k,n}t^{k}$ para todo $n\geq1$. Entonces
\begin{equation}\label{eqn_rel}
\Phi_{Q}(t,x,\Delta f(x))=x+\sum_{n=1}^{\infty}\textbf{B}_{n}[\mathfrak{A}(\Delta f)\cdot F]t^{n}
\end{equation}
adem\'as
\begin{equation}\label{eqn_rel_sum}
\rho_{q}\mathfrak{A}(\Delta f)\boxplus_{q}\rho_{q}\mathfrak{A}(\Delta g)=\sum_{n=1}^{\infty}\textbf{B}_{n}[\mathfrak{A}(\Delta f)\boxplus\mathfrak{A}(\Delta g)]\cdot Ft^{n}
\end{equation}
y
\begin{equation}\label{eqn_rel_prod}
\rho_{q}\mathfrak{A}(\Delta f)\circledast_{q}\rho_{q}\mathfrak{A}(\Delta g)=\sum_{n=1}^{\infty}\textbf{B}_{n}[\mathfrak{A}(\Delta f)\circledast\mathfrak{A}(\Delta g)]\cdot Ft^{n}
\end{equation}
donde $\textbf{B}=(\textbf{B}_{n})_{n\in\mathbb{N}}$ es la \textbf{matriz de conexi\'on} conectando
$(q_{n}(t))$ y $(t^{n})$
\begin{equation}
\textbf{B}_{n}=(b_{n,i})_{n\in\mathbb{N}_{1}}
\end{equation}
\begin{equation}
b_{n,i}=
\begin{cases}
0, & i<n,\\
\beta_{n,i}, & i\geq n
\end{cases}
\end{equation}
y $F=(1,\frac{1}{2!},\frac{1}{3!},...)$ y $\cdot$ es el producto de Hadamard de vectores, esto es,
producto componente a componente
\end{theorem}
\begin{proof}
Probaremos (\ref{eqn_rel}) solamente, pues una prueba para (\ref{eqn_rel_sum}) y (\ref{eqn_rel_prod})
son an\'aloga. Pongamos $A_{n}\equiv A_{n}([\delta^{n-1}f])$. Tenemos
\begin{eqnarray*}
\rho_{q}\mathfrak{A}(\Delta f)&=&\sum_{n=1}^{\infty}A_{n}(f)\frac{1}{n!}\sum_{k=1}^{n}\beta_{k,n}t^{k}\\
&=&A_{1}(f)\beta_{1,1}t+A_{2}(f)(\beta_{1,2}t+\beta_{2,2}t^{2})\frac{1}{2!}\\
&&+A_{3}(f)(\beta_{1,3}t+\beta_{2,3}t^{2}+\beta_{3,3}t^{3})\frac{1}{3!}+\cdots\\
&=&\left\{A_{1}(f)\beta_{1,1}+A_{2}(f)\beta_{1,2}\frac{1}{2!}+\cdots\right\}t\\
&&+\left\{A_{2}(f)\beta_{2,2}\frac{1}{2!}+A_{3}(f)\beta_{2,3}\frac{1}{3!}+\cdots\right\}t^{2}
+\cdots
\end{eqnarray*}
Por otro lado
\begin{eqnarray*}
\textbf{B}[\mathfrak{A}(\Delta f)\cdot F]&=&
\left(
\begin{array}{cccc}
\beta_{1,1}&\beta_{1,2}&\beta_{1,3}&\cdots\\
0&\beta_{2,2}&\beta_{2,3}&\cdots\\
0&0&\beta_{3,3}&\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{array}
\right)
\left(
\begin{array}{c}
\frac{A_{1}(f)}{1!}\\
\frac{A_{2}(f)}{2!}\\
\frac{A_{3}(f)}{3!}\\
\vdots
\end{array}
\right)\\
&=&(\textbf{B}_{1},\textbf{B}_{2},\textbf{B}_{3},...)[\mathfrak{A}(\Delta f)\cdot F].
\end{eqnarray*}
Juntando todo obtenemos el resultado deseado.
\end{proof}

La matriz de conexi\'on del flujo $\Phi(t,x,\Delta f(x))$ es la matriz identidad infinita
$\textbf{I}$. Si $\textbf{A}$ y $\textbf{B}$ son las matrices de conexi\'on de los flujos
$\Phi_{A}$ y $\Phi_{B}$, entonces $\textbf{BA}$ es la matriz de conexi\'on de
$\Phi_{A}\circledcirc\Phi_{B}$. Luego el conjunto de matrices de conexi\'on forman un grupo y
tenemos el siguiente

\begin{theorem}
Denote $\mathcal{C}$ el grupo de matrices de conexi\'on de conjuntos bases. Entonces
$\mathcal{F}$ es anti-isomorfo a $\mathcal{C}$.
\end{theorem}
\begin{proof}
Sean $A,B\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[\delta]]$ operadores delta con conjuntos base $(a_{n}(t))_{n\in\mathbb{N}}$ y
$(b_{n}(t))_{n\in\mathbb{N}}$, respectivamente y sea $\Phi_{C}=\Phi_{A}\circledcirc\Phi_{B}$
con $c_{n}(t)=a_{n}(\underline{b}(t))$. Ahora sean $a_{n}(t)=\sum_{k=0}^{n}\alpha_{k,n}t^{k}$,
$b_{n}(t)=\sum_{k=0}^{n}\beta_{k,n}t^{k}$ y $c_{n}(t)=\sum_{k=0}^{n}\alpha_{k,n}t^{k}$ con matrices
de conexi\'on $\textbf{A},\textbf{B}$ y $\textbf{C}$. Defina el mapa
$\mu:\mathcal{F}\rightarrow\mathcal{C}$ por $\mu(\Phi_{A})=\textbf{A}$. Entonces
$\mu(\Phi_{A}\circledcirc\Phi_{B})=\mu(\Phi_{C})=\textbf{B}\textbf{A}=\mu(\Phi_{B})\mu(\Phi_{A})$.
Adem\'as $\mu(\Phi)=\textbf{I}$ y $\mu(\Phi_{A}^{-1})=\mu(\Phi_{A})^{-1}$, en donde $\Phi_{A}^{-1}$
es el flujo inverso de $\Phi_{A}$. Luego $\mu$ es un anti-isomorfismo de grupos.
\end{proof}

Finalizamos este art\'iculo mostrando ejemplos de sistemas din\'amicos tipo delta. En particular
cuando $Q$ es alguno de los siguientes operadores: $\vartriangle$, $\triangledown$,
$E^{\alpha}\delta$ y $\log(I+\delta)$.

\section{Sistemas din\'amicos de diferencia forward}

Un sistema din\'amico discreto uni-dimensional es uno de la forma
\begin{equation}\label{eqn_sdd}
y_{n+1}=g(y_{n})
\end{equation}
en donde $n\in\mathbb{N}$. Este tipo de ecuaciones son de gran importancia por sus aplicaciones
en econom\'ia, ciencias sociales, biolog\'ia, f\'isica, ingenier\'ia, redes neurales, etc.
La ecuaci\'on (\ref{eqn_sdd}) generalmente no tiene soluci\'on exacta cuando
$g$ es una funci\'on no lineal. Sin embargo, los resultados anteriores permiten echar a un lado
esta limitaci\'on. Usamos la diferencia forward $\vartriangle$ de $y_{n}$, esto es,
$\vartriangle y_{n}=y_{n+1}-y_{n}$. Luego (\ref{eqn_sdd}) es equivalente a la ecuaci\'on
\begin{equation}
\vartriangle y_{n}=f(y_{n})
\end{equation}
en donde $f(x)=g(x)-x$.

El operador $\vartriangle$ es un operador delta con conjunto base $q=((t)_{n})_{n\in\mathbb{N}}$ donde
$(t)_{n}=t(t-1)(t-2)\cdots(t-n+1)$ y
\begin{enumerate}
\item $(t)_{0}=1$,
\item $(0)_{n}=0$,
\item $\vartriangle(t)_{n}=n(t)_{n-1}$
\end{enumerate}

Sea $L$ la representaci\'on umbral de $(t)_{n}$. Tenemos

\begin{definition}
Definimos un \textbf{sistema dinámico forward} sobre el anillo $R$ como la terna
$(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Phi_{\vartriangle})$ donde $R$ es el \textbf{conjunto de tiempos},
$L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el \textbf{espacio de fases} y $\Phi_{\vartriangle}$ es el \textbf{flujo forward}
\begin{equation}\label{flujo_diferencial}
\Phi_{\vartriangle}(t,x,\Delta f(x))=x+\sum_{k=1}^{\infty}A_{k}([\delta^{(k-1)}f(x)])\frac{(t)_{k}}{k!}.
\end{equation}
esto es, $\Phi_{\vartriangle}$ es el mapa definido
por $\Phi_{\vartriangle}:R\times L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ donde $S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
\end{definition}

Por definici\'on $\frac{(t)_{n}}{n!}=\binom{t}{n}$. Entonces
\begin{equation}\label{eqn_flujo_for}
\Phi_{\vartriangle}(t,x,\Delta f(x))=x+\sum_{k=1}^{\infty}A_{k}([\delta^{(n-1)}f(x)])\binom{t}{k}
\end{equation}
es la soluci\'on de (\ref{eqn_sdd}) para todo $t$ en el dominio $R$. Cuando $t\in\mathbb{N}$, entonces (\ref{eqn_flujo_for}) se reduce a
\begin{equation}
y_{n}=\Phi_{\vartriangle}(n,x,\Delta f(x))=x+\sum_{k=1}^{n}A_{k}([\delta^{(n-1)}f(x)])\binom{n}{k}
\end{equation}
y tenemos una soluci\'on exacta a la ecuaci\'on (\ref{eqn_sdd}) para todo
$f(x)=g(x)-x\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$. Es de inter\'es notar que $\Phi_{\vartriangle}$ es una funci\'on
de coeficientes binomiales. Luego es posible una interpretaci\'on combinatorial a estas soluciones.

Por un uso directo de la prueba de la Proposici\'on \ref{prop_pol_lin} obtenemos

\begin{equation}
\rho_{q}\mathfrak{A}(\Delta(ax+b))=\frac{1}{a}(ax+b)\left((1+a)^{t}-1\right)
\end{equation}

y si $f(x)=\prod_{k=1}^{m}(a_{k}x+b_{k})$ entonces

\begin{equation}\label{eqn_sdd_poli}
\Phi_{\vartriangle}(t,x,\Delta f(x))=x+\bcastq_{k=1}^{m}\frac{1}{a_{k}}(a_{k}x+b_{k})\left((1+a_{k})^{t}-1\right)
\end{equation}

es soluci\'on de la ecuaci\'on en diferencias

\begin{equation}
\vartriangle y_{n}=(a_{1}y_{n}+b_{1})(a_{2}y_{n}+b_{2})\cdots(a_{m}y_{n}+b_{m}).
\end{equation}

Como caso particular tenemos la ecuaci\'on logistica. Esta es la ecuaci\'on en diferencias de orden
uno no lineal m\'as simple y hasta ahora no exist\'ia una soluci\'on en forma cerrada para ella.
Esta ecuaci\'on es encontrada en cualquier libro de texto b\'asico sobre ecuaciones de
diferencias. Esta ecuaci\'on se escribe como
\begin{equation}\label{eqn_logis}
y_{n+1}=\mu y_{n}(1-y_{n}).
\end{equation}
Para resolverla debemos resolver la ecuaci\'on equivalente $\vartriangle y_{n}=y_{n}(\mu-1-\mu y_{n})=-y_{n}(\mu y_{n}-(\mu-1))$. Aqu\'i $f(x)=-x(\mu x-(\mu-1))$. Luego
\begin{eqnarray*}
\Phi_{\vartriangle}(t,x,\Delta f(x))&=&x+\Psi_{\vartriangle}(t,x,-\Delta x)\circledast_{\vartriangle}\Psi_{\vartriangle}(t,x,-\Delta(\mu x-(\mu-1)))\\
&=&x+(-x)\circledast_{\vartriangle}[\mu x-(\mu-1)](2^{t}-1)
\end{eqnarray*}

Con esta soluci\'on notamos que
\begin{eqnarray*}
\Phi_{\vartriangle}(t,0,\Delta f(x))&=&0\\
\Phi_{\vartriangle}\left(t,\frac{\mu-1}{\mu},\Delta f(x)\right)&=&\frac{\mu-1}{\mu}
\end{eqnarray*}
As\'i los ceros de $-x(\mu x-(\mu-1))$ corresponde con los puntos fijos de $\mu x(1-x)$.

Por otro lado, la ecuaci\'on log\'istica continua $y^{\prime}=-y(\mu y-(\mu-1))$ tiene flujo
\begin{equation}
\Phi(t,x,\Delta f(x))=\frac{(\mu-1)xe^{(\mu-1)t}}{\mu x(e^{(\mu-1)t}-1)+(\mu-1)}.
\end{equation}
Luego usando la representaci\'on umbral $L(t^{n})=(t)_{n}$ obtenemos otra forma de escribir la
soluci\'on de (\ref{eqn_logis})

\begin{equation}
L\left(\frac{(\mu-1)xe^{(\mu-1)t}}{\mu x(e^{(\mu-1)t}-1)+(\mu-1)}\right)=x+(-x)\circledast_{\vartriangle}[\mu x-(\mu-1)](2^{t}-1).
\end{equation}

Otra aplicaci\'on importante de las ecuaciones de diferencias la encontramos en la teor\'ia de
fractales, especialmente en mapas polinomiales iterativos. Uno de tales mapas son los mapas
cuadr\'aticos complejos $z_{n+1}=z_{n}^{2}+c$ \'utiles para construir conjuntos de Mandelbrot.
Por ser una ecuaci\'on en diferencias de orden uno podemos usar (\ref{eqn_sdd_poli})
para encontrar la soluci\'on exacta a dicho mapa.

Resolver
\begin{equation}
z_{n+1}=z_{n}^{2}+c
\end{equation}

es resolver la ecuaci\'on equivalente

\begin{equation}
\vartriangle z_{n}=z_{n}^{2}-z_{n}+c
\end{equation}
Para este caso tenemos que $f(x)=x^{2}-x+c=(x-\alpha)(x-\overline{\alpha})$ en donde
$\alpha=\frac{1+\sqrt{4c-1}i}{2}$.
Luego
\begin{eqnarray*}
z_{n}&=&\Phi_{\vartriangle}(n,x,\Delta f(x))\\
&=&x+\Psi_{\vartriangle}(n,x,\Delta(x-\alpha))\circledast_{\vartriangle}\Psi_{\vartriangle}(n,x,\Delta(x-\overline{\alpha}))\\
&=&x+[(x-\alpha)(2^{n}-1)]\circledast_{\vartriangle}[(x-\overline{\alpha})(2^{n}-1)]
\end{eqnarray*}

Igualmente como se hizo con la ecuaci\'on log\'istica tenemos
\begin{equation}
L\left(\frac{\overline{\alpha}e^{\sqrt{4c-1}it}(x-\alpha)-\alpha(x-\overline{\alpha})}{(x-\alpha)e^{\sqrt{4c-1}it}-(x-\overline{\alpha})}\right)=x+[(x-\alpha)(2^{n}-1)]\circledast_{\vartriangle}[(x-\overline{\alpha})(2^{n}-1)]
\end{equation}

\section{Sistemas din\'amicos de diferencia backward}

Si $Q$ es el operador backward $\triangledown$, entonces su conjunto base son los factoriales
rising $q=((t)^{n})_{n\in\mathbb{N}}$ donde $(t)^{n}=t(t+1)(t+2)\cdots(t+n-1)$ y
\begin{enumerate}
\item $(t)^{0}=1$,
\item $(0)^{n}=0$,
\item $\triangledown(t)^{n}=n(t)^{n-1}$
\end{enumerate}

Con el operador backward podemos resolver ecuaciones en diferencias de la forma

\begin{equation}
\triangledown y_{n}=f(y_{n})
\end{equation}
en donde $\triangledown y_{n}=y_{n}-y_{n-1}$. Las ecuaciones en diferencias backward se aplican en
las mismas \'areas en donde son aplicados las ecuaciones en diferencias forward.

Sea $L(t^{n})=(t)^{n}$ la representaci\'on umbral de $(t)^{n}$. Entonces tenemos la siguiente
definici\'on

\begin{definition}
Definimos un \textbf{sistema dinámico backward} sobre el anillo $R$ como la terna
$(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Phi_{\triangledown})$ donde $R$ es el \textbf{conjunto de tiempos},
$L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el \textbf{espacio de fases} y $\Phi_{\triangledown}$ es el \textbf{flujo backward}
\begin{equation}\label{flujo_diferencial}
\Phi_{\triangledown}(t,x,\Delta f(x))=x+\sum_{k=1}^{\infty}A_{k}([\delta^{(k-1)}f(x)])\frac{(t)^{k}}{k!}.
\end{equation}
esto es, $\Phi_{\triangledown}$ es el mapa definido
por $\Phi_{\triangledown}:R\times L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ donde $S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
\end{definition}

Por definici\'on $\frac{(t)^{k}}{k!}=\binom{t+k-1}{k}=\genfrac{\langle}{\rangle}{0pt}{}{t}{k}$. Entonces
\begin{equation}
\Phi_{\triangledown}(t,x,\Delta f(x))=x+\sum_{k=1}^{\infty}A_{k}([\delta^{(n-1)}f(x)])\genfrac{\langle}{\rangle}{0pt}{}{t}{k}
\end{equation}
El n\'umero $\genfrac{\langle}{\rangle}{0pt}{}{n}{k}$ es conocido como coeficiente binomial con repeticiones y es usado
para contar multiset de cardinalidad $k$ con elementos tomados de un conjunto de $n$ elementos.

Denote $\mathop{\mbox{\textsf{exp}}}\nolimits(a)$ la sucesi\'on $(a^{n})_{n\in\mathbb{N}_{1}}$ en $H_{R}$. En [3] fue mostrado que
\begin{equation}\label{eqn_R_prod}
\mathfrak{A}(a\Delta f(x))=\mathop{\mbox{\textsf{exp}}}\nolimits(a)\mathfrak{A}(\Delta f(x)),
\end{equation}

El siguiente resultado relaciona a los sistemas din\'amicos discretos de tipo forward y backward

\begin{theorem}
Para todo $f(x)\in\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$ se cumple
\begin{equation}
\Phi_{\triangledown}(t,x,\Delta f(x))=\Phi_{\vartriangle}(-t,x,-\Delta f(x)).
\end{equation}
\end{theorem}
\begin{proof}
Sabemos que $\genfrac{\langle}{\rangle}{0pt}{}{t}{k}=(-1)^{k}\binom{-t}{k}$. As\'i
\begin{eqnarray*}
\Phi_{\triangledown}(t,x,\Delta f(x))&=&x+\sum_{k=1}^{\infty}A_{k}([\delta^{(n-1)}f(x)])\genfrac{\langle}{\rangle}{0pt}{}{t}{k}\\
&=&x+\sum_{k=1}^{\infty}A_{k}([\delta^{(n-1)}f(x)])(-1)^{k}\binom{-t}{k}\\
&=&\Phi_{\vartriangle}(-t,x,-\Delta f(x))
\end{eqnarray*}
en donde hemos usado el (\ref{eqn_R_prod}) con $a=-1$.
\end{proof}

Por otro lado, podemos usar el anterior teorema y la ecuaci\'on (\ref{eqn_sdd_poli}) para mostrar que

\begin{equation}
\Phi_{\triangledown}(t,x,\Delta f(x))=x+\bcastq_{k=1}^{m}\frac{1}{a_{k}}(a_{k}x+b_{k})((1-a_{k})^{t}-1)
\end{equation}

es soluci\'on de la ecuaci\'on en diferencia backward

\begin{equation}
\triangledown y_{n}=(a_{1}y_{n}+b_{1})(a_{2}y_{n}+b_{2})\cdots(a_{m}y_{n}+b_{m}).
\end{equation}

El teorema anterior implica que es suficiente con resolver uno de los dos tipos de ecuaciones de
diferencias y luego por aplicaci\'on de este obtener la soluci\'on del otro tipo.

\section{Sistemas din\'amicos tipo Abel}

Denote $A(\alpha)$ el operador delta $E^{\alpha}\delta$. Este operador es conocido como operador Abel
donde sus polinomios bases son los polinomios de Abel
$(t(t-n\alpha)^{n-1})_{n\in\mathbb{N}}$ y
\begin{enumerate}
\item $t(t)^{-1}=1$,
\item $0(0-n\alpha)^{n-1}=0$,
\item $A(\alpha)t(t-n\alpha)^{n-1})=nt(t-(n-1)\alpha)^{n-2}$
\end{enumerate}
Los polinomios de Abel fueron estudiados en [1] en la representaci\'on de funciones.
Queremos resolver ecuaciones de la forma

\begin{equation}
A(\alpha)\Phi=f(\Phi).
\end{equation}

Sea $L$ la representaci\'on umbral de $(t(t-n\alpha)^{n-1})_{n\in\mathbb{N}}$. Tenemos la
siguiente definici\'on

\begin{definition}
Definimos un \textbf{sistema dinámico tipo Abel} sobre el anillo $R$ como la terna
$(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Phi_{A(\alpha)})$ donde $R$ es el \textbf{conjunto de tiempos},
$L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el
\textbf{espacio de fases} y $\Phi_{A(\alpha)}$ es el \textbf{flujo Abel}
\begin{equation}\label{flujo_diferencial}
\Phi_{A(\alpha)}(t,x,\Delta f(x))=x+\sum_{n=1}^{\infty}A_{n}([\delta^{(n-1)}f(x)])\frac{t(t-n\alpha)^{n-1}}{n!}.
\end{equation}
esto es, $\Phi_{A(\alpha)}$ es el mapa definido
por $\Phi_{A(\alpha)}:R\times L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ donde $S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
\end{definition}

Usando (\ref{eqn_sdd_poli}) encontramos que

\begin{equation}
\Phi_{A(\alpha)}(t,x,\Delta f(x))=x+\bcastq_{k=1}^{m}\frac{1}{a_{k}}(a_{k}x+b_{k})\left(e^{tW(\alpha)/\alpha}-1\right),
\end{equation}
en donde $W$ es la funci\'on $W$ de Lambert, es la soluci\'on de la ecuaci\'on
\begin{equation}
A(\alpha)\Phi=(a_{1}\Phi+b_{1})(a_{2}\Phi+b_{2})\cdots(a_{n}\Phi+b_{n})
\end{equation}

La siguiente identidad relaciona los flujos de las ecuaciones $A(\alpha)\Phi=af(\Phi)$ y
$A(a\alpha)\Phi=f(\Phi)$

\begin{theorem}
Para todo $a\in R$ se cumple
\begin{equation}
\Phi_{A(\alpha)}(t,x,a\Delta f(x))=\Phi_{A(a\alpha)}(at,x,\Delta f(x)).
\end{equation}
\end{theorem}
\begin{proof}
\begin{eqnarray*}
\Phi_{A(\alpha)}(t,x,a\Delta f(x))&=&x+\sum_{n=1}^{\infty}a^{n}A_{n}([\delta^{(n-1)}f(x)])\frac{t(t-n\alpha)^{n-1}}{n!}\\
&=&x+\sum_{n=1}^{\infty}A_{n}([\delta^{(n-1)}f(x)])\frac{(at)(at-na\alpha)^{n-1}}{n!}\\
&=&\Phi_{A(a\alpha)}(at,x,\Delta f(x)).
\end{eqnarray*}
\end{proof}

Con esta identidad podremos, quizas, darle una estructura m\'as que de grupo uni par\'ametrico
que la que se obtiene con el flujo de tipo Abel, pero no ser\'a hecho en este art\'iculo.

\section{Sistemas din\'amicos tipo Touchard}

Los polinomios de Touchard son polinomios de la forma
\begin{equation}
\phi_{n}(t)=\sum_{k=0}^{n}S(n,k)t^{k}
\end{equation}
en donde $S(n,k)$ son los n\'umeros de Stirling de segunda clase.
Los polinomios de Touchard son polinomios base del operador delta de Touchard $T=\log(I+\delta)$ y
satisfacen
\begin{enumerate}
\item $\phi_{0}(t)=1$,
\item $\phi_{n}(0)=0$,
\item $\log(I+\delta)\phi_{n}(t)=n\phi_{n-1}(t)$.
\end{enumerate}

Sea $L$ la representaci\'on umbral de $\phi_{n}(t)$. Entonces tenemos la definici\'on

\begin{definition}
Definimos un \textbf{sistema dinámico tipo Touchard} sobre el anillo $R$ como la terna
$(R,L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]],\Phi_{T})$ donde $R$ es el \textbf{conjunto de tiempos}, $L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ el
\textbf{espacio de fases} y $\Phi_{T}$ es el \textbf{flujo Touchard}
\begin{equation}\label{flujo_diferencial}
\Phi_{T}(t,x,\Delta f(x))=x+\sum_{n=1}^{\infty}A_{n}([\delta^{(n-1)}f(x)])\dfrac{\phi_{n}(t)}{n!}.
\end{equation}
esto es, $\Phi_{T}$ es el mapa definido
por $\Phi_{T}:R\times L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]\times\Delta\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]\rightarrow L\mathop{\mbox{\textsl{H}}}\nolimits_{S}[[t]]$ donde
$S=\mathop{\mbox{\textsl{H}}}\nolimits_{R}[[x]]$.
\end{definition}

Usando (\ref{eqn_sdd_poli}) encontramos que

\begin{equation}
\Phi_{T}(t,x,\Delta f(x))=x+\bcastq_{k=1}^{m}\frac{1}{a_{k}}(a_{k}x+b_{k})\left(e^{t(e-1)}-1\right)
\end{equation}

es el flujo de Touchard de la ecuaci\'on

\begin{equation}
\log(I+\delta)\Phi=(a_{1}\Phi+b_{1})(a_{2}\Phi+b_{2})\cdots(a_{n}\Phi+b_{n})
\end{equation}

\section{Conclusi\'on}
En este art\'iculo los sistemas din\'amicos tipo delta fueron definidos y algunas de sus propiedades
estudiadas. La teor\'ia aqu\'i desarrollada es muy \'util para resolver sistemas din\'amicos tipo
delta no lineales. Particularmente los sistemas din\'amicos de diferencia y de tipo Touchard pueden
ser \'util para resolver problemas en combinatoria enumerativa. Como los polinomios de Abel
son usados en la teor\'ia algebraica estad\'istica, entonces los sistemas din\'amicos tipo Abel
podr\'ian ser \'util para aplicaciones en estad\'istica. En este art\'iculo la derivaci\'on fue
variada. Estamos interesados en variar el dominio de integridad y dejar fijo el operador delta
y de esta forma estudiar los distintos resultados que se pueden obtener. Otra posibilidad es
variar el factorial, esto es, construir una teor\'ia para un factorial generalizado.


\section{Introduction}

Accurate theoretical predictions of properties of the ground and excited states of few-atomic molecules are of crucial importance for preparation and interpretation of different fundamental physical experiments such as laser-synthesis of diatomics~\cite{Ferber2009,Orban2015}, laser cooling of diatomic and polyatomic molecules, experiments aimed at the search for the fundamental symmetry violations effects such as the electron electric dipole moment, etc.
In particular, both diagonal and off-diagonal matrix elements of magnetic dipole hyperfine structure (HFS) and other operators, strongly localized in vicinities of atomic nuclei and sometimes referred to as core property operators, are required for the accurate treatment of non-adiabatic effects~\cite{Petrov:13,Skripnikov:2019a}, parity non-conserving amplitudes~\cite{Geddes:18} and some other transition properties in molecules. In this regard, magnetic dipole hyperfine structure calculations are of particular interest since they are conventionally used as a probe of the accuracy of theoretical predictions.

\begin{figure}
\center
\includegraphics[scale=0.3]{KCs-PEC}
\caption{Potential energy curves of the $X^1\Sigma^+$ and $a^3\Sigma^+$ states of the KCs molecule.}\label{fig:pec}
\end{figure}

An optimization of paths of laser synthesis of diatomic alkali molecules requires exceptionally complete information on both the energy and radiative properties of their rovibronic states, typically provided by nonadiabatic models including both \textit{ab initio} and spectroscopic data \cite{Pazyuk2015,Klincare2012}.

In this connection the most challenging problem is the accurate description of the complex formed by the ground and the first excited electronic states of alkali dimers, $X{}^1\Sigma^+ \sim a{}^3\Sigma^+$ (see Fig. \ref{fig:pec} for KCs), since rovibrational levels belonging to this complex are exploited as the working ones at the first step of any laser synthesis procedure (photoassociation or magnetoassociation). The $X{}^1\Sigma^+ - a{}^3\Sigma^+$ mixing due to the hyperfine interaction results in extremely complex compositions of rovibronic states near the dissociation limit \cite{Ferber2009,Orban2015}. Most nonadiabatic models of the $X{}^1\Sigma^+ \sim a{}^3\Sigma^+$ complex constructed to date assumed hyperfine interaction matrix elements to be constant with respect to the internuclear distance $R$, and thus their values were derived from the atomic ones (for example, see \cite{Ferber2009} and references therein). However, there are indications that the hyperfine interaction matrix elements actually depend on $R$ at least in the bound region of potential energy curves (PECs) \cite{Docenko2011,Lysebo2013,Temelkov2015}, but the character of this dependence is quite unclear since no highly reliable calculations of the HFS matrix elements in alkali diatomics were reported. The clarification of this question can not only justify the neglect of $R$-dependence for large values of $R$, but also bring the nonadiabatic models to a new level of accuracy which will be adequate to the current level of high-resolution molecular spectroscopy.

In the present paper, we report \textit{ab initio} calculations of the HFS matrix elements in the KCs molecule. The relativistic Fock space coupled cluster (FS-RCC) method \cite{Visscher2001} combined with the denominator-shift technique was employed throughout the paper. This combination has been proven to be a powerful tool for modeling both the energy and radiative properties of such molecules \cite{Zaitsevskii2017,Znotins2019}. Being able to treat relativistic and correlation effects simultaneously, which is of particular importance for really high-precision simulations of systems containing heavy nuclei \cite{Eliav2017}, the FS-RCC method is not well-suitable for analytic calculation of density matrices and calculation of properties due to difficulties with the extraction of explicit wave functions \cite{Szalay1995,Zaitsevskii2018}. This difficulty can be bypassed by using the finite-field (FF) approach. This approach is straightforward for diagonal matrix elements and is based on the Hellmann-Feynman theorem. However, its extension to the case of off-diagonal elements originally introduced in the multireference perturbation theory framework \cite{Zaitsevskii1998} and recently generalized for the case of the Fock space coupled cluster method \cite{Zaitsevskii2018} also allows one to avoid analytic construction of CC transition density matrices. Moreover, within this approach, all off-diagonal matrix elements of interest can be obtained simultaneously. Here we report the first application of the off-diagonal FF technique to calculation of core property matrix elements.
In Sec. \ref{sec:theory} we present the theory underlying the magnetic dipole hyperfine interaction in molecules and the extension of the finite-field technique to off-diagonal HFS matrix elements. Technical details describing the computational approach are given in Sec. \ref{sec:comput}. Sec. \ref{sec:results} reports the calculations of the HFS matrix elements mixing the $X{}^1\Sigma^+$ and $a{}^3\Sigma^+$ states of the KCs molecule, including the discussion of the corresponding selection rules, the role of core relaxation and correlation, the dependence of the HFS matrix elements on internuclear separation and the analysis of possible errors. The final Sec. \ref{sec:conclusion} concludes on our research.

\section{Theoretical considerations}\label{sec:theory}

\subsection{Molecular hyperfine structure}

The magnetic dipole hyperfine interaction of the nuclear magnetic moment $\bm{\mu} = \mu\mu_N \mathbf{I}/I$ (where $\mu$ is the value of the nuclear magnetic moment in nuclear magnetons $\mu_N$ and $\mathbf{I}$ is the nuclear spin) with electrons is given by the following operator (see, e.g.~\cite{Norman2018})
\begin{equation}\label{eq:hhfs}
H_{\rm hf}= \bm{\mu} \cdot\bm{T},
\end{equation}
where the electronic part of the operator is
\begin{equation}\label{eq:hfsEl}
\bm{T}=\sum_i \frac{[\mathbf{r}_i\times \bm{\alpha}_i]}{r_i^3},
\end{equation}
$\bm{\alpha}$ is the set of Dirac alpha matrices and $\mathbf{r}_i$ is a radius-vector of electron $i$ with respect to the center of the chosen nucleus with the magnetic moment $\mu$.

In the present paper we are interested in matrix elements of the operator (\ref{eq:hfsEl}) between two electronic states with wavefunctions $\psi_n$ and $\psi_m$:
\begin{equation}\label{eq:aperp}
A_{\eta}^{nm} = \frac{\mu}{I} \braket{ \psi_n | T_{\eta} | \psi_m},
\end{equation}
where $\eta=x,y,z$.
In these designations the usual ``pa\-ral\-lel'' component of the HFS constant, i.e. diagonal matrix element of the HFS operator, for the given electronic state $\psi_n$ of the diatomic molecule with the total electronic angular momentum projection $\Omega$ onto the internuclear axis $z$ is:
\begin{equation}\label{eq:apar}
A_{||} = \frac{1}{\Omega} A_{\eta}^{nn}=\frac{\mu}{I \Omega} \braket{ \psi_n | T_z | \psi_n}.
\end{equation}

The diagonal matrix elements (\ref{eq:apar}) are non-zero for all $\Omega > 0$; off-diagonal matrix elements are non-zero for electronic states with $\Delta \Omega = 1$ or $0^+$ and $0^-$ states. For the case of the $X^1\Sigma^+ \sim a^3\Sigma^+$ complex in KCs these selection rules result in non-zero matrix elements $1 - 1$, $0^+ - 0^-$, $0^- - 1$ and $0^+ - 1$ (we note that the ${}^3\Sigma^+$ state consists of the $0^-$ and $1$ spin-orbit-split states, while the ${}^1\Sigma^+$ state remains unsplit, $0^+$). At the K(${}^2S$)+Cs(${}^2S$) dissociation limit all these non-zero matrix elements converge to the atomic HFS constant (up to a factor 1/2, arising from the $\Omega = 1/2$ value for the atomic ${}^2S$ state).

It should be emphasized that in all the previous studies of HFS in alkali dimers the non-relativistic approximation to the Hamiltonian (\ref{eq:hhfs}) was employed; the Fermi-contact term was shown to be dominant in the HFS mixing of the $X{}^1\Sigma^+$ and $a{}^3\Sigma^+$ states \cite{Lysebo2013}. The relativistic form of the HFS Hamiltonian (\ref{eq:hhfs}) is more general and precise, including the paramagnetic spin-orbit, Fermi-contact, and spin-dipole hyperfine operators.

\subsection{Finite-field approach to hyperfine interaction matrix element calculations}

Since multireference coupled cluster models do not provide closed expressions for wavefunctions, the evaluation of property matrix elements through constructing pure-state and transition density matrices is normally very costly \cite{Szalay1995} or some linearization approximations are required. However,  most properties represented by one-electron operators can be calculated in a more convenient manner by numerical differentiation with respect to  the amplitude of some external perturbation. This approach usually called the finite-field method is widely used to evaluate expectation values. For example, the magnetic dipole hyperfine structure constant for the electronic state $\psi_n$ can be considered as the partial derivative of molecular energy with respect to the parameter $\lambda$ formally setting the strength of the hyperfine interaction, $H_{\rm hf} \to \lambda H_{\rm hf}$:
\begin{equation}
\braket{\psi_n|T_{\eta}|\psi_n} = \frac{\partial E_n}{\partial \lambda} \left| \begin{array}{l}\\_{\!\lambda=0}\end{array}\right. .
\label{ffdiag}
\end{equation}
The counterpart of Eq.~(\ref{ffdiag}) for off-diagonal matrix elements can be obtained from the approximate Hellmann-Feynman-like relation for effective operators formulated in Ref.~\cite{Zaitsevskii1998}. For the case of the magnetic dipole hyperfine structure, the $\braket{\psi_n|T_{\eta}|\psi_m}$ matrix element between two electronic states with wavefunctions $\psi_n$ and $\psi_m$ satisfy the approximate relation
\begin{equation}
\braket{\psi_n|T_{\eta}|\psi_m} \approx (E_m - E_n) \left< \tilde{\psi}_n^{\perp\perp}(\lambda_\eta) \left| \frac{\partial}{\partial \lambda_\eta} \tilde{\psi}_m(\lambda_\eta)\right. \right> \left| \begin{array}{l}\\_{\!\lambda_\eta=0}\end{array}\right.\!\! .
\label{ffhfs}
\end{equation}
Here $\eta=x,\,y,\,z$ and $\tilde{\psi}^{\perp\perp}$, $\tilde{\psi}$ denote left and right eigenvectors of the field-dependent non-Hermitian FS-RCC effective Hamiltonian acting in the field-independent (constructed for $\lambda=0$) model space, and $E_{m}$ and $E_{n}$ are the field-free energies of the involved states. The derivative in the r.h.s. of Eq. (\ref{ffhfs}) is readily estimated  by the finite difference method. Although the formula (\ref{ffhfs}) involves only the effective Hamiltonian eigenvectors (the model space projections of many-electron wavefunctions), the resulting matrix elements estimates implicitly incorporate the bulk of the contributions from the remainder parts of these wavefunctions~\cite{Zaitsevskii1998,Zaitsevskii2018}.

Note that the approximate matrix elements $\braket{\psi_n|T_{\eta}|\psi_m}$ and $\braket{\psi_m|T_{\eta}|\psi_n}^*$ defined by Eq.~(\ref{ffhfs}) generally do not coincide due to the non-Hermitian nature of the FS-RCC effective Hamiltonian. To prevent the appearance of unphysical differences between these quantities (which are normally within 1\%), the preliminary transformation of the non-Hermitian FSRCC effective Hamiltonians to the Hermitian form \emph{via} the symmetric orthogonalization of their eigenvectors can be performed.

For the case of the KCs molecule the differentiation (\ref{ffdiag}) allows one to obtain only the diagonal $1 - 1$ matrix element; the other three ($0^+ - 0^-$, $0^- - 1$ and $0^+ - 1$) can be calculated within the relation (\ref{ffhfs}). Furthermore, in alkali dimer molecules the lowest $0^-$ and $1$ electronic states corresponding to the $a^3\Sigma^+$ state are nearly not perturbed by the spin-orbit interaction with higher-lying states and hence are actually degenerate (see Fig. \ref{fig:pec}). This results in insurmountable numerical instabilities when using (\ref{ffhfs}) and the $0^- - 1$ matrix element cannot be accessed within the current formulation of the finite-field technique. However, this does not give rise to any problem since the $1 - 1$ and $0^- - 1$ matrix elements are expected to have nearly the same values (and the same dependence on internuclear distance, which is even more important here).

\section{Computational details}\label{sec:comput}

The relativistic Fock space coupled cluster method with single and double excitations (FS-RCCSD), recently demonstrated to be ideally suitable for high-precision modeling of energy and radiative properties of alkali diatomics \cite{Znotins2019}, was used throughout.
The coupled cluster calculations were carried out within the EXP-T program package \cite{Oleynichenko2020}
\footnote{See also http://www.qchem.pnpi.spb.ru/expt, accessed on 28 June 2020}.
Molecular spinors and molecular integrals, including the HFS ones, were calculated within the DIRAC code \cite{DIRAC17,Saue2020}.

Since the magnetic dipole hyperfine interaction is mainly determined by the behavior of electronic wavefunctions in the vicinity of the atomic nucleus, the all-electron four-component relativistic calculations with the Dirac-Coulomb Hamiltonian were performed in order to achieve high and controllable accuracy. An alternative two-step approach \cite{Titov2005} based on a restoration of spinors in the vicinity of a heavy nucleus is formally less precise (however, it can be used in combined schemes to treat high-order correlation effects, basis set corrections, etc. to increase the final accuracy ~\cite{Skripnikov:16b}).
In the present calculations, all electrons were correlated; the role of core relaxation and correlations is discussed below in Sec. \ref{sec:results}.

The ground state of KCs${}^{2+}$ was used as the reference state in the FS-RCC calculations and the electronic states of the neutral KCs molecule were constructed in the $(0h,2p)$ Fock space sector (two particles over the Fermi vacuum). The active space comprised 12 lowest Kramers pairs of spinors with the total angular momentum projection $|m_j| = 1/2$ , 6 pairs with $|m_j| = 3/2$ and 2 pairs with $|m_j| = 5/2$. At the dissociation limit this corresponds to the $ns$, $np$, $(n-1)d$ and $(n+1)$s spinors of the isolated K and Cs atoms. In order to suppress numerical instabilities due to intruder states in the FS-RCC calculations, the new modification of denominator shift technique \cite{Zaitsevskii2017} based on the simulation of an imaginary shift was used (see Appendix A).
We have used flexible uncontracted basis sets $(18s14p6d1f)$ for K and $(28s22p15d5f)$ for Cs (see Supplementary \cite{Supplementary}). These sets are based on the Dyall's basis set family \cite{Dyall2009}\footnote{Available from the Dirac web site, http://dirac.chem.sdu.dk, accessed on 28 June 2020.} and augmented with functions required to reproduce the atomic HFS constants and low-lying energy levels with the fairly good accuracy (see Table \ref{tab:athfs}).
The finite-field approach with the
numerical differentiation step $\Delta\lambda = 10^{-8}$ MHz${}^{-1}$ ($\approx 1.3 \cdot 10^{-4}$ a.u.) was employed to calculate hyperfine interaction expectation and transition matrix elements; this value was found to be optimal for ensuring numerical stability.

\begin{table}[h]
\center
\def1.2{1.2}
\caption{Atomic hyperfine structure constants and energy levels of $^{39}$K and $^{133}$Cs calculated within the FS-RCCSD method and compared with experimental data.}\label{tab:athfs}
\begin{tabular}{ccrrcrr}
\hline
\multicolumn{2}{c}{State} & \multicolumn{2}{c}{A, MHz} & & \multicolumn{2}{c}{
Transition energy, cm${}^{-1}$} \\
\cline{3-4} \cline{6-7}
& & FS-RCCSD & exptl\cite{Tiecke2011,Steck2010} & & FS-RCCSD & exptl\cite{Sansonetti2005} \\
\hline

\multicolumn{7}{c}{\it Potassium} \\

$4s$ & ${}^2S_{1/2}$ & 222.4   & 230.9 & & 0     & 0     \\
$4p$ & ${}^2P_{1/2}$ & 26.7     & 27.8 & & 12966 & 12985 \\
     & ${}^2P_{3/2}$ & 6.0      &  6.1 & & 13025 & 13043 \\

\multicolumn{7}{c}{\it Cesium} \\

$6s$ & ${}^2S_{1/2}$ & 2252.7 & 2298.2 & & 0     & 0     \\
$6p$ & ${}^2P_{1/2}$ & 280.3   & 291.9 & & 11184 & 11178 \\
     & ${}^2P_{3/2}$ & 48.6     & 50.3 & & 11732 & 11732 \\
\hline
\end{tabular}
\end{table}

We considered only the most abundant isotopomer ${}^{39}$K${}^{133}$Cs; the HFS constants for other isotopes can be obtained simply by scaling. The qualitative picture of energy levels at the K(${}^2S$)+Cs(${}^2S$) dissociation limit is strongly dominated by the HFS induced by the ${}^{133}$Cs nucleus ($I$ = 7/2, $A_{6s ^2S_{1/2}}$ = 2298.1579425 MHz \cite{Steck2010}); the effect induced by the ${}^{39}$K nucleus ($I$ = 3/2, $A_{4s ^2S_{1/2}}$ = 230.8598601(3) MHz \cite{Tiecke2011}) is much smaller. The picture is actually the same for the two other natural potassium isotopes (${}^{40}$K and ${}^{41}$K). Since the magnetic dipole hyperfine interaction is a one-electron property essentially defined by the behavior of the wavefunctions in the vicinity of nuclei, the HFS matrix element dependencies $A(R)$ can be calculated for the K and Cs nuclei separately. The contributions to the hyperfine structure induced by nuclear electric quadrupole moments are negligible for alkali diatomics \cite{Aldegunde2008,Lysebo2013} and hence are not considered here. The values of nuclear dipole magnetic moments, $\mu_{^{39}K}$ = 0.39147$\mu_N$ and $\mu_{^{133}Cs}$ = 2.582025$\mu_N$, were taken from Ref.~\cite{Stone2005}.

Special attention is to be paid to the accuracy limitations of the model employed. The most important factors restricting the overall accuracy of the calculated HFS matrix elements are the following ones:
(i) the approximation to the re\-la\-ti\-vis\-tic Hamiltonian (lack of Breit and QED terms),
(ii) restrictions of the point dipole model of nuclear magnetization distribution,
(iii) the approximation to the cluster operator (e.g. CCSD) used in the FS-RCC calculation and
(iv) errors caused by the basis set incompleteness.
The latter factor can be suppressed rather efficiently by using extended uncontracted basis sets designed specifically for HFS calculations. The contribution of the Breit two-electron interaction to the $A_{6s\ ^2S_{1/2}}$ constant of Cs atom was found to be less than 0.5\% and even smaller for lighter alkali atoms \cite{Safronova1999}, thus the Dirac-Coulomb approximation can be regarded as precise enough. Errors arising from the point-dipole approximation were estimated recently to be smaller than 1\% \cite{Ginges:2018,Prosnyak2020}. The effect of higher-order excitations in the cluster operator was found to be small for the Cs $6s\ ^2S_{1/2}$ HFS constant (ca.~$\sim$2\%) \cite{Haase2020}. Finally, errors due to non-compensated denominator shifts used to suppress intruder states in a wide range of $R$ are negligible for such low-lying electronic states \cite{Zaitsevskii2017,Znotins2019}. However, the clarification of the $R$-dependence of HFS matrix elements does not require extremely high accuracy; an error of about 5\% is to be regarded as quite satisfactory.

\section{Results and discussion}\label{sec:results}

The magnetic dipole hyperfine interaction matrix elements between the $X^1\Sigma^+_{0^+}$ and $a^3\Sigma^+_{0^-,1}$ states of ${}^{39}$K${}^{133}$Cs as functions of the internuclear distance are shown in Fig. \ref{fig:hfs}. The obtained $R$-dependencies of the $0^-$ -- $0^+$ coupling elements closely resemble those reported in \cite{Lysebo2013} for the Fermi-contact interaction in the $a^3\Sigma^+_u$ state of homonuclear alkali dimers.

\begin{figure}[h]
\center
\includegraphics[width=0.49\textwidth]{hfs-K-Cs.pdf}
\caption{Hyperfine coupling matrix elements mixing the $X^1\Sigma^+_{0^+}$ and $a^3\Sigma^+_{0^-,1}$ states of ${}^{39}$K${}^{133}$Cs for the HFS induced by the ${}^{133}$Cs nucleus (top panel) and the ${}^{39}$K nucleus (bottom panel). The black horizontal ticks show the corresponding experimental values for the atomic hyperfine structure constant $A_{at}/2$.}\label{fig:hfs}
\end{figure}

The widely used assumption that the hyperfine coupling matrix elements are constant is not correct for the ``bonding'' region of potential curves ($R <8$ \AA) where the deviations from the corresponding atomic values reach 15\% (see Fig. \ref{fig:pec}). However, these deviations are small for a larger $R$. Thus the variation of the matrix elements of the magnetic dipole hyperfine interaction with $R$ is not expected to affect significantly the hyperfine mixing near the dissociation limit.
This conclusion is in agreement with the experimentally observed hyperfine splittings which are nearly constant with respect to the vibrational quantum number for the $a^3\Sigma^+$ state (see Fig. 5 in \cite{Ferber2009}).
This fact legitimizes the use of atomic values for the modeling of laser synthesis \textit{via} Feshbach resonances or photoassociation.

The variations of HFS matrix elements functions are essentially confined to the ``bonding'' domain of $R$ values. It is thus natural to suppose that these variations originate mainly from the redistribution of valence-shell electrons during the formation of chemical bonds.
In other words, these variations can be considered as ``chemical shifts''. This feature of the HFS matrix elements can be used to reduce the number of adjustable parameters in non-adiabatic models considering heteronuclear alkali dimers.

The elucidation of the role of core relaxation and correlations is of great interest since the neglect of these effects can in principle lead to significant computational savings, making feasible calculations of core properties for large systems. Within the FS-RCCSD approach, core correlations are associated with the double excitations in the cluster operator which involve two core spinors. Omitting these excitations, one neglects the core correlations but retains the bulk of the core relaxation and core-valence correlation which are described by cluster terms involving only one core spinor. Fig.~\ref{fig:core} presents the $R$-dependencies of the $1 - 1$ and $0^+ - 0^-$ matrix elements of the Cs-induced HFS calculated with different sizes of the Cs core. It can be seen that electron correlations restricted to the core including all spinors with the principal quantum numbers $n \le 3$ practically do not contribute to the HFS matrix elements. On the contrary, the results are distorted significantly if the core relaxation and core-valence correlations are omitted as well (i.e. all core spinors are frozen at the post-SCF correlation calculation step): the deviation from the results with the full correlation treatment (black solid curve) reaches 5\% for the frozen $n \le 3$ core and 10\% for the frozen $n \le 4$ core. The latter error is of the same order of magnitude as that of the two-step non-variational core restoration scheme \cite{Titov2005}, which also completely neglects core relaxation and core-valence correlation contributions. However, since these contributions are actually independent from the chemical environment (see Fig. \ref{fig:core}), it might be reasonable to keep the core spinors frozen in molecular calculations and to extract the missing contributions from atomic calculations. A similar scheme can be used for evaluating other core properties. Furthermore, HFS matrix elements calculated as functions of molecular geometry with the frozen core can be adjusted to experimental atomic values simply by appropriate shift.

\begin{figure}[h]
\center
\includegraphics[width=0.49\textwidth]{KCs_HFS_Core.pdf}
\caption{Hyperfine coupling matrix elements for the HFS induced by the ${}^{133}$Cs nucleus obtained at the different levels of treatment of core electrons. Black solid line: all electrons are correlated; ``frozen'': core shells with the principal quantum numbers $n$ are completely excluded from correlation treatment; dotted line, ``$n \le 4$ cv only'': only core-valence correlations and relaxation effects for the $n \le 4$ core shells of Cs are incorporated (the corresponding results for the smaller Cs cores are practically indiscernible from those of the full correlation treatment). }
\label{fig:core}
\end{figure}

\section{Conclusions}\label{sec:conclusion}

The relativistic four-component Fock space coupled cluster method was used to calculate both diagonal and off-diagonal matrix elements of the magnetic dipole hyperfine interaction operator within the finite-field approach. The numerical stability of the solutions of the FS-RCC equations was ensured by the imaginary shift simulation technique. From the results for the dissociation limit, the overall error of our estimates should be about 3-5\%.

Our results demonstrate that the HFS coupling matrix elements in KCs are nearly constant for large and intermediate internuclear distances ($R > 8$ \AA, i.\ e. including the area of small separations of the $X^1\Sigma^+$ and $a^3\Sigma^+$ states in energy). This finding justifies the use of unscaled atomic hyperfine structure data for modeling the Feshbach resonances in alkali diatomics. For smaller $R$ the variation of HFS matrix elements become significant; the deviation of these entities from the corresponding atomic values reaches ca. 15\% in the vicinity of equilibrium distances for the ground and strongly bound excited states. A detailed study of the effect of the core relaxation, core-core and core-valence electron correlations on the computed matrix elements was performed. Core-core correlations practically do not contribute to the matrix elements under study. Rather significant (up to $\sim 10$\%) core relaxation and core-valence correlation contributions were found nearly independent on $R$. Thus the essential information on the variation of the HFS matrix elements upon chemical bond formation (``chemical shift'') can be obtained via relatively inexpensive calculations within the frozen core approximation.

\section{Acknowledgements}

Authors are grateful to E. A. Bormotova and A. V. Stolyarov for fruitful discussions.

This work has been carried out using computing resources of the federal collective usage centre Complex for Simulation and Data Processing for Mega-science Facilities at NRC ``Kurchatov Institute'', \url{http://ckp.nrcki.ru/}, and computers of Quantum Chemistry Lab at NRC ``Kurchatov Institute'' -- PNPI.

The reported study was funded by RFBR according to the research project \textnumero19-33-50122\textbackslash19.

\section{Appendix A}

Although the strategy of solving the relativistic many-electron problem generally followed that described in Refs.~\cite{Zaitsevskii2017,Zaitsevskii2018a}, we employed a somewhat different scheme to suppress numerical instabilities which is described below.

Due to the model space completeness, the straightforward application of the conventional FS-RCC formulation in wide ranges of internuclear separations is blocked by numerical instabilities caused by intruder states (see e.g. Refs.~\cite{Evangelisti1987,Zaitsevskii2017} for detailed discussions). The presence of intruder states manifests itself as the appearance of poles of the $D_K^{-1}$ factors arising from small or zero energy denominators in FS-RCC amplitude equations at certain nuclear geometries ($D_K$ stands for the energy denominator for the excitation $K$ \cite{Zaitsevskii2017}).

Similarly to the technique used in our previous work~\cite{Zaitsevskii2017}, we replace the original energy denominators $D_K$ in the FS-RCC amplitude equations by their shifted counterparts in such a way that numerical instabilities are suppressed without significant distortions of the description of low-lying states. However, the particular Ans\"{a}tze for the dependencies of real shifts on the original denominator values proposed in Refs.~\cite{Eliav2005,Zaitsevskii2017,Zaitsevskii2018a} can give rise to non-negligible systematic errors and rather strong dependencies of results on the particular choice of shift parameters. Though this problem is efficiently solved through extrapolating to the zero shift limit \cite{Eliav2005,Zaitsevskii2018a}, a less cumbersome approach is desirable for molecular excited state calculations in wide areas of nuclear configurations.

To construct a potentially better Ansatz, let us turn to the idea of intruder state avoidance by means of imaginary denominator shifts~\cite{Surjan1996,Forsberg1997}. The imaginary analog of the denominator shift proposed in Ref.~\cite{Zaitsevskii2017} would imply the replacement
\begin{equation}
D_K\;\Longrightarrow\;D^C_K(m)=D_K+s_K\,
\left(\displaystyle\frac{|s_K|}{\left| D_K+s_K i \right|}\right)^{m} i,
\label{imags}
\end{equation}
where the $s_K$ parameters (shift amplitudes) are real; these parameters are supposed to have the same value for all classes of excitations. The non-negative integer $m$ (attenuation parameter) is common for all excitations $K$. The imaginary part of $D_K^C$ tends to $s_K$ for nearly zero denominators $D_K$ and becomes small for ''well-defined'' (large negative or positive) $D_K$. The larger the $m$ value, the smaller the shifting of ``well-defined'' energy denominators. The essential advantage of imaginary shifts compared to their real counterpart introduced in Ref.~\cite{Zaitsevskii2017} arises from the fact that $D_K^C(m)^{-1}$ with any non-zero $s_K$ have no poles on the real axis regardless the behaviour of $D_K$ as functions of the nuclear geometry. Note, however, that the straightforward application of (\ref{imags}) could lead to non-physical imaginary contributions to effective Hamiltonian eigenvalues; furthermore, one would lose the possibility to perform FS-RCC calculations of diatomic (or other highly symmetric) molecules in real numbers. Therefore instead of using Eg.~(\ref{imags}), we replace the original energy denominators $D_K$ in the FS-RCC amplitude equations by real shifted denominators $D^{\prime}_K(m)$
\begin{equation}
D_K\;\Longrightarrow\;D^{\prime}_K(m)=D_K + \frac{s_K^2}{D_K}\left(\displaystyle\frac{s_K^2}{  D_K^2+s_K^2}\right)^{m}
\label{imagreal}
\end{equation}
satisfying the requirement
\begin{equation}
\frac{1}{D^{\prime}_K(m)}={\rm Re}\left(\frac{1}{D_K^C(m)}\right)
\end{equation}
and thus simulating in a sense the modification (\ref{imags}) without leaving the real axis.

The modification (\ref{imagreal}) is a straightforward generalization of the ``imaginary shifting'' widely used in second-order multireference many-body perturbation theories \cite{Forsberg1997,Witek2002}. Assuming $m=0$ in Eq.~(\ref{imagreal}), one immediately arrives at the expression for energy denominators introduced in the cited works. The increase of $m$ should reduce the errors introduced through the replacement of original denominators by shifted ones but can slow down (and finally ruin) the convergence of the iterative solution of the modified cluster equations in the presence of intruder states. A series of single-point calculations demonstrated that the use of a larger $m$ ($m\ge{}2$) is essential for suppressing the dependence of the resulting energies on the choice of shift amplitudes $s_K$. All results presented in Sec.~\ref{sec:results} were obtained with shift parameters $s=-0.2$ a.u. and $m=3$ in the $(0h,2p)$ Fock space sector; FS-RCC equations for the $(0h,0p)$ and $(0h,1p)$ were not modified.