Continuing from our chat in $\hslash$, I have to concur with the other participants that you are conflating very many things together, and that that is the real source of all your confusions regarding this topic.
In fact, even your title is conflating stuff and thus making your confusion apparent. However, it would go too far afield in linguistics to disentangle that.
I do, however, like your own answer to this, as it shows that you have learnt some nice little things.
It is clear that your end goal here is to learn a bit about the variational principle in QM as it pertains to finding approximate eigenfunctions of a system.
However, the basis expansion that you have tried to write down as your first equation, is already a massive impediment to the understanding of what is being done. I am perfectly aware that in your quantum chemistry work, STO and GTO and so forth are always used in the way you are talking about, and they definitely do work. My point, however, is that if you focus too much on these trees, you will miss the forest. You should take a step back and contemplate the general situation.
The kind of trial functions that you should be thinking about, should look like
$$\tag{1$^\prime$}\langle x\vert\psi^\prime_{(N,\{a_i\})}\rangle=N(a_1x^2+1)e^{-a_2x^2-a_3x^4}$$
The reason why this is helpful (and somewhat standard in physics textbooks covering the variational principle in QM), is that one is under no delusions that this is being expanded in any form of basis. We are making a general argument over any random kind of attempted trial function, and it does not have to span the Hilbert space of possible quantum states, and obviously if it is not having enough parameters, it would then be unable to actually get us to a Hamiltonian eigenstate, no matter how good we choose the parameters $\{a_i\}$.
Now, you might think that by assuming the wavefunction to be normalised, you would be simplifying. However, in the context of perturbations and approximations, this is rarely helpful, and it applies here too. I have separated out the parameter $N$, and it is obvious that it is meant to be the normalising constant, found by enforcing $\langle\psi^\prime_{(N,\{a_i\})}\vert\psi^\prime_{(N,\{a_i\})}\rangle=1$
However, this forces your optimisation procedure to be really convoluted, unnaturally breaking this into two steps, and increasing the number of free parameters by one. Instead, it is far nicer to omit this normalisation constant and consider the form of the variational principle that looks like
$$\tag2\frac{\langle\psi_{\{a_i^b\}}\vert\hat{\mathcal H}\vert\psi_{\{a_i^b\}}\rangle}{\langle\psi_{\{a_i^b\}}\vert\psi_{\{a_i^b\}}\rangle}\geqslant E_0$$
Note that this is a complicated situation. The theorem is true for any set of $\{a_i\}$, but we are only interested in the lowest expected energy. That is, considering the unnormalised trial wavefunctions
$$\tag1\langle x\vert\psi_{\{a_i\}}\rangle=(a_1x^2+1)e^{-a_2x^2-a_3x^4}$$
we are really interested in the stationary values
$$\tag3\forall\hat{\mathcal H}\exists\{a_i^b\}\in\mathbb C^n\ |\qquad\left.\frac{\partial\ }{\partial a_i}\frac{\langle\psi_{\{a_i\}}\vert\hat{\mathcal H}\vert\psi_{\{a_i\}}\rangle}{\langle\psi_{\{a_i\}}\vert\psi_{\{a_i\}}\rangle}\right|_{\{a_i\}=\{a_i^b\}}=0$$
these conditions net us a few sets of $\{a_i^b\}$ whereby the $E_t$ in
$$\tag4E_t\overset{\text{def}}=\frac{\langle\psi_{\{a_i^b\}}\vert\hat{\mathcal H}\vert\psi_{\{a_i^b\}}\rangle}{\langle\psi_{\{a_i^b\}}\vert\psi_{\{a_i^b\}}\rangle}$$
are stationary with respect to minor variations in $\{a_i\}$, and we really only care for the smallest value of $E_t$ and its trial wavefunction, i.e. its set of $\{a_i^b\}$.
Let us take another step back. Given any basis $\left|b_k\right>$, we can express any wavefunction in terms of the basis. i.e.
$$\begin{align}
\tag5\forall\left|\psi\right>\in\mathcal H\exists{c_k}\in\mathbb C^N\ |\qquad\left|\psi\right>&=\sum_{k=1}^N\left|b_k\right>c_k\\
\tag6&=\sum_{k=1}^N\left|b_k\right>\left<b_k\vert\psi\right>
\end {align}$$
where the equality of $c_k=\left<b_k\vert\psi\right>$ in Equation (6) is only true if $\left|b_k\right>$ is an orthonormal basis. If you are dealing with single atoms, the atomic orbitals can have this be true; however, it is much more common in molecules and condensed matter that the kind of STO and GTO that you are using, fail to be orthonormal, and in fact are usually overcomplete. Then it would not be deterministic, but we can still work with them. It represents but mild complications, and not major revisions to the conceptualisation.
The consequences of having a basis at all, i.e. Equation (5), means that, even without optimising, any trial function will have $c_k=c_k(\{a_i\})$
$$\tag7\forall\{a_i\}\in\mathbb C^n\exists\{c_k\}\in\mathbb C^N\ |\qquad\left|\psi_{\{a_i\}}\right>=\sum_{k=1}^N\left|b_k\right>c_k$$
It does not matter that, in the overcomplete case, there will be more than one set of $\{c_k\}$ that satisfies this, only that we can always find some set that works.
Now we specialise to the situation considered in the theorem of the variational principle. As you noted in your answer, the spectral theorem guarantees that $\hat{\mathcal H}$ has a complete orthonormal eigenbasis $\left|E_k\right>$ with ground state energy $E_0$ being its minimum energy. This means that we can substitute the version of Equation (7) that looks like Equation (6) into Equation (4) combined with Equation (2) to get
$$\tag8E_{t,\text{min}}=\frac{\sum_k\langle\psi_{\{a_i^b\}}\vert E_k\rangle E_k\langle E_k\vert\psi_{\{a_i^b\}}\rangle}{\sum_k\langle\psi_{\{a_i^b\}}\vert E_k\rangle\langle E_k\vert\psi_{\{a_i^b\}}\rangle}=E_0+\frac{\sum_{k\neq0}\langle\psi_{\{a_i^b\}}\vert E_k\rangle(E_k-E_0)\langle E_k\vert\psi_{\{a_i^b\}}\rangle}{\sum_k\langle\psi_{\{a_i^b\}}\vert E_k\rangle\langle E_k\vert\psi_{\{a_i^b\}}\rangle}\geqslant E_0$$
where the last sum is necessarily positive by definition that $E_0$ is the minimal energy; all the subtraction $E_k-E_0>0$ You might find Equation (8) a little easier to understand in the familiar notation
$$\tag{8$^\prime$}\frac{\sum_k E_k|c_k|^2}{\sum_k|c_k|^2}=E_0+\frac{\sum_{k\neq0}(E_k-E_0)|c_k|^2}{\sum_k|c_k|^2}\geqslant E_0$$
And how do we know that this combinations won't actually (for some unknown reason) be actually below the system's ground energy ?
Note that we have not assumed that the trial wavefunction looks like any basis in this derivation, or be expanded in any basis. It is clear that this is just the argument that the average of stuff that are necessarily bigger than some minimum, can thus only be bigger than or equal to said minimum. There is thus no possibility of accidentally getting the trial wavefunction to have energy lower than the system's ground energy $E_0$. i.e. $E_0\leqslant E_{t,\text{min}}$ is a strict inequality proved by this theorem to hold in general for any kind of variation.
When we do it in the real world, using atomic orbitals, or just any suitable set of basis functions,
Do the basis functions that we use (especially in linear variations) need to be eigenfunctions of the Hamiltonian of the system?
Manifestly not. In the case of trying atomic orbitals, what you are doing is choosing for $a_i$ to now just be $c_k$ of some basis of orbitals. In Equation (8), I have chosen painstakingly to stay generic, and thus it does not matter that your $\lvert\psi_{\{a_i\}}\rangle$ are now themselves expanded in some random basis, you just theoretically apply the unknown $\left<E_k\right|$ to them in accordance with Equation (8) and obtain the same result.
