  The modified HFB (MHFB) theory at finite temperature is derived, which
conserves the unitarity relation of the particle-density matrix. This is
achieved by constructing a modified quasiparticle-density matrix, where the
fluctuation of the quasiparticle number is microscopically built in. This
matrix can be directly obtained from the usual quasiparticle-density matrix by
applying the secondary Bogoliubov transformation, which includes the
quasiparticle occupation number. It is shown that, in the limit of constant
pairing parameter, the MHFB theory yields the previously obtained modified BCS
(MBCS) equations. It is also proved that the modified quasiparticle RPA, which
is based on the MBCS quasiparticle excitations, conserves the Ikeda sum rule.
The numerical calculations of the pairing gap, heat capacity, level density,
and level density parameter within the MBCS theory are carried out for
$^{120}$Sn. The results show that the superfluid - normal phase transition is
completely washed out. The applicability of the MBCS up to a temperature as
high as $T\sim$ 5 MeV is analyzed in detail.
