Here's my shot at it and my whole thought process so we are checking each other.
If we put in the Hamiltonian into the exponent and talk it through, we get the following:
$$
\int e^{-\beta( \sum_i p^{2}_i/2m_1 +\sum_j p^{2}_j/2m_2)}
$$
but what is the measure of integration? well they are non interacting, so we have to look at each individual particle, of each type, and then sum over all the possible momenta they could have
$$
\int e^{-\beta( \sum_i p^{2}_i/2m_1 +\sum_j p^{2}_j/2m_2)}\Pi_i d^3p_i \Pi_j d^3p_j
$$
where I have ignored anything to do with space recklessly. now my guess for the density of states part is
$$
\frac{d^3p d^3 x}{\hbar^3}=\frac{4\pi V p^2 dp}{\hbar^3}
$$
(there is no spatial dependence, so I just integrated over the whole volume) so really I had
$$
\int e^{-\beta H}\frac{d^{3n}p d^{3n} x}{\hbar^{3n}}
$$
and I need to put that in
$$
\frac{(4\pi)^{N_1 +N_2}V^{N_1+N_2}}{\hbar^{3(N_1+N_2)}} \int_{0}^{\infty} e^{-\beta( \sum_i p^{2}_i/2m_1 +\sum_j p^{2}_j/2m_2)}(\Pi_i p^{2}_i dp_i)( \Pi_j p^{2}_j dp_j)
$$
This is a ton of do-able individual integrals, all multiplied, because the sum in the exponents are multiplied down below. I have been really sloppy about numerical factors (like $e^{-\beta N_2 \Delta}$) and constants, sorry. I hope this helps/is correct.
There are a few misconceptions here. First, the difference between Maxwell-Boltzmann statistics and Fermi/Bose statistics is not centered on distinguishability. Even in classical thermodynamics, one has to deal with questions of indistinguishable particles.
Instead, you should think of Maxwell-Boltzmann as the classical limit of quantum statistics, corresponding to high temperature and low density. Quantum effects become important when the average intermolecular distance is on the order of the de Broglie wavelength of the particles. Alternatively, one could say that classical statistics is valid as long as the probability that any given quantum state is occupied is much, much less than 1.
The question of indistinguishability can be approached classically or from a quantum mechanical viewpoint. In the quantum picture, particles are identical if they share precisely the same quantum numbers (labels, essentially), which are determined by the Hamiltonian of the system. These quantum numbers refer to intrinsic properties such as mass, spin, and electric charge, as well as extrinsic properties which correspond to values of observables which commute with the Hamiltonian such as total angular momentum. In solid state physics, particles which are localized to distinct lattice sites can also be distinguished by spatial location.
Classically, the question is a bit more vague. From the standpoint of the Hamiltonian physics, the state of the system is specified by the precise position and momentum of each constituent particle. As no two particles have the same position and momentum, all particles are distinguishable in the sense that swapping their labels results in a noticeably different state.
However, when we move to actual thermodynamics, we no longer define the state of a system as a singular point in $2N$-dimensional phase space. Instead, we deal with macrostates of the system, which are characterized by macroscopic quantities like energy and volume. If two particles exchange their position and momenta and the system remains in the same macrostate, the particles can be considered to be identical, and in fact need to be in order to resolve the Gibbs paradox.
Also, to answer your question - yes, two hydrogen atoms are indistinguishable if they have the same quantum numbers. Ortho/para-hydrogen refer to hydrogen molecules, not hydrogen atoms, but nonetheless they are also indistinguishable if they have the same quantum numbers.
Best Answer
Bosonic $1D$ $N$ harmonic oscillators allow `exceptionally' a closed form of the canonical partition function:
$Z_{N}=\prod_{n=1}^{N}\frac{q^{1/2}}{1-q^{n}}$ where $q=e^{-\beta \hbar \omega}$
This expression resembles a sort of a grand canonical partition function for a system of bosonic ``phonons" having finite number of possible energy spectra with vanishing chemical potential.
One can derive this expression by noting,
Area of a Young tableux = sum of the length of columns = sum of the length of rows