text:
\section{Introduction}\label{sec1}\setcounter{equation}{0}

Probabilistic limit theorems for dynamical systems and Markov chains is a well studied topic.
One  way to derive such results is relying on some quasi-compactness (or spectral gap)
of an appropriate transfer or Markov operator, together with a suitable perturbation theorem
(see \cite{Neg1}, \cite{Neg2}, \cite{GH} and \cite{HH}). This  quasi-compactness  can often be verified only via
an appropriate Ruelle-Perron-Frobenius (RPF) theorem, which is the main key for thermodynamic formalism type constructions.
Probabilistic limit theorems for
random dynamical systems  and Markov chains in random dynamical environments were also studied in the literature (see, for instance, \cite{Kifer-1996}, \cite{Kifer-1998}, \cite{book}, \cite{Aimino} ,\cite{drag} and references therein). In these circumstances, the probabilistic behaviour of the appropriate process is determined by compositions of random operators, and not of a single operator, so no spectral theory can be exploited, and instead, many of these results rely on an appropriate version
of the RPF theorem for random operators. For instance, the large deviations theorem in \cite{Kifer-1996} rely on
such results for real operators, and
the conditions guaranteeing that the central limit theorem in \cite{Kifer-1998}  holds true
 can be often verified by an RPF theorem for random real operators. Relying on certain contraction properties of random complex transfer and Markov operators, with respect to a complex version of the Hilbert projective metric due to H.H. Rugh \cite{Rug} (see also \cite{Dub1} and \cite{Dub2}), we proved in \cite{book} an RPF theorem for random complex operators and presented the appropriate random complex thermodynamic formalism type constructions,
which was one of the main keys in the proof of  versions of the Berry-Esseen theorem and the local central limit
theorem (LCLT) for certain processes in random dynamical environment, and in the proof of some nonconventional LCLT for dynamical systems (see Chapters 2 and 7 of \cite{book}).
 In this paper we  use this RPF theorem in order to derive additional limit theorems, as described in the following paragraphs.

Let $({\Omega},{\mathcal F},P,{\theta})$ be an invertible measure preserving system, let ${\mathcal X}$ be a compact metric space, ${\mathcal E}_{\omega}\subset{\mathcal X}$ be a measurable family of compact subsets and $u_{\omega}:{\mathcal E}_{\omega}\to{\mathbb R}$ be a random function, where ${\omega}\in{\Omega}$. Let $T_{\omega}:{\mathcal E}_{\omega}\to{\mathcal E}_{{\theta}{\omega}}$
be a random distance expanding map  and let
$\xi_0^{{\omega}},\xi_1^{{\theta}{\omega}},\xi_2^{{\theta}^2{\omega}},...$ be a Markov chain in the random dynamical environment $({\Omega},{\mathcal F},P,{\theta})$. In this paper we consider sequences of random variables having either the form  $S_n^{\omega}=\sum_{j=0}^{n-1}u_{{\theta}^j{\omega}}\circ T_{\omega}^j(x)$ or the form ${\mathcal S}_n^{\omega} =\sum_{j=0}^{n-1}u_{{\theta}^j{\omega}}(\xi_j^{{\theta}^j{\omega}})$, where
$T_{\omega}^j=T_{{\theta}^{j-1}{\omega}}\circ\cdots\circ T_{{\theta}{\omega}}\circ T_{\omega}$ and $x$ is distributed according to a special (Gibbs) measure $\mu_{\omega}$.
Recall the following result. Let $X_1,X_2,X_3,...$ be a sequence of centered, independent and identically distributed random variables with finite moments of all orders, and set ${\sigma}^2={\mathbb E} X_1^2$ and $\zeta={\mathbb E} X_1^3$.
Then  for each integer  $k\geq 2$ the limit
\[
{\gamma}_k=\lim_{n\to\infty}n^{-[\frac k2]}{\mathbb E}\big(\sum_{j=1}^n X_j\big)^k
\]
of the $k$-th normalized moment
exists. Moreover, ${\gamma}_2={\sigma}^2$ and ${\gamma}_3=\zeta$, and when
 $k$ is even we have ${\gamma}_k=C_k{\sigma}^k$, where $C_k=2^{-\frac k2}(\frac{k}{2}!)^{-1}k!$, while for odd $k$'s larger than $2$ we have
${\gamma}_k=D_k{\sigma}^{k-3}\zeta$, where $D_k=\frac{k!}{3!}2^{-\frac12(k-3)}(\frac{k-3}{2}!)^{-1}$. This result is a direct consequence of the multinomial theorem, and we refer the readers to Lemma 4.1 and Corollary 4.2 in \cite{Liverani} for a generalization to several weakly dependent sequences of random variables.
In this paper, we will prove generalize these results for the sequences  $S_n^{\omega}$ and ${\mathcal S}_n^{\omega}$ defined above by showing that the $k$-th normalized moments of the random variables $S_n^{\omega}$ and ${\mathcal S}_n^{\omega}$
behave  like powers of ergodic averages of the second and third derivative at $0$ of some pressure function, which will imply that the asymptotic moments exist, satisfy the desired relations and that they do not depend on ${\omega}$ when $({\Omega},{\mathcal F},P,{\theta})$ is ergodic.
The fact that the second (normalzied) moments behave like an ergodic average was proved in \cite{Kifer-1998} in a more general setup, and here we show that this is true (in our context) also for all the higher moments.

When $({\Omega},{\mathcal F},P,{\theta})$ satisfies some mixing conditions then we obtain almost sure convergence rate of order $n^{-\frac12}\ln n$ towards the asymptotic moments, which allows us to obtain almost optimal convergence rate (of the same order) in the corresponding quenched central limit theorem for the sequences $n^{-\frac12}S_n^{\omega}$ and $n^{-\frac12}{\mathcal S}_n^{\omega}$. The arguments in the proof  of the above convergence rate also yields that for sufficiently well mixing base maps ${\theta}$ the partial sums generated by the random RPF triplets (corresponding to an appropriate parametric family of random transfer operators) and  all of their parametric derivatives satisfy several probabilistic limit theorems (see Remark \ref{RPF trip LimThms}), which also seems to be a new result.

Let $X_1,X_2,X_3,...$ be a sequence of centered  iid random variables with finite moments. Then,
under certain assumptions on the behavior of the characteristic function of $X_1$, a classical result (see \cite{Feller}) states that (Edgeworth) expansions of the form
\[
\sup_{s\in{\mathbb R}}\big|\sqrt{2\pi}P(\sum_{j=1}^n X_j\leq{\sigma}\sqrt n s)-\int_{-\infty}^s e^{-\frac{t^2}{2}}dt-\sum_{j=1}^d n^{-\frac{j}2}P_j(s)e^{-\frac{s^2}2}\big|=o(n^{-\frac{d}2})
\]
hold true for some polynomials $P_1,P_2,...,P_d$.
Edgeworth expansion have been obtained  in \cite{Neg2} for some classes of Markov chains and in \cite{Coelho} and \cite{Liverani}  for deterministic dynamical systems (and Markov chains), and in this paper
we will also obtain Edgeworth type expansions for the  processes $S_n^{\omega}$ and ${\mathcal S}_n^{\omega}$  described earlier,
with random polynomials $P_i=P_{i,{\omega},n}$ whose degree depends only on $i$, and the coefficients of $P_{i,{\theta}^{-n}{\omega},n}$
converge as $n\to\infty$.
When the dynamical environment $({\Omega},{\mathcal F},P,{\theta})$ satisfies certain mixing conditions we also obtain certain convergence rate in the later limits. In general, in order to obtain Edgeworth type expansions by analysing the asymptotic behaviour of the characteristic functions
$\varphi_n(t)={\mathbb E} e^{it S_n}$ of an underlying sequence $S_n,\,n\geq1$ of random variables, it is also necessary to have estimates on the decay of these functions (as $n\to\infty$) for large $t$'s, which is beyond the applications of the random complex
thermodynamic formalism (the decay estimates sufficient for the
existence of Edgeworth expansions in the deterministic case are discussed
in \cite{Coelho} and \cite{Liverani}). Here ${\mathbb E}$ denotes expectation with respect to $P$.
In the deterministic case, some of these estimates were obtained using the spectral theory of appropriate perturbations of the original transfer or Markov operator. In the circumstances of this paper there is no single operator, but a random family of operators, and we will use the  approach taken in \cite{book} in order to obtain the appropriate estimates in the case of iterates of random operators.

The paper is organized as follows. In Section \ref{sec2} we  formulate our results for random variables $S_n^{\omega}$
whose characteristic function is related to certain random operators ${\mathcal A}_{it}^{\omega},\,t\in{\mathbb R}$ which  satisfy some complex thermodynamic formalism type assumption, together with certain decay rates of the norms of the iterates
${\mathcal A}_{it}^{{\omega},n}={\mathcal A}_{it}^{{\theta}^{n-1}{\omega}}\circ{\mathcal A}_{it}^{{\theta}^{n-2}{\omega}}\circ\cdots\circ{\mathcal A}_{it}^{\omega}$ of these operators for large $t$'s.
In Section \ref{MomSec} we  prove the results concerning the normalized asymptotic moments and in Section \ref{EdgeSec} we  prove the results concerning Edgeworth expansions. Section \ref{secEXM} is devoted to examples of random dynamical systems $T_{\omega}$ and Markov chains $\xi_0^{\omega},\xi_1^{{\theta}{\omega}},\xi_2^{{\theta}^2{\omega}},...$ in random environments satisfying our assumptions.

