Approximating power series coefficients --- Why does a clearly illegitimate method (sometimes) work so well? For reasons that don't matter here,
I want to estimate the power series coefficients
$t_{ij}$ for the rational function
$$T(x,y)=  {(1+x)(1+y)\over 1- x y(2+x+y+x y)}=\sum_{i,j} t_{ij}x^iy^j$$
Using a method that I cannot justify, I  get
highly accurate estimates when $i=j$ and highly inaccurate estimates when
$|i-j|$ strays at all far from zero.
My questions are:

Q1)  Why does my apparently illegitimate method work so well when $i=j$?


Q2)  Why does the answer to Q1) not apply when $i\neq j$ ?

(Of course, once the answer to Q1) is known, the answer to Q2) might be
self-evident.)

I'll first present the method, then explain why I think it shouldn't work,
then present the evidence that it works anyway when $i=j$, and then present
the evidence that it rapidly goes haywire when $i\neq j$.
The Apparently Illegitimate Method:
Note that $t_{ij}=t_{ji}$, so we can limit ourselves to estimating
$t_{j+k,j}$ for $k\ge 0$.
I)  Define
$$T_k(y)=\sum_jt_{k+j,j}y^j$$
For example, a residue calculation gives
$$ T_0(y)= {1-y-\sqrt{1-4y+2y^2+y^4}\over y\sqrt{1-4y+2y^2+y^4}} $$
It turns out that all of the $T_k$ share a branch point at $\zeta\approx .2956$ and are analytic in the disc $r<\zeta$.
II)  Write
$$L_k=\lim_{y\mapsto \zeta} T_k(y)\sqrt{y-\zeta}$$.
Discover that $L_0\approx 1.44641$ and $L_k=L_0/\zeta^{k/2}$.
III)  Approximate
$$T_k(y)\approx L_k/\sqrt{y-\zeta}$$
IV) Expand the right hand side in a power series around $y=0$ and equate
coefficients to get
$$t_{ij}\approx \pm{L_0\over\sqrt{\zeta}}\pmatrix{-1/2\cr j\cr}\zeta^{-(i+j)/2}
\approx \pm  2.66036 \pmatrix{-1/2\cr j\cr}\zeta^{-(i+j)/2}\qquad(E1)$$
Remarks:

*

*Obviously one could try to improve this approximation
at Step III by using more terms in the power series for $T_k$ at $y=\zeta$.
This doesn't seem to help, except when $k=0$, in which case the original approximation is already quite good.


*For $k\ge 2$, $T_k(y)$ has a zero of
order $k-1$ at the origin.  Thus one could modify this method by approximating
$T_k(y)/(y^{k-1})$ instead of $T_k(y)$
This yields
$$t_{ij}\approx \pm{2.66036}\pmatrix{-1/2\cr 1-i+2j}\zeta^{-(i+j)/2}\qquad(E2)$$
(E2) is (much) better than (E1) in the range $i\ge 2j+1$, where it gets
exactly the correct value, namely zero.  Otherwise, it seems neither systematically better nor worse.
Why Nothing Like This Should Work:   The expansion of $T_k(y)$ at
$\zeta$ contains nonzero terms of the form
$A_{i,j}(\zeta-y)^j$ for all positive integers $j$.  (I'm writing $i=j+k$ to
match up with the earlier indexing.)  The truncation at Step III throws all
these terms away.  Therefore the expansion around the origin in Step IV
ignores (among other things) the contribution of $A_{ij}$ to the estimate
for $t_{ij}$.  So unless we can control the sizes of the $A_{ij}$, we
have absolutely no control over the quality of the estimate.
And in fact, even when $k=0$, the $A_{j,j}$ are not small.
For example, $t_{8,8}=8323$ and  my estimate for $t_{8,8}$ is a
respectable $8962.52$.  But $A_{8,8}$, which should have contributed to that
estimate and got truncated away, is equal to $58035$.  It seems remarkable
that I can throw away multiple terms of that size and have the effects nearly cancel.
I'd like a conceptual explanation for this.
But When $i=j$, It Works Anyway:

and these get even better if you truncate just slightly farther out.
Why any explanation can't be too general:

 A: The paper 'A New Method for Computing Asymptotics of Diagonal Coefficients of Multivariate Generating Functions' by A. Raichev and M. Wilson has the precise machinery that can solve this problem.  Get a copy and these brief notes correspond to their symbols, for the diagonal case
$$ f_{n,n} = [x^n \, y^n ] \frac{I(x)}{J(x)} = [x^n \, y^n ] \frac{(1+x)(1+y)}{1-xy(2+x+y+xy)}.$$
Solve the simultaneous system for zeros of the denominator $J$; in Mathematica,
Solve[ { x D[J,x] == y D[J,y], J==0 },{x,y} ].
The proper solution must have both $x$ and $y$ positive.  That set is
$$ \mathbf c=(\rho,\rho),\,\rho=(\tau-2/\tau-1)/3, \, \tau=(17+3\sqrt{33})^{1/3} \approx 0.543689.$$
This solution set, with identical $c_1 = c_2$, falls under the purview of a simplified calculation, in which it can be shown
$$  f_{n,n} \sim \rho^{-2n} \frac{1.5009481}{\sqrt{n}}.$$
(The true amplitude can be written in terms of $\rho$, and as many decimal places as wanted are possible, but I'm not going to bother to typeset it.)
For comparison:

*

*$n=40$, $\text{true}=3.4601\times 10^{20}$, $\text{asym} = 3.5261\times 10^{20}$ , $\text{absolute % err} = 1.91\%$.


*$n=200$, $\text{true}=7.6554\times 10^{104}$, $\text{asym} = 7.6847\times 10^{104}$ , $\text{absolute % err} = 0.38\%$.
For the non-diagonal case, you will be looking at $f_{an,bn}$.  The machinery should work, though it is more complicated. You'll get roots that depend on $(a,b)$ and have to solve a complicated determinant to get the amplitude, also dependent on $(a,b)$.  The question is, why does the non-diagonal case deviate rapidly from the diagonal?  Is there a way to understand this qualitatively?  I believe the answer lies in the modified form
$$ f_{an,bn} \sim c_1(a,b)^{-a n} c_2(a,b)^{-b n} \cdot \operatorname{amp}(a,b)/\sqrt{n}.$$
The amplitude will vary only like a polynomial upon changing $(a,b)$, but the first two factors have an exponential dependence.
