\section{Introduction}\label{Sect1}

The classical renewal theory is an area of applied probability dealing with nondecreasing standard random walks and various derived processes like {\it renewal process, first passage time, overshoot, undershoot} etc. A good overview of the renewal theory can be found in \cite{Asmussen:2003}, \cite{Resnick:2002} and more recent accounts \cite{Iksanov:2016} and \cite{Mitov+Omey:2014}.

Assume that a general branching process (a.k.a. Crump--Mode--Jagers branching process) is generated by a standard random walk $S^{(1)}$ with nonnegative steps. Clearly, the random sequence $S^{(j)}$ defined by the birth times in the $j$th generation of the process ($j\geq 2$) is much more complicated than the standard random walk $S^{(1)}$ defining the birth times in the $1$st generation. It is natural to call $(S^{(j)})_{j\geq 2}$ {\it iterated standard random walk on a general branching process tree}. This should not be confused with iterated renewal processes treated in \cite{Sen:1981}. In this paper we initiate a systematic study of $S^{(j)}$ for $j\geq 2$ and its derived processes, our primary purpose being obtaining counterparts of the classical renewal-theoretic results. Actually, our setting will be a bit more general than that outlined above. We shall develop elements of renewal theory for iterated {\it perturbed} random walks rather than standard random walks, thereby making our results more general.

Now it is time to set the scene. Let $(\xi_i, \eta_i)_{i\in\mathbb{N}}$ be independent copies of a $\mathbb{R}^2$-valued random vector $(\xi, \eta)$ with arbitrarily dependent components. Denote by $(S_i)_{i\geq 0}$ the zero-delayed standard random walk with increments $\xi_i$ for $i\in\mathbb{N}$, that is, $S_0:=0$ and $S_i:=\xi_1+\cdots+\xi_i$ for $i\in\mathbb{N}$. Define
\begin{equation*}
T_i:=S_{i-1}+\eta_i,\quad i\in \mathbb{N}.
\end{equation*}
The sequence $T:=(T_i)_{i\in\mathbb{N}}$ is called {\it perturbed random walk} (PRW, in short). A survey of various results for the so defined PRWs can be found in the book \cite{Iksanov:2016}. An incomplete list of more recent papers addressing various aspects of the PRWs includes \cite{Alsmeyer+Iksanov+Marynych:2017, Duchamps+Pitman+Tang:2019, Iksanov+Pilipenko+Samoilenko:2017, Pitman+Tang:2019, Pitman+Yakubovich:2019, Rashytov:2018}.

In what follows we assume that $\xi$ and $\eta$ are almost surely (a.s.) positive. Put $N(t):=\sum_{i\geq 1}\1_{\{T_i\leq t\}}$ and $V(t):=\mathbb{E} N(t)$ for $t\geq 0$.  It is clear that
\begin{equation}\label{equ}
V(t)=\mathbb{E} U((t-\eta)_{+})=(U\ast G)(t)=\int_{[0,\,t]}U(t-y){\rm d}G(y),\quad t\geq 0,
\end{equation}
where, for $t\geq 0$, $U(t):=\sum_{i\geq 0}\mathbb{P}\{S_i\leq t\}$ is the renewal function and $G(t)=\mathbb{P}\{\eta\leq t\}$. As usual, $x_{+}:=\max(x,0)$. Here and in what follows we denote by $u \ast v$ the Lebesgue--Stieltjes convolution of two functions $u,v$ of locally bounded variation. We also use the notation $u^{\ast(j)}$, $j\in\mathbb{N}$, for the $j$th convolution power of $u$.

Now we provide more details about the construction of a general branching process (already mentioned at the beginning of the section) in the special case it is generated by $T$. At time $0$ there is one individual, the ancestor. The ancestor produces offspring (the first generation) with birth times given by the points of $T$. The first generation
produces the second generation. The shifts of birth times of the second generation
individuals with respect to their mothers' birth times are distributed according to copies of $T$, and for different mothers these copies are independent. The second generation produces the
third one, and so on. All individuals act independently of each other.

For $t\geq 0$ and $j\in\mathbb{N}$, denote by $N_j(t)$ the number of the $j$th generation individuals with birth times $\leq t$ and put $V_j(t):=\mathbb{E} N_j(t)$, and $V(t):=0$ for $t<0$. Then $N_1(t)=N(t)$, $V_1(t)=V(t)$ and
$$
V_j(t)=(V_{j-1}\ast V)(t)=\int_{[0,\,t]} V_{j-1}(t-y){\rm d}V(y),\quad j\geq 2,\quad t\geq 0.
$$
The basic decomposition that sheds light on the properties of $N_j:=(N_j(t))_{t\geq 0}$
and also demonstrates its recursive structure is
\begin{equation}\label{basic1232}
N_j(t)=\sum_{r\geq 1}N^{(r)}_{j-1}(t-T_r)\1_{\{T_r\leq t\}}=\sum_{k\geq 1}N^{(k)}_1(t-T^{(j-1)}_k)\1_{\{T^{(j-1)}_k\leq t\}} ,\quad j\geq 2,\quad t\geq 0,
\end{equation}
where $N_{j-1}^{(r)}(t)$ is the number of successors in the $j$th generation with birth times within $[T_r,t+T_r]$ of the first generation individual with birth time $T_r$; $T^{(j-1)}:=(T^{(j-1)}_k)_{k\geq 1}$ is some enumeration of the birth times in the $(j-1)$th generation; $N_1^{(k)}(t)$ is the number of children in the $j$th generation with birth times within $[T^{(j-1)}_k,t+T^{(j-1)}_k]$ of the $(j-1)$th generation individual with birth time $T^{(j-1)}_k$. By the branching property, $(N_{j-1}^{(1)}(t))_{t\geq 0}$, $(N_{j-1}^{(2)}(t))_{t\geq 0},\ldots$ are independent copies of $N_{j-1}$ which are also independent of $T$, and $(N_1^{(1)}(t))_{t\geq 0}$, $(N_1^{(2)}(t))_{t\geq 0},\ldots$ are independent copies of $(N(t))_{t\geq 0}$ which are also independent of $T^{(j-1)}$.  Note that, for $j\geq 2$, $N_j$ is a particular instance of a random process with immigration at random times (the term was introduced in \cite{Dong+Iksanov:2020}, see also \cite{Iksanov+Rashytov:2020}).

Our motivation behind introducing the iterated perturbed random walks is at least three-fold.

\begin{enumerate}[leftmargin=*]
\item [1)] For each integer $j\geq 2$, the sequence $T^{(j)}$ and the process $N_j$ are a natural generalization of the perturbed random walk $T$ and the counting process $(N(t))_{t\geq 0}$. It is interesting to investigate to which extent the renewal-theoretic properties of $T$ and $(N(t))$ are inherited by $T^{(j)}$ and $N_j$.
Thus, the activity undertaken in the present article can be thought of as the development of renewal theory for the iterated perturbed random walks.

\item[2)] The sequence $(T^{(j)})_{j\in\mathbb{N}}$ is a particular instance of a branching random walk in which the first generation point process is $(N(t))_{t\geq 0}$, the counting process of a perturbed random walk. Alternatively, and this is our preferable viewpoint, for $j\in\mathbb{N}$, $T^{(j)}$ can be interpreted as the sequence of birth times in the $j$th generation of a general branching process.
Therefore, the results of the present article contribute towards better understanding of how the births occur within a particular generation.  Being of intrinsic interest for the theory of general branching processes, this information also sheds light on the organization of levels (the sets of vertices located at the same distance from the root) of some random trees (for instance, random recursive trees and binary search trees) that can be constructed as family trees of general branching processes stopped at suitable random times. We refer to \cite{Holmgren+Janson:2017} for more details and examples of the embeddable random trees.

\item[3)] Renewal theory for perturbed random walks is an inevitable ingredient for investigation of nested occupancy scheme in random environment generated by stick-breaking. Referring to  \cite{Buraczewski+Dovgay+Iksanov:2020, Iksanov+Marynych+Samoilenko:2020} for more details we only mention that the latter scheme is a generalization of the classical Karlin infinite balls-in-boxes occupancy scheme \cite{Gnedin+Hansen+Pitman:2007, Karlin:1967}. Unlike the Karlin scheme in which the collection of boxes is unique, there is a nested hierarchy of boxes, and the hitting probabilities of boxes are defined in terms of iterated stick-breaking. Assuming that $n$ balls have been thrown, denote by $K_n(j)$ the number of occupied boxes in the $j$th level which is the basic object of interest. It turns out that whenever $j=j_n=o((\log n)^{1/2})$ (the case of fixed $j$ is included) the distributional behavior of $K_n(j)$ as $n\to\infty$ is the same as that of $N_j(\log n)$, when the underlying perturbed random walk $T$ is appropriately chosen.

\end{enumerate}

We call the $j$th generation {\it early}, {\it intermediate} or {\it late} depending on whether $j$ is fixed, $j=j(t)\to\infty$ and $j(t)=o(t)$ as $t\to\infty$, or $j=j(t)$ is of order $t$. In view of Proposition \ref{bigg} given in Section \ref{sec:height} there are no other generations. Assume, for the time being, that $j$ is a late generation and that $T$ is a collection of random points, not necessarily the perturbed random walk. Nevertheless, we retain the notation $N_j$ and $V_j$. In this case the asymptotic behavior of $V_j$ and $N_j$ is well-understood. For instance, a delicate counterpart of the key renewal theorem for $V_j$ which includes both a version of the elementary renewal theorem and a version of Blackwell's theorem can be found in Theorem A of \cite{Biggins:1979}. For the corresponding a.s.\ result for $N_j$, see Theorem B of the same paper and Theorem 4 in \cite{Biggins:1992}. A strong law of large numbers for $N_j(bj)$ for appropriate $b>0$ is given in formula (1.1) of \cite{Biggins:1979}. From these and the other results of this flavor it follows that $N_j$ forgets what was happening in the early history and particularly in the $1$st generation. The behavior of $N_j$ is universal for a wide class of input processes (responsible for the $1$st generation). It is driven by limit theorems available for general branching processes like convergence of the Biggins martingales, large deviations etc.

While the present paper deals with some intermediate generations, the early generations which admit a much simpler analysis will be treated in a separate paper \cite{Iksanov etal:2021}. One may expect that the behavior of the iterated perturbed random walks in the early and intermediate generations is very different from that in the late generations. When $j$ is a non-late generation, the process $N_j$ should inherit, for the most part, the properties of $N$, possibly in a modified form. This statement is confirmed by counterparts of the elementary renewal theorem (Theorems \ref{thm:elem1} and \ref{thm:elem2}), the key renewal theorem (Theorem \ref{thm:keyren}) and Blackwell's theorem (Corollary \ref{black}) which are our main results.

The remainder of the paper is structured as follows. Our main findings are formulated in Section \ref{results} and then proved in Section \ref{proofs}. Also, Section \ref{results} contains two previously known results concerning $N_j$ and $V_j$. To our knowledge, all the results presented in this paper form the state-of-the-art as far as the intermediate generations of the iterated perturbed random walks are concerned. Finally, the appendix collects a rate of convergence result and counterparts of Blackwell's theorem and the key renewal theorem for the perturbed random walks.

\section{Results}\label{results}

\subsection{The height of a confined general branching process tree}\label{sec:height}

For $t>0$, put $$H(t):=\inf\{j\in\mathbb{N}: N_j(t)=0\}$$ and note that $N_j(t)=0$ a.s.\ for all $j\geq H(t)$. We call the variable $H(t)$ the {\it height of a general branching process tree} generated by a perturbed random walk $T$ and confined to the strip $[0,t]$. The result given below is of principal importance for our classification of generations (early, intermediate, late).
\begin{assertion}\label{bigg}
For each $t\geq 0$, $H(t)<\infty$ a.s. Furthermore,
\begin{equation}\label{height}
\lim_{t\to\infty}\frac{H(t)}{t}=\frac{1}{\gamma}\in (0,\infty)\quad\text{{\rm a.s.}},
\end{equation}
where $\gamma:=\sup\{z>0: \mu(z)<1\}$ and $\mu(z):=\inf_{s>0}(e^{zs}\frac{\mathbb{E} e^{-s\eta}}{1-\mathbb{E} e^{-s\xi}})$ for $z>0$.
\end{assertion}
\begin{proof}
By assumption, $\mathbb{P}\{\eta=0\}=0$. This entails $\lim_{s\to\infty}\frac{\mathbb{E} e^{-s\eta}}{1-\mathbb{E} e^{-s\xi}}=0$ and thereupon
$$
\lim_{z\to 0+}\mu(z)=0.
$$
Also, $\lim_{z\to\infty}\mu(z)=\lim_{s\to 0+}\frac{\mathbb{E} e^{-s\eta}}{1-\mathbb{E} e^{-s\xi}}=\infty$. This shows that $\gamma\in (0,\infty)$.

Recall that, for $n\in\mathbb{N}$, $(T_r^{(n)})_{r\in\mathbb{N}}$ denotes some enumeration of the birth times in the $n$th generation of the general branching process. Put $B(n):=\inf_{r\geq 1} T_r^{(n)}$. By the famous Biggins result (Corollary on p.~635 in \cite{Biggins:1977}),
\begin{equation}\label{leftmost}
\lim_{n\to\infty}\frac{B(n)}{n}=\gamma\quad\text{a.s.}
\end{equation}
Since, for $n\in\mathbb{N}$ and $t>0$, $\{H(t)>n\}=\{B(n)\leq t\}$ and, according to \eqref{leftmost}, $\lim_{n\to\infty}B(n)=+\infty$ a.s., we infer $H(t)<\infty$ a.s.

Finally, we have $B(H(t))>t\geq B(H(t)-1)$ a.s. The left-hand inequality ensures $\lim_{t\to\infty}H(t)=+\infty$ a.s.\ which together with \eqref{leftmost} proves \eqref{height} with the help of a standard sandwich argument.
\end{proof}

It is seldom possible to find the constant $\gamma$ explicitly. Here is one happy exception. Let $(\xi,\eta)=(|\log W|, |\log(1-W)|)$, where $W$ has a uniform distribution on $[0,\,1]$. The distribution of the sequence $(e^{-T_i})_{i\in\mathbb{N}}$ is known as the {\it Griffiths--Engen--McCloskey distribution} with parameter $1$. In this case, $\mu(z)=ez$ for $z>0$ which gives $\gamma=e^{-1}$.

\subsection{Counterparts of the elementary renewal theorem for intermediate generations}\label{subsec:elem}

