ClearAll["Global`*"]
(*Entropie systému neinteragujících kvantových oscilátorů*)
SLHO=kB M ϵ/(ℏ ω)Log[(1+ℏ ω/2/ϵ)/(1-ℏ ω/2/ϵ)]+kB M Log[ϵ/(ℏ ω)]+kB M/2Log[1-(ℏ ω/(2ϵ))^2]

FullSimplify[%]

(*inverzní teplota β=1/(kB T) jako funkce energie*)
βϵ=FullSimplify[D[SLHO,ϵ]/(kB M),Assumptions->{ϵ>0,ω>0,ℏ>0,M>0,kB>0}]

ϵT=FullSimplify[Solve[βϵ==1/(kB T),ϵ],Assumptions->{ϵ>0,ω>0,ℏ>0,M>0,T>0,kB>0}]

(*hustota energie jako funkce teploty*)
ϵ/.ϵT[[1]]

(*tepelná kapacita na jeden oscilátor*)
c=D[%,T]

Limit[c,T->0]

(*Klasická limita: ℏω<<ϵ*)
SLHO
Series[SLHO,{ℏ,0,1}]