Fermi-Dirac Distribution

Clear["Global`*"]

Use the Fermi-Dirac distribution.

nFD[ϵ_, μ_, T_] := (1/(Exp[(ϵ - μ)/(k T)] + 1)) // Quiet

The chemical potential can by found by the normalization condition, N = (∫_0^∞)n(E)dE .  For a typical free electron gas model of a conductor, a value of 5 eV is reasonable.  Use that for an illustrative graph.  For extreme ranges, Mathematica issues warnings about loss of precision.  We will ignore them here since they don’t impact our qualitative results here.

k = QuantityMagnitude[UnitConvert[Quantity[1., "BoltzmannConstant"], "Electronvolts"/"Kelvins"]] (* Boltzmann's Constant in eV/K *)

μ = 5.00; (* in eV *)

This plot shows the FD distribution for different temperatures, and keeps the very-low temperature graph on for comparison.

N[nFD[5,μ, 1 ]]

Manipulate[Plot[{nFD[ϵ, μ, 0.1], nFD[ϵ, μ, T]}, {ϵ,0, 10},
LabelStyle -> Larger, AxesLabel -> {"ϵ (eV)", "(Subscript[n, FD])"},  
PlotLegends->{"Low T", StringForm["T = `` (K)", T]}, 
PlotRange -> {{0, 10}, {0, 1}}, ImageSize -> Scaled[0.7]] // Quiet, {T, 10, 10000, Appearance ->"Open"}]