The simplest result of the renewal theory, called the elementary renewal theorem, tells us that
$$
U(t)=\sum_{i\geq 0}\mathbb{P}\{S_i\leq t\}~\sim~\frac{t}{{\tt m}},\quad t\to\infty,
$$
where ${\tt m}:=\mathbb{E} \xi<\infty$. Here and hereafter, the notation $f(t)\sim g(t)$ means that the ratio $f(t)/g(t)$ tends to $1$ as $t\to\infty$.

From \eqref{equ} it follows that, without any assumptions on $\eta$,
\begin{equation}\label{eq:elem_v_1}
V(t)~\sim~\frac{t}{{\tt m}},\quad t\to\infty.
\end{equation}
This is a counterpart of the elementary renewal theorem for the perturbed random walks.

In this section we state two results on the first-order behavior of the convolutions powers $V_j$ of $V$. Our first result, Theorem \ref{thm:elem1}, deals with `early intermediate' generations satisfying $j=j(t)\to\infty$ and $j(t)=o(t^{1/2})$ as $t\to\infty$ as well as early generations. At this point we stress that even though both Theorem \ref{thm:elem1} and Theorem \ref{thm:elem2} hold true for early generations, the assumptions of these theorems are too restrictive as far as early generations are concerned. We refer to the forthcoming article \cite{Iksanov etal:2021} for a proper version of the elementary renewal theorem in early generations. Recall the standard notation $x\wedge y=\min (x,y)$ for $x,y\in\mathbb{R}$.
\begin{thm}\label{thm:elem1}
Assume that either (i) $\mathbb{E}\xi^r<\infty$ for some $r\in (1,\,2]$ or (ii) $\mathbb{P}\{\xi>t\}~\sim~bt^{-r}$ for some $r\in (1,\,2)$ and some $b>0$. Suppose further that $\mathbb{E}(\eta\wedge t)=O(t^{2-r})$ as $t\to\infty$ with the same $r$ as in (i) or (ii). Then, for any integer-valued function $j=j(t)$ satisfying
$j(t)=o(t^{(r-1)/2})$ as $t\to\infty$,
\begin{equation}\label{eq:eq_elem1}
V_j(t)~\sim~\frac{t^j}{{\tt m}^j j!},\quad t\to\infty,
\end{equation}
where ${\tt m}=\mathbb{E}\xi<\infty$.
\end{thm}
\begin{rem}
The condition $\mathbb{E}\eta^{r-1}<\infty$ is sufficient for $\mathbb{E}(\eta\wedge t)=O(t^{2-r})$, $t\to\infty$. This follows from
\begin{multline*}
\mathbb{E}(\eta\wedge t)=\int_0^t \mathbb{P}\{\eta>y\}{\rm d}y\leq \int_0^t \Big(\frac{t}{y}\Big)^{2-r}\mathbb{P}\{\eta>y\}{\rm d}y\\
=t^{2-r}\int_0^\infty y^{r-2}\mathbb{P}\{\eta>y\}{\rm d}y=(r-1)^{-1}\mathbb{E}\eta^{r-1}t^{2-r}.
\end{multline*}
\end{rem}

Specializing Theorem \ref{thm:elem1} to $r=2$ gives the following corollary which has already been obtained
via a slightly different argument in formula (4.6) of \cite{Buraczewski+Dovgay+Iksanov:2020}.

\begin{cor}\label{cor:r=2}
Assume that $\mathbb{E}\xi^2<\infty$ and $\mathbb{E} \eta<\infty$. Then relation \eqref{eq:eq_elem1} holds for any integer-valued function $j=j(t)$ satisfying $j(t)=o(t^{1/2})$ as $t\to\infty$.
\end{cor}

Given next is a quite surprising result which shows that the convolution power $V_j$ exhibits a phase transition in the generations $j$ satisfying $j=j(t)\sim {\rm const}\,\cdot t^{1/2}$ as $t\to\infty$. Here, further moment and smoothness assumptions seem to be indispensable. In particular, we assume that the distribution of $\xi$ is spread-out, which means that some convolution power of the distribution function $t\mapsto \mathbb{P}\{\xi\leq t\}$ has an absolutely continuous component.

\begin{thm}\label{thm:elem2}
Assume that the distribution of $\xi$ is spread-out,
that $\mathbb{E} \xi^3<\infty$ and $\mathbb{E} \eta^2<\infty$. Then,
for any integer-valued function $j=j(t)$ satisfying
$j(t)=o(t^{2/3})$ as $t\to\infty$,
$$
V_j(t)~\sim~ \frac{t^j}{{\tt m}^j j!}\exp{\left(\frac{\gamma_0{\tt m}j^2}{t}\right)},\quad t\to\infty,
$$
where
\begin{equation}\label{eq:gamma_0_def}
\gamma_0:=\int_{[0,\,\infty)}{\rm d}(V(y)-{\tt m}^{-1}y)=\lim_{t\to\infty}(V(t)-{\tt m}^{-1}t)=\frac{\mathbb{E}\xi^2}{2{\tt m}^2}-\frac{\mathbb{E}\eta}{{\tt m}}
\end{equation}
may be positive, negative or zero.
\end{thm}
\begin{rem}
Assume that $(\xi, \eta)=(|\log W|, |\log(1-W)|)$, where $W$ is a random variable having a uniform distribution on $[0,1]$. Then
$$
V_j(t)=\frac{t^j}{j!},\quad j\in\mathbb{N},\quad t\geq 0.
$$
This is in line with the asymptotics provided by Theorem \ref{thm:elem2}, for, in this case, $\gamma_0=0$ and ${\tt m}=1$.
\end{rem}

\subsection{Counterparts of the key renewal theorem and Blackwell's theorem for intermediate generations}

In the renewal theory the key renewal theorem is usually obtained as a corollary to Blackwell's theorem. We proceed differently by first proving a counterpart of the key renewal theorem (Theorem \ref{thm:keyren}) and then obtain a counterpart of Blackwell's theorem (Corollary \ref{black}) as a corollary. Recall that the distribution of $\xi$  is nonlattice if it is not concentrated on any centered lattice of the form $(dn)_{n\in\mathbb{N}_0}$ for some $d>0$.
\begin{thm}\label{thm:keyren}
Let $f: [0,\infty)\to [0,\infty)$ be a directly Riemann integrable (dRi) function on $[0,\infty)$. Assume that either (a) or (b)
below holds true:
\begin{itemize}
\item[(a)] the distribution of $\xi$ is nonlattice, the conditions of Theorem \ref{thm:elem1} hold for some $r\in (1,\,2]$ (or some $r\in (1,\,2)$) and $j(t)=o(t^{(r-1)/2})$ as $t\to\infty$;
\item[(b)] the conditions of Theorem \ref{thm:elem2} hold and $j(t)=o(t^{2/3})$ as $t\to\infty$.
\end{itemize}
Then
\begin{equation}\label{eq:keyren}
(f\ast V_j)(t)=\int_{[0,\,t]}f(t-y){\rm d}V_j(y)~\sim~ \Big(\frac{1}{{\tt m}}\int_0^\infty f(y){\rm d}y\Big) V_{j-1}(t),\quad t\to\infty,
\end{equation}
where ${\tt m}=\mathbb{E}\xi<\infty$	, and $V_{j-1}(t)$ on the right-hand side can be replaced with $t^{j-1}/(m^{j-1} (j-1)!)$ in the case (a),
or with $t^{j-1}/(m^{j-1} (j-1)!)\exp{(\gamma_0{\tt m}j^2/t)}$ in the case (b).
\end{thm}

Upon taking $f(y)=\1_{[0,\,h)}(y)$ in Theorem \ref{thm:keyren} we immediately obtain the following.

\begin{cor}\label{black}
Let $h>0$ be fixed. Under the assumptions of Theorem \ref{thm:keyren}
\begin{equation}\label{bla_early}
V_j(t+h)-V_j(t)~\sim~\frac{h}{m}V_{j-1}(t),\quad t\to\infty.
\end{equation}
\end{cor}

\subsection{A couple of previously known results}

In this section we collect two previously known facts concerning the asymptotic behaviour of $N_j$ in the intermediate generations. They are borrowed from \cite{Iksanov+Marynych+Samoilenko:2020} and stated here for integrity and the reader's convenience. We write ${\overset{{\rm f.d.d.}}\longrightarrow}$ to denote weak convergence of finite-dimensional distributions.

\begin{thm}[Multivariate central limit theorem for $(N_j(t))_{t\geq 0}$]\label{main5}
Assume that ${\tt s}^2={\rm Var}\,\xi\in (0,\infty)$ and $\mathbb{E} \eta<\infty$. Let $j=j(t)$ be any positive integer-valued function satisfying $j(t)\to \infty$ and $j(t)=o(t^{1/2})$ as $t\to\infty$. Then, as $t\to\infty$,
\begin{equation}\label{clt22}
\left(\frac{\lfloor j(t)\rfloor^{1/2}(\lfloor j(t)u\rfloor-1)!}{({\tt s}^2{\tt m}^{-2\lfloor j(t)u\rfloor-1}t^{2\lfloor j(t)u\rfloor-1})^{1/2}}\bigg(N_{\lfloor j(t)u\rfloor}(t)-V_{\lfloor j(t)u \rfloor}(t)\bigg)\right)_{u>0}\\
{\overset{{\rm f.d.d.}}\longrightarrow}~ \Bigg(\int_{[0,\,\infty)}e^{-uy}{\rm d}B(y)\Bigg)_{u>0},
\end{equation}
where $(B(v))_{v\geq 0}$ is a standard Brownian motion.
\end{thm}

According to Proposition 3.1, Theorems 3.2 and 3.3 in \cite{Buraczewski+Dovgay+Iksanov:2020}, the centering $V_{\lfloor j(t)u \rfloor}(t)$ in \eqref{clt22} can be replaced by its leading term
$$
t^{\lfloor j(t)u \rfloor}/((\lfloor j(t)u \rfloor)!{\tt m}^{\lfloor j(t)u \rfloor}),
$$
provided that $j(t)=o(t^{1/3})$. For functions $t\mapsto j(t)$ which grow faster, this is not always the case. Plainly, the possibility/impossibility of such a replacement is justified by a second-order behavior of $V_j$. It should come as no surprise that second-order results for $V_j$ require more restrictive assumptions on the distributions of $\xi$ and $\eta$ than the corresponding first-order results. The following proposition which is concerned with the rate of convergence in the elementary renewal theorem for $V_j$ was proved in Proposition 8.1 of \cite{Iksanov+Marynych+Samoilenko:2020}.
\begin{assertion}\label{prop:sec_order_behaviour}
Assume that the distribution of $\xi$ has an absolutely continuous component, that $\mathbb{E} e^{\beta_1\xi}<\infty$, $\mathbb{E} e^{\beta_2\eta}<\infty$ for some $\beta_1,\beta_2>0$ and
$$
\gamma_0=\frac{\mathbb{E}\xi^2}{2{\tt m}^2}-\frac{\mathbb{E} \eta}{{\tt m}}>0.
$$ Then
\begin{equation}\label{prec}
V_j(t)-\frac{t^j}{j!{\tt m}^j}~\sim~\frac{\gamma_0 jt^{j-1}}{(j-1)!{\tt m}^{j-1}},\quad t\to\infty
\end{equation}
whenever $j=j(t)=o(t^{1/2})$ as $t\to\infty$ ($j$ is allowed to be fixed).
\end{assertion}
Formula \eqref{prec} can be thought of as a generalization of formulae \eqref{eq:gamma_0_def} and \eqref{eq:sec_order_U}. These provide the second-order behaviour of the functions $V$ and $U$, respectively.

\section{Proofs}\label{proofs}

\subsection{Preparatory results}\label{Prep}

Recall that $U$ denotes the renewal function for $(S_n)_{n\in\mathbb{N}_0}$. According to Lorden's inequality which holds whenever $\mathbb{E}\xi^2<\infty$,
\begin{equation}\label{lord}
U(t)-{\tt m}^{-1}t \leq c_0,\quad t\geq 0,
\end{equation}
where $c_0:=\mathbb{E} \xi^2/{\tt m}^2$ and ${\tt m}=\mathbb{E}\xi<\infty$. See \cite{Carlsson+Nerman:1986} for a nice proof under the assumption that the distribution of $\xi$ is nonlattice. Let $S_0^\ast$ be a random variable with distribution $\mathbb{P}\{S_0^\ast\in{\rm d}x\}={\tt m}^{-1}\mathbb{P}\{\xi>x\}\1_{(0,\infty)}(x){\rm d}x$. The basic formula (2) of the last cited paper which reads
\begin{equation}\label{eq:S0}
\mathbb{E} U(t-S_0^\ast)={\tt m}^{-1}t,\quad t\geq 0,
\end{equation}
holds true in the lattice case as well. Thus, the argument given in \cite{Carlsson+Nerman:1986} proves \eqref{lord} in general. We also note that
\begin{equation}\label{eq:sec_order_U}
\lim_{t\to\infty}\left(U(t)-{\tt m}^{-1}t\right)=\frac{\mathbb{E}\xi^2}{2{\tt m}^2}
\end{equation}
whenever the distribution of $\xi$ is nonlattice and $\mathbb{E}\xi^2<\infty$.

Since $V(t)\leq U(t)$ for $t\geq 0$ we infer
\begin{equation}\label{lord1}
V(t)-{\tt m}^{-1}t \leq c_0,\quad t\geq 0.
\end{equation}
On the other hand, assuming that $\mathbb{E}\eta<\infty$ (whereas the assumption $\mathbb{E}\xi^2<\infty$ is not needed here),
\begin{eqnarray*}
V(t)-{\tt m}^{-1}t&=&\int_{[0,\,t]}(U(t-y)-{\tt m}^{-1}(t-y)){\rm d}G(y)\\&\hphantom{==}-& {\tt m}^{-1} \int_0^t (1-G(y)){\rm d}y\geq-{\tt m}^{-1}\int_0^t (1-G(y)){\rm d}y\geq
-{\tt m}^{-1}\mathbb{E}\eta
\end{eqnarray*}
having utilized $U(t)\geq {\tt m}^{-1}t$ for $t\geq 0$ which is a consequence of Wald's identity $t\leq \mathbb{E} S_{\nu(t)}={\tt m} U(t)$, where $\nu(t):=\inf\{k\in\mathbb{N}: S_k>t\}$ for $t\geq 0$. Thus, we have shown that,  under the assumptions $\mathbb{E}\xi^2<\infty$ and $\mathbb{E}\eta<\infty$,
\begin{equation}\label{lord2}
|V(t)-{\tt m}^{-1}t|\leq c_L,\quad t\geq 0
\end{equation}
where $c_L=\max(c_0, {\tt m}^{-1}\mathbb{E}\eta)$.

\subsection{Results on convolution powers of functions of a linear growth and proofs of Theorems \ref{thm:elem1} and \ref{thm:elem2}}

The results presented here are concerned with the following purely analytic problem. Assume that a nondecreasing function $f$ exhibits a linear growth, that is, $f(t)~\sim~at$ as $t\to\infty$ for some $a>0$. Then, for fixed $j\in\mathbb{N}$,
$$
f^{\ast(j)}(t)~\sim~\frac{a^j t^j}{j!},\quad t\to\infty.
$$
Imposing various assumptions on the behavior of $f(t)-at$ we shall extend this asymptotics to the case when $j=j(t)$ diverges to infinity as $t\to\infty$.

\begin{assertion}\label{prop:convolutions1}
Let $f:\mathbb{R}\to [0,\,\infty)$ be a nondecreasing right-continuous function vanishing on the negative half-line and satisfying
\begin{equation}\label{eq:asymp_exp_assump}
f(t)=at+O(t^{\alpha}),\quad t\to\infty
\end{equation}
for some $a>0$ and $\alpha\in [0,1)$.
Then, for any integer-valued function $j=j(t)$ such that
$j(t)=o(t^{(1-\alpha)/2})$ as $t\to\infty$,
$$
f_j(t):=f^{\ast(j)}(t)~\sim~\frac{a^j t^j}{j!},\quad t\to\infty.
$$
\end{assertion}
\begin{proof}
According to \eqref{eq:asymp_exp_assump} there exists $C\geq 1$ such that
\begin{equation}\label{eq:uniform_bound_f}
-C(t+1)^{\alpha} \leq f(t)-at \leq C(t+1)^{\alpha},\quad t\geq 0.
\end{equation}
For $j\in\mathbb{N}$ and $t\geq 0$, put
$$
r_j(t):=\int_{[0,\,t]}f_j(t-y){\rm d}(f(y)-ay)=\int_{[0,\,t]}\left(f(t-y)-a(t-y)\right){\rm d}f_{j}(y)
$$
and note that
$$
f_j(t)=r_{j-1}(t)+a\int_0^t f_{j-1}(y){\rm d}y,\quad j\geq 2,~~ t\geq 0.
$$
By virtue of \eqref{eq:uniform_bound_f}, we conclude that
$$
|r_j(t)|\leq C(t+1)^{\alpha}f_j(t),\quad j\in\mathbb{N},~~ t\geq 0.
$$
Using this bound and the mathematical induction we obtain
\begin{equation}\label{eq:ineq}
W_j^{-}(t)\leq f_j(t)\leq W_j^{+}(t),\quad j\in\mathbb{N},~~t\geq 0,
\end{equation}
where $W_j^{\pm}$ is defined recursively by $W_0^{\pm}(t):=1$ and
$$
W_j^{\pm}(t)=\left(\pm C(t+1)^{\alpha}W^{\pm}_{j-1}(t)+a\int_{0}^t W_{j-1}^{\pm}(y){\rm d}y\right)_{+},\quad j\in\mathbb{N},~~ t\geq 0.
$$
Here, we recall that $x_{+}=\max(x,0)$ and note that taking the nonnegative part is only relevant for $W_j^{-}$ ensuring its nonnegativity, whereas it can be omitted for $W_j^+$.

It remains to show that
\begin{equation}\label{eq:asymp_for_bounds}
W_j^{\pm}(t)~\sim~\frac{a^j t^j}{j!},\quad t\to\infty.
\end{equation}
To this end, we first prove by induction that
\begin{equation}\label{eq:bound_w_j_plus}
W_j^{+}(t)\leq \frac{a^j t^j}{j!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(t+1)^{\alpha(j-i)+i}}{i!},\quad j\in\mathbb{N},~~ t\geq 0.
\end{equation}
While for $j=1$ this follows immediately because $W_1^{+}(t)=C(t+1)^{\alpha}+at$,
the induction step works as follows
\begin{align*}
W_{j+1}^+(t)&\leq C(t+1)^{\alpha}\left(\frac{a^j t^j}{j!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(t+1)^{\alpha(j-i)+i}}{i!}\right)\\
&+a\int_0^{t}\left(\frac{a^j y^j}{j!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(y+1)^{\alpha(j-i)+i}}{i!}\right){\rm d}y\\
&\leq \sum_{i=0}^{j}\binom{j}{i}\frac{a^i C^{j+1-i}(t+1)^{\alpha(j+1-i)+i}}{i!}+\frac{a^{j+1}t^{j+1}}{(j+1)!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^{i+1}C^{j-i}}{i!}\frac{(t+1)^{\alpha(j-i)+i+1}}{\alpha(j-i)+i+1}\\
&\leq \frac{a^{j+1}t^{j+1}}{(j+1)!}+\sum_{i=0}^{j}\binom{j}{i}\frac{a^i C^{j+1-i}(t+1)^{\alpha(j+1-i)+i}}{i!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^{i+1}C^{j-i}}{(i+1)!}(t+1)^{\alpha(j-i)+i+1}\\
&= \frac{a^{j+1}t^{j+1}}{(j+1)!}+\sum_{i=0}^{j}\binom{j}{i}\frac{a^i C^{j+1-i}(t+1)^{\alpha(j+1-i)+i}}{i!}+\sum_{i=1}^{j}\binom{j}{i-1}\frac{a^iC^{j+1-i}}{i!}(t+1)^{\alpha(j+1-i)+i}\\
&= \frac{a^{j+1}t^{j+1}}{(j+1)!}+\sum_{i=0}^{j}\binom{j+1}{i}\frac{a^i C^{j+1-i}(t+1)^{\alpha(j+1-i)+i}}{i!}
\end{align*}
having utilized the binomial identity $\binom{j}{i}+\binom{j}{i-1}=\binom{j+1}{i}$ for the last step. Further,
\begin{align*}
&\hspace{-1cm}\frac{j!}{a^j t^j}\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(t+1)^{\alpha(j-i)+i}}{i!}\sim\frac{j!}{a^j (t+1)^j}\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(t+1)^{\alpha(j-i)+i}}{i!}\\
&\leq \sum_{i=0}^{j-1}\left(\frac{j!}{i!}\right)^{2}\left(\frac{C}{a}\right)^{j-i}(t+1)^{(1-\alpha)(i-j)} \leq \sum_{i=0}^{j-1}(j^{j-i})^2(Ca^{-1})^{j-i}(t+1)^{(1-\alpha)(i-j)}\\
&\leq \sum_{i\geq 1} \left(\frac{Ca^{-1} j^2}{(t+1)^{1-\alpha}}\right)^i=\frac{Ca^{-1} j^2}{(t+1)^{1-\alpha}}\left(1-\frac{Ca^{-1} j^2}{(t+1)^{1-\alpha}}\right)^{-1}.
\end{align*}
Thus, in view of the assumption $j(t)=o(t^{(1-\alpha)/2})$ we have
\begin{equation}\label{eq:asymp_for_bounds1}
\limsup_{t\to\infty}\frac{j!}{a^j t^j}W_j^{+}(t)\leq 1.
\end{equation}
To prove that
\begin{equation}\label{eq:bound_w_j_minus}
W_j^{-}(t)\geq \left(\frac{a^j t^j}{j!}-\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(t+1)^{\alpha(j-i)+i}}{i!}\right)_{+},\quad j\in\mathbb{N},~~ t\geq 0
\end{equation}
we use a similar reasoning. While \eqref{eq:bound_w_j_minus} is obviously true for $j=1$, we obtain with the help of induction, for $j\geq 2$,
\begin{align*}
&\hspace{-0.3cm}W_{j+1}^-(t)\geq -C(t+1)^{\alpha}W_{j}^{-}(t)+a\int_0^t W_{j}^{-}(y){\rm d}y\geq -C (t+1)^\alpha
W_{j}^{+}(t)+a\int_0^t W_{j}^{-}(y){\rm d}y\\
&\geq -C(t+1)^{\alpha}W_{j}^{+}(t)+a\int_0^t \left(\frac{a^j y^j}{j!}-\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^iC^{j-i}(y+1)^{\alpha(j-i)+i}}{i!}\right){\rm d}y\\
&\geq -C(t+1)^{\alpha}\left(\frac{a^j t^j}{j!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^iC^{j-i}(t+1)^{\alpha(j-i)+i}}{i!}\right)\\
&+a\int_0^t \left(\frac{a^j y^j}{j!}-\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^i C^{j-i}(y+1)^{\alpha(j-i)+i}}{i!}\right){\rm d}y\\
&\geq \frac{a^{j+1}t^{j+1}}{(j+1)!}-\left(\sum_{i=0}^{j}\binom{j}{i}\frac{a^i C^{j+1-i}(t+1)^{\alpha(j+1-i)+i}}{i!}+\sum_{i=0}^{j-1}\binom{j}{i}\frac{a^{i+1}C^{j-i}}{i!}\int_0^t (y+1)^{\alpha(j-i)+i}{\rm d}y\right)\\
&\geq \frac{a^{j+1}t^{j+1}}{(j+1)!}-\sum_{i=0}^{j}\binom{j+1}{i}\frac{a^i C^{j+1-i}(t+1)^{\alpha(j+1-i)+i}}{i!}.
\end{align*}
We have used \eqref{eq:ineq} and \eqref{eq:bound_w_j_plus} for the second and the fourth inequality, respectively. Since $W_{j+1}^-$ is nonnegative we arrive at \eqref{eq:bound_w_j_minus}. Thus,
\begin{equation}\label{eq:asymp_for_bounds2}
\liminf_{t\to\infty}\frac{j!}{a^j t^j}W_j^{-}(t)\geq 1.
\end{equation}
Combining \eqref{eq:asymp_for_bounds1} and \eqref{eq:asymp_for_bounds2} yields \eqref{eq:asymp_for_bounds}, thereby finishing the proof of Proposition \ref{prop:convolutions1}.
\end{proof}

\begin{proof}[Proof of Theorem \ref{thm:elem1}]
Theorem \ref{thm:elem1} is an immediate consequence of Proposition \ref{prop:convolutions1} and formula \eqref{Vsecond} of Lemma \ref{lem:secondorder} given in the appendix.
\end{proof}

The next results provides asymptotics of convolution powers $f^{\ast(j)}$ for $j=j(t)$ which may grow faster then $t^{1/2}$ under the assumption that the function $|f(t)-at|$ has a finite total variation and satisfies an additional integrability assumption. We shall use a convention that, for a function $x:\mathbb{R}\to\mathbb{R}$, $x^{\ast(0)}(t)=\1_{[0,\infty)}(t)$, $t\in\mathbb{R}$. Also, we shall write $\mathcal{V}_I(x)$ for the total variation of $x$ on the (possibly infinite) interval $I$. Finally, if $x$ is a function of a finite total variation on $[a,\,b]$, $-\infty\leq a<b\leq\infty$ and $y$ is a measurable function on $I$, we stipulate that
$$
\int_{[a,\,b]} y(t)|{\rm d}x(t)|=\int_{[a,\,b]} y(t){\rm d}\left(\mathcal{V}_{[a,\,t]}(x)\right),
$$
where the integral on the right-hand side is understood in the Lebesgue--Stieltjes sense.
\begin{assertion}\label{prop:convolutions2}
Let $f:\mathbb{R}\mapsto [0,\,\infty)$ be a nondecreasing right-continuous function vanishing on the negative half-line. Assume that the function $\varepsilon$ defined by
\begin{equation}\label{eq:asymp_exp_assump2}
\varepsilon(t):=f(t)-at,\quad t\geq 0,
\end{equation}
for some $a>0$, satisfies
\begin{equation}\label{eq:asymp_inegrability1}
\int_{[0,\,\infty)}y|{\rm d}\varepsilon(y)|<\infty.
\end{equation}
Then, for any integer-valued function $j=j(t)$ such that $j(t)=o(t^{2/3})$ as $t\to\infty$,
\begin{equation}\label{eq:asymp_exp_res2}
f_j(t):=f^{\ast(j)}(t)~\sim~\frac{a^j t^j}{j!}\exp{\left(\frac{\gamma_0 j^2}{a t}\right)},\quad t\to\infty,
\end{equation}
where $\gamma_0:=\int_{[0,\,\infty)}{\rm d}\varepsilon(y)=\lim_{t\to\infty}(f(t)-at)$.
\end{assertion}
\begin{proof}
The function $\varepsilon$, as the difference of two nondecreasing functions, has a finite total variation on every finite interval. In particular, \eqref{eq:asymp_inegrability1} entails
$$
\int_{[0,\,\infty)}|{\rm d}\varepsilon(y)|\leq \int_{[0,\,1)}|{\rm d}\varepsilon(y)|+\int_{[1,\,\infty)}y|{\rm d}\varepsilon(y)|<\infty.
$$
Thus, $\varepsilon$ has a finite total variation on $[0,\infty)$. Write
$$
\int_0^{\infty}|\varepsilon(y)-\gamma_0|{\rm d}y=\int_0^{\infty}\left|\int_{(y,\,\infty)}{\rm d}\varepsilon(z)\right|{\rm d}y\leq \int_0^{\infty}\int_{(y,\,\infty)}\left|{\rm d}\varepsilon(z)\right|{\rm d}y=\int_{[0,\,\infty)}y\left|{\rm d}\varepsilon(y)\right|<\infty
$$
having utilized integration by parts for the last equality. Hence, \eqref{eq:asymp_inegrability1} implies that
\begin{equation}\label{eq:asymp_inegrability2}
\int_0^\infty |\varepsilon(y)-\gamma_0|{\rm d}y<\infty.
\end{equation}
Now we modify \eqref{eq:asymp_exp_assump2} in a neighborhood of the origin, so that the
essential properties of $\varepsilon$ given by  \eqref{eq:asymp_inegrability1} and \eqref{eq:asymp_inegrability2} are preserved. Put
\begin{equation}\label{eq:decomposition_modified}
f(t)=(at+\gamma_0)_++\widetilde{\varepsilon}(t)=:\ell(t)+\widetilde{\varepsilon}(t),\quad t\in\mathbb{R}.
\end{equation}
Note that both summands
can be non-zero in a bounded left neighborhood of the origin, yet
\begin{equation}\label{eq:asymp_inegrability3}
\int_{\mathbb{R}}|\widetilde{\varepsilon}(y)|{\rm d}y<\infty\quad\text{and}\quad \int_{\mathbb{R}}|y| |{\rm d}\widetilde{\varepsilon}(y)|<\infty
\end{equation}
because $t\mapsto \varepsilon(t)-\gamma_0-\widetilde{\varepsilon}(t)$  has a bounded support. The advantage of \eqref{eq:decomposition_modified} is justified by a
simple formula for the convolution powers of $\ell$, namely,
$$
\ell^{\ast(j)}(t)=\frac{(at+\gamma_0 j)^j_{+}}{j!},\quad j\in\mathbb{N},~~ t\in\mathbb{R}.
$$
To check this we use the mathematical induction. While, for $j=1$, the formula is trivial,
the induction step works as follows: for $t\geq -a^{-1}\gamma_0(j+1)$,
\begin{align*}
\ell^{\ast(j+1)}(t)=\int_{\mathbb{R}}\frac{(a(t-y)+\gamma_0 j)^j_{+}}{j!}{\rm d}\ell(y)=a\int_{-\gamma_0 a^{-1}}^{t+j\gamma_0 a^{-1}}\frac{(a(t-y)+\gamma_0 j)^j}{j!}{\rm d}y\\
=\int_0^{at+\gamma_0(j+1)}\frac{z^j}{j!}{\rm d}z=\frac{(at+\gamma_0(j+1))^{j+1}}{(j+1)!},
\end{align*}
and $\ell^{\ast(j+1)}(t)=0$ for $t<-a^{-1}\gamma_0(j+1)$.

We intend to prove \eqref{eq:asymp_exp_res2}. Using \eqref{eq:decomposition_modified} we obtain
$$
f^{\ast(j)}(t)=\ell^{\ast(j)}(t)+\sum_{k=0}^{j-1}\binom{j}{k}\left(\ell^{\ast(k)}\ast\widetilde{\varepsilon}^{\ast(j-k)}\right)(t),\quad t\in \mathbb{R}.
$$
We are going to show that the second summand is asymptotically negligible with respect to
$\ell^{\ast(j)}(t)$ whenever $j(t)=o(t^{2/3})$. Assume this has already been done. Then
\eqref{eq:asymp_exp_res2} follows immediately because, for large enough $t$,
$$
f^{\ast(j)}(t)~=~\ell^{\ast(j)}(t)~=~\frac{a^jt^j}{j!}\left(1+\frac{\gamma_0 j}{at}\right)^j~=~\frac{a^jt^j}{j!}\exp{\left(j\log \left(1+\frac{\gamma_0j}{at}\right)\right)}.
$$
The right-hand side is asymptotically equivalent to $\frac{a^jt^j}{j!}\exp{\left(\frac{\gamma_0 j^2}{at}\right)}$ whenever $j=j(t)=o(t^{2/3})$ as $t\to\infty$.

Passing to the analysis of
$$
R_j(t):=\sum_{k=0}^{j-1}\binom{j}{k}\left(\ell^{\ast(k)}\ast\widetilde{\varepsilon}^{\ast(j-k)}\right)(t),\quad t\geq 0
$$
we first check that
\begin{equation}\label{eq:finite_variation}
\mathcal{V}_{\mathbb{R}}(\ell\ast \widetilde{\varepsilon})\leq \widetilde{C}<\infty
\end{equation}
for an absolute constant $C>0$. For $t\in\mathbb{R}$,
$$
(\ell\ast\widetilde{\varepsilon})(t)=\int_{\mathbb{R}}\widetilde{\varepsilon}(t-y){\rm d}\ell(y)=a\int_{-a^{-1}\gamma_0}^\infty \widetilde{\varepsilon}(t-y){\rm d}y=a\int_{-\infty}^{t+a^{-1}\gamma_0}\widetilde{\varepsilon}(y){\rm d}y.
$$
Thus, \eqref{eq:finite_variation} holds with $\widetilde{C}:=a\int_{\mathbb{R}}|\widetilde{\varepsilon}(y)|{\rm d}y$. Put
$$
g_{i,j}(t):=\mathcal{V}_{(-\infty,\,t]}(\ell^{\ast(i)} \ast \widetilde{\varepsilon}^{\ast(j)}),\quad i,j\in\mathbb{N}_0,\quad t\in\mathbb{R}.
$$
Then, for $i,j\in\mathbb{N}$,
\begin{multline*}
g_{i,j}(t)=\mathcal{V}_{(-\infty,\,t]}\left((\ell^{\ast(i-1)} \ast \widetilde{\varepsilon}^{\ast(j-1)})\ast (\ell\ast \widetilde{\varepsilon})\right)\leq \mathcal{V}_{(-\infty,\,t]}(\ell^{\ast(i-1)} \ast \widetilde{\varepsilon}^{\ast(j-1)})\mathcal{V}_{(-\infty,\,t]} (\ell\ast \widetilde{\varepsilon})\\
\leq \mathcal{V}_{(-\infty,\,t]}(\ell^{\ast(i-1)} \ast \widetilde{\varepsilon}^{\ast(j-1)})\mathcal{V}_{\mathbb{R}} (\ell\ast \widetilde{\varepsilon})\leq \widetilde{C}g_{i-1,j-1}(t),\quad t\in\mathbb{R},
\end{multline*}
where we have used that the total variation of the convolution of two functions is bounded by the product of their total variations, see Theorem 1.3.2(c) in \cite{Rudin:1962}. Iterating this inequality we conclude that
$$
|R_j(t)|\leq \sum_{k=0}^{j-1}\binom{j}{k}g_{k,j-k}(t)\leq \sum_{k\leq j/2}\binom{j}{k}\widetilde{C}^k g_{0,j-2k}(t)+\sum_{j/2<k<j}\binom{j}{k}\widetilde{C}^{j-k} g_{2k-j,0}(t),\quad t\in\mathbb{R}.
$$
Note that $g_{0,j-2k}(t)\leq \mathcal{V}_{\mathbb{R}}(\widetilde{\varepsilon}^{\ast(j-2k)})\leq (\mathcal{V}_{\mathbb{R}}(\widetilde{\varepsilon}))^{j-2k}\leq \widetilde{C}_1^{j-2k}$ for $\widetilde{C}_1:=\int_{\mathbb{R}}|{\rm d}\widetilde{\varepsilon}(y)|<\infty$. Therefore,
$$
\sum_{k\leq j/2}\binom{j}{k}\widetilde{C}^k g_{0,j-2k}(t)\leq \sum_{k\leq j/2}\binom{j}{k}\widetilde{C}^k\widetilde{C}^{j-2k}_1\leq (\widetilde{C}\widetilde{C}_1^{-1}+\widetilde{C}_1)^j=o\left(\frac{a^jt^j}{j!}\left(1+\frac{\gamma_0 j}{at}\right)^j\right),\quad t\to\infty,
$$
for $b^j=b^{j(t)}$ grows slower than $\frac{a^jt^j}{j!}\left(1+\frac{\gamma_0 j}{at}\right)^j$ as $t\to\infty$ for an arbitrary finite constant $b>0$. Now we analyze the second sum
\begin{multline*}
\sum_{j/2<k<j}\binom{j}{k}\widetilde{C}^{j-k} g_{2k-j,0}(t)=\sum_{j/2<k<j}\binom{j}{k}\widetilde{C}^{j-k}\mathcal{V}_{(-\infty,\,t]}(\ell^{\ast(2k-j)})=\sum_{j/2<k<j}\binom{j}{k}\widetilde{C}^{j-k}\ell^{\ast(2k-j)}(t)\\
=\sum_{j/2<k<j}\binom{j}{k}\widetilde{C}^{j-k}\frac{(at +\gamma_0(2k-j))^{2k-j}}{(2k-j)!}=\sum_{1\leq k<j/2}\binom{j}{k}\widetilde{C}^{k}\frac{(at +\gamma_0(j-2k))^{j-2k}}{(j-2k)!}.
\end{multline*}
Here, the second equality follows from monotonicity of $\ell$ and the third equality holds for $t$ large enough. It is important for what follows that, for $k<j/2$ and $t>0$,
$$
\frac{(at +\gamma_0(j-2k))^{j-2k}}{(j-2k)!}\leq \frac{a^{j-2k}t^{j-2k}}{(j-2k)!}\exp{\left(\frac{\gamma_0 (j-2k)^2}{at}\right)}.
$$

\noindent
{\sc Case $\gamma_0\geq 0$.} We obtain, for $t>0$,
\begin{align*}
&\hspace{-0.3cm}\sum_{1\leq k<j/2}\binom{j}{k}\widetilde{C}^{k}\frac{(at +\gamma_0(j-2k))^{j-2k}}{(j-2k)!}\leq
\exp{\left(\frac{\gamma_0 j^2}{at}\right)}\sum_{1\leq k<j/2}\binom{j}{k}\widetilde{C}^{k}\frac{a^{j-2k}t^{j-2k}}{(j-2k)!}\\
&=
\frac{a^jt^j}{j!}\exp{\left(\frac{\gamma_0 j^2}{at}\right)}\sum_{1\leq k<j/2}\frac{(j!)^2}{(j-k)!(j-2k)!}\frac{1}{k!}\frac{\widetilde{C}^{k}}{a^{2k}t^{2k}}\leq
\frac{a^jt^j}{j!}\exp{\left(\frac{\gamma_0 j^2}{at}\right)}\sum_{k\geq 1}j^{3k}\frac{1}{k!}\frac{\widetilde{C}^{k}}{a^{2k}t^{2k}}\\
&=\frac{a^jt^j}{j!}\exp{\left(\frac{\gamma_0 j^2}{at}\right)}\left(\exp{\left(\frac{\widetilde{C}j^3}{a^2t^2}\right)}-1\right).
\end{align*}
The last factor converges to zero whenever $j=j(t)=o(t^{2/3})$, whence the claim.

\noindent
{\sc Case $\gamma_0<0$.} Arguing in the same vein it is enough to check that
$$
\sum_{1\leq k<j/2}\frac{1}{k!}\left(\frac{\widetilde{C}j^3}{a^2t^2}\right)^k\exp{\left(\frac{\gamma_0(j-2k)^2}{at}\right)}=o\left(\exp{\left(\frac{\gamma_0j^2}{at}\right)}\right),\quad t\to\infty
$$
which is equivalent to
$$
I_t:=\sum_{1\leq k<j/2}\frac{1}{k!}\left(\frac{\widetilde{C}j^3}{a^2t^2}\right)^k\exp{\left(\frac{4|\gamma_0|k(j-k)}{at}\right)}=o(1),\quad t\to\infty.
$$
Invoking the inequality $\exp{\left(\frac{4|\gamma_0|k(j-k)}{at}\right)}\leq \exp{(4|\gamma_0|a^{-1}k)}$ for $1\leq k<j$ and large enough $t$ we infer $$I_t\leq \sum_{1\leq k<j/2}\frac{1}{k!}\left(\frac{\widetilde{C}j^3\exp(4|\gamma_0|a^{-1})}{a^2t^2}\right)^k\leq\exp\left(\frac{\widetilde{C}j^3}{a^2t^2}\exp{(4|\gamma_0|a^{-1})}\right)-1~\to~0,\quad t\to\infty.$$
The proof of Proposition \ref{prop:convolutions2} is complete.
\end{proof}

\begin{proof}[Proof of Theorem \ref{thm:elem2}]
We intend to apply Proposition \ref{prop:convolutions2}. To this end, it is enough to check that, under the assumptions of Theorem \ref{thm:elem2},
$$\int_{[0,\,\infty)}y|{\rm d}(V(y)-{\tt m}^{-1}y)|<\infty.$$
Recall that $V=U\ast G$ and denote by ${\rm Id}$ the identity function on $[0,\infty)$, that is, ${\rm Id}(t):=t_{+}=t\1_{[0,\infty)}(t)$ for $t\in\mathbb{R}$. Then
$$
V-{\tt m}^{-1}{\rm Id}=(U-{\tt m}^{-1}{\rm Id})\ast G-{\tt m}^{-1}({\rm Id}\ast (1-G)).
$$
Using this and integration by parts yields
\begin{align*}
&\int_{[0,\,\infty)}y|{\rm d}(V(y)-{\tt m}^{-1}y)|=-\int_{[0,\,\infty)}y{\rm d}\mathcal{V}_{[y,\,\infty)}(V-{\tt m}^{-1}{\rm Id})=\int_{[0,\,\infty)}\mathcal{V}_{[y,\,\infty)}(V-{\tt m}^{-1}{\rm Id}){\rm d}y\\
&\hspace{2cm}\leq \int_{[0,\,\infty)}\mathcal{V}_{[y,\,\infty)}(U-{\tt m}^{-1}{\rm Id}){\rm d}y+{\tt m}^{-1}\int_{[0,\,\infty)}\mathcal{V}_{[y,\,\infty)}({\rm Id}\ast (1-G)){\rm d}y\\
&\hspace{2cm}=\int_{[0,\,\infty)}y|{\rm d}(U(y)-{\tt m}^{-1}y)|+{\tt m}^{-1}\int_0^{\infty}\int_y^{\infty} (1-G(z)){\rm d}z{\rm d}y.
\end{align*}
The first summand is finite by Remark 3.1.7(ii) on p.~121 in \cite{Frenk:1982} and the second is finite in view of
the assumption $\mathbb{E}\eta^2<\infty$.

The explicit form of $\gamma_0$ follows from the decomposition
$$
V(t)-{\tt m}^{-1}t=\int_{[0,\,t]}(U(t-y)-{\tt m}^{-1}(t-y)){\rm d}G(y)-{\tt m}^{-1}\int_{[0,\,t]}y{\rm d}G(y)-{\tt m}^{-1}t(1-G(t)),
$$
in which the first summand converges to $(2{\tt m}^2)^{-1}\mathbb{E}\xi^2$ by the dominated convergence theorem, \eqref{lord} and \eqref{eq:sec_order_U}; the second converges to $-{\tt m}^{-1}\mathbb{E}\eta$ and the third tends to zero as $t\to\infty$.
\end{proof}

Finally, we give a general result on the behavior of $f^{\ast(j)}$ for arbitrary $j=j(t)=o(t)$. Unfortunately, this result can seldom be applied to the counting function $V$ but is of independent interest and has at least two merits. On the one hand, it gives a probabilistic explanation of a rather mysterious appearance of the exponent in \eqref{eq:asymp_exp_res2}. On the other hand, it may be used for guessing the behaviour of $V_j$ for $j=j(t)$ growing at least as fast as $t^{2/3}$.
\begin{assertion}\label{prop:convolutions3}
Let $(\widetilde{S}_j)_{j\in\mathbb{N}_0}$ be a nondecreasing zero-delayed standard random walk with $K(t):=\mathbb{P}\{\widetilde{S}_1\leq t\}$ for $t\in\mathbb{R}$. Assume that, for some $a>0$,
$$
f(t)=at-\int_0^t(1-K(y)){\rm d}y,\quad t\geq 0.
$$
Then
$$
f^{\ast(j)}(t)=\frac{\mathbb{E} (at-\widetilde{S}_j)_{+}^j}{j!},\quad j\in\mathbb{N},~~ t\geq 0.
$$
In particular, if $\mathbb{E} \widetilde{S}_1^2<\infty$ and $j=j(t)=o(t^{2/3})$ as $t\to\infty$, then \eqref{eq:asymp_exp_res2} holds with $\gamma_0=-\mathbb{E} \widetilde{S}_1$.
\end{assertion}
\begin{proof}
Replacing $K$ with $t\mapsto K(at)$ we can and do assume that $a=1$, that is, $f(t)=\int_0^t K(y){\rm d}y$ or, in short, $f=K\ast {\rm Id}$. Then
$$
f^{\ast(j)}(t)=\left(({\rm Id})^{\ast(j)}\ast K^{\ast(j)}\right)(t)=\int_{[0,\,t]}\frac{(t-y)^j}{j!}{\rm d}K^{\ast(j)}(y)=\frac{\mathbb{E} (t-\widetilde{S}_j)_{+}^j}{j!},\quad t\geq 0.
$$
If $\mathbb{E} \widetilde{S}_1^2<\infty$, then $j=j(t)=o(t^{2/3})$ as $t\to\infty$ implies that
\begin{equation}\label{eq:derivation_rw_asymptotics}
\mathbb{E} (t-\widetilde{S}_j)_{+}^j~\sim~t^j\exp\left(\frac{\gamma_0j^2}{t}\right),\quad t\to\infty.
\end{equation}
This can be justified as follows. We first note that $\gamma_0<0$. Further, in the decomposition
\begin{equation}\label{eq:decomposition_gen_conv}
\mathbb{E}\left(1-\frac{\widetilde{S}_j}{t}\right)_{+}^j=\mathbb{E} (e^{j\log(1-\widetilde{S}_j/t)}\1_{\{\widetilde{S}_j\leq t/2\}})+\mathbb{E}\left(1-\frac{\widetilde{S}_j}{t}\right)_{+}^j\1_{\{\widetilde{S}_j\in (t/2,\,t)\}}
\end{equation}
the second summand is bounded by $2^{-j}$ and $2^{-j}=o\left(\exp\left(\frac{\gamma_0j^2}{t}\right)\right)$ as $t\to\infty$, for $j^2/t=o(j)$. The first summand in \eqref{eq:decomposition_gen_conv} can be bounded with the help of the inequalities
$$
-x-x^2 \leq \log(1-x)\leq -x,\quad x\in [0,\,1/2]\quad\text{and}\quad 1-x\leq e^{-x},\quad x\in\mathbb{R}.
$$
Indeed, we obtain, for $j\geq 4$,
\begin{multline*}
\mathbb{E} e^{-j\widetilde{S}_j/t}\left(1-\frac{j\widetilde{S}_j^2}{t^2}\right)\leq \mathbb{E} e^{-j\widetilde{S}_j/t}\left(1-\frac{j\widetilde{S}_j^2}{t^2}\right)\1_{\{\widetilde{S}_j\leq t/2\}}\leq \mathbb{E} e^{-j\widetilde{S}_j/t}e^{-j\widetilde{S}_j^2/t^2}\1_{\{\widetilde{S}_j\leq t/2\}}\\
\leq \mathbb{E} (e^{j\log(1-\widetilde{S}_j/t)}\1_{\{\widetilde{S}_j\leq t/2\}})\leq \mathbb{E} e^{-j\widetilde{S}_j/t}\1_{\{\widetilde{S}_j\leq t/2\}}\leq \mathbb{E} e^{-j\widetilde{S}_j/t}.
\end{multline*}
For $\lambda\geq 0$, put $\phi(\lambda):=\mathbb{E} e^{-\lambda\widetilde{S}_1}$. In view of $\mathbb{E}\widetilde{S}_1^2<\infty$ we infer
$$
\mathbb{E} e^{-j\widetilde{S}_j/t}=\phi^j(j/t)=\left(1+\frac{\gamma_0 j}{t}+O\left(\frac{j^2}{t^2}\right)\right)^j.
$$
The right-hand side is asymptotically equivalent to $\exp(\gamma_0 j^2/t)$ as $t\to\infty$ under the assumption $j=j(t)=o(t^{2/3})$. Finally, the relation
$$
\mathbb{E} e^{-j\widetilde{S}_j/t}\left(\frac{j\widetilde{S}_j^2}{t^2}\right)=o\left(\mathbb{E} e^{-j\widetilde{S}_j/t}\right),\quad t\to\infty
$$
can be checked using the equality $\mathbb{E} e^{-j\widetilde{S}_j/t}\widetilde{S}_j^2=\frac{\partial^2}{\partial \lambda^2}(\phi^j (\lambda))\Big|_{\lambda=j/t}$ in conjunction with the assumptions $j=j(t)=o(t^{2/3})$ and $\mathbb{E}\widetilde{S}_1^2<\infty$.
\end{proof}

\begin{rem}\label{rem:remove}
In the setting of Proposition \ref{prop:convolutions3}, assume that $\mathbb{E} \widetilde{S}_1^3<\infty$ and $j=j(t)=o(t^{3/4})$ as $t\to\infty$. We state, without going into details (which become rather technical), that
$$\mathbb{E}\big(at-\widetilde{S}_j\big)_+^j~\sim~a^jt^j\exp{\left(\gamma_0 j^2/t+(\gamma_1/2-\gamma_0^2)j^3/t^2\right)},\quad t\to\infty,
$$
where $\gamma_0=-\mathbb{E} \widetilde{S}_1$ and $\gamma_1:=\mathbb{E} \widetilde{S}_1^2$.
\end{rem}

\subsection{Proof of Theorem \ref{thm:keyren} }

For $t\geq 0$, put $g(t):=\int_{[0,\,t]}f(t-y){\rm d}V(y)$ and $I:={\tt m}^{-1}\int_0^\infty f(y){\rm d}y$. By Lemma \ref{key}(a), given $\varepsilon>0$ there exists $t_0>0$ such that $|g(t)-I|\leq \varepsilon$ whenever $t\geq t_0$. Also, by Lemma \ref{keylight}, $g(t)\leq J$ for some $J>0$ and all $t\geq 0$. Hence, for $t\geq t_0$,
\begin{multline}\label{eq:keyren_eq1}
(f\ast V_j)(t)=(g\ast V_{j-1})(t)=\int_{[0,\,t]}g(t-y){\rm d}V_{j-1}(y)=\int_{[0,\,t-t_0]}g(t-y){\rm d}V_{j-1}(y)\\
+\int_{(t-t_0,\,t]}g(t-y){\rm d}V_{j-1}(y) \leq (I+\varepsilon)V_{j-1}(t)+J (V_{j-1}(t)-V_{j-1}(t-t_0)).
\end{multline}
We claim that
\begin{equation}\label{eq:keyren_increment}
\lim_{t\to\infty}\frac{V_{j(t)-1}(t)-V_{j(t)-1}(t-t_0)}{V_{j(t)-1}(t)}=0.
\end{equation}
Note that \eqref{eq:keyren_increment} is not a direct consequence of the elementary renewal theorem, for the theorem provides the asymptotics of $V_{j(t-t_0)-1}(t-t_0)$ rather than $V_{j(t)-1}(t-t_0)$ which is actually needed for \eqref{eq:keyren_increment}. To prove \eqref{eq:keyren_increment} we write with the help of \eqref{subad}
$$
0\leq V_{j(t)-1}(t)-V_{j(t)-1}(t-t_0)=\int_{[0,\,t]}(V(t-y)-V(t-t_0-y)){\rm d}V_{j(t)-2}(y)\leq U(t_0)V_{j(t)-2}(t),
$$
for all $t\geq 0$. Thus, \eqref{eq:keyren_increment} follows from
$$
\lim_{t\to\infty}\frac{V_{j(t)-2}(t)}{V_{j(t)-1}(t)}=0,
$$
which is a consequence of Theorems \ref{thm:elem1} and \ref{thm:elem2} applied with $j=j(t)-1$ and $j=j(t)-2$.

Combining \eqref{eq:keyren_eq1} and \eqref{eq:keyren_increment} we obtain
$$
\limsup_{t\to\infty}\frac{(f\ast V_j)(t)}{V_{j-1}(t)}\leq I.
$$
The converse inequality for the limit inferior follows analogously. The remaining statements of Theorem \ref{thm:keyren} are secured by Theorems \ref{thm:elem1} and \ref{thm:elem2}.

\section{Appendix}

In this section we shall prove counterparts for perturbed random walks of some standard renewal-theoretic results. Recall that, under the sole assumption ${\tt m}=\mathbb{E}\xi<\infty$,
$$\lim_{t\to\infty}\frac{U(t)}{t}=\lim_{t\to\infty}\frac{V(t)}{t}=\frac{1}{{\tt m}}.$$ We start by discussing the rate of convergence in both limit relations.
\begin{lemma}\label{lem:secondorder}
Assume that either (i) $\mathbb{E}\xi^r<\infty$ for some $r\in (1,\,2]$, or (ii) $\mathbb{P}\{\xi>t\}~\sim~bt^{-r}$ for some $r\in (1,\,2)$ and some $b>0$ as $t\to\infty$. Then
\begin{equation}\label{eq:u_two_terms_exp}
U(t)=\frac{t}{{\tt m}}+O(t^{2-r}),\quad t\to\infty,
\end{equation}
where ${\tt m}=\mathbb{E}\xi<\infty$. In the case (i) when $r\in (1,2)$, the big $O$ can be replaced with a little $o$.

Under the additional assumption $\mathbb{E}(\eta\wedge t)=O(t^{2-r})$ as $t\to\infty$,
\begin{equation}\label{Vsecond}
V(t)=\frac{t}{{\tt m}}+O(t^{2-r}),\quad t\to\infty.
\end{equation}
\end{lemma}
\begin{proof}
First, we focus on \eqref{eq:u_two_terms_exp}.

\noindent
{\sc Case (i)}. If $r=2$, then \eqref{eq:u_two_terms_exp} follows from Lorden's inequality \eqref{lord}. Assume now that $r\in (1,2)$ and note that, for any $p>1$, $\mathbb{E} \xi^p<\infty$ is equivalent to $\mathbb{E} (S_0^\ast)^{p-1}<\infty$. The situation is not excluded that $\mathbb{E} S_0^\ast<\infty$ in which case $\mathbb{E}\xi^2<\infty$, so that $U(t)={\tt m}^{-1}t +O(1)={\tt m}^{-1}t+O(t^{2-r})$ as $t\to\infty$. Thus, in what follows we can and do assume that $\mathbb{E} S_0^\ast=\infty$. Then $$U(t)-{\tt m}^{-1}t=\int_{[0,\,t]}\mathbb{P}\{S_0^\ast>t-y\}{\rm d}U(y)~\sim~{\tt m}^{-1}\int_0^t \mathbb{P}\{S_0^\ast>y\}{\rm d}y,\quad t\to\infty,$$ where the equality is nothing else but \eqref{eq:S0}, and the asymptotic relation follows from Theorem 4 in \cite{Sgibnev:1981}. Now $\mathbb{E} (S_0^\ast)^{r-1}<\infty$ entails $\mathbb{P}\{S_0^\ast>t\}=o(t^{1-r})$, whence $\int_0^t \mathbb{P}\{S_0^\ast>y\}{\rm d}y=o(t^{2-r})$ as $t\to\infty$.

\noindent {\sc Case (ii)}. In this case $\mathbb{P}\{S_0^\ast>t\}\sim b({\tt m}(r-1))^{-1} t^{-(r-1)}$ as $t\to\infty$. This implies that
\begin{equation*}
U(t)-\frac{t}{\tt m}~\sim~\frac{1}{{\tt m}}\int_0^t \mathbb{P}\{S_0^\ast>y\}{\rm d}y~\sim~\frac{b}{{\tt m^2}(r-1)(2-r)}t^{2-r},\quad t\to\infty
\end{equation*}
and thereupon \eqref{eq:u_two_terms_exp}. Under the additional assumption that the distribution of $\xi$ is nonlattice relation \eqref{eq:u_two_terms_exp} also follows from Theorem 2.2 in \cite{Mohan:1976}.

Finally, relation \eqref{Vsecond} follows from the equality (which has already appeared in Section \ref{Prep})
$$V(t)-{\tt m}^{-1}t=\int_{[0,\,t]}(U(t-y)-{\tt m}^{-1}(t-y)){\rm d}G(y)-{\tt m}^{-1}\mathbb{E} (\eta\wedge t)$$
because each summand is $O(t^{2-r})$ by \eqref{eq:u_two_terms_exp} and the assumption of the theorem, respectively.
\end{proof}

We continue by noting that
\begin{equation}\label{subad}
V(x+y)-V(x)\leq U(y),\quad x,y\in\mathbb{R}.
\end{equation}
Indeed, for $x,y\geq 0$,
\begin{eqnarray}\label{eqV}
V(x+y)-V(x)&=&\mathbb{E} (U(x+y-\eta)-U(x-\eta))\1_{\{\eta\leq x\}}+\mathbb{E} U(x+y-\eta)\1_{\{x<\eta\leq x+y\}}\\&\leq& U(y)(\mathbb{P}\{\eta\leq x\}+\mathbb{P}\{x<\eta\leq x+y\})\leq U(y)\notag
\end{eqnarray}
having utilized subadditivity and monotonicity of $U$ for the penultimate inequality. If $x,y<0$, then both sides of \eqref{subad} are zero. Finally, we use monotonicity of $V$ to obtain: if $x<0$ and $y\geq 0$, then $V(x+y)-V(x)=V(x+y)\leq V(y)\leq U(y)$; and if $x\geq 0$ and $y<0$, then $V(x+y)-V(x)\leq 0=U(y)$.

Lemmas \ref{bla} and \ref{key} are counterparts of Blackwell's theorem and the key renewal theorem, respectively. Observe that the presence of the $\eta_k$ plays no role, and the results are of the same form as for renewal functions.

\begin{lemma}\label{bla}
Let $h>0$ be any fixed number.

\noindent (a) Assume that the distribution of $\xi$ is nonlattice and ${\tt m}=\mathbb{E}\xi<\infty$. Then $$\lim_{t\to\infty} (V(t+h)-V(t))={\tt m}^{-1}h.$$

\noindent (b) Assume that ${\tt m}=\infty$ (the assumption that the distribution of $\xi$ is nonlattice is not needed). Then
\begin{equation}\label{bla0}
\lim_{t\to\infty} (V(t+h)-V(t))=0.
\end{equation}
\end{lemma}
\begin{proof}
(a) According to Blackwell's theorem,
\begin{equation}\label{blackwell}
\lim_{t\to\infty}(U(t+h)-U(t))={\tt m}^{-1}h.
\end{equation}
In view of \eqref{blackwell}, $\lim_{t\to\infty}(U(t+h-\eta)-U(t-\eta))\1_{\{\eta\leq t-t^{1/2}\}}={\tt m}^{-1}h$ a.s. Recalling \eqref{subad} we infer
$$
\lim_{t\to\infty}\mathbb{E} (U(t+h-\eta)-U(t-\eta))\1_{\{\eta\leq t-t^{1/2}\}}={\tt m}^{-1}h
$$
by Lebesgue's dominated convergence theorem. Another appeal to \eqref{subad} yields $$\mathbb{E} (U(t+h-\eta)-U(t-\eta))\1_{\{t-t^{1/2}<\eta\leq t\}}\leq U(h)\mathbb{P}\{t-t^{1/2}<\eta\leq t\},$$ and the right-hand side converges to $0$ as $t\to\infty$. Finally, by monotonicity,
$$\mathbb{E} U(t+h-\eta)\1_{\{t<\eta\leq t+h\}}\leq U(h)\mathbb{P}\{t<\eta\leq t+h\},$$ and the right-hand side converges to $0$ as $t\to\infty$. Invoking the first equality in \eqref{eqV} with $x=t$ and $y=h$ completes the proof of part (a).

\noindent (b) If the distribution of $\xi$ is nonlattice, then, by Blackwell's theorem,
\begin{equation}\label{blackwell1}
\lim_{t\to\infty}(U(t+h)-U(t))=0.
\end{equation}
If the distribution of $\xi$ is $d$-lattice, then, by Blackwell's theorem, \eqref{blackwell1} holds for $h=jd$, $j\in\mathbb{N}$. However, using monotonicity of $U$ we can ensure that \eqref{blackwell1} holds for any fixed $h>0$ in both nonlattice and lattice cases. With this at hand, repeating verbatim the proof of part (a) we arrive at \eqref{bla0}.
\end{proof}

\begin{lemma}\label{key}
Let $f: \mathbb{R}\to \mathbb{R}$ be a directly Riemann integrable (dRi) function on $\mathbb{R}$.

\noindent (a) Assume that ${\tt m}<\infty$ and that the distribution of $\xi$ is nonlattice. Then
$$\lim_{t\to\infty} \int_{[0,\,\infty)} f(t-y){\rm d}V(y)= {\tt m}^{-1} \int_\mathbb{R} f(y){\rm d}y.$$

\noindent (b) Assume that ${\tt m}=\infty$ (the assumption that the distribution of $\xi$ is nonlattice is not needed). Then
$$\lim_{t\to\infty} \int_{[0,\,\infty)} f(t-y){\rm d}V(y)=0.$$

If $f$ is dRi on $[0,\infty)$ or $(-\infty, 0]$, then the ranges of integration $[0,\,\infty)$ and $\mathbb{R}$ should be replaced with $[0,\,t]$ and $[0,\infty)$ or $[t,\infty)$ and $(-\infty, 0]$, respectively.
\end{lemma}
\begin{proof}
(a) We only prove the claim under the assumption that $f$ is dRi on $\mathbb{R}$ which is equivalent to the fact that $f_{+}$ and $f_{-}$ (nonnegative and nonpositive parts of $f$) are dRi on $\mathbb{R}$. Thus, we can and do assume that $f\geq 0$ on $\mathbb{R}$. Obviously, it is enough to show that $$\lim_{t\to\infty} \int_{[0,\,t]} f(t-y){\rm d}V(y)= {\tt m}^{-1} \int_0^\infty f(y){\rm d}y$$ and that $$\lim_{t\to\infty} \int_{(t,\,\infty)} f(t-y){\rm d}V(y)= {\tt m}^{-1} \int_{-\infty}^0 f(y){\rm d}y.$$ The proof of the first relation with $U$ replacing $V$ can be found on p.~241--242 in \cite{Resnick:2002}. We only check the second limit relation by following closely aforementioned Resnick's proof.

We proceed via three steps complicating successively the structure of $f$.

\noindent {\sc Step 1}. Suppose first that
$$f(t)=\1_{[(n-1)h,\,nh)}(t),\quad t<0$$ for fixed nonpositive integer $n$ and $h>0$. Then $f(t-y)=1$ if, and only if, $y\in
(t-nh,\,t-(n-1)h]$ which entails $$\int_{(t,\,\infty)}f(t-y){\rm d}V(y)=V(t-(n-1)h)-V(t-nh).$$ By Lemma \ref{bla}(a), the last
difference tends to ${\tt m}^{-1}h$ as $t\to\infty$, thereby proving that $$\lim_{t\to\infty} \int_{(t,\,\infty)}f(t-y){\rm d}V(y)={\tt m}^{-1}h={\tt m}^{-1}\int_{-\infty}^0 f(y){\rm d}y.$$

\noindent {\sc Step 2}. Suppose now that $$f(t)=\sum_{n\leq 0} c_n\1_{[(n-1)h,\,nh)}(t),\quad t<0,$$ where
$\big(c_n\big)_{n\leq 0}$ is a sequence of nonnegative numbers satisfying $\sum_{n\leq 0}c_n<\infty$. An argument similar to that used in the previous step enables us to assert that $$\int_{(t,\,\infty)}f(t-y){\rm d}V(y)=\sum_{n\leq 0}c_n\big(V(t-(n-1)h)-V(t-nh)\big).$$ Using Lemma \ref{bla}(a) in combination with \eqref{subad} we infer with the help of Lebesgue's dominated convergence theorem, $$\lim_{t\to\infty} \int_{(t,\,\infty)} f(t-y){\rm d}V(y)={\tt m}^{-1}h \sum_{n\leq 0} c_n={\tt m}^{-1}\int_{-\infty}^ 0 f(y){\rm d}y.$$

\noindent {\sc Step 3}. Let now $f$ be an arbitrary nonnegative dRi function on $\mathbb{R}$ (actually, for the present proof it is enough it is dRi on $(-\infty, 0)$). For each $h>0$, put $$\overline{f}_h(t):=\sum_{n\leq 0}\underset{(n-1)h\leq
y<nh}{\sup}\,f(y)\1_{[(n-1)h,\,nh)}(t),\quad t<0$$ and $$\underline{f}_h(t):=\sum_{n\leq 0}\underset{(n-1)h\leq y<nh}{\inf}\,f(y)\1_{[(n-1)h,\,nh)}(t), \ \ t<0.$$ By the definition of direct Riemann integrability, $$\sum_{n\leq 0}\underset{(n-1)h\leq y<nh}{\sup}\,f(y)<\infty \quad \text{and}\quad \sum_{n\leq 0}\underset{(n-1)h\leq y<nh}{\inf}\,f(y)<\infty$$ for each $h>0$. Thus, the functions $\overline{f}_h$ and $\overline{f}_h$ have the same structure as the functions discussed in Step 2. According to the result of Step 2,
$$\lim_{t\to\infty} \int_{(t,\,\infty)}\overline{f}_h(t-y){\rm d}V(y)={\tt m}^{-1}h\sum_{n\leq 0}\underset{(n-1)h\leq y<nh}{\sup}\,f(y)=:{\tt m}^{-1}\overline{\sigma}(h)$$ and $$\lim_{t\to\infty} \int_{(t,\,\infty)}\underline{f}_h(t-y){\rm d}V(y)={\tt m}^{-1}h\sum_{n\leq 0}\underset{(n-1)h\leq y<nh}{\inf}\,f(y)=:\mu^{-1}\underline{\sigma}(h)$$ for all $h>0$.
Since, for each $h>0$, $$\underline{f}_h(t)\leq f(t)\leq \overline{f}_h(t),\quad t<0,$$ it follows that
\begin{eqnarray*}
{\tt m}^{-1}\underline{\sigma}(h)&=&\underset{t\to\infty}{\lim\inf}\,\int_{(t,\,\infty)}\underline{f}_h(t-y){\rm d}V(y)\leq \underset{t\to\infty}{\lim\inf}\,\int_{(t,\,\infty)} f(t-y){\rm
d}V(y)\\&\leq& \underset{t\to\infty}{\lim\sup}\,\int_{(t,\,\infty)}f(t-y){\rm d}V(y)\leq \underset{t\to\infty}{\lim\sup}\,\int_{(t,\,\infty)}\overline{f}_h(t-y){\rm d}V(y)\\&=&{\tt m}^{-1}\overline{\sigma}(h).
\end{eqnarray*}
We have $\lim_{h\to 0+}\,\big(\overline{\sigma}(h)-\underline{\sigma}(h)\big)=0$ by the definition of direct Riemann integrability. Also, it is known that
$\lim_{h\to 0+}\, \overline{\sigma}(h)=\int_{-\infty}^0 f(y){\rm d}y$. Letting $h\to 0+$ in the last chain of inequalities completes the proof of part (a).

\noindent (b) Use part (b) of Lemma \ref{bla} in place of part (a) and proceed as above.
\end{proof}

Sometimes it is the case that the precision of Lemma \ref{key} is not needed. In this situation the following `light' version, borrowed from Lemma 9.1 in \cite{Iksanov+Marynych+Samoilenko:2020}, may suffice.
\begin{lemma}\label{keylight}
Let $f: \mathbb{R}\to [0,\infty)$ be a dRi function on $\mathbb{R}$. Then for some $r>0$ and all $x\in\mathbb{R}$
\begin{equation}\label{bound}
\int_{[0,\,\infty)} f(x-y){\rm d}V(y)\leq r.
\end{equation}

