Context:
The present question is similar to [1]. But whereas in [1] there is a single constant, $p$, raised to a power of $x$, here we have the product of $(p)^{k}\,(1-p)^{n - k}$. I have looked in books of special functions to generalize the solution found in [1] for said purposes. I have sought a solution in terms of (i) hypergeometric functions or (ii) Legendre polynomials. I have not found any generalization. Thus, my path towards a solution is halted.
Is this a tractable problem?
What is a fruitful approach?
Introductory Formalism:
Allow a binomial distribution as in [2].
Provide that its support is
${\displaystyle k\in \{0,\dots ,n\}}$, and
$p ∈ [0,1]$; and
the probability mass function, g, is
$g\left(k\right) = {\displaystyle {\frac {n!}{k! (1-k)!}} \, p^{k}\,(1-p)^{1-k}} $.
Writing the expected value as $E(\cdot)$, we have
$E(a) = a$, for constant $a$;
$E(k) = n\,p$; and
$E(a^{k}) = (1 - p + a\,p)^n$; for constant $a$ (as per [3]).
Question:
What is the expected value of the probability mass function:
$E(g) = \sum_{k=0}^{n} \textstyle \left[{n \choose k}\,p^{k}(1-p)^{n-k}\right]^2 =$ ?
Bibliography:
[1] Squared binomial coefficient
[2] https://en.wikipedia.org/wiki/Binomial_distribution
[3] Covariance of product of two functions of two binomial distributions