Expected value of determinant of simple infinite random matrix

Very nice problem!

Let me recall you that the determinant of $n \times n$ matrices with entries in $\{0,1\}$ is related to the one of $n+1 \times n+1$ matrices with entries in $\{-1,+1\}$: replace the zeros by $-1$'s and add a row of $-1$'s on the top and a column of ones on the right (you may want to read this arXiv).

I can now tell you that, besides the trivial cases $p=0$ and $p=1,$ your problem is solved for $p=\frac{1}{2}$.

Indeed it was studied by T. Tao and V. Vu in arXiv. They proved that for the matrix $M_n$ of size $n \times n$ where the entries are i.i.d. Bernoulli random variables in $\{-1,+1\},$ the probability that $|\det(M_n)|$ is close to $\sqrt{n!}$ tends to one:

$$P \left(|\det (M_n)|\geq \sqrt{n!}\exp(-29n^{1/2}\ln^{1/2}n)\right)=1-o(1).$$

I agree with user39115!

I will give a heuristic from random matrix theory because we know the global behaviour of the eigenvalue. First $$A=p 1 +\sqrt{N(p-p^2)}\frac{B}{\sqrt{N}} $$ where $1$ is the matrix with only 1 entries and $$B_{i,j}=\begin{cases} \frac{1-p}{\sqrt{(p-p^2)}} &\text{with probability } p\\ \frac{-p}{\sqrt{(p-p^2)}} &\text{ with probability }1-p \end{cases} $$ $B$ is a random matrices independent entries with zero mean and variance 1 so we are exactly in set up of generalise Wigner Matrices https://terrytao.files.wordpress.com/2011/02/matrix-book.pdf. The eigenvalue of $B/\sqrt{N}$ converge to the Uniform law on the unit circle. $$\mu_n = \frac{1}{N}\sum_{\lambda\in \sigma(B/\sqrt{N})}\delta_\lambda \rightarrow \frac{1}{\pi}\mathcal{U}_{x^2+y^2\leq 1}$$ $A$ is just a rank one perturbation of $B$ so it does not change the global distribution of the eigenvalue of $B$ only add the larger eigenvalue $\lambda_1\sim pN$. $$\log(|det(A)|)\sim \log(pN)+N\log(\sqrt{N(p-p^2)})+\log(|det(B/\sqrt{N})|)$$ and we believe that $$\frac{1}{N}\log(|det(B/\sqrt{N})|=\frac{1}{N}\sum_{\lambda\in\sigma(B/\sqrt{N})} \log(|\lambda|)\sim \frac{1}{\pi}\int_{x^2+y^2\leq 1}\log(\sqrt{x^2+y^2})dxdy$$ and the last integral is equal to $-1/2$.Therefore we should get $$\log(|det(A)|)\sim N\log(\sqrt{\frac{N(p-p^2)}{e}}) $$ note that we recover $\sqrt{N!}$ of user39115.

So this is my heuristic (which rigourous proof is probably extremly hard.)

Here is a solution to Hipstpaka's question about $\det(A)^2$. I don't have enough space to answer in a comment. I don't know where a statement appears in the literature, but the proof uses a standard technique discussed for instance in Enumerative Combinatorics, vol. 2, Exercise 5.64.

Write $\varepsilon_w$ for the sign of the permutation $w\in S_n$. Then \begin{eqnarray*} \sum_A \det(A)^2 & = & \sum_{i,j=1}^n\sum_{a_{ij}=0,1} \left(\sum_{w\in S_n} \varepsilon_w a_{1,w(1)}\cdots a_{n,w(n)}\right)^2\\ & = & \sum_{i,j=1}^n\sum_{a_{ij}=0,1} \sum_{u,v\in S_n} \varepsilon_u\varepsilon_v a_{1,u(1)}\cdots a_{n,u(n)} a_{1,v(1)}\cdots a_{n,v(n)}. \end{eqnarray*} Let $f(u,v)$ be the number of distinct variables occurring among $a_{1,u(1)},\dots, a_{n,u(n)},a_{1,v(1)},\dots, a_{n,v(n)}$. The product $a_{1,u(1)}\cdots a_{n,u(n)}a_{1,v(1)}\cdots a_{n,v(n)}$ is 1 with probability $p^{f(u,v)}$ and is otherwise 0. Moreover, $f(u,v)= 2n-\mathrm{fix}(uv^{-1})$, where $\mathrm{fix}(uv^{-1})$ denotes the number of fixed points of $uv^{-1}$. If $E$ denotes expectation, then we get $$ E(\det(A)^2) = \sum_{u,v\in S_n}\varepsilon_u\varepsilon_v p^{2n-\mathrm{fix}(uv^{-1})}. $$ Setting $w=uv^{-1}$ and noting that $\varepsilon_u\varepsilon_{wu^{-1}} = \varepsilon_w$, we get \begin{eqnarray*} E(\det(A)^2) & = & \sum_{w\in S_n} p^{2n-\mathrm{fix(w)}} \sum_u\varepsilon_u\varepsilon_{wu^{-1}}\\ & = & n!p^{2n}\sum_{w\in S_n}p^{-\mathrm{fix(w)}}\varepsilon_w. \end{eqnarray*} Let $g(n)=\sum_{w\in S_n}p^{-\mathrm{fix(w)}}\varepsilon_w$. By standard generating function techniques (Enumerative Combinatorics, vol. 2, Section 5.2) we have \begin{eqnarray*} \sum_{n\geq 0} g(n)\frac{x^n}{n!} & = & \exp\left( \frac 1p x-\frac{x^2}{2}+\frac{x^3}{3} -\frac{x^4}{4}+\cdots\right)\\ & = & (1+x)\exp \left( \frac 1p-1\right)x. \end{eqnarray*} It is now routine to extract the coefficient of $x^n$ and then compute $$ E(\det(A)^2)= n!\,p^n(p-1)^{n-1}(1+(n-1)p). $$

In general we don't have $E(\det(A)^2)=E(|\det(A)|)^2$, so this does not directly answer the original question.

For a related question (and answer) see Expected determinant of a random NxN matrix.