I got the hamiltonians for n1 and n2 as (n1+1/2)hw and (n2+1/2)hw. Since the an ensemble is made by N distinguishable pairs of quantum oscillators, the general canonical partition function for the system is 1/N!((sum(-BH(n1))sum(-BH(n2))), where B is the thermodynamic beta.
I got to the step...