If $f$ is dRi on $[0,\infty)$ or $(-\infty, 0]$, then the range of integration $[0,\,\infty)$ should be replaced with $[0,\,x]$ or $[x,\infty)$ and then \eqref{bound} holds for all $x\geq 0$ or all $x\leq 0$, respectively.
\end{lemma}

\vspace{5mm}

\noindent {\bf Acknowledgement}. The present work was supported by National Research Foundation of Ukraine (project 2020.02/0014 ``Asymptotic regimes of perturbed random walks: on the edge of modern and classical probability'').


\section{Introduction}
\label{sec:I}

SU($N$) Mott insulators are representative examples of quantum systems
with a large local Hilbert space where quantum fluctuations can be
strongly enhanced and exotic quantum phases could be stabilized.
All the SU($N$) spin operators are present in the effective model
for the Mott insulators and would be able to shuffle all the spin states
rather actively in the local Hilbert space. The system can be quite
delocalized within the local Hilbert space, that is to enhance
quantum fluctuations and induce exotic quantum phases. This aspect
is fundamentally different from the SU(2) Mott insulators with
large-$S$ local moments that also has a large local Hilbert space.
For the large-$S$ SU(2) Mott insulators, the pairwise Heisenberg interaction
is quite ineffective to delocalize the spin states in the large-$S$
Hilbert space, and thus quantum fluctuations are strongly suppressed.
Thus, it is the conventional wisdom not to search for exotic quantum phases
among the large-$S$ SU(2) Mott insulators, but among the spin-1/2 quantum
magnets with a strong frustration. In contrast, the emergence of the SU($N$)
Mott insulators brings a new searching direction for exotic quantum phases.

