Tiling a Tilted Square

We will solve the problem in three steps, where each time we build up from the previous step and add layers of complexity.

Step 1: one quarter

First, consider the original problem, but with one quarter of the shape:

If we ask about tiling this shape using the procedure outlined in the problem statement, then we can write a recursion in the general case. To see why, start by picking one of the edge squares. This divides the shape into two similar but smaller shapes. For example, if there are $n=6$ edge squares and we pick the third one from the top, we obtain:

We are left with the cases $n=2$ and $n=3$, which we can further subdivide. If we call $c_n$ the number of tilings of this shape, we therefore have the recursion:
\[
c_0=1,\quad\text{and}\quad c_{n+1} = \sum_{k=0}^n c_k c_{n-k}\quad\text{for }n\geq 0
\]This is the well-known recurrence relation for the Catalan numbers. The first few Catalan numbers are (starting from $n=0$):
\[
\{c_n\} = \{1, 1, 2, 5, 14, 42, 132, 429, 1430,\dots\}
\]and a general formula is given by:
\[
c_n = \frac{1}{n+1}\binom{2n}{n}
\]

Step 2: one half

Now, consider the original problem, but with one half of the shape:

If we ask about tiling this shape, we can form a recurrence relation like we did in Step 1. This time, when we pick an edge square, we draw one large rectangle and we are left with a half-shape (up top) and identical quarter-shapes (on either side). For example, if there are $n=6$ edge squares and we pick the fourth from the top, we obtain:

If we call $d_n$ the number of tilings of this shape, we therefore have the recursion:
\[
d_0=1,\quad\text{and}\quad
d_{n+1} = \sum_{k=0}^n c_{n-k}^2 d_k\quad\text{for }n\geq 0
\]The first few numbers in this sequence are (starting from $n=0$):
\[
\{d_n\} = \{1, 1, 2, 7, 38, 274, 2350, 22531, 233292,\dots\}
\]

Step 3: the whole thing

Now consider the original problem. You can probably guess the pattern by now… By picking one of the edge squares, we subdivide the problem into four half-shapes in identical pairs (north-south and east-west). If we call $e_n$ the number of tilings of the whole shape, we therefore have the formula:
\[
e_0=1,\qquad\text{and}\quad
e_{n+1} = \sum_{k=0}^n d_{n-k}^2 d_k^2\quad\text{for }n\geq 0
\]The first few numbers in this sequence are (starting from $n=0$):
\[
\{e_n\} = \{1, 1, 2, 9, 106, 3002, 153432, 11209105, 1027079042\}
\]

Unfortunately, none of this is particularly satisfying since we do not have a closed-form solution for $e_n$. Let’s try to find one…

Attempt at a closed-form solution

A good place to start is to see how the closed-form solution for the Catalan numbers is derived. One way is to use generating functions. The idea is that we define an infinite polynomial (a power series) where the coefficients are the sequence we care about. For Catalan numbers, we have:
\[
C(x) = \sum_{n=0}^\infty c_n x^n
\]Now notice that the recurrence relation for Catalan numbers is a convolution, which we can obtain by squaring $C(x)$:
\begin{align*}
C(x)^2 &= \sum_{m=0}^\infty \sum_{k=0}^\infty c_m c_k x^{m+k} \\
&= \sum_{n=0}^\infty \left(\sum_{k=0}^n c_k c_{n-k} \right)x^n \\
&= \sum_{n=0}^\infty c_{n+1} x^n \\
&= \frac{1}{x}\left( C(x)-1\right)
\end{align*}Solving for $C(x)$, we obtain:
\[
C(x) = \frac{1-\sqrt{1-4x}}{2x}
\](The other root of the quadratic can be excluded since it does not satisfy $C(0)=c_0=1$) From here, we can perform a series expansion via the binomial theorem and extract the coefficient of $x^n$, which yields the formula $c_n = \frac{1}{n+1}\binom{2n}{n}$.

We can use a similar argument to obtain a generating function for $d_n$. To this effect, define:
\[
D(x) = \sum_{n=0}^\infty d_n x^n
\]Now, we can write:
\begin{align*}
D(x) &= 1 + \sum_{n=0}^\infty d_{n+1} x^{n+1} \\
&= 1 + \sum_{n=0}^\infty \sum_{k=0}^n c_{n-k}^2 d_k x^{n+1} \\
&= 1 + \sum_{k=0}^\infty \sum_{n=k}^\infty c_{n-k}^2 d_k x^{n+1} \\
&= 1 + \sum_{k=0}^\infty \sum_{n=0}^\infty c_{n}^2 d_k x^{n+k+1} \\
&= 1 + \sum_{k=0}^\infty d_k x^k \sum_{n=0}^\infty c_n^2 x^{n+1} \\
&= 1 + D(x)\sum_{n=0}^\infty c_n^2 x^{n+1}
\end{align*}Therefore, we can express the generating function for $d_n$ in terms of the generating function of squared Catalan numbers:
\[
\hat C(x) := \sum_{n=0}^\infty c_n^2 x^{n+1}\qquad\text{and}\qquad
D(x) = \frac{1}{1-\hat C(x)}
\]According to Mathematica, the series involving squared Catalan numbers can be evaluated in terms of a hypergeometric function. Namely:
\[
\hat C(x) =
\frac{1}{4} \bigl(\, _2F_1(-\tfrac{1}{2},-\tfrac{1}{2};1;16 x)-1\bigr)
\]Unfortunately, this isn’t particularly helpful as there does not appear to be any way to obtain a formula for the coefficient of $x^n$ in the series for $D(x)$.

