Number of certain (0,1)-matrices, Stanley's Enumerative Combinatorics

Stanley’s Enumerative Combinatorics (http://www-math.mit.edu/~rstan/ec/ec1.pdf) contains next fact:
1.1.3 Example. Let f(n) be the number of n × n matrices M of $0$’s and $1$’s such that every
row and column of M has three 1’s. For example, f(0) = 1, f(1) = f(2) = 0, f(3) = 1. The
most explicit formula known at present for f(n) is
$$f(n)=6^{-n}{(n!)}^2\sum\frac{(-1)^{\beta}(\beta+3\gamma)!2^\alpha 3^{\beta}}{\alpha!\beta!\gamma!^26^{\gamma}}$$
where the sum ranges over all (n + 2)(n + 1)/2 solutions to α + β + γ = n in nonnegative
integers.

I need proof of this fact. (i.e. reference to book or articles that contains proof this fact).

Solutions Collecting From Web of "Number of certain (0,1)-matrices, Stanley's Enumerative Combinatorics"

You want the coefficient of $x_1^3\cdots x_n^3y_1^3\cdots y_n^3$ in $\prod_{i=1}^n \prod_{j=1}^n (1+x_iy_j)$. Expand the logarithm, delete all terms with an exponent greater than 3, apply the exponential function exp, and simplify. For a generalization, see the paper by Musiker-Odama at https://www.mtholyoke.edu/courses/gcobb/REU_MCMC/papers.html.

For $n=4$, there are $4$ options for placing the single $0$ in the first row, then $3$ in the second row, $2$ in the third row and $1$ in the fourth row, for a total of $24$. Searching for $1,0,0,1,24$ at OEIS yields OEIS sequence A001501; the entry contains lots of references.

This problem was first solved by Ron Read in this 1958 PhD thesis. I didn’t check if exactly this formula appears, but there are several that look very similar. One combinatorial way to prove it is like this:

Let $A_1,\ldots,A_n$ and $B_1,\ldots,B_n$ be disjoint sets of size 3. Consider bijections from $A_1\cup\cdots\cup A_n$ to $B_1\cup\cdots\cup B_n$. There are $(3n)!$ in total. Use inclusion-exclusion to find the number such that each $A_i$ maps to 3 different $B_j$s. (There are two things that can go wrong: mapping two elements of $A_i$ to two elements of $B_j$, and mapping all of $A_i$ to all of $B_j$; the counting is easy since the involved $B_j$s are distinct.) Then make an $n\times n$ matrix whose $(i,j)$ entry is the number of elements of $A_i$ which map to $B_j$. Finally, note that each matrix comes from $(3!)^{2n}$ bijections.

An answer can also be found in the section 6.3 of the book “L. Comtet, Advanced Combinatorics, Springer 1974” and asymptotic formula $f(n)\sim e^{-2}36^{-n}(3n)!$ can be derived from a general asymptotic formula in the paper “C. J. Everett, P. R. Stein, The asymptotic number of integer stochastic matrices, Disc. Math. Vol. 1, No. 1 (1972) 55-72”.