SU($N$) Mott insulators are not a theoretical fantasy, but exist in nature.
It has been shown that the ultracold alkaline-earth atoms (AEA) on optical lattices
can simulate quantum many-body physics with an SU($N$) symmetry
without any fine-tuning~\cite{Gorshkov2010}.
The nuclear spin of fermionic AEA can be as large as ${I=9/2}$ for \ce{^87Sr}, while the outer shell
electrons give a total spin ${S=0}$ and makes the hyperfine coupling inactive.
This observation effectively extends the realization of SU($N$)
magnets up to ${N=2I+1}$.
Early efforts by Congjun Wu explored total spin-3/2 alkaline fermions on optical
lattices where the SU(4) symmetry can be achieved with fine
tuning~\cite{PhysRevLett.91.186402,PhysRevB.77.134449,PhysRevLett.95.266404,WU_2006}.
Quantum Monte Carlo simulations were introduced later to study magnetic properties of the SU$(2N)$ Hubbard model~\cite{PhysRevB.88.125108,PhysRevLett.112.156403}.
There is also some effort in searching for an emergent SU($N$)
symmetry in real materials particularly for the SU(4) case.
The two-orbital Kugel-Khomskii model can become SU(4) symmetric
after some fine tuning~\cite{Kugel__1982,PhysRevB.60.6584}.
Experimental and numerical evidence also suggests that \ce{Ba3CuSb2O9} could be a prominent candidate on a decorated honeycomb lattice~\cite{PhysRevB.100.205131}, though the Cu-Sb dumbbell is quenched rather than an active degree of freedom.
The SU(4) Heisenberg model has further been proposed for the spin-orbit-coupled
Mott insulator $\alpha$-\ce{ZrCl3} where the degree of freedom is the spin-orbit-entangled
${J=3/2}$ local moment~\cite{PhysRevB.82.174440} on a honeycomb lattice~\cite{PhysRevLett.121.097201}.
\begin{figure}[b]
    \centering
    \includegraphics[width=0.48\textwidth]{Fig01.pdf}
    \caption{(a) Triangular lattice with background U(1) gauge flux $\phi$ within each plaquette.
    (b) The spinon dispersions for chiral spin liquids with $\phi=5\pi/6$.}
    \label{fig:01}