Continuing in this fashion, we can also define $E(x)$ as the generating function for $e_n$, which we can express in terms of the square of the generating function for $d_n^2$. Namely:
\[
E(x) = \sum_{n=0}^\infty e_n x^n,\qquad
\hat D(x) = \sum_{n=0}^\infty d_n^2 x^{n+1},\qquad
E(x) = 1 + \tfrac{1}{x}\hat D(x)^2
\]But again, this doesn’t really seem helpful as we can’t evaluate $D(x)$ or much less $\hat D(x)$.

If anybody else can make progress on this problem I would love to hear your approach!

Asymptotics

It was pointed out by commenter MarkS that $c_n^2/d_n$ and $c_n^4/e_n$ appear to tend to finite limits as $n\to\infty$. Here is what we get when we plot these quantities up to $n=2000$:

Is there a way we can find the exact values of these limits? Maybe! We’ll make use of the following fact. If the sequences $a_n$ and $b_n$ have corresponding generating functions $A(x)$ and $B(x)$, and these have radius of convergence $R$, then we can write:
\[
\lim_{n\to\infty}\frac{a_n}{b_n} = \lim_{x\to R^-} \frac{A(x)}{B(x)}
\]This works so long as $A(x)$ and $B(x)$ go to $\infty$ as $x\to R^-$, because any finite truncation of the series will be dominated by its last term (largest power of $x$), so the ratio of truncated series just behaves like the ratio of its last terms, which is what we care about (apologies for the hand-waving; hopefully this makes sense!).

Applying this idea, we would like to write:
\[
\lim_{n\to\infty}\frac{c_n^2}{d_n} = \lim_{x\to R^-}\frac{\hat C(x)}{x D(x)} = \lim_{x\to R^-} x \hat C(x)\bigl(1-\hat C(x)\bigr)
\]But there is just one problem… We already established that
\[
\hat C(x) = \sum_{n=0}^\infty c_n^2 x^{n+1} =
\frac{1}{4} \bigl(\, _2F_1(-\tfrac{1}{2},-\tfrac{1}{2};1;16 x)-1\bigr)
\]and as it turns out, the series for this function has a radius of convergence of $R=\frac{1}{16}$. If you’re curious as to why, remember that this is a series where the coefficients are squared Catalan numbers. Catalan numbers have a well-known property that
\[
c_n \sim \frac{4^n}{n^{3/2}\pi},\qquad\text{i.e.,}\quad
\lim_{n\to\infty}\frac{c_n}{\left(\frac{4^n}{n^{3/2}\pi}\right)} = 1
\]Therefore, $c_n^2 \sim 16^n/n^3\pi$, so a necessary condition for $\hat C(x)$ to converge is that $|x|\leq \frac{1}{16}$ (the general term of the series must tend to zero).

So what’s the problem? Here is a 3D plot of $\hat C(z)$ for complex $z$ (vertical axis is magnitude, color-coded by argument) and plotted for $|z|\leq \frac{1}{16}$.

As we can see, nothing is going to infinity near the boundary. How is it possible? It turns out the radius of convergence is $\frac{1}{16}$ because there is a different kind of discontinuity at this point; the argument (rather than the magnitude) becomes discontinuous. You can see it more clearly when the plot is zoomed out to $|z|\leq 1$ (notice the color discontinuity along the positive real axis).

So how do we deal with this problem? We need to transform the series so that it goes to infinity as $x\to\frac{1}{16}$. One way to do this is to take derivatives. This will yield a series with terms like $n a_n x^{n-1}$ and $n b_n x^{n-1}$, and the ratio remains unchanged! Here is a plot of the first derivative, $\hat C'(z)$:

This is better, as now the function appears non-differentiable near $z=\frac{1}{16}$, but it still doesn’t go to infinity. Let’s differentiate again! Here is a plot of $\hat C'{}'(z)$:

Now we’re in business. So the limit we’re looking for is:
\[
\lim_{n\to\infty}\frac{c_n^2}{d_n}
= \lim_{x\to\frac{1}{16}} \frac{\frac{d^2}{dx^2}\Bigl(\frac{1}{x}\hat C(x)\Bigr)}{\frac{d^2}{dx^2}\Bigl(\frac{1}{1-\hat C(x)}\Bigr)}=\left(5-\frac{4}{\pi}\right)^2\approx 13.88874349
\]I used Mathematica in the last step, which can analytically differentiate and evaluate hypergeometric functions. In conclusion, we have the asymptotic formula:

$\displaystyle
d_n \sim \frac{c_n^2}{\left( 5-\tfrac{4}{\pi}\right)^2}
\sim \frac{16^n}{\pi\left( 5-\tfrac{4}{\pi}\right)^2 n^3}
$

Now as for $e_n$, I’m stuck again. In principle we could use the same technique, which would require evaluating:
\[
\lim_{n\to\infty}\frac{c_n^4}{e_n}
= \lim_{x\to R^-} \frac{\frac{d^k}{dx^k}\sum_{n=0}^\infty c_n^4 x^n}{\frac{d^k}{dx^k}E(x)}
= \lim_{x\to R^-} \frac{\frac{d^k}{dx^k}\sum_{n=0}^\infty c_n^4 x^n}{\frac{d^k}{dx^k}\left(1+\frac{1}{x}\hat D(x)^2\right)}
\]where we differentiate however many times necessary to ensure the series diverges as $x\to R^-$. Mathematica is able to evaluate the numerator (it’s hypergeometric functions again, but a different kind this time), so I was able to determine that $R = \frac{1}{256}$ and the smallest $k$ we can use is $5$ (yikes!). But the sticking point is the denominator. Although we found an asymptotic expansion for $d_n$, we don’t have an actual formula. Without this, we can’t evaluate the general term of $\hat D(x)^2$, which depends on all the $d_n$.

Again, if anybody has ideas, let me know in the comments!