Deriving Generating function for centered trinomial coefficients

Here is a variation based upon G. P. Egorychev's Classic: Integral Representation and the Computation of Combinatorial Sums. We start with the central trinomial coefficients: \begin{align*} [x^n](1+x+x^2)^n\qquad\qquad n\geq 0 \end{align*} We consider the function \begin{align*} f(x)=1+x+x^2\tag{1} \end{align*} and derive a function $y=y(x)$: \begin{align*} y(x)=\frac{x}{f(x)}=\frac{x}{1+x+x^2}\qquad\qquad y^{\prime}(x)=\frac{1-x^2}{(1+x+x^2)^2 }\tag{2} \end{align*}

With $f(x)$ and $y(x)=\frac{x}{f(x)}$ we can now apply the substitution rule (Rule 5, one-dimensional case) from section 1.2.2 in G. P. Egorychev's book as follows: \begin{align*} \color{blue}{[x^n](f(x))^n=[y^n]\left.\left(\frac{1}{f(x)y^{\prime}(x)}\right)\right|_{x=g(y)}}\tag{3} \end{align*} with $g(y)$ the inverted function given by $y=y(x)$ in (2).

We obtain from (1) - (3): \begin{align*} \color{blue}{[x^n]}&\color{blue}{\left(1+x+x^2\right)^n}\\ &=[y^n]\left.\left(\frac{1}{\left(1+x+x^2\right)\frac{d}{dx}\left(\frac{x}{1+x+x^2}\right)}\right)\right|_{x=g(y)}\\ &=[y^n]\left.\frac{1+x+x^2}{1-x^2}\right|_{x=g(y)}\\ &\,\,\color{blue}{=[y^n]\frac{1}{\sqrt{1-2y-3y^2}}}\tag{4} \end{align*} and the claim follows.

In (4) we use the identity \begin{align*} 2y=\frac{2x}{1+x+x^2}&=1-3\left(\frac{x}{1+x+x^2}\right)^2-\left(\frac{1-x^2}{1+x+x^2}\right)^2\\ &=1-3y^2-\left(\frac{1-x^2}{1+x+x^2}\right)^2\\ \frac{1+x+x^2}{1-x^2}&=\left(1-2y-3y^2\right)^{-\frac{1}{2}} \end{align*}

$c_n$ is the coefficient of $x^n$ in $(1 + x + x^2)^n$. It follows that its generating function is the diagonal of the rational generating function

$$F(x, y) = \frac{1}{1 - y(1 + x + x^2)} = \sum_{n \ge 0} y^n (1 + x + x^2)^n = \sum f_{n, m} x^n y^m$$

in the sense that $c_n = f_{n, n}$. It's a general fact (which you can find stated, for example, as Theorem 6.3.3 in Stanley's Enumerative Combinatorics, Vol. II) that the diagonal of a bivariate rational generating function is algebraic and can be computed using contour integration, as explained in Stanley, and you can also see my blog post Extracting the diagonal. We can do the computation as follows. Write $C(r) = \sum c_n r^n$. Then for sufficiently small $r$ we have

$$\frac{1}{2 \pi i} \int_{\gamma} \frac{F(rz, rz^{-1})}{z} \, dz = C(r^2)$$

where $\gamma$ is the contour given by the unit circle. In our case the integrand is

$$\frac{F(rz, rz^{-1})}{z} = \frac{1}{z - r - r^2 z - r^3 z^2}$$

which, as a meromorphic function of $z$, has poles given by the zeroes of the denominator. These are the zeroes of a quadratic $r^3 z^2 + (r^2 - 1) z + r$, which are then

$$z_0, z_1 = \frac{(1 - r^2) \pm \sqrt{1 - 2r^2 - 3r^4}}{2r^3}$$

by the quadratic formula. We only need to consider the residue at a pole inside our contour for small $r$, and as $r \to 0$ the $+$ zero goes to infinity so we only need to consider the $-$ zero

$$z_0 = \frac{(1 - r^2) - \sqrt{1 - 2r^2 - 3r^4}}{2r^3}.$$

The residue at this pole is

$$\lim_{z \to z_0} \frac{z - z_0}{-r^3(z - z_0)(z - z_1)} = \frac{1}{-r^3(z_0 - z_1)} = \frac{1}{\sqrt{1 - 2r^2 - 3r^4}}$$

so the residue theorem gives

$$C(r^2) = \frac{1}{\sqrt{1 - 2r^2 - 3r^4}}$$

as desired.

Now some more general facts can be used to deduce asymptotics. The dominant singularity of $C(z) = \frac{1}{\sqrt{1 - 2z - 3z^2}} = \frac{1}{\sqrt{(1 - 3z)(1 + z)}}$ occurs at $z = \frac{1}{3}$. Around this singularity $C(z)$ looks like $\frac{1}{\sqrt{\frac{4}{3}(1 - 3z)}}$ which gives (using e.g. binomial expansion together with Stirling's formula) that the leading order asymptotic of $c_n$ is

$$\boxed{ c_n \sim \sqrt{\frac{3}{4 \pi n}} \, 3^n }.$$

This is in agreement with the comment left by Vaclav Kotesovec on the OEIS page, and in particular implies that the true value of $\lim_{n \to \infty} \frac{c_{n+1}}{c_n}$ is $3$ exactly. For much more on this topic see Chapter VI.1 of Flajolet and Sedgewick's Analytic Combinatorics.