\end{figure}
More recently, it has been proposed that the Mott insulating and
superconducting behaviors in twisted bilayer graphenes can be captured by a two-orbital
Hubbard model with an emergent SU(4) symmetry on a moir\'{e} triangular
lattice~\cite{PhysRevLett.121.087001}. Other work even suggested that double moir\'{e} layers,
built from transition metal dichalcogenides or graphenes, separated from one another
by a thin insulating layer would be a natural platform to realize Hubbard models on the
triangular lattice with SU(4) or SU(8) symmetries~\cite{Ashvin:2020,zhang2021su4}.

Owing to enhanced quantum fluctuations for SU($N$) Mott insulators,
the pioneering theoretical works~\cite{PhysRevLett.103.135301,PhysRevB.84.174441}
by Hermele \textit{et al}, have obtained
the topological chiral spin liquid (CSL) ground states with
intrinsic topological orders even for an unfrustrated square lattice
when ${N \geq 5}$~\cite{PhysRevLett.103.135301,PhysRevB.84.174441,PhysRevA.93.061601}.
Depending on the atom occupation numbers, the system can support both Abelian
and non-Abelian statistics for the anyonic excitations. Since unfrustrated
lattices such as square and honeycomb lattices already bring exotic and interesting
physics~\cite{PhysRevLett.107.215301,PhysRevLett.105.265301,PhysRevX.2.041013},
the frustrated lattices could further harbor nontrivial quantum phases, for example the triangular lattice with SU(3) and SU(4) spins \cite{10.21468/SciPostPhys.8.5.076,PhysRevB.101.245159,PhysRevResearch.2.023098}.
In this work, we focus on the SU($N$) Heisenberg model
on a triangular lattice where each lattice site comprises the fundamental
representation of the SU($N$) group. For AEA,
this corresponds to one atom per lattice site and is known to be most stable
against the three-body loss. Because the $1/N$ filling is kept throughout,
the large-$N$ limit differs fundamentally from the large-$N$ extension
of the spin-1/2 SU(2) models where the 1/2 is kept and two sites can form
a SU(2) spin singlet~\cite{Subir}. Here, as $N$ sites are needed to form a singlet, the
valence cluster solid (VCS) state is generically disfavored in the large-$N$ limit.
Instead, by our large-$N$ calculation, two types of CSL states with background U(1) gauge flux $\phi = \pi-\pi/N$ and $\pi-\pi/(2N)$ (see Fig.~\ref{fig:01}(a)) are identified as the ground and lowest competing states for $6 \le N \le 9$, respectively.
For smaller $N$'s, various symmery-broken cluster/stripe states are obtained.
We expect the large-$N$ results are more reliable
when $N$ is large.

