I assume that in the question it is intended to ask for the expectation value $\mathbb{E}\{\det(AA^T)\}$, not $\mathbb{E}(AA^T)$.
I have no good idea yet how to attack the general Toeplitz case, but for the special case of square circulant matrices one can use the fact that their eigenvalues can explicitly be expressed in terms of the matrix elements by a discrete Fourier transform. Below I will show that in this case the answer to the question is yes, and I will also present some evidence for the validity stronger statements concerning the growth rate of
Denoting as $a_i$ (with $i=0,...,n-1$) the elements in the first row of a $n \times n$ circulant matrix $A$, the eigenvalues $\lambda_k$ of $A$ are given by:
$$
\lambda_k = \sum_{i=0}^{n-1} a_i\, \omega^{ik} \quad (k=0,...,n-1) \ ,
\hspace{2cm} (1)
$$
with the primitive $n^{\rm th}$ root of unity $\omega={\rm e}^{\frac{2\pi{\rm i}}n}$.
The determinant of $A$ can then be expressed as
$$
\det A = \prod_{k=0}^{n-1} \left( \sum_{i=0}^{n-1} a_i\, \omega^{ik} \right) \ ,
\hspace{2cm} (2)
$$
which I will use as starting point for this (partial) answer to your question.
For later convenience, let us make a variable change to new random variables $\varepsilon_i = 2a_i-1 \in \{-1,1\}$. Using that $\sum_{i=0}^{n-1} \omega^{ik} = n\, \delta_{k0}$, and splitting off the $k=0$ factor in $(2)$, we obtain:
$$
\det A = 2^{-n} \left( n + \sum_{i=0}^{n-1} \varepsilon_i \right) \prod_{k=1}^{n-1} \left( \sum_{i=0}^{n-1} \varepsilon_i\, \omega^{ik} \right) \ .
\hspace{2cm} (3)
$$
Since $A$ is a real matrix and thus $\det(AA^{\rm T}) = (\det A)^2$, the moments of the $\det(AA^{\rm T})$ distribution are equal to the even moments of the $\det A$ distribution. To calculate these moments, powers of $(3)$ have to be averaged over the i.u.d. variables $\varepsilon_i \in \{-1,1\}$. This may be done in a convenient way with help of a generating function $F$ defined as follows:
\begin{align}
F(x_0,...,x_{n-1})
&= \mathbb{E}\left\{ \exp \sum_{i,k=0}^{n-1} x_k\, \varepsilon_i\, \omega^{ik} \right\}
= \mathbb{E}\left\{ \prod_{i=0}^{n-1} \exp\left( \varepsilon_i \sum_{k=0}^{n-1} x_k\, \omega^{ik} \right) \right\} \\
&= \prod_{i=0}^{n-1} \cosh\left( \sum_{k=0}^{n-1} x_k\, \omega^{ik} \right) \ .
\hspace{2cm} (4)
\end{align}
Let us describe the $p^{\rm th}$ moment of the $\det A$ distribution, for our $n \times n$ circulant matrices, by $N_{\rm cir}(p,n) = 2^{pn} {\big\langle} (\det A)^p {\big\rangle}$. With $(3)$ and $(4)$ we can express this quantity as follows by derivatives of the generating function (using the shorthand notation $\partial_k = \frac\partial{\partial x_k}$):
$$
N_{\rm cir}(p,n)
= \left[ \left( {\big(} n + \partial_0 {\big)} \prod_{k=1}^{n-1} \partial_k \right)^p F \right]_{x_0=x_1=...=x_{n-1}=0} \ .
\hspace{1cm} (5)
$$
In order to further evaluate the right hand side of $(5)$, we first note that $(4)$ implies that any multiple derivative of $F$ with respect to the $x_k$ can, according to the product rule, be written as a sum over many terms, each of which is a product of (i) exactly $n$ factors of either $\cosh_i = $ $\cosh\left( \sum_{k=0}^{n-1} x_k\, \omega^{ik} \right)$ or $\sinh_i =$ $\sinh\left( \sum_{k=0}^{n-1} x_k\, \omega^{ik} \right)$ (one for each index $i=1,...,n$), and moreover (ii) some factors $\omega^{ik}$ (whose number is equal to the order of the derivative), each of which may have different values of $i$ and $k$. Each new derivative $\partial_k$ on one of these terms "pulls out" one new factor $\omega^{ik}$ together with a summation over $i$, and at the same time converts one of the $\cosh_i$ factors into a $\sinh_i$ or vice versa. In the end, terms having any $\sinh_i$ factors left will give a vanishing contribution this sum, since $\sinh(0)=0$; so in particular, the entire sum can be nonvanishing only if the total number of derivatives is even. This is certainly the case for any even $p$, and in particular for $p=2$ which is the case of interest here.
Furthermore, carrying out the summations over $i$ leads to certain constraints between the different $k$ values in each term, since for instance $\sum_{i=0}^{n-1} \omega^{ik} \omega^{il} = n\delta_{kl}$ etc.. This again eliminates a large number of terms, and in the end only a certain narrow class of terms remains, each of which is real and contributes with a combinatorial factor and one factor of $n$ for each remaining $i$ sum. The combinatorial factors can be estimated using a diagrammatic approach. (I could explain it in some detail later, if there is enough interest and if I find the time). For $p=2$ this procedure leads to
$$
N_{\rm cir-}(2,n) < N_{\rm cir}(2,n) < N_{\rm cir+}(2,n) \ ,
\hspace{2cm} (6)
$$
with the following explicit bounds (here $\lfloor \frac{n-1}2 \rfloor$ denotes the integer part of $\frac{n-1}2$):
\begin{align}
N_{\rm cir-}(2,n) &= 2^{\lfloor \frac{n-1}2 \rfloor} (n+1)! \ \ , \\
N_{\rm cir+}(2,n) &= 2^{\lfloor \frac{n-1}2 \rfloor} (n+1) \, n^n \ .
\hspace{2cm} (7)
\end{align}
For general $\{0,1\}$-matrices we have, as mentionend in the question,
$$
N_{\rm gen}(2,n) = (n+1)! \ ,
\hspace{2cm} (8)
$$
so that $(6)$ and $(7)$ in particular imply that always $N_{\rm cir}(2,n) > N_{\rm gen}(2,n)$. This answers the question for circulant matrices.
The results above can be further visualized by considering the following expression for what we may call the "asymptotic prefactor of $n$ in the exponent", or a little sloppy "relative growth rate":
$$
f_{\alpha,c}(n) = \ln N_\alpha(2,n)/n - \ln(n{+}c) \ ,
$$
where $\alpha \in \{{\rm cir{-}}, {\rm cir{+}}, {\rm cir},{\rm gen}\}$, and $c$ is a positive constant of our choice (which won't affect the limit $n \to \infty$ of $f_{\alpha,c}(n)$ but can be adjusted to make the curves look nice). For the three expressions in $(7)$ and $(8)$,
$$
\lim_{n \to \infty} f_{\alpha,c}(n) =
\begin{cases} \frac1 2 \ln2 -1 \approx -0.65, \quad \alpha = {\rm cir{-}} \\ \quad \frac1 2 \ln2 \approx 0.35, \quad \alpha = {\rm cir{+}} \\ \qquad -1 , \qquad \alpha = {\rm gen} \end{cases} \ .
\hspace{2cm} (9)
$$
(Note that the Hadamard absolute upper limit (for arbitrary $n \times n$ matrices $A$ with all elements lying in the closed unit disk) corresponds to a $n \to \infty$ limit of $2\ln2 \approx 1.39$).
A plot of numerical results for $f_{\alpha,c}(n)$ (with $c=4$ and $n = 1,...,20$) is shown below. It suggests that also $\lim_{n \to \infty} f_{{\rm cir},c}(n)$ exists, and thet it has a value ${\approx}{-0.3}$ (red curve).
The following questions then arise naturally:
- Can the existence or non-existence of $\lim_{n \to \infty} f_{{\rm cir},c}(n)$ be proved?
- If the limit exists, can its value be calculated, or can better bounds be provided than those in the first two lines of $(9)$?
