How to find the asymptotics of a linear two-dimensional recurrence relation Let $d$ be a positive number.
There is a two dimensional recurrence relation as follow:
$$R(n,m) = R(n-1,m-1) + R(n,m-d)$$
where
$R(0,m) = 1$ and $R(n,0) = R(n,1) = \cdots = R(n, d-1) = 1$ for all $n,m>0$.
How to analyze the asymptotics of $R(n, kn)$ for fixed $k$?
It is easy to see that
$$R(n, kn) = O\left( c_{k,d}^{n} \cdot (n+k+d)^{O(1)} \right)$$
Is there a way (or an algorithm) to find $c_{k,d}$ given $k$ and $d$?
PS: I have calcuated the bivariate generating function of $R(\cdot, \cdot)$:
\begin{align}
f(x,y) 
&= \frac{1 - xy - y^{d} + xy^{d}}{(1 - x)(1 - y)(1 - xy - y^{d})} \\
&= \frac{1}{(1 - x)(1 - y)} + \frac{xy^{d}}{(1 - x)(1 - y)(1 - xy - y^{d})} \\
\end{align}
 A: Your generating function $f(x,y)$ is convergent on a polydisk $|x|<\epsilon$, $|y|< \delta$, for some $\epsilon, \delta < 1$. We can reduce $\epsilon$ so that $\epsilon \ll \delta^k$ and the domain of $f(x,y)$ to the product of the polydisk $|x|<\epsilon$ and the annulus $\delta' < |y| < \delta$, with $\delta' \sim \delta$. Then $g(z,y) = f(z/y^k, y)$ is well-defined and analytic on the product domain of $|z|\lesssim \epsilon/\delta^k$ and $\delta' < |y| < \delta$. The generating function $g(z,y)$ has a unique expansion that is a power series in $z$ and a Laurent series in $y$. Since $x^n y^m = (z/y^k)^n y^m = z^n y^{m-kn}$, the coefficients $R(n,kn)$ are the coefficients of $z^n y^0$ in this double expansion. Hence, using the Cauchy integral formula, the generating function for $R(n,kn)$ is now given by the contour integral
$$
  \sum_{n=0}^\infty R(n,kn) z^n = \frac{1}{2\pi i} \oint_{|y|\sim \delta} g(z,y) \frac{dy}{y} .
$$
Of course, the integration contour is now allowed to deform, as long as it doesn't hit any singularities of the integrand.
The above integrand has a convenient partial fraction expansion expansion with respect to $z$,
$$
  \frac{g(z,y)}{y} = \frac{y^k}{(1-y-y^d)(y^k-z)}
    - \frac{y^{d+k-1}(1-y^d)}{(1-y)(1-y-y^d) (y^{k-1} (1-y^d)-z)} .
$$
Expanding in $z$ under the integral gives the formula
$$
  R(n,kn) = \frac{1}{2\pi i} \oint \frac{dy}{y^{kn} (1-y-y^d)}
    - \frac{1}{2\pi i} \oint \frac{y^d dy}{(1-y)(1-y-y^d) y^{n(k-1)} (1-y^d)^n} ,
$$
whose leading asymptotics can be estimated by residues or steepest descent.
Using residues on the first term, all the contributions come from the roots of $1-y-y^d=0$. The root with the smallest magnitude $y_*$ will give the leading contribution. Some experimentation shows that $y_*$ is the unique positive real root. The leading asymptotic term then looks like
$$
  C_* y_*^{-kn} .
$$
In the special case $d=1$, $y_* = 1/2$. In the special case $k=1$, the second term has no poles inside its contour and hence evaluates to zero. So in the sequel we can assume that $k>1$.
Applying steepest descent to the second integral, we find the stationary phase points at the roots $y_\star$ of $y^d = (k-1)/(k+d-1) + O(1/n)$. Let $y_\star$ be that root (actually we just need its $n\to \infty$ limit) which minimizes the magnitude of $y_\star^{k-1} (1-y_\star^d)$. It is easy too see that $y_\star = \sqrt[d]{(k-1)/(d+k-1)}$. Then the leading asymptotic term is
$$
  C_\star n^{-1/2} y_\star^{-n(k-1)} (1-y_\star^d)^{-n}
  = C_\star n^{-1/2} w_\star^{-kn} ,