\begin{table}[t]
    \caption{
    The SCM results of the ground and
    lowest competing states for ${2 \le N \le 9}$. The average background fluxes of doubled unit-cell stripe states for $N=4$ and $5$ satisfy $\phi_{\text{avg}}=\pi-\pi/N$.}
    \label{tab:A}
    \begin{ruledtabular}
        \begin{tabular}{ccc}
            $N$  & Ground state & Lowest competing state \\
            \midrule
            2 & dimer state & CSL $\phi = \pi/2$ \\
            3 & three-site complex VCS & CSL $\phi = 2\pi/3$ \\
            4 & four-site VCS &  doubled unit-cell stripe\\
            5 & doubled unit-cell stripe & CSL  $\phi = 4\pi/5$\\

            $6\le N \le 9$ & CSL $\phi = \pi-\pi/N$ & CSL $\phi = \pi-\pi/(2N)$
        \end{tabular}
        \end{ruledtabular}
\end{table}

The rest of this paper is organized as follows.
The general Heisenberg model of SU($N$) spins is introduced and simplified at the large-$N$ saddle point in Sec.~\ref{sec:II}.
The self-consistent minimization algorithm is implemented to solve the reduced mean-field spinon Hamiltonian.
The technique details of this algorithm are described in Sec.~\ref{sec:III}.
It is emphasized that the optimized solutions strictly satisfy the local constraints.
In Sec.~\ref{sec:IV}, both the ground states and lowest competing states for $2 \le N \le 9$ are reported and analysed, especially for two types of CSL states.
Finally, in Sec.~\ref{sec:V}, this work conclude with a discussion about relevant systems and the experimental probes.

