Let $G$ be a finite abelian group, $n$ a positive integer and let $G^n$ denote the direct product of $n$ copies of $G$. We say an element of $G^n$ is *full* if it acts as a nonidentity element of $G$ in each of the factors of $G^n$.


Now consider the following random process. Sample a full group element $(g\_1,g\_2,\ldots,g\_n)$ uniformly at random from $G^n$. Now generate a subgroup $H$ of $G^n$ consisting of all elements of the form $\bigl(g\_1^{a\_1},\ldots,g\_n^{a\_n}\bigr)$ for integers $a\_i$. If we do this $k$ times (sampling with or without replacement, more on this below), then we can let $H\_j$ denote the subgroup generated on the $j$th iteration. Now we take the union of these subgroups and define 
$$ N\_k = \left| \bigcup\_{j=1}^k H\_j \right| \ ,$$
where $|\cdot|$ denotes the cardinality of the set. Finally, let $\mu\_k = \mathbb{E}(N\_k)$.


When $G$ is also a simple group, this is easy to calculate. So let's specialize to the simplest nontrivial case, $G = \mathbb{Z}/2 \times \mathbb{Z}/2$. (The other simplest case is $\mathbb{Z}/4$, but I'm not as interested in that one.) Then my question is,



> 
> How does $\mu\_k$ grow as a function of $k$? 
> 
> 
> 


I am principally interested in lower bounds on $\mu\_k$ in the case where the sampling is done uniformly *without* replacement on the set of full elements of $G^n$. Clearly sampling *with* replacement gives a lower bound, and it's much easier to work with. If you can say something about the variance of $N\_k$ too, that would be outstanding.



Generating a group by randomly sampling generators
Let's get a lower bound by considering sampling *with* replacement, for which we can even get the exact answer. If we sample $k$ full elements, what is the probability that they agree at exactly $l$ locations? For $G = \mathbb{Z}/2 \times \mathbb{Z}/2$, the answer is
$$
p(k,l) = {n \choose l} \left(\left(\frac{1}{3}\right)^{k-1}\right)^l \left(1-\left(\frac{1}{3}\right)^{k-1}\right)^{n-l}.
$$
When two full elements agree in exactly $l$ places, they jointly generate $2\times2^n - 2^l$ elements of $G^n$. Continuing this logic, an inclusion-exclusion principle lets us count the expected number of elements after $k$ steps, namely,
$$
\mu\_k = \sum\_{j=1}^k {n \choose j} (-1)^{j-1} \sum\_{l=0}^n p(j,l) 2^l .
$$
This formula isn't very nice because when $k$ is large we have to sum over as many as $3^n$ terms. Fortunately, we can use the binomial theorem twice to get rid of the sum on $j$. We first expand the binomial term in $p(j,l)$, then collapse the sum on $j$. The result is
$$
\mu\_k = 4^n - \sum\_{l=0}^n \sum\_{m=0}^{n-l} {n \choose l} {n-l \choose m} (-3)^m 6^l \bigl(1-3^{-(l+m)}\bigr)^k .
$$
This is much easier to evaluate than the previous formula. 


While this bound is really good for small $k$, it is not good for very large $k$. Since $\mu\_k$ is concave in $k$ (see the comments in Will Sawin's answer), the convex hull of two lower bounds is again a lower bound. So we can add in the fact that the point $(3^n,4^n)$ is the answer for sampling without replacement and find the line which connects this point back to the previous bound in an optimal way. This is hard to do analytically, but can be done for specific values of $n$. This gives a very good lower bound. 