$$
where $w_\star = \left(\frac{k-1}{d+k-1}\right)^{\frac{1-1/k}{d}} \left(\frac{d}{d+k-1}\right)^{1/k} = \frac{s^{s/k}}{(s+1)^{(s+1)/k}}$ with $s=(k-1)/d$, as pointed out in the comments.
The conclusion is that
$$
  R(n,kn) = O\left((\min(y_*,w_\star)^{-k})^n\right) .
$$
Experimentation suggests that $y_* < w_\star$ in all cases, except for $(k,d)=(2,1)$, when $y_* = w_\star = 1/2$. Otherwise the $y_*$ contribution always seems to dominate over the $w_\star$ contribution. At least that's how it seemed for not too large values of $d$.
Doing some rudimentary asymptotic calculations for large $d$, it seems that $y_* \sim 1- \frac{\log d}{d} + o(\log(d)/d)$. A better approximation is
$$
  y_* \sim 1 - \frac{W(d)}{d} + \frac{W(d)^3}{2(W(d)+1) d^2} + O(W(d)^3/d^3) ,
$$
where $W(d)$ is the Lambert W function. On the other hand, $w_\star$ has a minimum as a function of $k$ (treating it as a continuous variable) around $k_\star \sim W(d) + 1 + O(W(d)^2/d)$, the minimum reaches roughly
$$
   w_\star \sim 1 - \frac{W(d)}{d-1}  \lesssim y_* .
$$
One should push the asymptotics of $k_\star$ one more order to get a better estimate of the difference, but experiments do show that the minimum of $w_\star$ does dip below $y_*$, for instance this happens near $(k,d) = (5, 205)$ or $(6,700)$. So it looks like the $w_\star$ asymptotic contribution will dominate in small ranges of $k \sim k_\star(d)$ if it happens to fall near an integer. To get the size of that window, one could go to a quadratic Taylor approximation of $w_\star$.
With a bit of extra work, one could also extract the coefficients $C_*$ and $C_\star$ from the integral formula.
A: This is to complement Blanco's answer by showing that
\begin{equation*}
    R(n,kn)=\exp\{(C_{k,d}+o(1))\,n\} \tag{0}\label{0}
\end{equation*}
(as $n\to\infty$), where $k\ge1$ and $d\ge1$ (are fixed),
\begin{equation*}
    C_{k,d}:=\frac k{1+y_{k,d}\,d}\,\big(\ln(1+y_{k,d})+y_{k,d}\,\ln(1+1/y_{k,d})\big),
\end{equation*}
\begin{equation*}
    y_{k,d}:=\max\Big(x_d,\frac{k-1}d\Big), 
\end{equation*}
and $x_d$ is the unique positive root of the equation
\begin{equation*}
    x_d(1+x_d)^{d-1}=1. \tag{0.5}\label{0.5}
\end{equation*}
In view of Blanco's answer, it is enough to show that
\begin{equation*}
    B(n):=B(n,kn)=\exp\{(C_{k,d}+o(1))\,n\}, \tag{1}\label{1}
\end{equation*}
where
\begin{equation*}
    B(n):=\max\{c_{a,b}\colon (a,b)\in E_{n,k,d}\},
\end{equation*}
\begin{equation*}
    E_{n,k,d}:=\{(a,b)\colon 0\le a\le n-1,b\ge0,bd+a\le kn-d,\ a,b \text{ are integers}\},
\end{equation*}
\begin{equation*}
    c_{a,b}:=\binom{a+b}b=\binom{a+b}a. 
\end{equation*}
Note that $c_{a,b}$ is increasing in $a$ and in $b$.
Note also that $c_{a,b}=(a+b)^{O(1)}=n^{O(1)}$ for $(a,b)\in E_{n,k,d}$ if $a=O(1)$ or $b=O(1)$.
So, it remains to consider the case when $a\to\infty$ and $b\to\infty$. Then, by Stirling's formula,
\begin{equation*}
\begin{aligned}
    \ln c_{a,b}\sim a\ln(1+b/a)+b\ln(1+a/b). 
\end{aligned}
\tag{2}\label{2}
\end{equation*}
Also,
\begin{equation*}
\begin{aligned}
    \frac{c_{a+d,b-1}}{c_{a,b}}&=b\frac{(a+d+1)\cdots(a+d+b-1)}{(a+1)\cdots(a+b)} \\ 
    &= b\frac{(a+b+1)\cdots(a+b+d-1)}{(a+1)\cdots(a+d)} \\ 
    &\sim b\frac{(a+b)^{d-1}}{a^d} =\frac ba\Big(1+\frac ba\Big)^{d-1}. 