\section{The SU($N$) Heisenberg model in the large-$N$ approximation}
\label{sec:II}

We begin with an SU($N$) Heisenberg model
on the triangular lattice where each site comprises
the fundamental representation of the SU($N$) group.
This model can be obtained from the strong coupling limit of an
SU($N$) Hubbard model with $1/N$ filling or one particle per site.
The SU($N$) Heisenberg model is given as
\begin{equation}
\label{eq:SpinHam}
    {\mathcal H}=J \sum_{\langle {\boldsymbol r} {\boldsymbol r}' \rangle}
    S_{\alpha\beta}({\boldsymbol r})S_{\beta\alpha}({\boldsymbol r}'),
\end{equation}
where $J$ is the antiferromagnetic exchange interaction and the sum is
taken over the nearest neighbor bonds and spin flavors.
The SU($N$) spin operators
can be simply expressed with the Abrikosov fermion representation
${S_{\alpha\beta}({\boldsymbol r})=f_{{\boldsymbol r}\alpha}^{\dagger}
f_{{\boldsymbol r}\beta} }$, and ${\alpha,\beta=1,\ldots, N}$.
Hereafther, a summation over repeated indices in the form of Greek letters is supposed unless otherwise specified.
A local constraint on the fermion occupation
${f_{{\boldsymbol r}\alpha}^{\dagger}f_{{\boldsymbol r}\alpha}=1}$
is imposed to reduce the enlarged Hilbert space.
In principle, an SU$(N$) singlet could be formed by $N$ sites. But the SU($N$) exchange rapidly transforms one $N$-site singlet to other sets of $N$-site singlets, resulting in a failure of the conventional understanding.
Instead, the spin Hamiltonian
Eq.~\eqref{eq:SpinHam} is solvable in the limit $N \rightarrow \infty$
via a large-$N$ saddle point approximation in the imaginary-time
path integral formulation.
The corresponding partition function can be expressed as
\begin{eqnarray}
{\mathcal Z} = \int  {\mathcal D} \chi^{\dagger}  {\mathcal D} \chi {\mathcal D} \mu
{\mathcal D} f^{\dagger} {\mathcal D} f e^{- {\mathcal{S}}}.
\end{eqnarray}
The action ${\mathcal S}$ is given as
\begin{widetext}
\begin{eqnarray}
{\mathcal S} &=& \int_0^{\beta} d\tau \Big[
\sum_{\boldsymbol r} f^{\dagger}_{{\boldsymbol r}\alpha}  \partial_\tau^{}  f^{}_{{\boldsymbol r}\alpha}
+ \sum_{\langle  {\boldsymbol r} {\boldsymbol r}'\rangle }
\big(\chi_{  {\boldsymbol r} {\boldsymbol r}' }^{}
f^{\dagger}_{{\boldsymbol r}\alpha} f^{}_{{\boldsymbol r}'\alpha} + \mathrm{ h.c.}
\big)
+ N \sum_{\langle  {\boldsymbol r} {\boldsymbol r}'\rangle }
\frac{|\chi_{ {\boldsymbol r}{\boldsymbol r}'}^{}|^2}{\mathcal J} + \sum_{\boldsymbol r} \mu_{\boldsymbol r}^{}
( f^\dagger_{{\boldsymbol r}\alpha} f^{}_{\boldsymbol r \alpha}- 1)
\Big],
\end{eqnarray}
where $\mu_{\boldsymbol r}$ is the Lagrangian multiplier to enforce the Hilbert
space constraint, $\chi_{ {\boldsymbol r}{\boldsymbol r}'}$ is the auxiliary field
to decouple the fermion operators, and ${{\mathcal J} \equiv NJ}$. As the action
${\mathcal S}$ scales linearly with $N$, the large $N$ limit leads to a saddle point
approximation that results in the saddle point equations,
${\chi_{ {\boldsymbol r}{\boldsymbol r}'}^{}=-{\mathcal{J}}\langle {f_{ {\boldsymbol r}'\alpha}^{\dagger} f^{}_{{\boldsymbol r}\alpha}} \rangle/N}$, ${\langle {f_{ {\boldsymbol r}\alpha}^{\dagger} f_{ {\boldsymbol r}\alpha}^{}} \rangle =1}$, and
the saddle point or mean-field Hamiltonian for the fermionic spinons is
\begin{eqnarray}
\label{eq:MFTHam}
    \mathcal{H}_{\text{MF}}& =&
    \frac{N}{\mathcal{J}}\sum_{\langle {\boldsymbol{rr'}} \rangle}
    |\chi_{\boldsymbol{rr}'}^{}|^2
    +\sum_{\langle \boldsymbol{rr}' \rangle}
     (\chi_{\boldsymbol{rr}'}^{}f_{\boldsymbol{r}\alpha}^{\dagger}f_{\boldsymbol{r}'\alpha}^{}+ \mathrm{h.c.})
+\sum_{\boldsymbol{r}}\mu_{\boldsymbol{r}}^{}(1-f_{\boldsymbol{r}\alpha}^{\dagger}
f^{}_{\boldsymbol{r}\alpha}).
\end{eqnarray}
\end{widetext}
In the following, we will search for the saddle point with the
lowest mean-field energy $E_{\text{MF}}$ numerically and discuss the properties
of ground states.
Before that, we first discuss the lower bound of
$E_{\text{MF}}$ and the bound saturation conditions for reference.
An exact lower bound on $E_{\text{MF}}$
for generic lattices~\cite{PhysRevB.42.2526} was first obtained by
Rokhsar for the half-filling, and was shown to
be saturated by valence bond states with various spin singlet coverings.
These states break the lattice translation but preserve
the spin-rotation symmetry, and fluctuations beyond the mean field
can break the high degeneracy among them~\cite{Subir}.
The lower bound was generalized to
the $1/N$ filling with~\cite{PhysRevB.84.174441}
\begin{equation}
\label{eq:simplexbound}
    E_{\text{MF}} \ge - N_s \frac{N-1}{2 N}\mathcal{J}_{\text{max}},
\end{equation}
where $N_s$ is the number of lattice sites.

We set ${\mathcal{J}_{\text{max}}=\mathcal{J}}$
for each bond. The saturation for Eq.~\eqref{eq:simplexbound}
is reached by a $N$-simplex VCS state composed by $N$-site simplices.
That is, every site on the lattice is directly connected to the other ${N-1}$ sites
within the same cluster by a single bond with an exchange coupling
$\mathcal{J}_{\text{max}}$.
On a $d$-dimensional lattice, the $N$-simplex VCS with $N>d+1$ is prohibited without fine-tuning of the exchange.
Thus, there are only two-simplex
and three-simplex VCS's on the triangular lattice. For ${N \ge 4}$,
possible cluster states are general $N$-site ones satisfying a stricter
energy bound~\cite{PhysRevB.84.174441}
$E_{\text{MF}} \ge -{ \mathcal{N}_b\mathcal{J}_{\text{max}}}/{N}$,
where $\mathcal{N}_b$ is the total number of isolated bonds in the lattice.

\begin{table}[b]
    \caption{\label{tab:B}
    The SCM results of the ground state and the lowest competing state energies for
    ${2 \le N \le 9}$.
    The energy is in units of $N\mathcal{J}N_s=N^2J N_s$.}
    \begin{ruledtabular}
        \begin{tabular}{c c c}
        $N$ & Ground state & Lowest competing state \\
        \midrule
        2 & $-0.1250000$ & $-0.1202034$ \\
        3 & $-0.1111111$ & $-0.0999171$ \\
        4 & $-0.0781250$ & $-0.0760440$ \\
        5 & $-0.0581877$ & $-0.0581046$ \\
        6 & $-0.0455285$ & $-0.0443810$ \\
        7 & $-0.0364651$ & $-0.0357994$ \\
        8 & $-0.0297864$ & $-0.0293744$ \\
        9 & $-0.0247509$ & $-0.0244826$ \\
        \end{tabular}
        \end{ruledtabular}
\end{table}

\section{The SCM algorithm}
\label{sec:III}

To determine the saddle-point solutions, we closely follow the numerical self-consistent minimization
(SCM) algorithm developed in Refs.~\onlinecite{PhysRevLett.103.135301,PhysRevB.84.174441}.
The technical details of the algorithm are described in the following.
During one energy optimization process for a given geometry and periodic boundary conditions, the algorithm starts from initializing the fields $\chi_{\boldsymbol{rr}'}$
at each bond randomly with ${\chi_{\boldsymbol{rr'}}=|\chi_{\boldsymbol{rr}'}|e^{i\phi_{\boldsymbol{rr}'}}}$,
with a uniform distribution of amplitudes ${|\chi_{\boldsymbol{rr}'}| \in [0.02,0.20]}$ and
phases ${\phi_{\boldsymbol{rr}'}\in[0,2\pi]}$.
The chemical potentials are set to be the default value $\mu_{\bm{r}}=0$ in the beginning. Obviously, the local constraints are violated here in general and the deviation of the local fermion density can be denoted as $\delta n_{\bm{r}} = 1 - \braket{f_{\bm{r}\alpha}^{\dagger}f_{\bm{r}\alpha}}$. The expectation value is taken using the ground state of $\mathcal{H}_{\text{MF}}$ with the unchanged $\mu_{\bm{r}}$ at this stage. To obtain the correct density, one must adjust the chemical potential $\mu_{\bm{r}}$ by $\delta \mu_{\bm{r}}$ at each site. It is found that to the lowest order, the desired adjustment of chemical potentials $\delta\mu_{\bm{r}}$ can be expressed as
\begin{equation}\label{eq:response}
    \delta \mu_{\bm{r}} = \sum_{\bm{r}'}G^{-1}(\bm{r}-\bm{r}',0)\delta n_{\bm{r}'},
\end{equation}
where $G^{-1}(\bm{r}-\bm{r}',0)$ is nothing but the inverse of the density-density correlation at zero frequency.
The elements of the correlation $G(\bm{r}-\bm{r}',0)$ can be calculated at the single particle level from the mean-field Hamiltonian $\mathcal{H}_{\text{MF}}$, with the help of the standard linear response theory in principle.
However, it follows that the inversion of $G(\bm{r}-\bm{r}')$ naively diverges.
Physically, this is because a uniform adjustment of the chemical potential at every site is trivial and the fermion density will keep intact.
This obstacle is overcome by following the treatment suggested in Ref.~\onlinecite{PhysRevB.84.174441}; that is diagonalizing $G(\bm{r}-\bm{r}',0)$ and only focusing on its non-zero eigenvalues $g_{i}$.
Note that the index $i$ refers to the basis where $G(\bm{r}-\bm{r}',0)$ is diagonalized.
In such a basis, the relationship Eq.~\eqref{eq:response} becomes
\begin{equation}
    \delta \mu_{i} = \delta n_{i} / g_{i}.
\end{equation}
Here, the index $i$ should not be summed.
For all indices with vanishing eigenvalues $g_i$, $\mu_{i}$ is taken to be zero concurrently.
Then, the adjustment of the chemical potential $\delta \mu_{\bm{r}}$ is well-defined.
A simple replacement of $\mu_{\bm{r}}\rightarrow\mu_{\bm{r}}+\delta \mu_{\bm{r}}$ gives a new mean-field Hamiltonian $\mathcal{H}_{\text{MF}}$ and related ground state.
In consequence, it affects the local fermion density conversely, resulting in a new deviation $\delta n_{\bm{r}}$.
These processes construct a self-consistent procedure and the problem of searching for an appropriate set of chemical potential deviation $\delta \mu_{\bm{r}}$ can be solved by iterating the procedure until the density is uniform.
This is the core of the algorithm to preserve the local single-occupation constraints $n_{\bm{r}}=\braket{f_{\bm{r}\alpha}^{\dagger}f_{\bm{r}\alpha}}=1$.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.48\textwidth]{Fig02.pdf}
    \caption{Ground states in the large-$N$ limit for (a) ${N=2}$, (b) ${N=3}$, (c) ${N=4}$ and (d) ${N=5}$.
    For $2 \le N \le 4$, any state covering all sites by the same dimer/cluster unit has the same energy in each case.
    We merely show one of such configurations. For ${N=5}$, there is a stripe pattern with doubled
    unit-cell along one of lattice vectors. The grayness of bonds indicates the relative magnitude of expectation values $|\chi_{\boldsymbol{rr}'}|$. Black (white) bonds refer to maximal (zero) $|\chi_{\boldsymbol{rr}'}|$.}
    \label{fig:02}
\end{figure}

With the modified chemical potentials calculated in the previous step, an update set of auxiliary fields $\chi_{\bm{rr}'}$ can be determined via ${\chi_{\boldsymbol{rr}'}'=-{\mathcal{J}}\langle {f_{\boldsymbol{r}'\alpha}^{\dagger}f_{\boldsymbol{r}\alpha}} \rangle/N}$, but the local constraints may be violated once again.
The amended auxiliary fields and chemical potentials can be treated as a new and better starting point.
The procedure described in the previous step is thus implemented iteratively until reaching a converge of the ground state energy within a given numerical error.
It has been shown in Ref.~\onlinecite{PhysRevB.84.174441} that the energy of the final state must be less or equal to that of the initial state after the optimization process.
Therefore, we eventually obtain a local energy minimum.
In order to reach the global minimum as much as possible, we start from at least 50 randomly initialized fields, and reap a collection of local minima satisfying the single-occupation constraints.
The lowest two are accepted as the best results of the ground state and lowest competing state energies.

\section{The mean-field results}
\label{sec:IV}

We describe the ground states and lowest competing states of the mean-field Hamiltonian
Eq.~\eqref{eq:MFTHam} from the SCM algorithm. Because different geometries can
accommodate different candidate ground states especially cluster states with unknown dimensions,
in the calculation, we consider all unit cells of a parallelogram geometry $\ell_1 \times \ell_2$ with
${\ell_{1,2} \le 2N}$ for ${2 \le N \le 5}$ and ${\ell_{1,2} \le N}$ for ${N \ge 6}$, respectively.
In some ordered cases for $N=4,5$, larger unit-cell sizes are also
considered for confirmation. Each unit cell is repeated by $L_{1,2}$ $({\ge 20})$ times
along the triangular Bravais lattice vectors to form a superlattice.
The superlattice itself has periodic boundary conditions.
In practice, for each case, we ran the SCM procedure at least 50 times with different random seeds to avoid any missing of ground states caused by numerical problem. The results are listed in Table~\ref{tab:A} and Table~\ref{tab:B}.

\begin{figure}[t]
    \centering
    \includegraphics[]{Fig03.pdf}
    \caption{The lowest competing state in the large-$N$ limit for ${N=4}$.
    This is a stripe pattern with doubled unit-cell along one of lattice vectors.
    The grayness of bonds indicates the relative magnitude of expectation value $|\chi_{\boldsymbol{rr}'}|$.}
    \label{fig:03}
\end{figure}

For $N=2$ and $3$, the lowest energies we found saturate the bound Eq.~\eqref{eq:simplexbound}, meaning that the ground states are dimer and three-site simplex VCS, respectively.
Actually, the obtained ground states are highly degenerate because any
state with each site covered by one dimer/cluster unit has the
same energy in the large-$N$ limit. An ordered dimer or three-site simplex VCS state, as illustrated in Figs.~\ref{fig:02}(a-b), is expected to be selected beyond the mean field.
Note that, the non-zero expectation values $|\chi_{\bm{rr}'}|$ in Fig.~\ref{fig:01}(a) and (b) can only take $1/2$ and $1/3$, respectively.
The true ground state for $N=2$ is the well-known
$120^{\circ}$ N\'{e}el state and differs from the mean-field here.
The $N=3$ case was also shown to have a three-sublattice magnetic order in previous numerics
~\cite{PhysRevLett.97.087205,PhysRevB.85.125116}.
This reflects the deficiency of the
large-$N$ approximation for small $N$ cases. In Table~\ref{tab:A}, we further find that
the lowest competing states are CSL states with $\phi = \pi-\pi/N$.

The $N=4$ case is more compelling due to the rapid growth of experimental proposals \cite{PhysRevLett.117.167202,PhysRevLett.121.225303,Ashvin:2020,zhang2021su4}
and large-$N$ approximation could provide useful insight here. We find the ground state
energy do not saturate the bound Eq.~\eqref{eq:simplexbound} as excepted, but saturate
the stricter bound discussed previously. The ground states are four-site VCS depicted in
Fig.~\ref{fig:02}(c), accompanied by a large degeneracy for
the same reason as $N=2$ and $3$.
The non-zero expectation values are $|\chi_{\bm{rr}'}|=1/4$.
A similar plaquette order was also reported in a recent work~\cite{zhang2021su4}.
In Fig.~\ref{fig:03}, we depict the lowest competing state for ${N=4}$.
One can see that the lattice is covered by stripes with three different bond expectation values $|\chi_{\boldsymbol{rr}'}|\approx 0.030$, $0.158$, and $0.224$.
Along one of lattice vectors, there is a unit-cell doubling
that breaks the lattice translation.
The background fluxes through each plaquette are inhomogeneous as well and manifest the same unit-cell doubling pattern.
Specifically, the average flux is a constant value $\phi_{\text{avg}}=\pi-\pi/N$ where $N=4$.
The same ordered pattern has also been obtained through the DMRG study from the Kugel-Khomskii
model~\cite{Chaoming:2019}. They attribute the symmetry breaking to symmetry allowed
Umklapp interactions in certain finite geometries. In our results, however, such a stripy state is not the ground state even in the two-dimensional limit.

From ${N=5}$, neither the bound Eq.~\eqref{eq:simplexbound} nor
the stricter one can be saturated. Thus, the ground states are no longer VCS.
Our numerical calculation finds that the ground state for ${N=5}$ is very similar
to the lowest competing state for $N=4$ (see Fig.~\ref{fig:02}(d)), except that the bond expectation values can only take $|\chi_{\boldsymbol{rr}'}|\approx 0.072$, $0.145$, and $0.179$, and the average background flux $\phi_{\text{avg}}=\pi-\pi/N$ shifts accordingly.
It also breaks the lattice translation symmetry along one of lattice vectors
and manifests itself as a stripe pattern with a doubled unit-cell.
The lowest competing state for $N=5$ is a CSL with $\phi = 4\pi/5$.
With further increasing $N$, the frustration is enhanced. Eventually, the ground states
become CSL states with $\phi = \pi-\pi/N$ when ${6\le N \le 9}$.
Correspondingly, the lowest competing states we found share the identical form.
They are CSL states as well except that the background magnetic fluxes shift
to $\phi = \pi -\pi/(2N)$.

With the two types of CSL states identified, we now discuss the properties of these topological liquid states. The CSL is characterized by the mean-field saddle point
\begin{gather}
    \chi_{\boldsymbol{rr}'}=|\chi| e^{i a_{\boldsymbol{rr}'}},\\
    \mu_{\boldsymbol{r}}=0.
\end{gather}
All bonds on the lattice have a uniform expectation value for $|\chi|$ but modulated by a
U(1) gauge field $a_{\boldsymbol{rr}'}$ so that the flux $\phi$ on each plaquette is a constant.
The CSL breaks both parity and time-reversal. The bond phase field $a_{\boldsymbol{rr}'}$ is
treated as a fluctuating $\mathrm{U}(1)$ gauge field coupled to the fractionalized spinons.
By checking the behavior of the mean-field Hamiltonian Eq.~\eqref{eq:MFTHam} at the CSL saddle points,
we find that both types of CSL states have a fermion band structure with $N$ bands where
only the lowest is filled (see Fig.~\ref{fig:01}(b)). The Fermi level lies in the gap between the lowest two bands,
thus all discussion can be applied to both CSL states.
Furthermore, the first type of CSL with $\phi_{u,d}=\pi-\pi/N$ on the triangular lattice can be mapped to the counterpart on a square lattice up to a time-reversal. If we regard the adjacent up and down triangles as a unit shown in Fig.~\ref{fig:04}, the phase of the shared bond has no contribution to the total flux, and we have the relationship
\begin{equation}
    \phi_{u}+\phi_{d}=-\phi_{s} \mod 2\pi,
\end{equation}
where $\phi_s$ is the background U(1) flux through each square plaquette. In Ref.~\cite{PhysRevLett.103.135301}, Hermele \textit{et al}, found that, the CSL states are ground states for $5 \le N \le 10$ on the square lattice in the mean-field calculation.
As the spinon is gapped out by the emergent U(1) gauge flux pattern,
the Chern-Simons term enters the theory for U(1) gauge fluctuations.
After integrating the gapped spinons, one obtains a topological quantum field theory
with a Chern-Simons term corresponding to the chiral Abelian topological order
and anyonic statistics. The spinon is converted in anyons with a statistical angle ${\pi \pm \pi/N}$.
Gapless chiral states carrying spin degree of freedom are also supported by the CSL as edge modes and the low-energy theory is described by the SU(N)$_1$ WZW model.

\section{Discussion}
\label{sec:V}

To summarize, we study the Heisenberg antiferromagnets with SU($N$) symmetry on the triangular lattice.
In the large-$N$ approximation, a variety of ground states and lowest competing states are identified for different parameter $N$.
At the mean-field level, the ground state for $2 \le N \le 4$ is a $N$-site VCS state
with a large degeneracy. Specifically, ordering patterns with doubled unit cell and average background flux $\phi_{\text{avg}}=\pi-\pi/N$ are found
in the $N=4$ and $5$ cases. These ordered states break the lattice translation symmetry
along one of lattice vectors and become the lowest competing state and ground state for
$N=4$ and $5$, respectively. The frustration from SU($N$) exchange interaction
is enhanced when $N > 5$, resulting two types of CSL states as the lowest two states for
$6 \le N \le 9$. Among them, the CSL states with $\phi=\pi-\pi/N$ have a lower energy,
and have its counterpart on the square lattice. Although the true ground states
are not what we found for $N<4$, our calculation can provide useful insight for
cases $N \ge 4$ where large-$N$ approximation becomes more reliable.
In a very recent DMRG study of an SU(4) spin model on the triangular lattice, phase diagrams for integer fillings were obtained and compared with conventional MFT ones \cite{zhang2021su4}.
The quite good agreement of the phase boundaries determeted by two methods suggests that $N = 4$ is perhaps large enough for the mean-field analysis we performed in this work.

\begin{figure}[t]
    \centering
    \includegraphics[]{Fig04.pdf}
    \caption{The definition of the emergent U(1) gauge flux $\phi$ for CSL on the (a) square and (b) triangular lattices.
    The subscripts $s$, $u$ and $d$ refer to the square, up triangular and down triangular plaquettes, respectively.}
    \label{fig:04}
\end{figure}

Thanks to the development of ultracold experimental techniques,
the SU($N$) Mott insulators have been realized with AEA
on various optical lattices using the Pomeranchuk cooling~\cite{Taie_2012},
even the Mott crossover and SU($N$) antiferromagnetic spin correlations were recently observed
with \ce{^173Yb} atoms~\cite{PhysRevX.6.021030,taie2020observation}.
Nontrivial physics of multicomponent fermions with broken SU$(N)$ symmetry were also proposed on this platform~\cite{PhysRevA.98.063628}.
The emergency of the widely concerned SU(4) or even SU(8) symmetric
interaction has been proposed in the twisted bilayer graphene
and double moir\'{e} layers systems recently~\cite{Ashvin:2020,PhysRevLett.121.097201}.
This moir\'{e} lattice system could be an novel platform
to detect possible phases in this work.
The ultracold atom systems may have many restrictions
in the detection of anyonic spinon excitations and edge states.
 Nevertheless, spin-dependent Bragg
spectroscopy may be used to detect the spinon continuum~\cite{PhysRevB.84.174441,
PhysRevA.93.061601,PhysRevLett.103.135301},
 singlet-triplet oscillation technique can detect the nearest-neighbor spin correlation~\cite{taie2020observation} for CSLs,
 and the lattice potential could be adjusted to
 localize and manipulate the anyonic quasiparticles~\cite{PhysRevA.93.061601,PhysRevLett.103.135301,PhysRevB.84.174441}.
 For solid-state systems, more techniques are available.
 These include but are not restricted to quantized thermal Hall transport~\cite{Gao_2020}
 for the edge modes, scanning tunneling microscope of anyons at defects~\cite{PhysRevX.8.011037},
 or even angle-resolved photon emission measurement for the
spinon signatures~\cite{PhysRevB.87.045119}.

\begin{acknowledgments}
We thank Shizhong Zhang for some discussion
and the proofreading of the manuscript. This work is supported by the
research funds from the Ministry of Science and Technology of China
with grant No.2016YFA0301001, No.2018YFE0103200 and No.2016YFA0300500,
by Shanghai Municipal Science and Technology Major Project with Grant
No.2019SHZDZX04, and from the Research Grants Council of Hong Kong
with General Research Fund Grant No.17303819 and No.17306520.
\end{acknowledgments}