We adopt a new approach, thermo entangled representation, to study time evolution of density operator in thermal environment. We then investigate the analytical expressions of Wigner function (WF) evolution of arbitrary number excited coherent states (ECSs) and excited even (odd) coherent states (EECSs, EOCSs) in thermal environment, respectively. In addition, their nonclassicality is numerically discussed by exploring the negativity of WF with decay time in thermal channel, respectively. It is found that WF loses its non-Gaussian nature and becomes Gaussian after long times.