\section{Preliminaries and main results}\setcounter{equation}{0}\label{sec2}
Our setup consists of a complete probability space
$({\Omega},{\mathcal F},P)$ together with an invertible $P$-preserving transformation ${\theta}:{\Omega}\to{\Omega}$,
of a compact metric space $({\mathcal X},\rho)$ normalized in size so that
$\text{diam}{\mathcal X}\leq 1$ together with the Borel ${\sigma}$-algebra ${\mathcal B}$, and
of a set ${\mathcal E}\subset{\Omega}\times {\mathcal X}$ measurable with respect to the product ${\sigma}$-algebra ${\mathcal F}\times{\mathcal B}$ such that the fibers ${\mathcal E}_{\omega}=\{x\in {\mathcal X}:\,({\omega},x)\in{\mathcal E}\},\,{\omega}\in{\Omega}$ are compact. The latter yields (see \cite{CV} Chapter III) that the mapping ${\omega}\to{\mathcal E}_{\omega}$ is measurable with respect to the Borel ${\sigma}$-algebra induced by the Hausdorff topology on the space ${\mathcal K}({\mathcal X})$ of compact subspaces of ${\mathcal X}$ and the distance function $\rho(x,{\mathcal E}_{\omega})$ is measurable in ${\omega}$ for each $x\in {\mathcal X}$.
Furthermore, the projection map $\pi_{\Omega}({\omega},x)={\omega}$ is
measurable and it maps any ${\mathcal F}\times{\mathcal B}$-measurable set to an
${\mathcal F}$-measurable set (see ``measurable projection" Theorem III.23 in \cite{CV}).

Let $(H_1^{\omega},\|\cdot\|_1)$ be a Banach space of complex valued functions on ${\mathcal E}_{\omega}$, containing the constant functions, so that $\|g\|_1\geq\sup|g|$ for any $g\in H_1^{\omega}$, and $\|\textbf{1}\|_1=1$,
where $\textbf{1}$ denotes the function taking the constant value $1$ (regardless of its domain).
Denote by $(H_1^{\omega})^*$ the space of all continuous linear functionals on $(H_1^{\omega},\|\cdot\|_1)$, equipped with the operator norm $\|\cdot\|_1$.
Let ${\mathcal A}_z^{\omega}:H_1^{\omega}\to H_1^{{\theta}{\omega}},\,z\in{\mathbb C}$ be a family of continuous linear operators, $S_n^{\omega}:{\mathcal E}_{\omega}\to{\mathbb R},\,n\geq1$ be a sequence of Borel measurable
functions and $\mu_{\omega}$ be a random probability measure on ${\mathcal E}_{\omega}$
so that $P$-a.s. we have $\mu_{\omega}(S_n^{\omega})=0$, $\|S_n^{\omega}\|_{L^k({\mathcal E}_{\omega},\mu_{\omega})}<\infty$ for any $k\geq1$
and  for any $z\in{\mathbb C}$ and $n\geq1$,
\[
\mu_{\omega}(e^{zS_n^{\omega}})=\mu_{{\theta}^n{\omega}}({\mathcal A}_z^{{\omega},n}\textbf{1})
\]
where ${\mathcal A}_z^{{\omega},n}={\mathcal A}_z^{{\theta}^{n-1}{\omega}}\circ{\mathcal A}_z^{{\theta}^{n-2}{\omega}}\circ\cdots\circ{\mathcal A}_z^{\omega}$.
 In this paper we will study certain (asymptotic) properties the distribution of
$S_n^{\omega}(x)$, when $x$ is distributed according to $\mu_{\omega}$,  and ${\omega}$ ranges over a set of full $P$-probability.
Our basic requirements from the operators are described in the following
\begin{assumption}\label{BaseAss}
There exist constants $\rho,C>0$ and ${\delta}\in(0,1)$  so that $P$-a.s. for any
$z\in B(0,\rho):=\{\zeta\in{\mathbb C}: |\zeta|<\rho\}$
there is a
triplet  consisting of a nonzero complex number ${\lambda}_{\omega}(z)$,  a function $h_{\omega}(z)\in H_1^{\omega}$
and a linear functional $\nu_{\omega}(z)\in(H_1^{\omega})^*$ such that $\nu_{\omega}(z)\textbf{1}=1$ and for any
 $g\in H_1^{\omega}$ and $n\geq1$,
\begin{equation}\label{ExpConvAss}
\left\|({\lambda}_{{\omega},n}(z))^{-1}{\mathcal A}_z^{{\omega},n}g-\nu_{\omega}(z)(g)h_{{\theta}^n{\omega}}(z)\right\|_1\leq C\|g\|_1{\delta}^n
\end{equation}
where  ${\lambda}_{{\omega},n}(z)=\prod_{j=0}^{n-1}{\lambda}_{{\theta}^j{\omega}}(z)$.
Moreover, the above triplet is measurable in ${\omega}$, analytic in $z$ and the random variables $\sup_{z: |z|<\rho}|{\lambda}_{\omega}(z)|$,  $\sup_{z: |z|<\rho}\|h_{\omega}(z)\|_1$ and $\sup_{z: |z|<\rho}\|\nu_{\omega}(z)\|_1$ are bounded. Furthermore,
$\nu_{\omega}(0)=\mu_{\omega}$, $h_{\omega}(0)=\textbf{1}$ and ${\lambda}_{\omega}(0)=1$.
\end{assumption}
We refer the readers' to Section \ref{secEXM} for several examples in which Assumption \ref{BaseAss} holds true.
Under Assumption \ref{BaseAss} we have $\nu_{\omega}(z)(h_{\omega}(z))=1$, ${\mathcal A}_z^{\omega} h_{\omega}(z)={\lambda}_{\omega}(z)h_{{\theta}{\omega}}(z)$
and $({\mathcal A}_z^{\omega})^*\nu_{{\theta}{\omega}}(z)={\lambda}_{\omega}(z)\nu_{\omega}(z)$ (see the beginning of Section \ref{MomSec}),
where $({\mathcal A}_z^{\omega})^*: (H_1^{{\theta}{\omega}})^*\to (H_1^{\omega})^*$ is the dual operator of ${\mathcal A}_z^{\omega}$ (i.e. the triplet from the assumption is a random complex Ruelle-Perron-Frobenius triplet, see Chapter 4 in \cite{book}).
In particular ${\mathcal A}_0^{\omega} \textbf{1}=\textbf{1}$.

Our first result is the following

\begin{theorem}\label{MomThm}

(i) By possibly decreasing $\rho$, $P$-a.s. we can define an analytic function $\Pi_{\omega}:B(0,\rho)\to{\mathbb C}$
so that $\Pi_{\omega}(0)=0$, ${\lambda}_{\omega}(z)=e^{\Pi_{\omega}(z)}$ and
$|\Pi_{\omega}(z)|\leq c_0$ for any $z\in B(0,\rho)$, for some constant $c_0$ which does not depend on ${\omega}$
and $z$.

(ii) Suppose that $\int S_n^{\omega}\mu_{\omega}=0$ for $P$-almost any ${\omega}$. For any $k\geq 2$ and $n\geq1$, let ${\gamma}_{k,n}$ be the random variable defined by
${\gamma}_{k,n}({\omega})=n^{-[\frac k2]}\int_{{\mathcal E}_{\omega}} (S_n^{\omega}(x))^kd\mu_{\omega}(x)$. Then, the limits
\[
{\gamma}_k({\omega})=\lim_{n\to\infty}{\gamma}_{k,n}({\omega})=\lim_{n\to\infty}{\gamma}_{k,n}({\theta}^{-n}{\omega})
\]
exist $P$-a.s. and in $L^p$, for any $1\leq p<\infty$. Moreover,
with ${\sigma}^2={\gamma}_2$, $\zeta={\gamma}_3$, $C_k=2^{-\frac k2}(\frac{k}{2}!)^{-1}k!$ and
$D_k=\frac{k!}{3!}2^{-\frac12(k-3)}(\frac{k-3}{2}!)^{-1}$,
when $k$ is even we have ${\gamma}_k=C_k{\sigma}^k$, while for odd $k$'s we have
${\gamma}_k=D_k{\sigma}^{k-3}\zeta$. Furthermore,
${\sigma}^2={\gamma}_2={\mathbb E}[\Pi_{\omega}''(0)|{\mathcal I}]$ and
$\zeta={\gamma}_3={\mathbb E}[\Pi_{\omega}'''(0)|{\mathcal I}]$, where ${\mathcal I}$ is the sub-$\sigma$ algebra of ${\mathcal F}$ containing only the ${\theta}$-invariant sets.
In particular,
when the measure preserving system $({\Omega},{\mathcal F},P,{\theta})$ is ergodic then all ${\gamma}_k({\omega})$'s do not depend on ${\omega}$.

(iii) In fact, we can write
\[
\int_{{\mathcal E}_{\omega}} (S_n^{\omega}(x))^kd\mu_{\omega}(x)=\mu_{{\theta}^n{\omega}}(h_{{\theta}^n{\omega}}^{(j)}(0))+\sum_{s=1}^{[\frac{k}{2}]}C_{{\omega},n}^{(k)}(s)n^s+{\varepsilon}_{{\omega},n}
\]
where $|{\varepsilon}_{{\omega},n}|\leq c_1\eta^n$ for some constants $c_1>0$ and $0<\eta<1$,
\[
\lim_{n\to\infty}C_{{\omega},n}^{(k)}([\frac{k}{2}])=\lim_{n\to\infty}C_{{\theta}^{-n}{\omega},n}^{(k)}([\frac{k}{2}])={\gamma}_k({\omega})
\]
 $C_{{\omega},n}^{(k)}([\frac{k}{2}]-1)=C_{{\omega},n}^{(k)}(1)=0$ and for all other $s$'s,
\begin{eqnarray*}
\lim_{n\to\infty}C_{{\theta}^{-n}{\omega},n}^{(k)}(s)=\\
\sum_{j=0}^k j!\binom{k}{j}\mu_{\omega}(h_{\omega}^{(k-j)}(0))\sum_{m_2,...,m_j}\big(\prod_{l=2}^j(i!)^{m_l}m_l!\big)^{-1}\prod_{l=2}^j({\mathbb E}_P[\Pi_{\omega}^{(l)}|{\mathcal I}])^{m_l}.
\end{eqnarray*}
Here the sum ranges over all possible choices of nonnegative integers $m_2,...,m_j$ so that $\sum_{l=2}^j m_l=s$
and $\sum_{l=2}^j lm_l=j$, and $f^{(l)}$ stands for the $l$-th derivatives of a function $f$ on $B(0,\rho)$ which takes values at some Banach space. Moreover, closed formulas for the $C_{{\omega},n}^{(k)}$'s can be recovered from the proof.
\end{theorem}
As in \cite{Liverani}, the method used in the proof of Theorem \ref{MomThm} yields that the limits ${\gamma}_k$ exist also without the assumption that $\int S_n^{\omega} d\mu_{\omega}=0$, but then the relations between the ${\gamma}_k$'s are more complicated and will not be discussed in this paper.

Now we will describe our results concerning  Edgeworth type expansions.
For any $r>0$ consider the following assumptions

\begin{assumption}\label{A3}
(i) There exists a random variable $a_{\omega}$ so that  $P$-a.s. we have
$\|{\mathcal A}_{it}^{{\omega},n}\|_1\leq a_{\omega}$ for any $t\in(-\rho,\rho)$, where $\rho$ comes from Assumption \ref{BaseAss}.

(ii) For any compact set $J\subset{\mathbb R}\setminus\{0\}$, $P$-.a.s for any $n\geq1$ we have
\[
\lim_{n\to\infty}n^{r+\frac12}\sup_{t\in J}\|{\mathcal A}^{{\omega},n}_{it}\|_1=0.
\]
\end{assumption}

\begin{assumption}\label{A4}
There exists a Banach space $(H_2^{\omega}, \|\cdot\|_2)$ containing the constant functions, such that $\|g\|_2\geq\sup|g|$ for any $g\in H_2^{\omega}$, and ${\mathcal A}_{it}^{\omega}$, $t\in{\mathbb R}$, can be (possibly)
extended to $H_2^{\omega}$ so that ${\mathcal A}_{it}$ is a continuous map between $H_2^{\omega}$ and $H_2^{{\theta}{\omega}}$. Moreover:

(i) There exists a random variable $R_{\omega}$ so that $P$-a.s.
 for any $t\in{\mathbb R}$ and $n\geq1$ we have $\|{\mathcal A}_{it}^{{\omega},n}\|_2\leq R_{\omega}(1+|t|)$.

(ii) There exist random variables $K_{\omega},D_{\omega},C_{\omega}>0$, $r_2({\omega})>0$  and $N({\omega})$ so that $P$-a.s. for any  $n\geq N({\omega})$,
\begin{equation}\label{EqA4}
\sup_{t: K_{\omega}\leq |t|\leq D_{\omega} n^{r}}\|{\mathcal A}_{it}^{{\omega},n}\|_2\leq C_{\omega} n^{-r_2({\omega})}.
\end{equation}
\end{assumption}

\begin{theorem}\label{EdgeThm}
Suppose that Assumptions \ref{BaseAss} and \ref{A3} hold true with $r=\frac12$. Set
$\Pi_{{\omega},n,2}=n^{-1}\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}''(0)$, and assume that ${\sigma}^2=\lim_{n\to\infty}\Pi_{{\omega},n,2}=\lim_{n\to\infty}\Pi_{{\theta}^{-n}{\omega},n,2}={\mathbb E}[\Pi_{\omega}''(0)]|{\mathcal I}]>0$.

(i) There exists a sequence of polynomials
$P_{{\omega},n,1}(s)=\sum_{j=0}^{m_1}a_{{\omega},n,j,1}s^j,\,n\geq1$ with random coefficients, whose degree $m_1$ does not depend on ${\omega}$ and $n$, so that $P$-almost surely we have
\begin{eqnarray*}
\lim_{n\to\infty}n^{\frac12}\sup_{s\in{\mathbb R}}\Big|\sqrt{2\pi}\mu_{\omega}\{x\in{\mathcal E}_{\omega}: S_n^{\omega}(x)\leq\sqrt n s\}-\\\frac{1}{\sqrt{\Pi_{{\omega},n,2}}}\int_{-\infty}^s e^{-\frac{t^2}{2\Pi_{{\omega},n,2}}}dt- n^{-\frac{1}2}P_{{\omega},n,1}(s)e^{-\frac{s^2}2}\Big|=0.
\end{eqnarray*}

(ii) Let $d\in{\mathbb N}$. If Assumptions \ref{A3} and  \ref{A4} hold true with $r=\frac{d}2$ and $({\Omega},{\mathcal F},P,{\theta})$ is ergodic, then for each $k\geq2$
there exists a sequences of polynomials
$P_{{\omega},n,k}(s)=\sum_{j=0}^{m_k}a_{{\omega},n,j,k}s^j,\,n\geq1$ with random coefficients and degree $m_k$ which depends only on $k$, so that $P$-almost surely we have
\begin{eqnarray*}
\lim_{n\to\infty}n^{\frac d2}\sup_{s\in{\mathbb R}}\Big|\sqrt{2\pi}\mu_{\omega}\{x\in{\mathcal E}_{\omega}: S_n^{\omega}(x)\leq\sqrt n s\}-\\\frac{1}{\sqrt{\Pi_{{\omega},n,2}}}\int_{-\infty}^s e^{-\frac{t^2}{2\Pi_{{\omega},n,2}}}dt-\sum_{j=1}^d n^{-\frac{j}2}P_{{\omega},n,j}(s)e^{-\frac{s^2}2}\Big|=0
\end{eqnarray*}

(iii) The coefficients $a_{{\omega},n,j,k}$ of the above
 polynomials are algebraic combinations of the derivatives of the functions $\mu_{{\theta}^n{\omega}}(h_{\omega}(z))$ and $n^{-1}\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}(z)$ at $z=0$. In particular, they
are uniformly bounded in ${\omega}$ and $n$, where ${\omega}$ ranges over a set of probability $1$,  and
$a_{{\theta}^{-n}{\omega},n,j,k}$ converges $P$-a.s. and in $L^p$, for any $p\in[1,\infty)$, as $n\to\infty$ towards a limit
$a_{{\omega},j,k}$
 The coefficients $a_{{\omega},n,j,k}$ and the
corresponding limits can be recovered from the proof.
\end{theorem}
Note that when (\ref{MomRateTemp}) below holds true  then  we obtain a converge rate of
order $n^{-\frac12}\ln n$ for all of the coefficients of
the $P_{{\theta}^{-n}{\omega},n,k}$'s.

\begin{remark}
A natural question is whether Theorem \ref{EdgeThm} (or a version of it with $S_n^{{\theta}^{-n}{\omega}}$ instead of $S_n^{\omega}$) can be obtained with polynomials which depend only on ${\omega}$ and not on $n$ and with ${\sigma}^2$ in place of $\Pi_{{\omega},n,2}$. Similar arguments to the ones in Remark 1.1 from \cite{Liverani} show that replacing $P_{{\omega},n,k}$ with a polynomial  $P_{{\omega},k}$ which does not depend on $k$ is possible only if $n^{\frac k2}|P_{{\omega},k}(x)-P_{{\theta},n,k}(x)|$ converges to $0$ as $n\to\infty$ uniformly in $x$, which means that the only candidate for such polynomials are the polynomials $P_{{\omega},k}$ whose coefficients are the limits of the coefficients of the $P_{{\omega},n,k}$'s (or the $P_{{\theta}^{-n}{\omega},n,k}$'s in the ``reverse version").
Proving such a result does not seem possible using only the converge of the coefficients of $P_{{\theta}^{-n}{\omega},n,k}$ even when ${\theta}$ satisfies some mixing conditions, since this convergence  is derived from the convergence of certain ergodic averages, which, in general, is slower than $n^{-\frac12}$ (even for independent summands). Note that even for Edgeworth expansions of order one (as well as for a Berry-Esseen theorem) it is not clear how to obtain the results with ${\sigma}^2$ in place of $\Pi_{{\omega},n,2}$, since this essentially requires better convergence rate of the variance of $n^{-\frac12}S_n^{\omega}$ towards ${\sigma}^2$. The difference between this variances and ${\sigma}^2$ behaves as an ergodic average, which leads again to the same problem discussed above (concerning the rate of converges of ergodic averages).
\end{remark}

\subsection{Almost sure convergence rates under mixing conditions}\label{MixSec}
When  $({\Omega},{\mathcal F},P,{\theta})$ satisfies certain  mixing assumptions, then, in Section \ref{MomRateSec}
we obtain  almost sure converges rate of the form
\begin{equation}\label{MomRateTemp}
\left|n^{-[\frac k2]}\int_{{\mathcal E}_{\omega}} (S_n^{\omega}(x))^kd\mu_{\omega}(x)-{\gamma}_k\right|\leq R_{{\omega},k} n^{-\frac12}\ln n
\end{equation}
where $R_{{\omega},k}$ is some random variable.
Under certain mixing condition the authors of \cite{Ste} obtained  converges rates of order
$(\ln n)^{\frac 32+{\delta}_0}n^{-\frac12},\,{\delta}_0>0$ for $k=2$, and in our setup we obtain such rates for all $k$'s, with
${\delta}_0=-\frac12$.
 Using (\ref{MomRateTemp}) with $k=2$, when ${\sigma}^2>0$ we
also derive in Section \ref{MomRateSec} almost optimal convergence rate in the central limit theorem of the form
\[
\sup_{s\in{\mathbb R}}\left|\mu_{\omega}\{x\in{\mathcal E}_{\omega}: S_n^{\omega}(x)\leq s\sqrt n\}-\frac1{\sqrt{2\pi{\sigma}^2}}\int_{-\infty}^s e^{-\frac{t^2}{2{\sigma}^2}}dt\right|\leq
c_{\omega} n^{-\frac12}\ln n.
\]

\section{Asymptotic moments: Proof of Theorem \ref{MomThm}}\label{MomSec}\setcounter{equation}{0}
First, since ${\lambda}_{\omega}(0)=1$, ${\lambda}_{\omega}(\cdot)$ is an analytic function and $|{\lambda}_{\omega}(z)|$ is bounded
uniformly in ${\omega}$ and $z$, where ${\omega}$ ranges over a set of probability $1$ and $z\in B(0,\rho)$, it is indeed possible to construct a function $\Pi_{\omega}(\cdot)$ which satisfies the
conditions stated in Theorem \ref{MomThm} (i), namely, an appropriate branch of the logarithm of ${\lambda}_{\omega}(z)$ can be defined on some deterministic neighborhood of $0$.

Next, we will show that  $({\lambda}_{\omega}(z),h_{\omega}(z),\nu_{\omega}(z))$ is an RPF triplet, namely that $P$-a.s. for any $z\in B(0,\rho)$,
\begin{equation}\label{RPF}
\nu_{\omega}(z)(h_{\omega}(z))=1,\,{\mathcal A}_z^{\omega} h_{\omega}(z)={\lambda}_{\omega}(z)h_{{\theta}{\omega}}(z)\,\,\text{ and }\,\,
({\mathcal A}_z^{\omega})^*\nu_{{\theta}{\omega}}(z)={\lambda}_{\omega}(z)\nu_{\omega}(z).
\end{equation}
Indeed,
in order to prove the second equality in (\ref{RPF}),
plug in ${\theta}^{-n}{\omega}$ in place of ${\omega}$ in (\ref{ExpConvAss}), and then choose $g=\textbf{1}$, in order to deduce that
\begin{eqnarray*}
{\mathcal A}_z^{\omega} h_{\omega}(z)=\lim_{n\to\infty}({\lambda}_{{\theta}^{-n}{\omega},n}(z))^{-1}{\mathcal A}_z^{\omega}({\mathcal A}_z^{{\theta}^{-n}{\omega},n}\textbf{1})\\=
{\lambda}_{\omega}(z)\lim_{n\to\infty}({\lambda}_{{\theta}^{-(n+1)}{\theta}{\omega},n+1}(z))^{-1}
{\mathcal A}_z^{{\theta}^{-(n+1)}{\theta}{\omega},n+1}={\lambda}_{\omega}(z)h_{{\theta}{\omega}}(z).
\end{eqnarray*}
Next, applying  (\ref{ExpConvAss}) with ${\theta}{\omega}$ in place of ${\omega}$ and with ${\mathcal A}_z^{\omega} g$ in place of $g$ we derive that for any $g\in H_1^{\omega}$,
\begin{eqnarray*}
\nu_{{\theta}{\omega}}(z)({\mathcal A}_z^{\omega} g)=\lim_{n\to\infty}\frac{{\mathcal A}_z^{{\theta}{\omega},n}({\mathcal A}_z^{\omega} g)}{{\lambda}_{{\theta}{\omega},n}(z)h_{{\theta}^{n+1}{\omega}}(z)}\\={\lambda}_{\omega}(z)\lim_{n\to\infty}\frac{{\mathcal A}_z^{{\omega},n+1}g}{{\lambda}_{{\omega},n+1}(z)h_{{\theta}^{n+1}{\omega}}(z)}=
{\lambda}_{\omega}(z)\nu_{\omega}(z)(g)
\end{eqnarray*}
which implies that $({\mathcal A}_z^{\omega})^*\nu_{{\theta}{\omega}}(z)={\lambda}_{\omega}(z)\nu_{\omega}(z)$. Note that the function $h_{\omega}(z)$ does not vanishes
when $z$ lies in a deterministic neighborhood of $0$ since $h_{\omega}(0)\equiv1$, the random variable $\sup_{|z|<\rho}\|h_{\omega}(z)\|_1$ is bounded and the map $z\to h_{\omega}(z)$ is analytic.
Finally, since
\[
h_{\omega}(z)=\lim_{n\to\infty}\frac{{\mathcal A}_z^{{\theta}^{-n}{\omega},n}\textbf{1}}{{\lambda}_{{\theta}^{-n}{\omega},n}(z)}
\]
and
\[
{\lambda}_{{\theta}^{-n}{\omega},n}(z)=
{\lambda}_{{\theta}^{-n}{\omega},n}(z)\nu_{{\theta}^{-n}{\omega}}(z)\textbf{1}=
\nu_{\omega}(z)({\mathcal A}_z^{{\theta}^{-n}{\omega},n}\textbf{1})
\]
we obtain that $\nu_{\omega}(z)(h_{\omega}(z))=1$. \qed

Now we will start proving Theorem \ref{MomThm}.
First, it follows from (\ref{RPF}) that ${\lambda}_{{\omega},n}(z)=\nu_{{\theta}^n{\omega}}(z)({\mathcal A}_z^{{\omega},n}\textbf{1})$, and
so, since $\nu_{\omega}(z)\textbf{1}=1$ and $\mu_{\omega}(e^{zS_n^{\omega}})=\mu_{{\theta}^n{\omega}}({\mathcal A}_z^{{\omega},n}\textbf{1})$,
\[
{\lambda}_{{\omega},n}'(0)=\nu_{{\theta}^n{\omega}}'(0)\textbf{1}+
\frac{ d\nu_{{\theta}^n{\omega}}(0)({\mathcal A}_z^{{\omega},n}\textbf{1})}{dz}\Big|_{z=0}=\frac{d\mu_{\omega}(e^{zS_n^{\omega}})}{dz}\Big|_{z=0}=\mu_{\omega}(S_n^{\omega})=0.
\]
Next, using Assumption \ref{BaseAss} for any $n\geq1$ we can write
\begin{equation}\label{I}
\mu_{\omega}(e^{zS_n^{\omega} })=\mu_{{\theta}^n{\omega}}({\mathcal A}_{z}^{{\omega},n}\textbf{1})=
{\lambda}_{{\omega},n}(z)\big(\mu_{{\theta}^n{\omega}}(h_{{\theta}^n{\omega}}(z))+{\delta}_{{\omega},n}(z)\big)
\end{equation}
where $|{\delta}_{{\omega},n}(z)|\leq C{\delta}^n$. Since ${\delta}_{{\omega},n}(\cdot)$ is an analytic function, it follows from the Cauchy integral formula that
for any $k\geq 1$ there exists a constant $Q_k$ so that
$|{\delta}_{{\omega},n}^{(k)}(z)|\leq Q_k{\delta}^n$ for any $z\in B(0,\frac \rho 2)$. Since ${\lambda}_{{\omega},n}(0)=1$ and $\sup_{z\in B(0,\rho)}|{\lambda}_{\omega}(z)|$ is a bounded random variable, the analyticity of ${\lambda}_{\omega}(\cdot)$ implies that for any $k\geq1$ we have
$|{\lambda}_{{\omega},n}^{(k)}(0)|\leq H_k n^k$, for some constant $H_k$ which does not depend on $n$ and ${\omega}$. Therefore, differentiating $k$ times both sides of (\ref{I})  at $z=0$ yields that
\begin{equation}\label{II}
\mu_{\omega}(S_n^{\omega} )^k=\sum_{j=0}^k\binom{k}{j}{\lambda}_{{\omega},n}^{(j)}(0)\mu_{{\theta}^n{\omega}}(h_{{\theta}^n{\omega}}^{(k-j)}(0))+d_{{\omega},n}
\end{equation}
for some random variable $d_{{\omega},n}$ so that
$|d_{{\omega},n}|\leq Ac^n$ for some $A\geq1$ and $c\in(0,1)$, which do not depend on ${\omega}$.
Since ${\lambda}_{{\omega},n}(z)=e^{\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}(z)}$ and $\Pi_{\omega}(0)=1$, it follows from the Fa\'a di Bruno formula that for any $j$,
\begin{equation}\label{***}
{\lambda}_{{\omega},n}^{(j)}(0)=j!\sum_{s=1}^{[\frac j2]}n^s\sum_{(m_2,...,m_k)\in {\Gamma}_{j,s}}
\big(\prod_{l=2}^j(l!)^{m_l}m_l!\big)^{-1}\prod_{l=2}^{j}\Big(n^{-1}\sum_{i=0}^{n-1}\Pi_{{\theta}^i{\omega}}^{(l)}(0)\Big)^{m_l}
\end{equation}
where ${\Gamma}_{j,s}$ is the set of all $j-1$-tuples $(m_2,...,m_j)$ of nonnegative integers
so that $\sum_l lm_l=j$ and $\sum_{l}m_l=s$,
and we took into account that
\[
\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}'(0)={\lambda}_{{\omega},n}'(0)=\mu_{{\omega}}(S_n^{\omega})=0.
\]
Observe that when $s=[\frac j2]$
and $j$ is even  we have ${\Gamma}_{j,s}=\{(\frac j2,0,0,...,0)\}$, while for odd $j$'s we have  ${\Gamma}_{j,s}=\{(\frac{j-3}2,1,0,...,0)\}$. Since $\Pi_{\omega}(z)$ is analytic in $z$ and uniformly bounded in ${\omega}$ and $z$
(where ${\omega}$ ranges over a set of probability $1$ and $z\in B(0,\rho)$),
for each $j$ there exists a constant $L_j$ so that
$\sup_{z\in B(0,\frac \rho 2)}|\Pi_{\omega}^{(j)}(z)|\leq L_j$ for $P$-almost any ${\omega}$.
When $k$ is even,
by considering the case when $j=k$, we conclude that
 there exist constants $R_k, k\geq2$ so that for even $k$'s we have
\begin{equation}\label{Rate even}
\left|\frac{\mu_{\omega}(S_n^{\omega} )^k}{n^{[\frac k2]}}-C_k\Big(n^{-1}\sum_{i=0}^{n-1}\Pi_{{\theta}^i{\omega}}^{(2)}(0)\Big)^{\frac k2}\right|\leq R_kn^{-1}
\end{equation}
where $C_k=2^{-\frac k2}(\frac k2!)^{-1}k!$. Now, suppose that $k$ is odd. Then the dominating terms come from considering $j=k$ and $j=k-1$. Differentiating both sides of the equality $\nu_{\omega}(z)(h_{\omega}(z))=1$ and plugging in $z=0$ we derive that
\[
\mu_{\omega}(h_{\omega}'(0))=-\nu_{\omega}'(0)\textbf{1}=0
\]
where we used that $h_{\omega}(0)\equiv\textbf{1}$ and that
$\nu_{\omega}(z)\textbf{1}=1$ for any $z$. Therefore, we actually only need to consider the case when $j=k$ also for odd $k$'s, and so, with $D_k=\frac{k!}{3!}2^{-\frac12(k-3)}(\frac{k-3}{2}!)^{-1}$
 for odd $k$'s large than $2$ we have
\begin{equation}\label{Rate odd}
\left|\frac{\mu_{\omega}(S_n^{\omega} )^k}{n^{[\frac k2]}}-D_k\Big(n^{-1}\sum_{i=0}^{n-1}\Pi_{{\theta}^i{\omega}}^{(2)}(0)\Big)^{\frac{k-3}2}\cdot
\Big(n^{-1}\sum_{i=0}^{n-1}\Pi_{{\theta}^i{\omega}}^{(3)}(0)\Big)\right|\leq R_kn^{-1}.
\end{equation}
Since  the $L^\infty$-norms of
$n^{-1}\sum_{i=0}^{n-1}\Pi_{{\theta}^{i}{\omega}}^{(l)}(0)$ are uniformly bounded in $n$ (for each $l$),
all the positive integer powers of these averages
converge almost surely and in $L^p$, for any $p\in [1,\infty)$, towards the appropriate power of
$E[\Pi_{\omega}^{(l)}(0)|{\mathcal I}]$.
Combining this with the above estimates we complete the proof of Theorem \ref{MomThm}.

\subsection{Convergence rates towards the asymptotic moments under mixing conditions}\label{MomRateSec}

We assume here that ${\Omega}={\mathcal Y}^z$ for some measurable space ${\mathcal Y}$, and that $({\Omega},{\mathcal F},P,{\theta})$ is the shift system
generated by a ${\mathcal Y}$-valued stationary process $\{\xi_n:\,n\in{\mathbb Z}\}$. Consider the situation when
${\mathcal A}_z^{\omega}$ depends only on
the $0$-th coordinate of ${\omega}$. Then  by (\ref{ExpConvAss}), when $z$ lies in some deterministic neighborhood of the origin, the functional
$\nu_{\omega}(z)$ can be
approximated in $L^\infty$ exponentially fast in $n$ by functions of the first $n$-th coordinates. Indeed, since $h_{\omega}(0)=\textbf{1} $, $h_{\omega}(z)$ is analytic and $\sup_{z\in B(0,\rho)}\|h_{\omega}(z)\|_1$ is bounded, we derive that there exist constants $r_0>0$ and ${\delta}_0>0$ so that $P$-a.s. for any $z\in B(0,r_0)$ we have $\inf|h_{\omega}(z)|>{\delta}_0$.
Since  $\nu_{\omega}(z)\textbf{1}=1$, it follows  now from (\ref{ExpConvAss}) that for any $z\in B(0,r_0)$ and $g\in H_1^{\omega}$ so that $\|g\|_1\leq1$,
\[
\Big\|\nu_{{\omega}}(z)(g)-\frac{{\mathcal A}^{{\omega},n}_zg}{{\mathcal A}^{{\omega},n}_z \textbf{1}}\Big\|_1\leq A{\delta}^n
\]
where we used that $\nu_{{\omega}}(z)(g)=\frac{\nu_{{\omega}}(z)(g)}{\nu_{{\omega}}(z)\textbf{1}}$. Here $A>0$ is some constant and ${\delta}$ comes from (\ref{ExpConvAss}).
By (\ref{RPF}) we have ${\lambda}_{\omega}(z)=\nu_{\omega}(z)({\mathcal A}_z^{\omega} \textbf{1})$, and therefore $P$-a.s. for any $z\in B(0,r_0)$
and $n\geq 1$ we have
\[
|{\lambda}_{\omega}(z)-q_{n,z}({\omega}_0,{\omega}_1,...,{\omega}_{n-1})|\leq C{\delta}^n
\]
where $C>0$ is some constant, ${\omega}=\{{\omega}_s:\,s\in{\mathbb Z}\}$ and $q_{n,z}:{\mathcal Y}^n\to{\mathbb C}$ is a family of measurable functions which is uniformly bounded  in ${\omega}_0,{\omega}_1,...,{\omega}_{n-1}$,  $z$ and $n$.
Since ${\lambda}_{\omega}(0)=1$ and $\sup_{z\in B(0,\rho)}|{\lambda}_{{\omega}}(z)|$ is a bounded random variable, there exist constants $a,b,\rho>0$ so that $P$-a.s. for any $z\in B(0,\rho)$ we have
$|{\lambda}_{\omega}(z)|\in(a,b)$. Therefore, also $\Pi_{{\omega}}(z)$ can be approximated exponentially fast in the above sense
when $z$ belongs to a deterministic neighborhood of $0$.
Thus, using the identifications
\[
{\theta}^j{\omega}=\{\xi_{s+j}:\,s\in{\mathbb Z}\},\,j\in{\mathbb Z}
\]
we derive that
for any $k$ there exists a constant $A_k$ and  a uniformly bounded family of functions $g_{k,n}:{\mathcal Y}^n\to{\mathbb C}$ so that for any $j\geq0$ and $n\geq1$,
\begin{equation}\label{ExpAppPress}
\left\|\Pi_{{\theta}^j{\omega}}^{(k)}(0)-g_{k,n}(\xi_j,\xi_{j+1},...,\xi_{j+n-1})\right\|_{L^\infty({\Omega},{\mathcal F},P)}\leq A_k{\delta}^n.
\end{equation}

Next, recall that the $\phi$-mixing (dependence) coefficients associated with the sequence $\{\xi_n: n\in{\mathbb Z}\}$ are given by
\[
\phi(n)=\sup\big\{|P(B|A)-P(B)|: A\in{\mathcal F}_{-\infty,k},\, B\in{{\mathcal F}_{k+n,\infty}},\, k\in{\mathbb Z},\, P(A)>0\big\}
\]
where ${\mathcal F}_{-\infty,s}$ is the ${\sigma}$-algebra generated by $\{\xi_l,\,l\leq s\}$ and ${\mathcal F}_{s,\infty}$
is the ${\sigma}$-algebra generated by $\{\xi_l: l\geq s\}$ (for any integer $s$).
Recall  (see \cite{Br}, Ch. 4) also that $\phi(n)$ can be written as
\begin{equation}\label{Phi-Rel-StPaper}
\phi(n)=\frac 12\sup\big\{\|{\mathbb E}[g|{\mathcal F}_{-\infty,s}]-{\mathbb E} g\|_{\infty}\,:  g\in L^\infty({\Omega},{\mathcal F}_{s+n,\infty},P),\,
\|g\|_{\infty}\leq1,\,s\in{\mathbb Z}\big\}.
\end{equation}
We will obtain convergence rates towards the asymptotic moments under the assumption that
\[
S_\phi:=\sum_{n=1}^\infty\phi(n)<\infty.
\]
Under this assumption, the next step of the proof of (\ref{MomRateTemp}) is to approximate $\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}^{(k)}(0)$ by martingales in the $L^\infty$ norm. The following arguments are classical and are given here for readers' convenience.
For each $m\geq0$, let ${\mathcal F}_{0,m}$ be the ${\sigma}$-algebra generated by $\xi_0,\xi_1,...,\xi_m$.
For any $k\geq 2$ and $r\geq 1$ set
\[
Y_j^{(k,r)}=g_{j,k}(\xi_j,...,\xi_{j+r-1})-{\mathbb E} g_{j,k}(\xi_j,...,\xi_{j+r-1}),\,j\geq0
\]
and let  $M^{(r,k)}_n=\sum_{j=0}^{n-1}X_j^{(k,r)}$ be the martingale
(w.r.t the filtration $\{{\mathcal F}_{0,n-1+r}: n\geq 0\}$) whose differences are given by
\[
X_j^{(r,k)}=Y_j^{(k,r)}+\sum_{s\geq j+1}{\mathbb E}[Y_s^{(k,r)}|\xi_0,\xi_1,...,\xi_{j+r}]-
\sum_{s\geq j}{\mathbb E}[Y_s^{(k,r)}|\xi_0,\xi_1,...,\xi_{j-1+r}].
\]
Then by (\ref{Phi-Rel-StPaper}), similarly to the proof of Theorem 2.8 in \cite{Haf-MD},
we have $\|X_j^{(r,k)}\|_{L^\infty}\leq C(1+r+S_\phi)$ for any $j$, where $C$ is some constant, and, using also (\ref{ExpAppPress}) we derive that
\begin{eqnarray*}
\left\|\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}^{(k)}(0)-n\int\Pi_{{\omega}'}^{(k)}(0)dP({\omega}')-M_n^{(k,r)}\right\|_{L^\infty}\\
\leq\left\|\sum_{j=0}^{n-1}Y_j^{(k,r)}-M_n^{(k,r)}\right\|_{L^\infty}+nA_k{\delta}^r\leq
Q_k(1+r+S_\phi+n{\delta}^r):=R
\end{eqnarray*}
where $Q_k$ is some constant.
Applying the Azuma-Hoeffding inequality (see, for instance, page 33 in \cite{Mil})  with the martingales $M_n^{(k,r)}$ and $-M_n^{(k,r)}$, we obtain that
that for any $t\geq 0$,
\begin{eqnarray*}
P\big(\big|\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}^{(k)}(0)-n\int\Pi_{{\omega}'}^{(k)}(0)dP({\omega}')\big|\geq t+R\big)
\leq P(|M_n^{(k,r)}|\geq t)\leq P(M_n^{(k,r)}\geq t)\\+P(-M_n^{(k,r)}\geq t)
=P(e^{{\lambda}_t M_n^{(k,r)}}\geq e^{t{\lambda}_t})+(e^{{\lambda}_t\cdot(-M_n^{(k,r)})}\geq e^{t{\lambda}_t})\\
\leq e^{-t{\lambda}_t}\big({\mathbb E} e^{{\lambda}_t M_n^{(k,r)}}+{\mathbb E}  e^{-{\lambda}_t M_n^{(k,r)}}\big)
\leq 2 e^{-t{\lambda}_t}e^{{\lambda}_t^2\sum_{j=0}^{m-1}\|X_j^{(k,r)}\|_{L^\infty}}
\leq 2e^{-\frac{t^2}{4nC(1+r+S_\phi)}}
\end{eqnarray*}
where  ${\lambda}_t=\frac{t}{2nC(1+r+S_\phi)}$, and in the third inequality we used the Markov inequality.
Taking $r=r_n$ of logarithmic order in $n$, we derive  that there exist constants $a_k,C_k,c_k,d_k>0$ so that for any $n\geq2$ and ${\varepsilon}>0$,
\[
P\big\{{\omega}: |\Pi_{{\omega},n}^{(k)}(0)-n\int \Pi_{{\omega}'}^{(k)}(0)dP({\omega}')|\geq n{\varepsilon}+a_k\ln n \big\}\leq C_ke^{-c_k\frac{{\varepsilon}^2 n}{\ln n}}.
\]
By taking ${\varepsilon}=u n^{-\frac12}\ln n$ for a sufficiently large $u$, we derive from the Borel Cantelli Lemma that for any $k$, $P$-a.s for any sufficiently large $n$ we have
\begin{equation}\label{P rate}
\left|\frac{1}{n}\Pi_{{\omega},n}^{(k)}(0)-\int \Pi_{{\omega}'}^{(k)}(0)dP({\omega}')\right|\leq b_k n^{-\frac12}\ln n
\end{equation}
where $b_k>0$ is some constant which depends only on $k$.
It follows from (\ref{P rate}), (\ref{Rate even}) and (\ref{Rate odd}) that the convergence rate towards the
asymptotic moments is at most of order $n^{-\frac12}\ln n$ (i.e. that (\ref{MomRateTemp}) holds true).

Next, the arguments in the proof of Theorem 7.1.1 in \cite{book}  show  that when ${\sigma}^2>0$ then there exist random variables $c_{\omega}$ and $d_{\omega}$ so that $P$-a.s. for any $n\geq1$,
\begin{eqnarray*}
\sup_{s\in{\mathbb R}}\big|\mu_{\omega}\{x\in{\mathcal E}_{\omega}: S_n^{\omega}(x)\leq s\sqrt n\}-\frac1{\sqrt{2\pi{\sigma}^2}}\int_{-\infty}^s e^{-\frac{t^2}{2{\sigma}^2}}dt\big|\\\leq
c_{\omega} n^{-\frac12}+d_{\omega}\big|\frac 1n\mu_{\omega}(S_n^{\omega})^2-{\sigma}^2\big|.
\end{eqnarray*}
Therefore, under the above mixing conditions we derive almost optimal convergence rate in the central limit theorem of the form
\[
\sup_{s\in{\mathbb R}}\left|\mu_{\omega}\{x\in{\mathcal E}_{\omega}: S_n^{\omega}(x)\leq s\sqrt n\}-\frac1{\sqrt{2\pi{\sigma}^2}}\int_{-\infty}^s e^{-\frac{t^2}{2{\sigma}^2}}dt\right|\leq
q_{\omega} n^{-\frac12}\ln n
\]
where $q_{\omega}$ is some random variable (which in general may not even be integrable).

\begin{remark}
When $\{\xi_n:\,n\geq1\}$ is a stationary geometrically ergodic Markov chain then $\phi(n)$ converges exponentially fast to $0$ (see \cite{Br}, Theorem 21.1), and in particular $S_\phi<\infty$. Such exponential convergence occurs also when $\xi_n$ has form $\xi_n(x)=\xi_0\circ S^n(x)$, where $S$ is a two sided topologically mixing subshift of finite type, $x$ is distributed according to some  Gibbs measure (see \cite{Bow}) and $\xi_0$ is a function of a finite number of coordinates.
If $\{\zeta_n: n\geq0\}$ is a one sided stationary process then we can define $\xi_n=(\xi_{n,j})_{j=-\infty}^{\infty}=(...,\zeta_{n-1},\zeta_n,\zeta_{n+1},...)$, where $\zeta_n$ appears in the $0$-the coordinate. In these circumstances, we can write ${\omega}=({\omega}_j)_{j\in{\mathbb Z}}$ and ${\omega}_j=({\omega}_{j,k})_{k=-\infty}^\infty$.
When ${\mathcal A}_z^{\omega}$ depends only on ${\omega}_{0,0}$ then all the arguments from this section hold true if we assume that the sequence $\phi_1(n),\,n\geq1$ given by
\[
\phi_1(n)=\sup\big\{|P(B|A)-P(B)|: A\in{\mathcal G}_{-\infty,k},\, B\in{{\mathcal G}_{k+n,\infty}},\, k\in{\mathbb N},\, P(A)>0\big\}
\]
is summable, where ${\mathcal G}_{-\infty,s}$ is the ${\sigma}$-algebra generated by $\{\zeta_0,\zeta_1,...,\zeta_s\}$ and ${\mathcal G}_{s,\infty}$ is the ${\sigma}$-algebra generated by $\{\zeta_s,\zeta_{s+1},...\}$. Note also that the above results hold true also when
 $\xi_n$ has the form $\xi_n(x)=\xi_0\circ T^n(x)$, where $T$ is a Young tower (see \cite{Young1} and \cite{Young2}) whose tails decay to $0$ sufficiently fast,  $x$ is distributed according to an appropriate Gibbs measure and $\xi_0$ is a bounded random variables which is measurable with respect to the partition which defines the tower. Young towers do not seem to be $\phi$-mixing with respect to the natural filtrations generated by cylinder sets, but they  satisfy certain type of ``graded $\phi$-mixing" conditions which allows use to obtain the results above (see Section 7 in \cite{Haf-Arrays} for several mixing properties of Young towers which involve certain $\phi$-mixing coefficients).
\end{remark}

\begin{remark}\label{RPF trip LimThms}
The argument above show that the RPF triplets ${\lambda}_{\omega}(z), h_{\omega}(z)$ and $\nu_{\omega}(z)$ and all their derivatives with respect to $z$ can be approximated by functions of the random variables $\xi_0,\xi_1,...,\xi_n$ (or $\xi_{-n},...,\xi_0$ in the case of $ h_{\omega}(z)$) with an error term which converges to $0$ exponentially fast when $n\to\infty$. Therefore, applying classical results from probability theory (for stationary and mixing processes) we obtain that for any real $t$ with a sufficiently small absolute value, the partial sums $\sum_{j=0}^{n-1}X_{i,k}({\theta}^j{\omega})$ generated by the random variables $X_{1,k}({\omega})={\lambda}_{\omega}^{(k)}(t)$,  $X_{2,k}({\omega})=\big(h_{\omega}^{(k)}(t)\big)(x)$ or $X_{3,k}({\omega})=\big(\nu^{(k)}_{\omega}(t)\big)(g)$ satisfy exponential concentration inequalities, as well as various limit theorems such as the CLT. Here we assume that ${\mathcal E}_{\omega}={\mathcal X}$ does not depend on ${\omega}$, $x\in{\mathcal X}$ and $g$ is a H\"older continuous function on ${\mathcal X}$.
\end{remark}

\section{Edgeworth type expansions: Proof of Theorem \ref{EdgeThm}}\label{EdgeSec}\setcounter{equation}{0}
The proof of Theorem \ref{EdgeThm} relies on classical arguments involving Fourier transforms which were used
 successfully in the i.i.d. case and in the deterministic case (i.e. when $|{\Omega}|=1$),
see, for instance, \cite{Feller}, \cite{Liverani} and references therein. Most of the arguments in the proof of Theorem \ref{EdgeThm} are modifications of the arguments in \cite{Liverani}. When the arguments are exactly as in \cite{Liverani}, we will just refer the reader's  there. The main difficulty is to verify Assumptions \ref{BaseAss}, \ref{A3} and \ref{A4}, and this is postponed to Section \ref{secEXM}.

For reader's convenience, we first will describe the main idea behind using Fourier transforms in order to derive Edgeworth expansions.
Recall first that by the
Esseen inequality (see \cite{Lin}, Ch. 4.1), there exists an absolute constant $A$ so that for any  $T>0$, a distribution function $F$ and an integrable function $G:{\mathbb R}\to{\mathbb R}$ with bounded first derivative such that $\lim_{x\to\infty} G(x)=0$,
\begin{equation}\label{Esseen}
\sup_{x\in{\mathbb R}}|F(x)-G(x)|\leq \int_{-T}^T\Big|\frac{f(t)-g(t)}{t}\Big|dt+\frac{A}{T}\sup_{x\in{\mathbb R}}|G'(x)|
\end{equation}
where $f(t)=\int e^{itx}dF(x)$ and $g(t)=\int e^{itx}dG(x)=\hat{G'}(-t)$ (here $\hat g$ is the Fourier transform of  a function $g$).  Let $F(s)=\mu_{\omega}\{x:\in{\mathcal E}_{\omega}: S_n^{\omega}\leq\sqrt n s\}$. The idea behind the proof is to find a
function $H$ so that $\frac{|\hat F(t)-H(t)|}{|t|}$
will be sufficiently small on some interval $[-T,T]$, where
$T$ is of order $n^{-\frac{d}{2}}$. In this case we could take $G(s)=\int_{-\infty}^s\hat H(x)dx$ (so $\hat G'(-t)=H(t)$) and it will remain to verify that $G$ has the desired form.
On intervals of the form $I_n({\delta})=[-{\delta}\sqrt n,{\delta}\sqrt n]$ for small ${\delta}$'s this approximation is done using Assumption \ref{BaseAss}, while on
$J_n({\delta},a)=[-a\sqrt n, a\sqrt n]\setminus I_n({\delta})$ (for large $a$'s) we will use Assumption \ref{A3}. When $d>1$ we will use Assumption \ref{A4} in order to obtain an appropriate estimate on $[-T,T]\setminus J_n(a,{\delta})$.

Let $d\geq1$. The first step of the proof of Theorem \ref{EdgeThm} is as follows. The arguments preceding (\ref{II}) show that for any $z\in B(0,\frac \rho 2)$ we can write
\begin{equation}\label{II.1}
\mu_{\omega}(e^{zS_n^{\omega}})=e^{\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}(z)}W_{{\theta}^n{\omega}}(z)+{\lambda}_{{\omega},n}(z){\delta}_{{\omega},n}(z)
\end{equation}
where $W_{{\omega}}(z)=\nu_{{\omega}}(z)(h_{{\omega}}(z))$, $d_{{\omega},n}(0)=0$
and $|{\delta}_{{\omega},n}^{(k)}(z)|\leq C_k{\delta}^n$. Since ${\delta}_{{\omega},n}$ is analytic and ${\delta}_{{\omega},n}(0)=0$, it follows that $|{\delta}_{{\omega},n}(z)|\leq C'_0|z|{\delta}^n$ for $z$'s in the above domain, where $C_0'$ is some constant. Moreover, it follows from Assumption \ref{A3} (i) together with (\ref{ExpConvAss}) that $P$-a.s. we have
$\sup\{|{\lambda}_{{\omega},n}(it)|:\,t\in(-r,r),n\geq1\}\leq C_0''c_{\omega}$, where $c_{\omega}$ is some random variable. Therefore,
$P$-a.s. for any $|t|<\frac \rho 2$ and $n\geq1$,
\begin{equation}\label{First}
|{\lambda}_{{\omega},n}(it){\delta}_{{\omega},n}(it)|\leq C_{\omega}|t|{\delta}^n
\end{equation}
for some random variable $C_{\omega}$.
Next, since $\Pi_{{\omega}}(0)=\Pi_{\omega}'(0)=0$ and all the derivatives of $\Pi_{\omega}(z)$ are uniformly bounded in ${\omega}$, where ${\omega}$ ranges over a set of probability $1$,
we can write
\[
\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}(\frac {it}{\sqrt n})=-\frac{t^2\Pi_{{\omega},n,2}}{2}+n\psi_{{\omega},n}(\frac {t}{\sqrt n})
\]
where
\[
\Pi_{{\omega},n,2}=n^{-1}\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}^{(2)}(0)\,\,\text{ and }\,\,\psi_{{\omega},n}(z)=\frac 1n\sum_{j=0}^{n-1}(\Pi_{{\theta}^j{\omega}}(iz)-\frac{(iz)^2}2\Pi_{{\theta}^j{\omega}}''(0)).
\]
By (\ref{II.1}) and (\ref{First}) we have for $P$-a.a. ${\omega}$,
\begin{equation}\label{II.2}
\mu_{\omega}(e^{\frac{itS_n^{\omega}}{\sqrt n}})=\exp\big(-\frac{t^2\Pi_{{\omega},n,2}}{2}+n\psi_{{\omega},n}(\frac{t}{\sqrt n})\big)W_{{\theta}^n{\omega}}(\frac{it}{\sqrt n})
+|t|o(n^{-\frac{d}{2}}).
\end{equation}
Let $k\geq0$. Then, by Theorem \ref{EdgeThm} (i),  there exists a constant $B_k>0$ so that $P$-a.s. we have $\sup_{z\in B(0,\rho)}|\psi_{{\omega},n}(z)|\leq B_k$ for any $n\geq1$. Observe also that $\psi_{{\omega},n}(0)=\psi_{{\omega},n}'(0)=\psi_{{\omega},n}''(0)=0$. Therefore we can write $\psi_{{\omega},n}(z)=z^2\psi_{{\omega},n,d}(z)+|z|^{d+2}\tilde\psi_{{\omega},n,d}(z)$, where $\psi_{{\omega},n,d}$ is a polynomial of degree $d$ and $\tilde\psi_{{\omega},n,d}$ is an analytic function which vanishes at $z=0$ and is bounded in some deterministic neighborhood of the origin by some constant which does not depend on ${\omega}$ and $n$. Using the latter notations, we can write
\[
\exp\big(n\psi_{{\omega},n}(\frac{it}{\sqrt n})\big)=\exp\left(t^2\psi_{{\omega},n,r}(\frac{t}{\sqrt n})+
n^{-\frac d2}t^{d+2}\tilde\psi_{{\omega},n,d}(\frac{t}{\sqrt n})\right).
\]
Since $\mu_{\omega}(h_{\omega}(0))$, $\mu_{\omega}$ is a probability measure, $h_{\omega}(z)$ is analytic in $z$ and $\sup_{z\in B(0,\rho)}\|h_{\omega}(z)\|_1$ is a bounded random variable, there exists a constant $0<r_0<\frac \rho 2$ so that $P$-a.s. we can develop a branch of $\ln W_{{\omega}}(z)$ in $B(0,r_0)$ which is analytic, uniformly bounded in ${\omega}$ (when ${\omega}$ ranges over a set of probability $1$) and takes the value $\ln W_{\omega}(0)=0$ at $z=0$.
Next, we can write $W_{{\omega}}(it)=1+W_{{\omega},d}(t)+t^d\tilde W_{{\omega},d}(t)$, where $W_{{\omega},d}(t)$ is a polynomial of degree $d$ which vanishes at $t=0$, whose coefficients are bounded random variables,
and $\tilde W_{{\omega},d}(t)$ is a $C^\infty$ function which vanishes at $t=0$,
whose derivatives are uniformly bounded in ${\omega}$ (around $0$). Similarly to the derivation of equality  (3.7) in \cite{Liverani},
and using the Taylor expansions of the functions $\exp(\cdot)$ and $\ln(\cdot)$ we derive that $P$-a.s. for any
$t\in(- n^{\frac 12}r_0,n^{\frac 12}r_0)$,

\begin{eqnarray}\label{(3.7)}
\exp\big(\frac{t^2\Pi_{{\omega},n,2}}{2}+n\psi_{{\omega},n}(\frac{it}{\sqrt n})\big)W_{{\theta}^n{\omega}}(\frac{it}{\sqrt n})\\=
\exp\big(\frac{t^2\Pi_{{\omega},n,2}}{2}+n\psi_{{\omega},n}(\frac{it}{\sqrt n})+\ln e^{W_{{\theta}^n{\omega}}(\frac{it}{\sqrt n})}\big)\nonumber\\=
\exp\Big(t^2\psi_{{\omega},n,d}(\frac{it}{\sqrt n})+\frac1{n^{d/2}}t^{d+2}\tilde\psi_{{\omega},n,d}(\frac{it}{\sqrt n})-\nonumber\\\sum_{k=1}^{d}\frac{(-1)^{k+1}}{k}\big(W_{{\omega},d}(\frac{it}{\sqrt n})\big)^k-\frac1{n^{d/2}}t^{d}\overline {W}_{{\omega},d}(\frac{it}{\sqrt n})\Big)\nonumber\\=
1+\sum_{m=1}^d\frac1{m!}\Big(t^2\psi_{{\omega},n,d}(\frac{it}{\sqrt n})-\sum_{k=1}^d \frac{(-1)^{k+1}}{k}\big(W_{{\omega},d}(\frac{it}{\sqrt n})\big)^k\Big)^m \nonumber\\+\frac1{n^{d/2}}t^{d+2}\tilde\psi_{{\omega},n,d}(\frac{it}{\sqrt n})-\frac1{n^{d/2}}t^{d}\overline{W}_{{\omega},d}(\frac{it}{\sqrt n})+t^{d+1}{\mathcal O}(n^{-\frac{d+1}{2}})\nonumber\\=
\sum_{k=0}^d n^{-\frac k2}A_{{\omega},n,k}(t)+
t^d n^{-\frac d2}\varphi_{{\omega},n}(\frac{t}{\sqrt n})+t^{d+1}{\mathcal O}(n^{-\frac{d+1}{2}}).\nonumber
\end{eqnarray}
Here $A_0\equiv1$ and all the other $A_{{\omega},n,k}$'s are polynomials of degree $s_k$ which does not depend on ${\omega}$ and $n$, whose coefficients are algebraic combinations of the  derivatives of $W_{{\omega},n}$ at $0$ and the derivatives of $\psi_{{\omega},n}$ at $0$, and $\varphi_{{\omega},n}(z)=z^2\tilde\psi_{{\omega},n,d}(z)-\overline{W}_{{\omega},n,d}(z)$, where $\overline{W}_{{\omega},n,d}(z)$ is the reminder of $\ln W_{{\omega},n}(iz)$ when approximated by powers of $W_{{\omega},n,d}(iz)$.
Note that $\varphi_{{\omega},n}(z)$ is analytic, that it vanishes at $z=0$ and that $\sup_{a\in B(0,r_1)}|\varphi_{{\omega},n}(z)|\leq C_1$ ($P$-a.s.) for some positive constants $r_1$ and $C_1$ which do not depend on ${\omega}$.
Set
\[
Q_{{\omega},n}(t)=\sum_{k=1}^d n^{-\frac k2}A_{{\omega},n,k}(t).
\]
Since $\Pi_{{\omega},n,2}$ converges to ${\sigma}^2>0$ as $n\to\infty$,
it follows from (\ref{(3.7)}) that for any sufficiently small ${\delta}_0$, $P$-a.s. we have
\begin{eqnarray}\label{(3.10)}
\int_{-{\delta}_0\sqrt n}^{{\delta}_0\sqrt n}\Big|\frac{{\lambda}_{{\omega},n}(\frac{it}{\sqrt n})W_{{\theta}^n{\omega}}(\frac{it}{\sqrt n})-e^{-\frac{t^2\Pi_{{\omega},n,2}}{2}}(1+Q_{{\omega},n})(t)}t\Big|dt\\=
\int_{-{\delta}_0\sqrt n}^{{\delta}_0\sqrt n}e^{-\frac{t^2\Pi_{{\omega},n,2}}{2}}\Big|\frac{\exp\big(n\psi_{{\omega},n}(\frac{it}{\sqrt n})+\ln W_{{\theta}^n{\omega}}(\frac{it}{\sqrt n})\big)-1-Q_{{\omega},n}(t)}{t}\Big|dt=o(n^{-\frac d2}).\nonumber
\end{eqnarray}
Combining this with (\ref{II.2}), we derive that for any sufficiently small ${\delta}_0>0$ for $P$-a.a. ${\omega}$,
\begin{equation}\label{3.12}
\int_{-{\delta}_0\sqrt n}^{{\delta}_0\sqrt n}\Big|\frac{\mu_{\omega}(e^{\frac{it S_n^{\omega}}{\sqrt n}})-e^{-\frac{t^2\Pi_{{\omega},n,2}}{2}}(1+Q_{{\omega},n}(t))}{t}\Big| dt=o(n^{-\frac d2}).
\end{equation}
Consider the functions $g_{{\omega},n}(t)=e^{-\frac{t^2\Pi_{{\omega},n,2}}{2}}(1+Q_{{\omega},n}(t))$.
Using the fact that for any $c>0$, the function $\sqrt{2\pi c}(it)^ke^{-\frac{t^2c}{2}}$ is the  Fourier transform of the function $e^{-\frac{t^2}{2c}}$, and using the integration by parts formula with integrals of the form
$\int x^je^{-cx^2}dx$,
we derive that
\begin{equation}\label{Finding Poly1}
\sqrt{2\pi \Pi_{{\omega},n,2}}g_{{\omega},n}=\textbf{F}\Big(\int_{-\infty}^te^{-\frac{x^2}{2\Pi_{{\omega},n,2}}}dx+\sum_{k=1}^d n^{-\frac{k}{2}}R_{{\omega},n}(t)e^{-\frac{t^2}{2\Pi_{{\omega},n,2}}}\Big)
\end{equation}
where $\textbf{F}$ is the Fourier transform operator and
 all the $R_{{\omega},n,k}$'s are polynomials whose coefficients
are linear combinations of the coefficients of the $A_{{\omega},n,k}$'s, and the degree of $R_{{\omega},n,k}$ depends only on $k$. Define $G_{{\omega},n}$ by $G_{{\omega},n}(t)=\int_{-\infty}^t\textbf F\big(g_{{\omega},n}\big)(-\xi)d\xi$. Then by the Fourier inversion theorem, $\textbf{F}(G'_{{\omega},n})=g_{{\omega},n}$ and $\xi\to\sqrt{2\pi \Pi_{{\omega},n,2}}\textbf{F}(g_{{\omega},n})(-\xi)$ is the function inside the Fourier transform on the right hand side of (\ref{Finding Poly1}).
Using also the integration by parts formula with integrals of the form $\int t^je^{-ct^2}dt=\int t^{j-1}te^{-ct^2}dt$, we derive that
\begin{equation}\label{Finding Poly1.1}
\sqrt{2\pi \Pi_{{\omega},n,2}}G_{{\omega},n}(t)=\int_{-\infty}^te^{-\frac{x^2}{2\Pi_{{\omega},n,2}}}dx+\sum_{k=1}^d n^{-\frac{k}{2}}P_{{\omega},k,n}(t)e^{-\frac{t^2}{2\Pi_{{\omega},n,2}}}
\end{equation}
where all the $P_{{\omega},n,k}$'s are polynomials whose coefficients
are linear combinations of the coefficients of the $A_{{\omega},n,k}$'s, and the degree of $P_{{\omega},n,k}$ depends only on $k$.

Next, applying (\ref{Esseen}) with  $F=F_{{\omega},n}(s)=\mu_{\omega}\{x\in{\mathcal E}_{\omega}: S_n^{\omega}\leq\sqrt n s\}$ and $G=G_{{\omega},n}$,  taking into account (\ref{3.12}), we obtain that for any ${\varepsilon}>0$ and $B>\frac{A}{\varepsilon}$, where $A$ comes from (\ref{Esseen}),
\begin{eqnarray*}
\sup_{s\in{\mathbb R}}|F_{{\omega},n}(s)-G_{{\omega},s}(s)|\leq \int_{-Bn^{\frac d2}}^{Bn^{\frac d2}}\Big|\frac{\mu_{\omega}(e^{\frac{it S_n^{\omega}}{\sqrt n}})-e^{-\frac{t^2\Pi_{{\omega},n,2}}{2}}(1+Q_{{\omega},n}(t))}t\Big|+\frac{A}{Bn^{\frac d2}}\\\leq I_1+I_2+I_3+{\varepsilon} n^{-\frac d2}
\end{eqnarray*}
where
\begin{eqnarray*}
I_1=\int_{-{\delta}_0\sqrt n}^{{\delta}_0\sqrt n}\Big|\frac{\mu_{\omega}(e^{\frac{it S_n^{\omega}}{\sqrt n}})-e^{-\frac{t^2\Pi_{{\omega},n,2}}{2}}(1+Q_{{\omega},n}(t))}t\Big|dt=o(n^{-\frac d2}),\\
I_2=\int_{{\delta}_0\sqrt n<|t|<Bn^{\frac d2}}\Big|\frac{\mu_{\omega}(e^{\frac{itS_n^{\omega}}{\sqrt n}})}t\Big|dt\,\,\text{ and }\,\,\\
I_3=\int_{|t|>\sqrt n{\delta}_0}e^{-\frac{t^2\Pi_{{\omega},n,2}}2}\big|t^{-1}(1+Q_{{\omega},n}(t))\big|dt.
\end{eqnarray*}
We note that the polynomials $P_{{\omega},k,n}$ appearing in the definition of $G$ depend on $n$.
Therefore, in principle, one cannot apply (\ref{Esseen}) directly. One can do it when $G'_{{\omega},n}$ is uniformly bounded in $n$ and $\lim_{x\to\infty}\sup_{n\geq1}|G_{{\omega},n}(x)|=0$ (these conditions are clearly satisfied in our circumstances).
In (\ref{3.12}) we have shown that $I_1=o(n^{-\frac d2})$. Since $\Pi_{{\omega},n,2}$ converges to ${\sigma}^2>0$ and the coefficients of all $A_{{\omega},k}$'s in the definition of $Q_{{\omega},n}$ are bounded in $n$, it follows that $I_3={\mathcal O}(e^{-cn})$ for some $c>0$. In order to complete the proof of Theorem \ref{EdgeThm}, it is sufficient to show that $I_2=o(n^{-\frac d2})$. Recall that
\[
\mu_{\omega}(e^{it S_n^{\omega}})=\mu_{{\theta}^n{\omega}}({\mathcal A}_{it}^{{\omega},n}\textbf{1})
\]
for any $t\in{\mathbb R}$.
In the case when $d=1$, using Assumption \ref{A3} with $r=\frac d2$ and the compact set
$J=\{t\in{\mathbb R}: {\delta}_0\leq |t|\leq B\}$ we obtain that
\[
I_2\leq \frac{1}{\sqrt n{\delta}_0}\sup_{t\in J}\|{\mathcal A}_{it}^{{\omega},n}\|_1={\mathcal O}(n^{-\frac d2})
\]
for some $c>0$, which completes the proof of Theorem \ref{EdgeThm} (i). When $d>1$, in order to complete the proof of Theorem \ref{EdgeThm} (ii), it remains to
show that, $P$-a.s. for any sufficiently large $C>0$
\[
\int_{C\sqrt n<|t|<Bn^{\frac d2}}\Big|\frac{\mu_{\omega}(e^{\frac{itS_n^{\omega}}{\sqrt n}})}t\Big|dt=o(n^{-\frac d2}).
\]
This follows from the following lemma, together with the estimate
$|\mu_{\omega}(e^{itS_n^{\omega}})|\leq\|{\mathcal A}_{it}^{{\omega},n}\|_2$, where $\|\cdot\|_2$ is the norm specified in Assumption \ref{A4}.
\begin{lemma}
Suppose that Assumption \ref{A4} holds true with $r=\frac{d}2$ and that $({\Omega},{\mathcal F},P,{\theta})$ is ergodic. Then for any $q>0$ there exist positive constants $K,D$ and $L$ so that $P$-a.s. for any sufficiently large $n$, and $t$ such that $K\leq |t|\leq Dn^{r}$,
\[
\|{\mathcal A}_{it}^{{\omega},n}\|_2\leq LR_{\omega} n^{-q}
\]
where $R_{\omega}$ is the random variable specified in Assumption \ref{A4}.
\end{lemma}

\begin{proof}
Set $r=\frac d2>0$ let $K_{\omega},D_{\omega},C_{\omega}$, $r_2({\omega})$, $R_{\omega}$  and $N({\omega})$ be the random variables specified in Assumption \ref{A4} with this specific $r$. Let $K,D,r_2,C$ and $N$ be positive numbers so that the set
\[
{\Gamma}=\{{\omega}\in{\Omega}:\,K_{\omega}\leq K,\,D_{\omega}\geq D,\,r_2({\omega})\geq r_2,\,C_{\omega}\leq C\,\text{ and }\,N({\omega})\leq N\}
\]
has positive probability. Set $p_0=P({\Gamma})$. By Birkhoff's ergodic theorem and Kac's lemma applied with the transformation induced on $\Gamma$ by ${\theta}$, we obtain that  $P$-a.s. there exists a sequence $m_1<m_2<m_3<...$ of positive integers so that
\[
{\theta}^{m_i}{\omega}\in{\Gamma}\,\text{ for all }\,i's\,\text{ and }\,\lim_{k\to\infty}\frac{m_k}{k}=\frac1{p_0}.
\]
It follows that for any $0<a<b<1$ and a sufficiently large $n$ there exists an index $k=k_{{\omega},n,a,b}$ so that
$an<m_k<bn$. Indeed, let $1>{\varepsilon}>0$. Then for any sufficiently large $k$ we have
$\frac{1-{\varepsilon}}{p_0}k<m_k<\frac{1+{\varepsilon}}{p_0}k$. Suppose that ${\varepsilon}$ is sufficiently small and let
 ${\delta}>0$ be so that $ap_0<{\delta}(1-{\varepsilon})<{\delta}(1+{\varepsilon})<bp_0$. Let
$(k_n)_{n=1}^\infty$ be a strictly increasing sequence of positive integers so that $k_n/n$ converges to ${\delta}$ as $n\to\infty$. Then for any sufficiently large $n$ we have
\[
an<\frac{1-{\varepsilon}}{p_0}k_n<m_{k_n}<\frac{1+{\varepsilon}}{p_0}k_n<bn.
\]

Next, let $q>0$ and let $m$ be a positive integer such that $mr_2>r+q$. For any sufficiently large $n$ we can find $k_1(n),k_2(n),...,k_m(n)$ so that for each $i$,
\[
\big(\sum_{j=1}^{i}2^{-j}+2^{-i-2}\big)n>m_{k_i(n)}>\big(\sum_{j=1}^{i}2^{-j}\big)n.
\]
Set ${\omega}_i={\theta}^{m_i}{\omega}$, ${\Delta}_i(n)=k_{i+1}(n)-k_{i}(n)$, $i=1,2,...,m-2$, ${\Delta}_0(n)=k_1(n)-m_1$,
${\Delta}_{m-1}(n)=n-k_{m-1}(n)$ and $k_0(n)=1$. Then there exists a constant $c>0$ which depends only
on $m$ so that
\[
{\Delta}_i(n)\geq cn\,\,\text{ for any }\,i=0,1,...,m-1.
\]
Now, $P$-a.s for
any large enough $n$ we can write for all $t$'s,
\[
 {\mathcal A}_{it}^{{\omega},n}={\mathcal A}_{it}^{{\omega}_{k_{m-1}(n)},{\Delta}_{m-1}(n)}\circ\cdots\circ{\mathcal A}_{it}^{{\omega}_{k_1(n)},{\Delta}_1(n)}\circ{\mathcal A}_{it}^{{\omega}_{k_n(0)},{\Delta}_0(n)}\circ{\mathcal A}_{it}^{{\omega},m_1}.
\]
If $K\leq |t|\leq Dn^d$ then, since ${\omega}_i\in{\Gamma}$ for all $i$'s, for any sufficiently large $n$ and $0\leq i\leq m-1$ we have
\[
\|{\mathcal A}_{it}^{{\omega}_{k_i(n)},{\Delta}_i(n)}\|_2\leq C\big({\Delta}_i(n)\big)^{-r_2}\leq
Cc^{-r_2}n^{-r_2}.
\]
By Assumption \ref{A4} we have $\|{\mathcal A}_{it}^{{\omega},m_1}\|_2\leq R_{\omega}(1+|t|)$ which is less than
$R_{\omega}(1+ Dn^r)$ for $t$'s in the above range. We conclude from submultiplicity of norms of operators
that $P$-a.s. for any sufficiently large $n$
and $t$ so that $K\leq |t|\leq Dn^r$ we have
\[
\|{\mathcal A}_{it}^{{\omega},n}\|_2\leq C^mc^{-mr_2}(1+D)R_{\omega} n^{-mr_2+r}\leq
 C^mc^{-mr_2}(1+D)R_{\omega} n^{-q}:=LR_{\omega} n^{-q}
\]
and the proof of the lemma is complete.
\end{proof}

Now we will explain how to derive Theorem \ref{EdgeThm} (iii).
It follows from the definition of the $P_{{\omega},n,k}$'s that their coefficients are algebraic combinations of the derivatives at $z=0$
of the functions $\mu_{\omega}(h_{\omega}(z))$ and $n^{-1}\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}(z)$.
Since the derivatives of these functions $0$ are uniformly bounded in ${\omega}$ and $n$ (when ${\omega}$ ranges over a set of full probability), these coefficient are uniformly bounded in ${\omega}$ and $n$, and the coefficients of $P_{{\theta}^{-n}{\omega},n,k}$ converge $P$-a.s. and in $L^p$ for any $p\in[1,\infty)$. When
 (\ref{P rate}) holds true then all the coefficients of the polynomials $P_{{\theta}^{-n}{\omega},n,k}, k=1,2,...$ converge towards their limits with rate of order $n^{-\frac12}\ln n$.

\section{Examples}\setcounter{equation}{0}\label{secEXM}

\subsection{Random distance expanding maps}\label{SkewProdSec}
Let $({\Omega},{\mathcal F},P,{\theta})$, $({\mathcal X},\rho)$ and ${\mathcal E}_{\omega}$ be as described at the beginning of Section \ref{sec2},
and let
$
\{T_{\omega}: {\mathcal E}_{\omega}\to {\mathcal E}_{{\theta}{\omega}},\, {\omega}\in{\Omega}\}
$
be a collection of  maps between the metric spaces
${\mathcal E}_{\omega}$ and ${\mathcal E}_{{\theta}{\omega}}$, so that
the map $({\omega},x)\to T_{\omega} x$ is measurable with respect to the $\sigma$-algebra ${\mathcal P}$
which is the restriction of ${\mathcal F}\times{\mathcal B}$ on ${\mathcal E}$.
Consider the
skew product transformation $T:{\mathcal E}\to{\mathcal E}$ given by
\begin{equation}\label{Skew product}
T({\omega},x)=({\theta}{\omega},T_{\omega} x).
\end{equation}
For any ${\omega}\in{\Omega}$ and $n\in{\mathbb N}$ consider the $n$-th step iterates $T_{\omega}^n$ given by
\begin{equation}\label{T om n}
T_{\omega}^n=T_{{\theta}^{n-1}{\omega}}\circ\cdots\circ T_{{\theta}{\omega}}\circ T_{\omega}: {\mathcal E}_{\omega}\to{\mathcal E}_{{\theta}^n{\omega}}.
\end{equation}
Our additional requirements concerning the family of maps $\{T_{\omega}:{\omega}\in{\Omega}\}$
are collected in the following assumptions.
\begin{assumption}[Topological exactness]\label{TopExAssRand}
There exist a constant $\xi>0$ and a random variable $n_{\omega}\in{\mathbb N}$ such that
$P$-a.s.,
\begin{equation}\label{TopExRand}
T_{\omega}^{n_{\omega}}(B_{\omega}(x,\xi))=
{\mathcal E}_{{\theta}^{n_{\omega}}{\omega}}\,\text{ for any } x\in{\mathcal E}_{\omega}
\end{equation}
where for any ${\omega}\in{\Omega}$, $x\in{\mathcal E}_{\omega}$ and $r>0$,
$\,B_{\omega}(x,r)$ denotes a ball in ${\mathcal E}_{\omega}$ around $x$ with radius $r$.
\end{assumption}

\begin{assumption}[The pairing property]\label{Ass pair prop}
There exist random variables ${\gamma}_{\omega}>1$ and $D_{\omega}\in{\mathbb N}$ such that $P$-a.s.
for any $x,x'\in{\mathcal E}_{{\theta}{\omega}}$ with $\rho(x,x')<\xi$
we can write
\begin{equation}\label{Pair1.0}
T_{\omega}^{-1}\{x\}=\{y_1,...,y_k\}\,\,\text{ and }\,\,T_{\omega}^{-1}\{x'\}=\{y_1',...,y_k'\}
\end{equation}
where $\xi$ is specified in Assumption \ref{TopExAssRand},
\[
k=k_{{\omega},x}=|T_{\omega}^{-1}\{x\}|\leq D_{\omega}
\]
and
\begin{equation}\label{Pair2.0}
\rho(y_i,y_i')\leq ({\gamma}_{\omega})^{-1}\rho(x,x')
\end{equation}
for any $1\leq i\leq k$.
\end{assumption}
These assumptions hold true for the random distance expanding maps considered in \cite{MSU} since in the setup there the inverse branches are contracting, see Section 2.5 in \cite{MSU}.
We refer the readers to Section 2.1 in \cite{MSU} for a  description of several examples of such maps $T_{\omega}$.

According to Lemma 4.11 in \cite{MSU} (applied with $r=\xi$),
there exist an integer valued random variable $L_{\omega}\geq1$ and ${\mathcal F}$-measurable
functions ${\omega}\to x_{{\omega},i}\in{\mathcal X},\,i=1,2,3,...$ so that
$x_{{\omega},i}\in{\mathcal E}_{\omega}$ for each $i$ and
\begin{eqnarray}\label{cover}
\bigcup_{k=1}^{L_{\omega}} B_{\omega}(x_{{\omega},k},\xi)={\mathcal E}_{\omega},
\,\,\,P\text{-a.s.}
\end{eqnarray}
Note that $L_{\omega}$ is constant in ${\omega}$ when
${\mathcal E}_{\omega}$ does not depend on ${\omega}$ (i.e. when ${\mathcal E}={\Omega}\times{\mathcal Y}$ for an appropriate
${\mathcal Y}\subset{\mathcal X}$).

Next, for any $g:{\mathcal E}\to{\mathbb C}$ and ${\omega}\in{\Omega}$ consider
the function $g_{\omega}:{\mathcal E}_{\omega}\to{\mathbb C}$ given by $g_{\omega}(x)=g({\omega},x)$.
Let $\phi, u:{\mathcal E}\to{\mathbb R}$ be measurable functions so that for $P$-a.a. ${\omega}$ the functions
$\phi_{\omega}$ and $u_{\omega}$ are H\"older continuous
with exponent ${\alpha}\in(0,1]$ which is
independent of ${\omega}$.
Let $z\in{\mathbb C}$ and consider the transfer operators ${\mathcal L}_z^{\omega},\,{\omega}\in{\Omega}$ which
map functions on ${\mathcal E}_{\omega}$ to functions on ${\mathcal E}_{{\theta}{\omega}}$ by
the formula\index{transfer operator!random!complex}
\begin{eqnarray}\label{TraOp}
{\mathcal L}^{\omega}_z g(x)=\sum_{y\in T_{\omega}^{-1}\{x\}}e^{\phi_{\omega}(y)+zu_{\omega}(y)}g(y).
\end{eqnarray}
Note that under Assumption \ref{Ass pair prop} the operators
${\mathcal L}^{\omega}_z,\,z\in{\mathbb C}$ are well defined and since $\phi_{\omega}$ and $u_{\omega}$ are H\"older continuous
they map a continuous function on ${\mathcal E}_{\omega}$
to a continuous function on ${\mathcal E}_{{\theta}{\omega}}$.
For any $n\in{\mathbb N}$ and $z\in{\mathbb C}$ consider the $n$-th
step iterates ${\mathcal L}_z^{{\omega},n}$ of the transfer operator given by
\begin{equation}\label{iter}
{\mathcal L}_z^{{\omega},n}={\mathcal L}_z^{{\theta}^{n-1}{\omega}}\circ\cdots\circ{\mathcal L}_z^{{\theta}{\omega}}\circ{\mathcal L}_z^{\omega}.
\end{equation}
Then
\begin{equation}\label{iter2}
{\mathcal L}_z^{{\omega},n}g(x)=\sum_{y\in (T_{\omega}^n)^{-1}\{x\}}e^{S_n^{\omega} \phi(y)+zS_n^{\omega} u(y)}g(y)
\end{equation}
where
$
S_n^{\omega} \psi=\sum_{i=0}^{n-1}\psi_{{\theta}^i{\omega}}\circ T_{\omega}^i
$
for any function $\psi:{\mathcal E}\to{\mathbb C}$, ${\omega}\in{\Omega}$ and $n\geq1$.
We also consider the (global) transfer operator ${\mathcal L}_z$
acting on functions  $g:{\mathcal E}\to{\mathbb C}$ by the formula
\begin{equation}\label{TraOp1}
{\mathcal L}_z g(s)=
\sum_{s'\in T^{-1}\{s\}=({\theta}^{-1}{\omega},y)}e^{\phi_{{\theta}^{-1}{\omega}}(y)+zu_{{\theta}^{-1}{\omega}}(y)}g(s')=
{\mathcal L}^{{\theta}^{-1}{\omega}}_z g_{{\theta}^{-1}{\omega}}(x),\,\,\,s=({\omega},x),
\end{equation}
namely ${\mathcal L}_z$  is generated by the skew product map $T$ and the function
$\phi+zu$.
Next, let $g:{\mathcal E}\to{\mathbb C}$ be a measurable function. Let
${\alpha}$ as described before (\ref{TraOp}) and set
\begin{eqnarray*}
v_{{\alpha},\xi}(g_{\omega})=\inf\{R: |g_{\omega}(x)-g_{\omega}(x')|\leq R\rho^{\alpha}(x,x')\,\text{ if }\,
\rho(x,x')<\xi\}\\
\text{and }\,\,\,\|g_{\omega}\|_1=\|g_{\omega}\|_{{\alpha},\xi}=\|g_{\omega}\|_\infty+v_{{\alpha},\xi}(g_{\omega})\hskip1cm
\end{eqnarray*}
where $\|\cdot\|_\infty$ is the supremum norm and $\rho^{\alpha}(x,x')=\big(\rho(x,x')\big)^{\alpha}$.
These norms are ${\mathcal F}$-measurable as a consequence
of Lemma 5.1.3 in \cite{book}. We denote by ${\mathcal H}_{{\alpha},\xi}^{\omega}=H_1^{\omega}$ the space of all functions
$f:{\mathcal E}_{\omega}\to{\mathbb C}$ so that $\|f\|_{{\alpha},\xi}<\infty$.

Our additional requirements from $T_{\omega},\phi_{\omega}$ and $u_{\omega}$ are specified in the following
\begin{assumption}\label{bound ass}
(i) The random variables $n_{\omega},D_{\omega},L_{\omega}$
are bounded and ${\gamma}_{\omega}-1$ is bounded from below by some positive
constant.

(ii) For any measurable function $g:{\mathcal E}\to{\mathbb C}$  so that $\|g_{\omega}\|_{{\alpha},\xi}\leq 1$ the map  $({\omega},x)\to{\mathcal L}_0 g({\omega},x)$ is measurable.

(iii) For $P$-a.a. ${\omega}$ we have $\phi_{\omega},u_{\omega}\in {\mathcal H}_{\omega}^{{\alpha},\xi}$  and
the random variables $\|\phi_{\omega}\|_{{\alpha},\xi}$ and $\|u_{\omega}\|_{{\alpha},\xi}$ are bounded.
\end{assumption}
Assumption \ref{bound ass} (i) holds true, for instance when ${\mathcal E}_{\omega}={\mathcal X}=[0,1)$ (and then $L_{\omega}$ does not depend on ${\omega}$), and each  $T_{\omega}$ is a piecewise distance expanding map on the unit interval with a bounded in ${\omega}$ number of monotonicity intervals (on which $T_{\omega}$ is onto $[0,1)$). Several multidimensional examples such as the case when $T_{\omega} x=(m_1({\omega}) x_1,...,m_d({\omega})x_d)\,\text{mod } 1$ can be also be given (here $m_i({\omega})\geq 2$ are integer valued and are bounded in ${\omega}$).
 We refer again the readers to Section 2.1 in \cite{MSU} for more examples of this nature. Different type of examples are random topologically mixing subshifts of finite type with bounded number of symbols, see \cite{Kifer-1998} for a description of such systems.
Note that under Assumption  \ref{bound ass} we have that ${\mathcal L}_z^{{\omega},n}({\mathcal H}_{\omega}^{{\alpha},\xi})\subset{\mathcal H}_{{\theta}^n{\omega}}^{{\alpha},\xi}$
and the corresponding operator norm satisfies
$\|{\mathcal L}_{it}^{{\omega},n}\|_{{\alpha},\xi}\leq B(1+|t|)$ where $B$ is some constant (see Lemma  5.6.1 in \cite{book}).

When the random variable $\|\phi_{\omega}\|_{{\alpha},\xi}$  is bounded then Assumption \ref{bound ass} (ii)
holds true if for any set $A\subset{\mathcal X}$ with a sufficiently small diameter the set
\[
\textbf{A}=\{({\omega},x): x\in T_{\omega}(A\cap{\mathcal E}_{\omega})\}
\]
is measurable (note that when ${\mathcal E}_{\omega}={\mathcal X}$ then $\textbf{A}=T({\Omega}\times A)$ and that $T_{\omega}$ is injective on ${\mathcal E}_{\omega}\cap A$ if the diameter of $A$ is smaller than $\xi$).
Indeed,  in these circumstances we can approximate (uniformly in ${\omega}$) the function $g_{\omega}\cdot e^{\phi_{\omega}}$, by sums of the form $\sum_{i=1}^m a_i({\omega}){\mathbb I}_{A_i\cap {\mathcal E}_{\omega}}$, where $\{A_i\}$ is a partition of ${\mathcal X}$ into sets with small diameter, ${\mathbb I}_{A_i\cap {\mathcal E}_{\omega}}$ is the indicator set of $A_i\cap{\mathcal E}_{\omega}$ and $a_i({\omega})$ are some random variables whose absolute value is bounded by $1$. Then ${\mathcal L}_0$ can be approximated by
\[
\sum_{i=1}^m a_i({\omega}){\mathbb I}(x\in T_{\omega}(A_i\cap{\mathcal E}_{\omega}))
\]
which are measureable functions of $({\omega},x)$. Observe that it is enough to assume that the $\textbf{A}$'s are measurbale for sets $A$ which are generated by finite collections of open balls using only the basic Boolean operations. For such sets  $\textbf{A}$ is measurable if ${\Omega}$ is a topological space, ${\mathcal F}$ is the Borel ${\sigma}$-algebra, ${\mathcal E}_{\omega}={\mathcal X}$ and $T$ is continuous. Another example is the case when $T_{\omega}$ takes only a countable number of ``values" $\{T_1,T_2,...\}$. We refer the readers to Chapter 4 in \cite{MSU} for more aspects of such measurability issues.

\subsubsection{Asymptotic moments: verification of Assumption \ref{BaseAss}}
We recall first the following random complex RPF theorem, which was stated in \cite{book} as Corollary 5.4.2.

\begin{theorem}\label{RPF rand T.O. general}
Suppose that Assumptions \ref{TopExAssRand}, \ref{Ass pair prop} and \ref{bound ass} hold true. Then there exists $\rho>0$ so that $P$-a.s.
for any $z\in U:=\{\zeta\in{\mathbb C}: |\zeta|<\rho\}$ there is a unique triplet consisting of a
 nonzero complex number $\boldsymbol{{\lambda}}_{\omega}(z)$
a function $\textbf{h}_{\omega}(z)\in{\mathcal H}^{{\alpha},\xi}_{\omega}$ and a continuous linear functional $\boldsymbol{\nu}_{\omega}(z)$
on ${\mathcal H}^{{\alpha},\xi}_{\omega}$ so that $\boldsymbol{\nu}_{\omega}(z)\textbf{1}=1$,
\begin{equation}\label{RPF1}
{\mathcal L}_z^{\omega} \textbf{h}_{\omega}(z)=\boldsymbol{{\lambda}}_{\omega}(z)\textbf{h}_{{\theta}{\omega}}(z),\,\,
({\mathcal L}_z^{\omega})^*\boldsymbol{\nu}_{{\theta}{\omega}}(z)=\boldsymbol{{\lambda}}_{\omega}(z)\boldsymbol{\nu}_{\omega}(z)\,\text{ and }\,
\boldsymbol{\nu}_{\omega}(z)\big(\textbf{h}_{\omega}(z)\big)=1.
\end{equation}
For any $z\in U$ the maps ${\omega}\to\boldsymbol{{\lambda}}(z)$ and $({\omega},x)\to \textbf{h}_{\omega}(z)(x),\,({\omega},x)\in{\mathcal E}$
are measurable and the family $\boldsymbol{\nu}_{\omega}(z)$ is measurable in ${\omega}$.
When $z=t\in{\mathbb R}$ then $\boldsymbol{{\lambda}}_{\omega}(t)>0$, the function $\textbf{h}_{\omega}(t)$ is strictly
positive, $\boldsymbol{\nu}_{\omega}(t)$ is a probability measure and the equality
$\boldsymbol{\nu}_{{\theta}{\omega}}(t)\big({\mathcal L}_t^{\omega} g)=\boldsymbol{{\lambda}}_{\omega}(t)\boldsymbol\nu_{{\omega}}(t)(g)$ holds true for any
bounded Borel function $g:{\mathcal E}_{\omega}\to{\mathbb C}$ and when the maps $T_{\omega}$ satisfy the assumptions from
\cite{MSU} then the triplet $(\boldsymbol{{\lambda}}_{\omega}(t),\textbf{h}_{\omega}(t),\boldsymbol{\nu}_{\omega}(t))$ coincides with the one
constructed there with the function $\phi+tu$.

Moreover, this triplet is analytic and uniformly bounded around $0$.
Namely, the maps
\[
\boldsymbol{\lambda}_{\omega}(\cdot):U\to{\mathbb C}, \textbf{h}_{\omega}(\cdot):U\to {\mathcal H}_{\omega}^{{\alpha},\xi}\,\text{ and }
\boldsymbol\nu_{\omega}(\cdot):U\to \big({\mathcal H}_{\omega}^{{\alpha},\xi}\big)^*
\]
are analytic, where $({\mathcal H}_{\omega}^{{\alpha},\xi})^*$ is the dual space of ${\mathcal H}_{\omega}^{{\alpha},\xi}$,
and there is a constant $C_0>0$ so that
\begin{equation}\label{UnifBound}
\max\Big(\sup_{z\in U}|\boldsymbol{\lambda}_{\omega}(z)|,\,
\sup_{z\in U}\|\textbf{h}_{\omega}(z)\|_{{\alpha},\xi},\, \sup_{z\in U}
\|\boldsymbol\nu_{\omega}(z)\|_{{\alpha},\xi}\Big)\leq C_0,\,\,P\text{-a.s.}
\end{equation}
where  $\|\nu\|_{{\alpha},\xi}$ is the
operator norm of a linear functional $\nu:{\mathcal H}_{\omega}^{{\alpha},\xi}\to{\mathbb C}$.

Furthermore, there exists constants $C>0$ and ${\delta}\in(0,1)$ so that $P$-a.s. for any $n\geq1$ and $g\in{\mathcal H}^{{\alpha},\xi}_{\omega}$ we have
\begin{equation}\label{ExpConv}
\left\|(\boldsymbol{\lambda}_{{\omega},n}(z))^{-1}{\mathcal L}_z^{{\omega},n}g-\boldsymbol\nu_{\omega}(z)(g)\textbf{h}_{{\theta}^n{\omega}}(z)\right\|_{{\alpha},\xi}\leq C\|g\|_{{\alpha},\xi}{\delta}^n.
\end{equation}
\end{theorem}
Note that in \cite{book} we have assumed that ${\mathcal L}_0g$ is measurable  for any measurable function $g$, which was only needed in order to insure that the RPF triplets are measurable in ${\omega}$, but taking a careful look at the arguments in Chapter 4 of \cite{book} we see that only Assumption \ref{bound ass} (ii) is needed.

We remark that the random  (Gibbs) measure $\mu_{\omega}=\textbf{h}_{\omega}(0)\boldsymbol{\nu}_{\omega}(0)$ satisfies
$(T_{\omega})_*\mu_{\omega}=\mu_{{\theta}{\omega}}$, or, equivalently,
$\mu=\int_{{\Omega}}\mu_{\omega} dP({\omega})$ is $T$-invariant, where $T({\omega},x)=({\theta}{\omega},T_{\omega} x)$ is the skew product map.
Consider now the transfer operator
operator ${\mathcal A}_z^{\omega}$ generated by the map $T_{\omega}$ and the function $\tilde\phi_{\omega}+zu_{\omega}$, where $\tilde\phi_{\omega}=\phi_{\omega}+\ln \textbf{h}_{\omega}(0)-\ln \textbf{h}_{{\theta}{\omega}}(0)\circ T_{\omega}-\ln\boldsymbol{{\lambda}}_{{\omega}}(0)$. Then ${\mathcal A}_z^{\omega}\textbf{1}=\textbf{1}$,
\[
\mu_{\omega}(e^{it S_n^{\omega} u})=\mu_{{\theta}^n{\omega}}({\mathcal A}_{it}^{{\omega},n}\textbf{1}),
\]
 ${\mathcal A}_z^{{\omega},n}({\mathcal H}_{\omega}^{{\alpha},\xi})\subset{\mathcal H}_{{\theta}^n{\omega}}^{{\alpha},\xi}$
and ${\mathcal A}_z^{\omega}$ also satisfies the above random complex RPF theorem. More precisely,
Assumption \ref{A3} is satisfied with the norm $\|\cdot\|_1=\|\cdot\|_{{\alpha},\xi}$ and triplets
$({\lambda}_{\omega}(z),h_{\omega}(z),\nu_{\omega}(z))$ which also satisfy ${\lambda}_{\omega}(0)=1$, $h_{\omega}(0)\equiv\textbf{1}$ and $\nu_{\omega}(0)=\mu_{\omega}$.   We refer the readers to the begging of Section 4 in \cite{SeqRPF} for the precise details.
Since $\|{\mathcal A}_{it}^{{\omega},n}\textbf{1}\|_{{\alpha},\xi} \leq C\|{\mathcal L}_{it}^{{\omega},n}\textbf{1}\|_{{\alpha},\xi}$ for some constant $C>0$ (see \cite{book}, Ch. 5), it is sufficient to show that ${\mathcal L}_{z}^{\omega}$ satisfy Assumption \ref{A3} with the above norm.
Henceforth, we will always assume that $\mu_{\omega}(u_{\omega})=0$ for $P$-a.a. ${\omega}$, which, due to $T$-invariance of $\mu$, implies that $\mu_{\omega}(S_n^{\omega} u)=0$ for each $n\geq1$. This is not really a restriction since we can always replace $u_{\omega}$ by $u_{\omega}-\int u_{\omega}(x)d\mu(x)$.

\subsubsection{Edgeworth expansions of order one: verification of Assumption \ref{A3}}

The first part of Assumption \ref{A3} is satisfied due to the random Lasota-Yorke type inequality stated in
Lemma 5.6.1 in \cite{book}.
Next, under the following assumption we showed in Chapter 7 of \cite{book} that the rest of Assumption \ref{A3} holds true for the operators ${\mathcal L}_z^{\omega}$.

\begin{assumption}\label{PerPointAss
(i)  ${\Omega}$ is a topological space, ${\mathcal F}$ is the corresponding
Borel ${\sigma}$-algebra, $({\Omega},{\mathcal F},P,{\theta})$ is ergodic and ${\theta}$ has a periodic point, namely there exist ${\omega}_0\in{\Omega}$ and $n_0\in{\mathbb N}$ so that ${\theta}^{n_0}{\omega}_0={\omega}_0$. Moreover, $P(U)>0$ for any open set $U$ which intersects the finite orbit $\{{\theta}^i{\omega}_0:\,0\leq i<n_0\}$ and
the spaces ${\mathcal E}_{\omega}$ are locally independent of ${\omega}$ around the points ${\theta}^{j}{\omega}_0,\, 0\leq j<n_0$.

(ii) The map ${\theta}$ is continuous at the points ${\omega}_0,{\theta}{\omega}_0,...,{\theta}^{n_0-1}{\omega}_0$, where ${\omega}_0$ and $n_0$ come from (i), and for any compact set $J\subset{\mathbb R}\setminus\{0\}$ the family of maps
${\omega}\to{\mathcal L}_{it}^{\omega},\,t\in J$ is equicontinuous continuous (with respect to the operator norm $\|\cdot\|_{{\alpha},\xi}$) at the points ${\theta}^{j}{\omega}_0,\,0\leq j<n_0$.

(iii) The function $S_{n_0}^{{\omega}_0}u$ is non-arithmetic
(aperiodic) with respect to the map $T^{n_0}_{{\omega}_0}$ in the classical sense of \cite{HH}, namely
if for any $t\in{\mathbb R}\setminus\{0\}$ there exist no
nonzero $g\in{\mathcal H}_{{\omega}_0}^{{\alpha},\xi}$ and ${\lambda}\in{\mathbb C}$, $|{\lambda}|=1$
such that
\begin{equation}\label{lattice condition 1.}
e^{itS_{n_0}^{{\omega}_0}u}g={\lambda} g\circ T_{{\omega}_0}^{n_0}.
\end{equation}
\end{assumption}
Note that under Assumption \ref{PerPointAss} we showed in Chapter 7 of \cite{book} that, in fact,
 for any compact set $J\subset{\mathbb R}\setminus\{0\}$, $P$-.a.s for any $n\geq1$ we have
\begin{equation}\label{ExpR}
\sup_{t\in J}\|{\mathcal A}^{{\omega},n}_{it}\|_1\leq d_{\omega} 2^{-cn}
\end{equation}
for some random variable $d_{\omega}>0$ and a constant $c>0$.
In Chapter 7 of \cite{book} we stated that (\ref{ExpR}) holds true under the additional assumption that $P(U)>0$ for any open set $U$, but in the proof we relied on that only for $U$'s which intersect the orbit $\{{\omega}_0,{\theta}{\omega}_0,...,{\theta}^{n_0-1}{\omega}_0\}$.
The condition about  continuity of ${\omega}\to{\mathcal L}_{it}^{\omega}$ holds true, for instance,
when the maps ${\omega}\to \phi_{\omega}, u_{\omega}\in{\mathcal H}^{{\alpha},\xi}_{\omega}$ are continuous at the
points ${\theta}^{j}{\omega}_0,\, 0\leq j<n_0$ and the  branches of $T_{\omega}$ are either locally constant as functions of ${\omega}$ there, or they are continuous
at these points with respect to an appropriate topology. For instance, we can just assume that $T_{\omega}$ does not depend on ${\omega}$ around the points ${\theta}^i{\omega}_0$. This is the case when $({\Omega},{\mathcal F},{\sigma})$  is a shift system built over a countable alphabet $\mathscr A$ so that $T_{\omega}=T_{{\omega}_0}$ depends only on the $0$-th coordinate of ${\omega}=({\omega}_i)$. We can consider any sub-shift of ${\Omega}$ instead, assuming it contains a point of the form $(...a,a,a,...)$  where $a\in\mathscr A^{n_0}$ for some $n_0$. Of course, we we also have to assume that $P(U)>0$ for the appropriate $U$'s.
When considering interval maps $T_{\omega}$ with deterministic monotonicity intervals $I_1,...,I_D$ (or similar multidimensional examples), then the condition about the continuity holds true if the restriction of the inverse of each map $T_{\omega}|_{I_i}$ is continuous in ${\omega}$ at ${\omega}={\theta}^j{\omega}_0$, $i=0,...,n_0-1$.

Note that when ${\mathcal X}$ is a smooth connected Riemannian manifold then we can also consider the norms $\|g\|_{C^1}=\sup|g|+\sup|Dg|$ instead of $\|\cdot\|_{{\alpha},\xi}$. In this case when each $T_{\omega}$ is also smooth and the map ${\omega}\to T_{\omega}$ is continuous in the $C^1$ topology then the transfer operators ${\omega}\to{\mathcal L}_{\omega}$ are continuous in the operator-norm topology (corresponding to the norm $\|\cdot\|_{C^1}$), see Proposition 5.3 in \cite{castro}. Therefore, when considering smooth functions $\phi_{\omega}$ and $u_{\omega}$ which are continuous in the $C^1$-topology as functions of ${\omega}$ we get that the conditions in item (ii) above are satisfied with $\|\cdot\|_{C^1}$ instead of $\|\cdot\|_{{\alpha},\xi}$, and all our results hold true.

\subsubsection{Edgeworth expansions of high orders: verification of Assumption \ref{A4}}
Now we will explain in which circumstances Assumption \ref{A4} is satisfied with any  $r$. We first recall the following results from  \cite{[3]}. Let $\mathscr{L}_b,\,b\in{\mathbb R}$ be the deteministic transfer operators  which are generated by two  piecewise smooth expanding functions $f:{\mathbb T}\to{\mathbb T}$ and $\tau:{\mathbb T}\to{\mathbb T}$ and is given by the formula
\[
\mathscr{L}_bh(x)=\sum_{y\in f^{-1}\{x\}}|f'(y)|^{-1}h(y)e^{ib\tau(y)}.
\]
Let $\|\cdot\|_2=\|\cdot\|_{{\textbf{BV}}}=\|\cdot\|_{L^1({\mathbb T})}+\text{Var}(\cdot)$ be the norm in the space of functions with bounded variation.
When $\tau$ is not cohomologous to a piecewise constant function (with respect to the base map $f$, see \cite{[3]}),
then by Proposition 1 in \cite{[3]}, there exist positive constants $b_0,\rho$ and ${\gamma}_2$ such that
\begin{equation}\label{Prop1}
\|\mathscr{L}_b^{n(b)}\|_{(b)}\leq e^{-n(b){\gamma}_2}\,\,\text{ for any }\,\,|b|\geq b_0,\,n(b):=[\rho\ln|b|]
\end{equation}
where $\|h\|_{(b)}=(1+|b|)^{-1}\|h\|_{\textbf{BV}}+\int_{{\mathbb T}}|h(x)|dx$. In \cite{Liverani} the authors showed that this implies
that Assumption \ref{A4} holds true in this deterministic case.

Next, we assume that ${\mathcal E}_{\omega}={\mathbb T}=[0,1)$ for any ${\omega}$ and that, in addition to the requirements at the beginning of Section \ref{SkewProdSec}, the functions $T_{\omega}$ and $u_{\omega}$ are piecewise smooth expanding functions and that $u_{\omega}$ takes values in ${\mathbb T}$. Then Assumptions \ref{BaseAss} holds true.
In order to keep our notations in line with the ones in \cite{[3]}, we set  $f_{\omega}=T_{\omega}$ and $\tau_{\omega}=u_{\omega}$. We also set
$\phi_{\omega}=-\ln |f_{\omega}'|$.
We assume here that there exists no random family $\eta_{\omega}$ of functions $\eta_{\omega}:{\mathbb T}\to{\mathbb T}$ so that $\tau_{\omega}-\eta_{{\theta}{\omega}}\circ f_{\omega}+\eta_{\omega}$ is piecewise constant ($P$-a.s.).
Then, when $({\Omega},{\mathcal F},P,{\theta})$ is ergodic, an appropriate quenched version of Proposition  1 in \cite{[3]}
 holds true when we replace $f^n$ and $\sum_{j=0}^{n-1}\tau\circ f^j$ with $f_{\omega}^n:=f_{{\theta}^{n-1}{\omega}}\circ f_{{\theta}^{n-2}{\omega}}\circ\cdots\circ f_{\omega}$ and $\sum_{j=0}^{n-1}\tau_{{\theta}^j{\omega}}\circ f_{\omega}^j$, respectively, when all the conditions in \cite{[3]} hold true uniformly in ${\omega}$
 in the following sense:
 there exist random variables $\rho_{\omega},{\gamma}_2({\omega}),b_{\omega}>0$ so that $P$-a.s., for any $b\in{\mathbb R}$,
\begin{equation}\label{RandomProp1}
\|{\mathcal L}^{{\omega},n_{\omega}(b)}_{ib}\|_{(b)}\leq e^{-n_{\omega}(b){\gamma}_2({\omega})}\,\,\text{ for any }\,\,|b|\geq b_{\omega},\,n_{\omega}(b):=[\rho_{\omega}\ln|b|].
\end{equation}

Before explaining how the proof of (\ref{RandomProp1}) differs from the proof of (\ref{Prop1}), we will rely on (\ref{RandomProp1}) in order to show that Assumption \ref{A4} holds true. Indeed, let
$b_0,\rho_1,\rho_2,{\gamma}_2>0$  be so that
\[
P\{{\omega}:\,\rho_1\leq \rho_{\omega}\leq\rho_2,\,{\gamma}_2({\omega})\geq{\gamma}_2\,\text{ and }
b_{\omega}\leq b_0\}:=a>0.
\]
Then by ergodicity of
$({\Omega},{\mathcal F},P,{\theta})$, for $P$-almost any ${\omega}$ there exists a sequence of positive integers $n_1<n_2<n_3<...$ so that for each $i$,
\[
\rho_1\leq \rho_{{\theta}^{n_i}{\omega}}\leq\rho_2,\,\, b_{{\theta}^{n_i}{\omega}}\leq b_0\,\text{ and }\, {\gamma}_2({\theta}^{n_i}{\omega})\geq {\gamma}_2
\]
and $\lim_{k\to\infty}k^{-1}n_k=a^{-1}$. It follows that for any  $c>0$ there exists a constant $c_1>0$ so that
for any sufficiently large $n$, there are indexes $i_1<i_2<...<i_s$, $s=s_n=[c_1\frac{n}{\ln n}]$, $i_j=i_j(n)$ with the property that $n_{i_j}+c\ln n<n_{i_{j+1}}<n-c\ln n$ for any $j=1,2,...,s-1$. In fact, since $n^{-1}\{\max k: n_k\leq n\}$ converges as $n\to\infty$ to $a$, we can always take $i_1=1$, $i_s<\frac n2-c\ln n$ and $c_1=\frac{a}{4(c+1)}$.

Next, it follows exactly as in Lemma 2 in \cite{[3]} that $\|{\mathcal L}_{ib}^{{\omega},n}\|_{(b)}\leq C$ for any $b$ and $n$, where $C\geq1$ is some constant (in fact, it is possible to take $C=1$, but it does not matter for our purpose). Let $d>0$
and let $b_1>\max(b_0,1)$ be so that $b_1^{\frac12\rho_1{\gamma}_2}>2C$.
 Then for any real $b$ so that $b_1\leq |b|\leq n^{d}$ we have
$\rho_2\ln|b|\leq\rho_2d\ln:=c\ln n$. Let $c_1$ and $i_1<i_2<...<i_s$ (for a sufficiently large $n$) be as describe in the previous paragraph  with $c=\rho_2d$.
For $n$ large enough, set
\begin{eqnarray*}
{\omega}_j={\omega}_j(n)={\theta}^{n_{i_j}}{\omega},\, m_j(b)=[\rho_{{\omega}_j}\ln|b|],\,
d_j(b)=d_j(b,n)=n_{i_{j+1}}-n_{i_{j}}-m_j(b),\\ \text{ where } n_{s+1}:=n, \text{ and }
{\omega}'_{j}(b)={\omega}'_j(b,n)={\theta}^{[\rho_{{\omega}_{j}}\ln|b|]}{\omega}_j.
\end{eqnarray*}
By writing
\[
{\mathcal L}_{ib}^{{\omega},n}={\mathcal L}_{it}^{{\omega}'_{s}(b),d_s(b)}\circ{\mathcal L}_{ib}^{{\omega}_s,m_s(b)}\circ\cdots\circ{\mathcal L}_{ib}^{{\omega}_2,m_2(b)}\circ{\mathcal L}^{{\omega}'_{1}(b),d_j(b)}_{ib}\circ{\mathcal L}_{ib}^{{\omega}_1,m_1(b)}\circ{\mathcal L}_{ib}^{{\omega},n_{i_1}}
\]
and using submultiplicativity of norms of operators, we deduce that $P$-a.s.
for any sufficiently large $n$ and real $b$ so that $b_1\leq |b|\leq n^{d}$ we have
\[
\|{\mathcal L}_{ib}^{{\omega},n}\|_{(b)}\leq C^{1+\frac{c_1 n}{\ln n}}e^{-\frac12c_1\rho_1{\gamma}_2\ln|b|\frac{n}{\ln n}}
=C\big(|b|^{\frac12\rho_1{\gamma}_2}C^{-1}\big)^{-\frac{c_1 n}{\ln n}}\leq
C2^{-\frac{c_1 n}{\ln n}}
\]
which implies that (\ref{EqA4}) from Assumption \ref{A4} holds true with $d_1=d$  and $K_{\omega}=b_1$ (since $|b|\leq n^d\leq 2^{2d\ln n}$).

Now we will explain the differences between the derivations of (\ref{Prop1}) and (\ref{RandomProp1}),
and show that a random variable $R_{\omega}$ with the properties required in Assumption \ref{A4} exists.
 First, in the above circumstances by Theorem 2.2 in \cite{Kifer-Thermo} we have $\boldsymbol{{\lambda}}_{\omega}(0)=1$
and $\boldsymbol{\nu}_{\omega}(0)=\text{Lebesgue}$.
By taking ${\alpha}=1$, the arguments in Chapter 5 of \cite{book} show that
$\|h_{\omega}(0)\|_{\textbf{BV}}$ and $\|1/h_{\omega}(0)\|_{\textbf{BV}}$ are both bounded random variables. Therefore, in order to show that Assumption \ref{A4} hold true with the operators ${\mathcal A}_{it}^{{\omega},n}$, it is sufficient to show that it holds true with the operators ${\mathcal L}_{it}^{{\omega},n}$ (and the norm $\|\cdot\|_2=\|\cdot\|_{\textbf{BV}}$).
Let
$\Lambda$ and $\tilde{\lambda}>1$ be positive constants so that $P$-a.s.
\[
\Lambda\geq\sup|f_{\omega}'|\geq\inf|f_{\omega}'|\geq\tilde{\lambda}.
\]
Then all the arguments in the proofs  in Section 2 of \cite{[3]} proceed similarly with the transfer operators ${\mathcal L}_{ib}^{{\omega},n}$ in place of $\mathscr L_{ib}^{n}$, since they only rely on the expansion properties of $f$ and $\tau$. In particular, $\|{\mathcal L}_{it}^{{\omega},n}\|_{\textbf{BV}}\leq R(1+|t|)$ for some constant $R>0$. Moreover, it is possible to construct  partitions ${\Omega}_n={\Omega}_n^{\omega},\,n\geq0$ almost exactly as in \cite{[3]} (see the paragraph preceding  Lemma 4 in \cite{[3]}), which in our case will  also depend on ${\omega}$: their recursive
construction in our nondeterministic setup proceeds exactly as in \cite{[3]}, expect that in the $n$-th step of the construction we consider the function $f_{{\theta}^{n}{\omega}}$ in place of $f$. The estimates derived in Lemma 4 in \cite{[3]} from there  hold true also for the ${\Omega}_n^{\omega}$'s, and the proof goes exactly in the same way as in \cite{[3]}.

Now we will explain how to generalize the ideas from Section 3 in \cite{[3]} to our non-deterministic setting.
Set $F_{\omega}(x,u)=(f_{\omega}(x),u+\tau_{\omega}(x))$, where $x,y\in{\mathbb T}$.
Then
\[
F^n_{\omega}(x,u):=F_{{\theta}^{n-1}{\omega}}\circ F_{{\theta}^{n-2}{\omega}}\circ\cdots\circ
F_{\omega}(x,u)=\big(f_{\omega}^n(x),u+\sum_{j=0}^{n-1}\tau_{{\theta}^j{\omega}}\circ f_{\omega}^j(x)\big)
\]
and the differential $DF_{\omega}(x,u)$ of $F_{\omega}$ does not depend on $u$.
 The fact that $f_{\omega}$ and $\tau_{\omega}$ are expanding uniformly in ${\omega}$ implies the  cones ${\mathcal K}=\{({\alpha},{\beta}): |\beta||{\alpha}|^{-1}\leq C_1\}$ are invariant under each $DF_{\omega}$. Here $C_1$ is a constant so that $P$-a.s.,
\[
C_1\geq\sup|\tau_{\omega}'|(\tilde{\lambda}-1)^{-1}.
\]
Thus, the notation of transversality from Section 3 in \cite{[3]} is naturally extended to our setup, using $DF^n_{\omega}(x)$ in place of the differential of $F^n(x,u)=\big(f^n(x),u+\sum_{j=0}^{n-1}\tau\circ f^j(x)\big)$.
That is, for any $y\in{\mathbb T}$, we will say that $x_1,x_2\in (f_{\omega}^n)^{-1}\{y\}$ are transversal if $DF^n_{\omega}(x_1){\mathcal K}\bigcap DF^n_{\omega}(x_2){\mathcal K}=\{0\}$.
Set $J_{{\omega},n}=|(f_{\omega}^n)'|^{-1}$, and
define the quantity $\varphi_{\omega}(n)$ by
\[
\varphi_{\omega}(n)=\sup_{y\in{\mathbb T}}\sup_{x_1\in (f_{\omega}^n)^{-1}\{y\}}
\sum_{x_2\in A_{\omega}(x_1,n,y)}J_{{\omega},n}(x_2)
\]
where $A_{\omega}(x_1,n,y)$ is the set of all preimages of $y$ by $f_{\omega}^n$ which are not transversal to $x_1$. We also define
the quantity $\tilde\varphi_{\omega}(n)$
by
\[
\tilde\varphi_{\omega}(n)=\sup_{y}\sup_{L}\sum_{\substack{
x\in (f_{\omega}^n)^{-1}\{y\}\\L\subset DF_{\omega}^n(x){\mathcal K}}}J_{{\omega},n}(x)\frac{\textbf{h}_{\omega}(x)}{\textbf{h}_{{\theta}^n{\omega}}(y)}
\]
where $L\subset{\mathbb R}^2$ is a line which passes through the origin
and $\textbf{h}_{\omega}=\textbf{h}_{\omega}(0)$ is the  random function from the RPF triplet, i.e. the one satsifying that ${\mathcal L}_0^{\omega} \textbf{h}_{\omega}=\textbf{h}_{{\theta}{\omega}}$ (recall that $\boldsymbol{{\lambda}}_{\omega}(0)=1$). Note that (see Chapter 5 in \cite{book}) the function $\textbf{h}_{\omega}$ is strictly positive and is bounded and bounded away from $0$ uniformly in ${\omega}$.
 Then $\tilde\varphi_{\omega}(n)\leq 1$ and it satisfies that
\[
\tilde\varphi_{\omega}(n+m)\leq \tilde\varphi_{\omega}(n)\cdot \tilde\varphi_{{\theta}^n{\omega}}(m).
\]
By Kingman's subadditive ergodic theorem, the limit $\lim_{n\to\infty}(\tilde\varphi_{\omega}(n))^{\frac1n}$ exists, $P$-a.s., and it does not depend on ${\omega}$. When this limit equals $1$ then, by considering the subadditive sequence $g_n=\int \ln\tilde\varphi_{\omega}(n)dP({\omega})$, it follows that $\tilde\varphi_{\omega}(n)=1$ for $P$-a.a. ${\omega}$ and all $n$'s.
Relying on that, the arguments from the proof of Lemma 7 in \cite{[3]} show that
\[
\limsup_{n\to\infty}\varphi_{\omega}(n)<1\,\,\text{for }\,P\text{-a.a. } {\omega}
\]
if there exists no random family $\eta_{\omega}$ of functions $\eta_{\omega}:{\mathbb T}\to{\mathbb T}$ so that $\tau_{\omega}-\eta_{{\theta}{\omega}}\circ f_{\omega}+\eta_{\omega}$ is piecewise constant ($P$-a.s.). The rest of the derivation of (\ref{RandomProp1}) proceeds now as in Section 4 of  \cite{[3]}, relying on the uniform (in ${\omega}$) expansion properties of the functions $f_{\omega}$ and $\tau_{\omega}$.

\begin{remark}
The assumption that $u_{\omega}$ is continuous can be dropped since all the results from \cite{book} hold true also when $u_{\omega}$ is piecewise H\"older continuous. We will not present here all the details and instead refer the readers to Assumption 5.2.1 in  \cite{book}, and to the results proceeding it.
\end{remark}

\subsubsection{Convergence rates towards the asymptotic moments}
In the circumstances of Section \ref{MomRateSec},
the operators ${\mathcal A}_z^{{\omega}}$ do not depend only on the $0$-th coordinate even when ${\mathcal L}_z^{\omega}$ depends on it (only). Still, in order to obtain the convergence rates stated in Section \ref{MixSec}, it is enough to obtain the appropriate almost sure rates in the convergence of the ergodic averages $n^{-1}\sum_{j=0}\Pi_{{\theta}^j{\omega}}^{(k)}(0)$ (see the arguments in Section \ref{MomSec}. Let us  denote by $\boldsymbol{\Pi}_{\omega}(z)$ the pressure function corresponding to $\boldsymbol{\lambda}_{{\omega}}(z)$. In (4.2) in \cite{SeqRPF} we have proved that there exist constants $c_1$ and $\rho_1$ so that  $P$-a.s. for any $z\in{\mathbb C}$ with $|z|<\rho_1$ and any $n\geq 1$ we have
\[
\left|\sum_{j=0}^{n-1}\Pi_{{\theta}^j{\omega}}(z)-\sum_{j=0}^{n-1}\boldsymbol\Pi_{{\theta}^j{\omega}}(z)\right|\leq c_1.
\]
Note that it is indeed possible to bound from below $\min |h_{\omega}(z)|$ when $z$ belongs to some deterministic neighborhood of the origin since $h_{\omega}(x)\geq C>0$ for some constant $C>0$ which does not depend on ${\omega}$ and $x$ (see (5.9.4) and (4.3.32) in \cite{book}).
Therefore, when ${\mathcal L}_z^{\omega}$ depend only on the $0$-th coordinate, we can just repeat the argument from Section \ref{MomRateSec}  with $\boldsymbol\Pi$ in place of $\Pi$ we obtain all the results stated in Section \ref{MixSec} under the mixing conditions from Section \ref{MomRateSec}.

\subsection{Markov chains with transition densities}\label{Densities}
For any ${\omega}\in{\Omega}$ denote by $B_{\omega}=H_1^{\omega}$ the Banach space of all bounded
Borel functions $g:{\mathcal E}_{\omega}\to{\mathbb C}$ together with the supremum norm $\|\cdot\|_\infty$.
For any $g:{\mathcal E}\to{\mathbb C}$ consider the functions $g_{\omega}:{\mathcal E}_{\omega}\to{\mathbb C}$ given by
$g_{\omega}(x)=g({\omega},x)$.
Then by Lemma 5.1.3 in \cite{book}, the norm ${\omega}\to\|g_{\omega}\|_\infty$ is a ${\mathcal F}$-measurable
function of ${\omega}$, for any measurable $g:{\mathcal E}\to{\mathbb C}$.

Let $r_{\omega}=r_{\omega}(x,y):{\mathcal E}_{\omega}\times{\mathcal E}_{{\theta}{\omega}}\to[0,\infty),\,{\omega}\in{\Omega}$ be a family of
integrable in $y$ Borel measurable functions, $m_{\omega},\,{\omega}\in{\Omega}$ be a family of Borel probability measures
on ${\mathcal E}_{\omega}$ and $u:{\mathcal E}\to{\mathbb R}$ be a measurable function so that $u_{\omega}\in B_{\omega},\,$ $P$-a.s.
and that the random variable $\sup|u_{\omega}|=\|u_{\omega}\|_\infty$ is bounded.
Consider the family of random operators $R_z^{\omega},\,z\in{\mathbb C}$ which map (bounded)
Borel functions $g$ on ${\mathcal E}_{{\theta}{\omega}}$ to Borel measurable functions on ${\mathcal E}_{\omega}$
by the formula
\begin{equation}\label{IntOp}
R^{\omega}_zg(x)=\int_{{\mathcal E}_{{\theta}{\omega}}} r_{\omega}(x,y)e^{zu_{{\theta}{\omega}}(y)}g(y)dm_{{\theta}{\omega}}(y).
\end{equation}
We will assume that
$R_0^{\omega}$ are Markov operators, namely that $R_0^{\omega} \textbf{1}=\textbf{1}$
where $\textbf{1}$ is the function which takes the constant value $1$ on ${\mathcal E}_{{\theta}{\omega}}$.
Observe that
\[
\|R_0^{\omega}\|_\infty:=
\sup_{g\in B_{{\theta}{\omega}}:\|g\|_\infty\leq 1}\|R_0^{\omega} g\|_\infty=\|R_0^{\omega}\textbf{1}\|_\infty
\]
and therefore for $P$-a.a. ${\omega}$
we have $\|R_z^{\omega}\|_\infty<\infty$ for any
$z\in{\mathbb C}$, namely,  $R_z^{\omega}$ is a continuous linear operator between the Banach spaces
$B_{{\theta}{\omega}}$ and $B_{\omega}$.

\begin{assumption}\label{Meas ass}
The functions $r_{\omega}$ are measurable in ${\omega}$ and the
the measures $m_{\omega}$ are measurable in the sense that
the map ${\omega}\to\int_{{\mathcal E}_{{\omega}}} g_{{\omega}}(x)dm_{{\theta}{\omega}}(x)$ is measurable for any
measurable function $g:{\mathcal E}\to{\mathbb C}$ so that $g_{\omega}(\cdot)$ is bounded $P$-a.s.
\end{assumption}
Under this assumption the map $({\omega},x)\to R_0^{\omega} g_{{\theta}{\omega}}(x)$ is measurable for any measurable function $g=g({\omega},x)=g_{\omega}(x)$.

For any ${\omega}\in{\Omega}$, $n\in{\mathbb N}$ and $z\in{\mathbb C}$ consider the
$n$-th order iterates $R_z^{{\omega},n}:B_{{\theta}^n{\omega}}\to B_{\omega}$
given by
\begin{equation}\label{Int Op oter 1}
R_z^{{\omega},n}=R_z^{\omega}\circ R_z^{{\theta}{\omega}}\circ\cdots\circ R_z^{{\theta}^{n-1}{\omega}}.
\end{equation}
Then we can write
\[
R_0^{{\omega},n}g(x)=\int _{{\mathcal E}_{{\theta}^n{\omega}}}r_{\omega}(n,x,y)g(y)dm_{{\theta}^n{\omega}}(y)
\]
for some family $r_{\omega}(n,\cdot,\cdot)=r_{\omega}(n,x,y):{\mathcal E}_{\omega}\times{\mathcal E}_{{\theta}^n{\omega}}\to[0,\infty)$
of integrable in $y$ Borel
measurable functions.
We will assume that the following random version of the two sided Doeblin condition holds true.
\begin{assumption}\label{Doeb}
There exist a bounded random variable $j_{\omega}\in{\mathbb N}$
and random variables $0<{\alpha}_m({\omega})\leq 1$, $m\in{\mathbb N}$ such that $P$-a.s.,
\begin{equation}\label{al j om}
{\alpha}_m({\omega})\leq r_{\omega}(m,x,y)\leq\big({\alpha}_m({\omega})\big)^{-1},
\end{equation}
for any $m\geq j_{\omega}\,$, $x\in{\mathcal E}_{\omega}$ and $y\in{\mathcal E}_{{\theta}^m{\omega}}$.
Moreover, let $j_0$ be so that $j_{\omega}\leq j_0$, $P$-a.s. Then there exists ${\alpha}>0$
so that ${\alpha}_{j_0}({\omega})\geq{\alpha}$.
\end{assumption}
Under the above assumptions, we showed in Chapter 6 of \cite{book} that the family of operators ${\mathcal A}_z^{\omega}=R_z^{{\theta}^{-1}{\omega}}$ satisfies Assumption \ref{BaseAss} with the measure preserving system $({\Omega},{\mathcal F},P,{\theta}^{-1})$.

 Let $\mu_{\omega}$ be any probability measure on ${\mathcal E}_{\omega}$ and let $\xi_n^{{\theta}^n{\omega}},\,n\geq1$
be the Markov chain with initial distribution $\mu_{\omega}$ whose $n$-th step operator is given by
$R_0^{{\omega},n}$. Set
\[
S_n^{\omega}=\sum_{j=0}^{n-1}u_{{\theta}^j{\omega}}(\xi_j^{{\theta}^j{\omega}}).
\]
Then,
\[
{\mathbb E} e^{itS_n^{\omega}}=\int\mu_{\omega}(R_{it}^{{\omega},n}\textbf{1})dP({\omega})
\]
and therefore all the results stated in Section \ref{sec2} hold true with the random variables $S_n^{{\theta}^{-n}{\omega}}$.

\subsubsection{Edgeworth expansions of order one: verification of Assumptions \ref{A3}}
The first part of Assumption \ref{A3} holds true with the supremum norm since $\|R_{it}^{\omega}\|_\infty\leq\|R_0^{\omega}\|_\infty=1$.
Let ${\omega}_0\in{\Omega}$ and $n_0\in{\mathbb N}$ be so that ${\theta}^{n_0}{\omega}_0={\omega}_0$. Suppose that for any $t\in{\mathbb R}\setminus\{0\}$ the spectral radius $r(t)$ of the operator $R_{it}^{{\omega}_0,n_0}$ is strictly less than $1$.
We refer the reader to Section 6.3 in \cite{Liverani} for conditions guaranteeing that $r(t)<1, t\not=0$ (see the verification of Assumption (A3) there). Then by Lemma 2.10.4 in  \cite{book}
 the second part Assumption \ref{A3} holds true when Assumption \ref{PerPointAss} is satisfied with $R_z^{\omega}$ in place of ${\mathcal L}_z^{\omega}$ and with the norm $\|\cdot\|_\infty$ instead of $\|\cdot\|_{{\alpha},\xi}$.
 Note that in Chapter 7 we have also verified Assumption \ref{A3} under certain type of non-lattice condition. These results were obtained under the assumption that $P(j_{\omega}=j_{{\theta}{\omega}}=1)>0$ (see Assumption 7.1.4 there), and it is also possible to extend them to more general cases (see Section 7.4.4 there).

\subsubsection{Edgeworth expansions of high orders: verification of Assumptions \ref{A4}}
We will show that Assumption \ref{A4} holds true with $H_2^{\omega}=B_{\omega}$  under the following
\begin{assumption}\label{AboveAss}
The system $({\Omega},{\mathcal F},P,{\theta})$ is ergodic and
the probability that
\begin{equation}\label{Cond}
\limsup_{|t|\to\infty}\|R_{it}^{{\omega},2}\|_\infty<1
\end{equation}
is positive.
\end{assumption}
Before verifying Assumption \ref{A4} we will provide conditions guaranteeing it holds true.
The following result is proved exactly as in \cite{Liverani}.
\begin{lemma}\label{Lem}
Let ${\omega}\in{\Omega}$ and suppose that ${\mathcal E}_{\omega}={\mathcal E}_{{\theta}{\omega}}$ is a compact connected smooth
manifold, that $u_{{\theta}{\omega}}$ and $u_{{\theta}^2{\omega}}$ are piecewise smooth and that $r_{\omega}$ is a $C^1$-function. Moreover, assume that
$j_{\omega}=j_{{\theta}{\omega}}=1$ and that the function $(y_1,y_2)\to u_{{\theta}{\omega}}(y_1)+u_{{\theta}^2{\omega}}(y_2)$ is not piecewise
constant. Then condition (\ref{Cond}) holds true with the above ${\omega}$.
\end{lemma}
Assumption \ref{Cond} holds true  when
the probability that ${\omega}$ satisfies the conditions of Lemma \ref{Lem} is positive.

Now we will verify Assumption \ref{A4}.
Under Assumption \ref{AboveAss}, there exist ${\varepsilon}_0,r_0>0$ so that $p_0=P({\Gamma}_0)>0$, where ${\Gamma}_0$ is the set of all ${\omega}$'s so that
\begin{equation*}
\sup_{t: |t|\geq r_0}\|R_{it}^{{\omega},2}\|_\infty\leq 1-{\varepsilon}_0.
\end{equation*}
Since ${\theta}$ is ergodic $P$-a.s. there exists a sequence $n_1<n_2<n_3<...$ of
positive integers so that
\[
\lim_{k\to\infty}\frac{n_k}{k}=\frac1{p_0}\,\,\text{ and }{\theta}^{-n_i}{\omega}\in{\Gamma}_0\,\,\forall\,i\in{\mathbb N}.
\]
Set $m_i=n_{2i}$, ${\omega}_i={\theta}^{-m_i}{\omega}$ and write
\begin{eqnarray*}
{\mathcal A}_{it}^{{\theta}{\omega},n}=R^{{\theta}^{-n+1}{\omega},n}_{it}=R_{it}^{{\theta}^{-n+1}{\omega}}\circ R_{it}^{{\theta}^{-n+2}{\omega}}\circ\cdots\circ
R_{it}^{{\omega}}\\
=R_{it}^{{\theta}^{2}{\omega}_{k_{n}},n-k_{n}-2}\circ R_{it}^{{\omega}_{k_{n}},2}\circ\cdots\circ R_{it}^{{\omega}_2,2}\circ R_{it}^{{\theta}^{2}{\omega}_1,m_2-m_1-2}R_{it}^{{\omega}_1,2}\circ R_{it}^{{\theta}^{2}{\omega}_2,m_1-2}
\end{eqnarray*}
where $k_n=\max\{i: m_i\leq n-3\}$. Since $\|R_{it}^{\omega}\|_\infty\leq 1$ and
$k_n/n$ converges to $\frac 12p_0$ as $n\to\infty$, we obtain from submultiplicativity of norms of operators that
$P$-a.s., for any $t\in{\mathbb R}$ so that $|t|\geq r_0$ and any sufficiently large $n$,
\[
\|{\mathcal A}_{it}^{{\omega},n}\|_\infty\leq (1-{\varepsilon}_0)^{k_{n-2}}\leq (1-{\varepsilon}_0)^{\frac14p_0n}
\]
which implies that Assumption \ref{A4} holds true. We note that Assumption \ref{A4} also
holds true with the norms $\|\cdot\|_\infty$ when appropriate random versions of conditions $(C)$ and $(B_k)$ from \cite{Neg2} are satisfied
(with positive probability).

redpajama_set_name: RedPajamaArXiv