No I believe it wants you to do this for the 3d model. However, note it asks for the average energy of a phonon not the average energy. This means you want U/N where U is the average energy of the system (proportional to T^4) and N here is the average number of phonons present (which will be proportional to T^3). The integrals in the 1996 paper will be exactly what you need for this calculation.
In 2006, yes I suspect you can leave result in terms of integrals. I don't think they expect you to know values of the Riemann zeta function etc.