\end{aligned}
\end{equation*}
So, for a fixed value of $bd+a$, the maximum of $c_{a,b}$ occurs when $b/a\to x_d$ (recall \eqref{0.5}).
Let now $(a,b)\in E_{n,k,d}$ be a maximizer of $c_{a,b}$ (such that $a\to\infty$ and $b\to\infty$). Then, since $c_{a,b}$ is increasing in $a$ and in $b$, we have
\begin{equation*}
    bd+a\sim kn. 
\end{equation*}
The conditions $b/a\to x_d$ and $bd+a\sim kn$ imply
\begin{equation*}
a\sim k\frac1{1+x_d\,d}\,n, \tag{3}\label{3}
\end{equation*}
and the latter condition is compatible with condition $a\le n-1$ only if
\begin{equation*}
    k\frac1{1+x_d\,d}\le1. \tag{4}\label{4} 
\end{equation*}
If this is the case, then
\begin{equation*}
b\sim k\frac{x_d}{1+x_d\,d}\,n, 
\end{equation*}
so that, by \eqref{2},
\begin{equation*}
\begin{aligned}
    \frac{\ln c_{a,b}}n\to \frac k{1+x_d\,d}\,(\ln(1+x_d)+x_d\,\ln(1+1/x_d)). 
\end{aligned}
\tag{5}\label{5}
\end{equation*}
Also, \eqref{4} implies $x_d\ge(k-1)/d$, so that $y_{k,d}=x_d$.
So, we have proved \eqref{1}, and thus \eqref{0} -- in the case when condition (3) is compatible with condition $a\le n-1$.
Otherwise, we have $a=n-1\sim n$ and still $bd+a\sim kn$, whence $b\sim(k-1)n/d$. So, by \eqref{2}, here
\begin{equation*}
\begin{aligned}
    \frac{\ln c_{a,b}}n\to \ln\Big(1+\frac{k-1}d\Big)+\frac{k-1}d\,\ln\Big(1+\frac d{k-1}\Big). 
\end{aligned}
\tag{6}\label{6}
\end{equation*}
Also, in this "incompatibility" case, we have $k\frac1{1+x_d\,d}\ge1$ -- cf. \eqref{4}. So, here $x_d\le(k-1)/d$ and hence $y_{k,d}=\frac{k-1}d$.
So, we have proved \eqref{1}, and thus \eqref{0} -- in the case when condition (3) is incompatible with condition $a\le n-1$.
Thus, in either case, \eqref{0} is proved. $\quad\Box$
A: The book Analytic Combinatorics in Several Variables by Pemantle and Wilson covers problems like this extensively.
A: This is a part of my approach, and I have not finished.
Firstly,
$$
\frac{1}{(1-x)(1-y)} = \sum_{i,j \geq 0} x^{i}y^{j}
$$
and
$$
\frac{1}{1 - xy - y^{d}} = \sum_{i \geq 0} \frac{x^{i}y^{i}}{(1 - y^{d})^{i + 1}} = 
\sum_{i, j \geq 0} \binom{i + j}{j} x^{i}y^{jd + i}$$
Thus
\begin{align}
\frac{1}{(1 - x)(1 - y)(1 - xy - y^{d})}
= &\left( \sum_{i,j \geq 0} x^{i}y^{j} \right) \cdot \left( \sum_{i, j \geq 0} \binom{i + j}{j} x^{i}y^{jd + i} \right) \\
= &\sum_{i,j \geq 0} \sum_{a \leq i, bd + a \leq j} \binom{a + b}{b} x^{i}y^{j}  \\
\end{align}
We have
$$f(x, y) = \sum_{i,j \geq 0} \left( 1 + \sum_{a \leq i - 1, bd + a \leq j - d} \binom{a + b}{b} \right) x^{i}y^{j}$$
and
$$R(n, m) = 1 + \sum_{a + 1 \leq n, (b + 1)d + a \leq m} \binom{a + b}{b}$$

Let
$$B(n,m) = \max_{a + 1 \leq n, (b + 1)d + a \leq m} \binom{a + b}{b}$$
then we have
$$B(n,m) \leq R(n,m) \leq nmB(n,m)$$
Thus $R(n, m) = O(B(n,m))\cdot (nm)^{O(1)}$
