I can do your warm up, but the more general problem will be harder because you can't use the same state space reduction idea.
Write for our state space the number of dimensions spanned so far, $\{0,1,2,3,\ldots,n\}$. One we have $k$ dimensions spanned, $2^k$ of the $2^n$ possible vectors will not give us any new dimensions, while $2^n-2^k$ of them will give us one new dimension. Therefore, on our state space evolution is according to the following stochastic matrix
$$P=\begin{pmatrix}
\frac{1}{2^n} & 1- \frac{1}{2^n} \\
& \frac{2}{2^n} & 1- \frac{2}{2^n} \\
&& \frac{4}{2^n} & 1- \frac{4}{2^n} \\
&&& \ddots\\
&&&& \frac{2^{n-1}}{2^n} & 1- \frac{2^{n-1}}{2^n} \\
&&&&& 1\\
\end{pmatrix}.$$
This defines a discrete discrete phase-type distribution, which we partition the matrix as
$$ P = \begin{pmatrix}
T & T^0 \\
0 & 1
\end{pmatrix}$$
and can compute moments using results in Lecture notes on phase-type distributions for 02407 Stochastic Processes by Bo Friis Nielsen (page 9) to be
$$\begin{align} \mathbb E(X) &= \begin{pmatrix}1 & 0 & \cdots & 0\end{pmatrix} (I-T)^{-1} \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix}\\
&= \begin{pmatrix}1 & 0 & \cdots & 0\end{pmatrix} \begin{pmatrix}
\frac{2^n}{2^n-1} & \frac{2^{n-1}}{2^{n-1}-1} & \cdots\\
0 & \frac{2^{n-1}}{2^{n-1}-1} & \cdots\\
0 & 0 & \cdots\\
\end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix}\\
&=\sum_{i=1}^n \frac{2^i}{2^i-1} =n+\sum_{i=1}^n \frac{1}{2^i-1}
\end{align}
$$
which converges to $n+E$, as you state with $E$ the Erdős–Borwein constant.
Using the formula in the notes for second factorial moments we can compute
$$\begin{align}\mathbb E(X(X-1)) &= 2! \begin{pmatrix}1 & 0 & \cdots & 0\end{pmatrix} (I-T)^{-2} T \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix}\\
&=\frac{2^n}{(2^n-1)^2} + \sum_{i=1}^{n-2} \left( \left(1-\frac{2^i}{2^n}\right) \left(\sum_{j=0}^j \frac{2^{n-j}}{2^{n-j}-1} \right) \frac{2^{n-i}}{2^{n-i}-1} + \frac{2^{i+1}}{2^n}\left( \sum_{j=0}^{i+1} \frac{2^{n-j}}{2^{n-j}-1}\right) \frac{2^{n-i-1}}{2^{n-i-1}-1} \right)
\end{align}$$
from which can can compute the variance. Numerically we see the convergence you describe:
n Var(X)
1 2.00000
2 2.44444
3 2.60771
4 2.67882
5 2.71212
6 2.72824
7 2.73618
8 2.74012
9 2.74208
10 2.74306
11 2.74355
12 2.74379
13 2.74391
14 2.74397
15 2.74400
16 2.74402
17 2.74403
18 2.74403
19 2.74403
20 2.74403