My intuition tells me that the ground state energy of bosons should
always be lower than the ground state of fermions -- no matter what
kind of interactions or other external properties we've chosen.
I think your intuition is usually correct; but it's possible to define systems where the ground state fermions would have less energy than the ground state bosons. First, a note on why fermions tend to have higher energies, than a note on how one can arrange for the bosons to have more energy than the fermions.
(1) Assume no particle interactions.
Under the assumption that the Hamiltonian for the boson and fermion are identical, the energies for single particles (i.e. $N=1$) will also be identical. This follows from the fact that the Schroedinger wave equation applies equally to a boson or fermion. In particular, under this assumption, the ground state energies are identical, call that energy $E_1$.
The simplest assumption for particle interactions is that there is none. In this case, the ground state energy for the boson is easy, it's just $N\; E_1$ as all the bosons are in the same state.
The ground state for the fermion will be a higher energy (due to the Pauli exclusion principle) except in the case that the ground state is $N$-fold degenerate in which case the bosons and fermions will have the same energy.
To some people, the above might be obvious in itself. For others, they might want a little less hand waving and a little more mathematics. So let the $N$ lowest energy eigenstates be $\psi_n(x)$ for $n=1,2,3,...,N$ with no degeneracy so that $H\psi_n = E_n\psi_n$. In the boson case the ground state has all the particles in this state so the wave function is a symmetrization of $\psi(x_1,x_2,...,x_n) = \psi_1(x_1)\psi_1(x_2)...\psi_1(x_n)$. But no matter how the positions are permuted, the energy of this state is $E_1+E_1+...+E_1 = N\;E_1$. Similarly, the energy for any permutation of the fermi ground state $\psi_1(x_1)\psi_2(x_2)...\psi_n(x_n)$ is $E_1+E_2+...E_N > N\;E_1$.
(1) Assume arbitrary particle interactions.
With arbitrary particle interactions, it's easy to create a situation where bosons have the same ground state energy as fermions. One simply adds energy to the bose wave functions without adding energy to the fermi states. Let's do this explicitly for a 2-particle wave function. To do this, we need to define the wave functions for $\psi_1 \times \psi_2$, that is, we need to define the tensor product. Let's use the simplest possible wave functions, spinors, and define the combined wave function $\psi_1\otimes \psi_2$ by:
$$|1\rangle=\left(\begin{array}{c}a_1\\b_1\end{array}\right)$$
$$|2\rangle=\left(\begin{array}{c}a_2\\b_2\end{array}\right)$$
$$|1\rangle\otimes |2\rangle=\left(\begin{array}{c}a_1a_2\\b_1a_2\\a_1b_2\\b_1b_2\end{array}\right).$$
Now (ignoring an unimportant factor of $\sqrt{1/2}$ from here on) a fermi symmetrization looks like this:
$$|1\rangle\otimes|2\rangle-|2\rangle\otimes|1\rangle=\left(\begin{array}{c}a_1a_2-a_2a_1\\b_1a_2-b_2a_1\\a_1b_2-a_2b_1\\b_1b_2-b_2b_1\end{array}\right)= \left(\begin{array}{c}0\\b_1a_2-b_2a_1\\a_1b_2-a_2b_1\\0\end{array}\right),$$
while a bose symmetrization looks like this:
$$|1\rangle\otimes|2\rangle+|2\rangle\otimes|1\rangle=\left(\begin{array}{c}a_1a_2+a_2a_1\\b_1a_2+b_2a_1\\a_1b_2+a_2b_1\\b_1b_2+b_2b_1\end{array}\right)= \left(\begin{array}{c}2a_1a_2\\b_1a_2+b_2a_1\\a_1b_2+a_2b_1\\2b_1b_2\end{array}\right).$$
So to add energy only to the bose symmetry case, all we have to do is to include a potential term that looks like:
$$V = \left(\begin{array}{cccc}E&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&E\end{array}\right).$$
Note that the above matrix has eigenvalues of 0 and $E$, but the $E$ eigenvalues are only accessible to a boson, the fermions automatically have an energy 0 for this potential. For a bose state to avoid the energy addition, it needs to have $a_1a_2=b_1b_2=0$. There are two ways to do this; either have $a_1=0$ or $a_2=0$. Since you can't have both $a_1$ and $a_2$ zero, or have both $b_1$ and $b_2$ zero, a little algebra will show that there is only one bose state that avoids the $E$:
$$\left(\begin{array}{c}0\\1\\1\\0\end{array}\right)$$
The above is the bose symmetrization of (1,0) x (0,1). On the other hand, the fermi symmetrization of these two states is:
$$\left(\begin{array}{c}0\\+1\\-1\\0\end{array}\right).$$
So to add an energy to the bose symmetrization but not the fermi symmetrization, use a potential of:
$$V = \left(\begin{array}{cccc}E&0&0&0\\0&E/2&E/2&0\\0&E/2&E/2&0\\0&0&0&E\end{array}\right).$$
The above has three eigenvectors with eigenvalue $E$ and only a single eigenvector with eigenvalue 0; this eigenvalue is accessible only to fermions.
You ask
how can you formally go from the summation over $s$ to the double sum over $s_1$ and $s_2$?
As we'll see in a moment, passing to the double sum relies on the mathematical fact (about tensor products of Hilbert spaces) that if $s_1$ labels a basis of states for system $1$, and if $s_2$ labels a basis of states for system $2$, then the set of all pairs $(s_1, s_2)$ labels a basis of states for the composite system.
Let's assume that the system at hand is a quantum system described by a state space (Hilbert space) $\mathcal H$. Let the (state) labels $s$ index an orthonormal basis of $\mathcal H$ consisting of eigenvectors $|s\rangle$ of the total system's Hamiltonian $H$;
\begin{align}
H|s\rangle = E(s)|s\rangle,
\end{align}
then the partition function is obtained by summing over these states;
\begin{align}
Z = \sum_s e^{-\beta E(s)}.
\end{align}
Notice that the state labels $s$ don't have to be numbers. They could, for example, be pairs of numbers or triples of numbers, whatever is most convenient to label the possible states for the system at hand.
Now, suppose that the system consists of a pair of subsystems $1$ and $2$, then the Hilbert space of the combined system can be written as a tensor product of the Hilbert spaces for the individual subsystems;
\begin{align}
\mathcal H = \mathcal H_1\otimes \mathcal H_2.
\end{align}
Let $H_1$ denote the hamiltonian for subsystem $1$, and let $H_2$ denote the hamiltonian for subsystem $2$. Let $\{|1, s_1\rangle\}$ be an orthonormal basis for $\mathcal H_1$ consisting of eigenvectors of $H_1$, and let $\{|2, s_2\rangle\}$ be an orthonormal basis for $\mathcal H_2$ consisting of eigenvectors of $H_2$, then we have the following basic fact;
The tensor product states
\begin{align}
|s_1, s_2\rangle := |1,s_1\rangle\otimes|2,s_2\rangle
\end{align}
yield an orthonormal basis for the full state space $\mathcal H = \mathcal H_1\otimes \mathcal H_2$.
In particular, the set of states one needs to sum over in the partition function can be enumerated by pairs $s=(s_1, s_2)$. Moreover, if the systems are non-interacting, then the Hamiltonian of the full system is essentially the sum of the Hamiltonians of the individual subsystems (with appropriate identity operator factors);
\begin{align}
H = H_1\otimes I_2 + I_1\otimes H_2,
\end{align}
so that if $E_1(s_1)$ and $E_2(s_2)$ denote the energy eigenvalues of $H_1$ and $H_2$ corresponding to the states $|1,s_1\rangle$ and $|2, s_2\rangle$, then the energy $E(s_1, s_2)$ of a tensor product basis state is their sum;
\begin{align}
H|s_1,s_2\rangle = E(s_1, s_2)|s_1, s_2\rangle = \big(E_1(s_1) + E_2(s_2)\big)|s_1,s_2\rangle,
\end{align}
and the partition function can be written as a sum over the tensor product basis states;
\begin{align}
Z = \sum_{(s_1,s_2)} e^{-\beta E(s_1, s_2)} = \sum_{s_1}\sum_{s_2} e^{-\beta\big(E_1(s_2) + E_2(s_2)\big)} = \sum_{s_1} e^{-\beta E_1(s_1)} \sum_{s_2} e^{-\beta E(s_2)} = Z_1Z_2
\end{align}
where in the second equality we have used the fact that a single sum over all possible pairs $(s_1, s_2)$ is equivalent to iterated sums over the possible values of $s_1$ and $s_2$.
Best Answer
For the partition sum, you have so sum $e^{-E}$ ($T=1$) over all possible eigenstates of the system where $E$ is the energy of the corresponding state.
Two bosons can be in the 10 states $|kl\rangle$, with $1\leq k \leq l \leq 4$ where we accounted for the degeneracy by introducing an additional state with $E_4 =2E$. The corresponding partition sum reads (we assume the particles to be noninteracting) $$ Z_B = \sum_{k\leq l} e^{-E_k- E_l} = 1+ e^{-E} + 3 e^{-2E} +2 e^{-3 E} +3 e^{-4E}.$$
Similarly, for fermions we have 6 states $|kl\rangle$, with $1\leq k < l \leq 4$ with the partition sum $$ Z_F = e^{-E} + 2 e^{-2E} +2 e^{-3 E} + e^{-4E}.$$
So the difference of the partition functions of a pair of two bosons and that of a pair of two fermions is ;-) $$ Z_B - Z_F = 1 + e^{-2E} +2 e^{-4E}.$$