黑體中不同頻率的光子數分布
光子（photon）屬於Bose-Einstein粒子，各狀態的光子個數不受包例不相容原理限制（關於受包立不相容原理限制的Fermi-Dirac粒子的統計物理性質，可參考 自由電子海的能量分布與Fermi能階。在黑體中假如有M個狀態、N個相同無法辨明的光子，則系統總共有\(W=H^M_N = C^{N+M-1}_{N}\)種可能的組態，此時系統的亂度為
$$ S = k_B \ln(W) = k_B \ln \left( \frac{(N+M-1)!}{N!(M-1)!}\right) $$補充說明：相同粒子在系統中的三種狀態數計算方式
在統計力學中，假如系統中有M個狀態、N個相同粒子，則計算系統組態數的方式有三種。
- Bose-Einstein粒子：各個狀態中的粒子數無限制，但即使在微觀下也沒有辦法辨明各個粒子（因粒子物質波會彼此干擾）－因此只有各個狀態的粒子數為可觀測之資訊。系統之組態數為\(W_{BE} = C^{N+M-1}_{N}\)。
- Fermi-Dirac粒子：各個狀態中的粒子數最多為1個，微觀下也沒有辦法辨明各個粒子。系統之組態數為\(W_{FD} = C^{M}_{N}\)。
- Maxwell-Boltzmann粒子：各個狀態中的粒子數無限制，且微觀下可辨明各個粒子；不過由於各粒子彼此相同，在巨觀下仍只有各個狀態的粒子數為可觀測之資訊。因此系統之組態數為\(W_{MB} = \frac{M^N}{N!}\)。
相同的M、N下，\(W_{FD} < W_{MB} < W_{BE}\)。
利用Stirling近似，當\(M\)和\(N\)皆趨近於無限時，系統亂度有底下的漸近關係
$$ S \rightarrow k_B \left( (N+M) \ln(N+M) - N \ln(N) - M \ln(M) \right) $$$$ = k_B M \left( (1+n) \ln(N+M) - n \ln(N) - \ln(M) \right) $$
$$ = k_B M \left( (1+n) \ln(1+n) - n \ln(n) \right) $$
其中\(n=\frac{N}{M}\)為每個狀態平均光子數。
定溫、定容，但粒子數可改變的系統中，我們定義總位能（grand potential）\(\Phi_G\)為
$$ \Phi_G = U - TS - \mu N $$其中\(U\)為系統內能、\(T\)為系統溫度、\(\mu\)為系統化學能。根據熱力學第二定律，在熱力學平衡的情況下，光子在各個狀態的分布\(N_i\)（\(i\)表不同狀態）會使得總位能有最小值。定義\(n_i = \frac{N_i}{M}\)，則系統總位能可以寫成
$$ \Phi_G = M \sum_i (\varepsilon_i - \mu) n_i - k_B T \left( (1+n_i) \ln(1+n_i) - n_i \ln(n_i) \right) $$其中\(\varepsilon_i\)為光子在狀態\(i\)時的內能。
當系統總位能有最小值時，其對各個\(n_i\)的偏微分應當為零：
$$ \left. \frac{\partial \Phi_G}{\partial n_i} \right|_{T,V, \mu} = M \left( \varepsilon_i - \mu - k_B T \left( \ln(1+n_i) - \ln(n_i) \right) \right) = 0 $$$$ \ln \left( \frac{1 + n_i}{n_i} \right) = \frac{\varepsilon_i - \mu}{k_B T} $$$$ \frac{1 + n_i}{n_i} = \exp \left(\frac{\varepsilon_i - \mu}{k_B T} \right) $$$$ n_i = \frac{1}{\exp \left(\frac{\varepsilon_i - \mu}{k_B T} \right) - 1} $$此即Bose-Einstein分布。而黑體中光子可以自由地被表面吸收獲釋放，因此該系統內並不需要考慮化學能，即\(\mu_i=0\)。代入單一光子所具有的內能\(\varepsilon_i = h|\vec{\nu}_i| = c |\vec{p}_i|\)（這裡\(h\)為普朗克常數、\(\vec{\nu}_i\)和\(\vec{p}_i\)各為狀態\(i\)之光波頻率與光子動量），則光子的在不同狀態的分布為
$$ f_i = \frac{1}{Z} \left( \exp \left(\frac{h\nu_i}{k_B T} \right) - 1 \right)^{-1} $$其中\(Z = \sum_i \left( \exp \left(\frac{h \nu_i}{k_B T} \right) - 1 \right)^{-1}\)。
為了簡化計算，我們試著用連續函數來近似\(n_i\)和\(n\)（當單位頻率下的狀態密度趨近無限，致使相鄰狀態對應到的光子頻率無限接近時，這樣的近似是成立的）。
在動量空間中，一無限小之動量向量量值所佔的體積為
$$ 4 \pi |\vec{p}|^2 d|\vec{p}| = \frac{4 \pi h^3}{c^3} \nu^2 d\nu $$我們在上式以及接下來把\(|\vec{\nu}|\)簡寫成\(\nu\)。
由於光波為電磁場構成的平面沿著平面法向量傳遞之波，每個\(\vec{\nu}\)對應到法向量，都可以找到兩組線性獨立的偏振光束做基底，構成該\(\vec{\nu}\)對應到的光波。因此在由位置和頻率向量組成的組態空間中，一無限小之動量向量內包含的狀態數為
$$ dM = 2 \int_{\mathbf{R}^3} \frac{4 \pi h^3}{c^3 h^3} \nu^2 d\nu d\vec{x} = \frac{8 \pi V}{c^3} \nu^2 d\nu $$補充說明：連續與離散組態空間的關係
根據萊曼積分的定義，單一方向的連續組態空間中某函數的積分為
$$ \int_{\mathbf{R} \times \mathbf{R}} g(x,p)dx dp = \lim_{(\Delta x, \Delta p) \rightarrow 0} \sum_{i,j} g(x_i, p_j) \Delta x \Delta p $$而現實中，目前物理理論能有效解釋到的\(\Delta x \Delta p\)最小值即為\(h\)，因此以積分式近似離散組態空間的加總時，有底下的關係式
$$ \sum_{i,j} g(x_i, p_j)=\frac{1}{\Delta x \Delta p}\sum_{i,j} g(x_i, p_j) \Delta x \Delta p \approx \frac{1}{h} \int_{\mathbf{R} \times \mathbf{R}} g(x,p)dx dp $$其中\(V\)為黑體體積。而一無限小之動量向量內包含的光子數則為
$$ dN = n(\nu) dM = \frac{8 \pi V \nu^2}{c^3} \left( \exp\left(\frac{h\nu_i}{k_B T} \right) - 1 \right)^{-1} d\nu $$定義\(n(\nu, T)= \frac{1}{V} \frac{dN}{d \nu} \)為單位體積黑體中，特定光波頻率所具有的光子數密度，則
$$ n(\nu, T) = \frac{8 \pi \nu^2}{c^3} \left( \exp\left(\frac{h\nu}{k_B T} \right) - 1 \right)^{-1} $$因此單位體積黑體中的光子總數\(\bar{N}\)為
$$ \bar{N}(T) = \int^{\infty}_{0} n(\nu, T) d \nu = \frac{8 \pi }{c^3}\int^{\infty}_{0} \frac{\nu^2}{\exp\left(\frac{h\nu}{k_B T} \right) - 1} d \nu $$利用變數變換\(u=\frac{h\nu}{k_B T}\)，上式可改寫成
$$ \bar{N}(T) = 8 \pi \left( \frac{k_B T}{c h} \right)^3 \int^{\infty}_{0} \frac{u^2}{e^u - 1} d u = 8 \pi \left( \frac{k_B T}{c h} \right)^3 \int^{\infty}_{0} u^2 e^{-u} \sum^{\infty}_{k=0} e^{-ku} du $$$$ \bar{N}(T) = 8 \pi \left( \frac{k_B T}{c h} \right)^3 \int^{\infty}_{0} u^2 \sum^{\infty}_{k=1} e^{-ku} du = 8 \pi \left( \frac{k_B T}{c h} \right)^3 \sum^{\infty}_{k=1} \int^{\infty}_{0} u^2 e^{-ku} du $$此時再做一次變數變換\(t = ku\)，則上述積分式可改寫成
$$ \bar{N}(T) = 8 \pi \left( \frac{k_B T}{c h} \right)^3 \sum^{\infty}_{k=1} \frac{1}{k^3} \int^{\infty}_{0} t^2 e^{-t} dt = 8 \pi \left( \frac{k_B T}{c h} \right)^3 \Gamma(3) \zeta(3) $$其中\(\Gamma(\cdot)\)為Gamma函數，當自變數為正整數\(m\)時有值\((m-1)!\)；\(\zeta(\cdot)\)為Zeta函數，當自變數為奇數時沒有辨法用有理數或基本常數（\(e\)、\(\pi\)等）表示。因此可把\(\bar{N}(T)\)寫成
$$ \bar{N}(T) = 16 \pi \left( \frac{k_B T}{c h} \right)^3 \zeta(3) $$而特定頻率的光子數機率密度函數\(f_n(\nu,T)\)則為
$$ f_n(\nu,T) = \frac{n(\nu, T)}{\bar{N}(T)} = \frac{1}{2 \zeta(3)} \left( \frac{h}{k_B T} \right)^3 \frac{\nu^2}{\exp\left(\frac{h\nu}{k_B T} \right) - 1} $$黑體中不同頻率的光子能量分布和韋恩定律
單位體積黑體中，特定光波頻率所具有的光子能量密度\(\varepsilon(\nu, T)\)為
$$ \varepsilon(\nu, T) = \varepsilon(\nu) n(\nu, T) = \frac{8 \pi h \nu^3}{c^3} \left( \exp\left(\frac{h\nu}{k_B T} \right) - 1 \right)^{-1} $$則單位體積黑體中的光子總能量\(\bar{\Epsilon}(T)\)為
$$ \bar{\Epsilon}(T) = \int^{\infty}_{0} \varepsilon(\nu, T) d \nu = \frac{8 \pi h}{c^3}\int^{\infty}_{0} \frac{\nu^3}{\exp\left(\frac{h\nu}{k_B T} \right) - 1} d \nu $$和前一節相同，利用變數變換\(u=\frac{h\nu}{k_B T}\)，上式可改寫成
$$ \bar{\Epsilon}(T) = \frac{8 \pi (k_B T)^4}{(ch)^3}\int^{\infty}_{0} \frac{u^3}{e^u - 1} d u = \frac{8 \pi (k_B T)^4}{(ch)^3} \Gamma(4) \zeta(4) = \frac{48 \pi (k_B T)^4}{(ch)^3} \zeta(4) $$當\(\zeta(\cdot)\)的應變數為偶數時，可利用傅立葉級數等方式，以有理數和基本常數表達其值。在本例中，\(\zeta(4) = \frac{\pi^4}{90}\)。因此
$$ \bar{\Epsilon}(T) =\frac{48 \pi (k_B T)^4}{(ch)^3} \zeta(4) = \frac{8 \pi^5 (k_B T)^4}{15 (ch)^3} = \frac{ \pi^2 (k_B T)^4}{15 (c \hbar)^3} $$補充說明：\(\zeta(4)\)的求法
\(\zeta(4)\)的求法有很多種；底下僅介紹利用傅立葉級數的作法。我們先定義一週期性函數：
$$ f(x) = \left\{ \begin{array}{lcl} x^5 &,& 0 < x < 1 \\ \\ \frac{1}{2} &,& x = 1 \\ \\ f(x + 1) &,& \forall x \end{array} \right. $$則上式可用一傅立葉級數表示：
$$ f(x) = \sum^{\infty}_{k= -\infty} c_k \exp (j 2 \pi k x) $$其中各係數\(c_k\)可以積分的方式求出：
$$ c_k = \int^{1}_{0} f(x) \exp (-j 2 \pi k x) dx $$$$ \int^{1}_{0} x^5 \exp (-j 2 \pi k x) dx = \left.\frac{x^5}{-j 2 \pi k} \right|^{1}_{0} + \frac{5}{j 2 \pi k} \int^{1}_{0} x^4 \exp (-j 2 \pi k x) dx $$$$ = \frac{j}{2 \pi k} + \frac{5}{j2 \pi k} \left( \left.\frac{x^4}{-j 2 \pi k} \right|^{1}_{0} + \frac{4}{j 2 \pi k}\int^{1}_{0} x^3 \exp (-j 2 \pi k x) dx \right) $$$$ = \frac{j}{2 \pi k} + \frac{5}{4 \pi^2 k^2} - \frac{5}{\pi^2 k^2} \left( \left.\frac{x^3}{-j 2 \pi k} \right|^{1}_{0} + \frac{3}{j 2 \pi k} \int^{1}_{0} x^2 \exp (-j 2 \pi k x) dx \right) $$$$ = \frac{j}{2 \pi k} + \frac{5}{4 \pi^2 k^2} - \frac{j5}{2 \pi^3 k^3} - \frac{15}{j2 \pi^3 k^3} \left( \left.\frac{x^2}{-j 2 \pi k} \right|^{1}_{0} + \frac{2}{j 2 \pi k} \int^{1}_{0} x \exp (-j 2 \pi k x) dx \right) $$$$ = \frac{j}{2 \pi k} + \frac{5}{4 \pi^2 k^2} - \frac{j5}{2 \pi^3 k^3} - \frac{15}{4 \pi^4 k^4} + \frac{15}{2 \pi^4 k^4} \left( \left.\frac{x}{-j 2 \pi k} \right|^{1}_{0}\right) $$$$ = \frac{j}{2 \pi k} + \frac{5}{4 \pi^2 k^2} - \frac{j5}{2 \pi^3 k^3} - \frac{15}{4 \pi^4 k^4} + \frac{j15}{4 \pi^4 k^5} $$現在代入\(f(1) = \frac{1}{2}\)：
$$ \sum^{\infty}_{k= -\infty} c_k = c_0 + \sum^{\infty}_{k= 0} (c_k + c_{-k}) = \frac{1}{2} $$已知\(c_0 = \int^{1}_{0} f(x) dx = \frac{1}{6}\)且\(c_{k}\)和\(c_{-k}\)互為共軛複數，因此
$$ \sum^{\infty}_{k= -\infty} c_k = \frac{1}{6} + 2\sum^{\infty}_{k= 0} \left( \frac{5}{4 \pi^2 k^2} - \frac{15}{4 \pi^4 k^4} \right) = \frac{1}{2} $$$$ \sum^{\infty}_{k= 0} \left( \frac{5}{4 \pi^2 k^2} - \frac{15}{4 \pi^4 k^4} \right) = \frac{1}{6} $$已知\(\zeta(2) = \frac{\pi^2}{6}\)的情況下，則可以導出
$$ \sum^{\infty}_{k= 0} \frac{1}{k^4} = \frac{4 \pi^4}{15} \left( \frac{5}{4 \pi^2} \sum^{\infty}_{k= 0} \frac{1}{k^2} - \frac{1}{6} \right) = \frac{\pi^4}{90} $$事實上，利用同樣的傅立葉級數法求出函數
$$ g(x) = \left\{ \begin{array}{lcl} x^3 &,& 0 < x < 1 \\ \\ \frac{1}{2} &,& x = 1 \\ \\ f(x + 1) &,& \forall x \end{array} \right. $$的傅立葉級數各項係數後，即可以得出\(\zeta(2)\)。這要在求\(\zeta(4)\)之前完成。之後，利用這個遞迴方法，理論上便可以得到任意偶數時Zeta函數之值。
而特定頻率的光子能量機率密度函數\(f_\varepsilon(\nu,T)\)則為
$$ f_\varepsilon(\nu,T) = \frac{\varepsilon(\nu, T)}{\bar{\Epsilon}(T)} = 15 \left( \frac{h}{\pi k_B T} \right)^4 \frac{\nu^3}{\exp\left(\frac{h\nu}{k_B T} \right) - 1} $$光子能量最強的\(f_\varepsilon(\nu,T)\)對頻率微分為零：
$$ \left. \frac{\partial f_\varepsilon(\nu,T)}{\partial \nu} \right|_T = 15 \left( \frac{h}{\pi k_B T} \right)^4 \frac{3 \nu^2 \left( \exp\left(\frac{h\nu}{k_B T} \right) - 1 \right) - \frac{h\nu^3}{k_B T}\exp\left(\frac{h\nu}{k_B T} \right)}{\left( \exp\left(\frac{h\nu}{k_B T} \right) - 1 \right)^2} = 0 $$$$ 3 \left( \exp\left(\frac{h\nu}{k_B T} \right) - 1 \right) - \frac{h\nu}{k_B T}\exp\left(\frac{h\nu}{k_B T} \right) = 0 $$$$ \left( \frac{h\nu^*}{k_B T} - 3 \right) \exp\left( \frac{h\nu^*}{k_B T} - 3 \right) = -3 e^{-3} $$$$ \frac{h\nu^*}{k_B T} - 3 = W_0(-3 e^{-3}) $$$$ \frac{\nu^*}{T}= \frac{k_B}{h} (3 +W_0(-3 e^{-3})) $$其中\(W_0(\cdot)\)為第零型Lambert W函數。由此我們便可得到韋恩定律（Wien’s displacement law）。
黑體表面輻射能量（史蒂芬－波茲曼定律）
由前面的推導中可以看出光子能量的機率分布僅和動量向量／頻率向量的量值有關，而和其方向無關。我們可以因而推論出在一定時間內、黑體一極小表面積\(dA\)上，接收到來自各個立體角\(d \Omega\)的光子能量皆應相同。這代表\(dA\)接收到來自各個立體角的光能功率\(dW\)為
$$ dW(\Omega) = \frac{c \cos(\phi)}{4 \pi} \bar{\Epsilon}(T) dA d \Omega $$其中\(\cos(\phi)\)為某立體角對應的方向向量和極小表面積法向量（方向朝內）的內積值，\(\phi\)為某立體角和極小表面積法向量的夾角。
定義\(\bar{W}\) 為黑體單位表面積接收到來自內部光子的功率，則可以透過下列積分式求出其：
$$ \bar{W} = \frac{c \bar{\Epsilon}(T)}{4 \pi} \int^{\frac{\pi}{2}}_{0} \int^{2 \pi}_{0} \cos(\phi) \sin(\phi) d \theta d \phi = \frac{c \bar{\Epsilon}(T)}{4} \int^{\frac{\pi}{2}}_{0} \sin(2 \phi) d \phi $$$$ \bar{W} = \frac{c \bar{\Epsilon}(T)}{4} = \frac{ \pi^2 {k_B}^4 T^4}{60 c^2 \hbar^3} = \sigma_{SB} T^4 $$此即史蒂芬－波茲曼定律（Stefan-Boltzmann law）。其中\(\sigma_{SB} = \frac{ \pi^2 {k_B}^4 }{60 c^2 \hbar^3} = 5.670 \times 10^{-8} W m^{-2} K^{-4}\)為史蒂芬－波茲曼常數。
一般而言，黑體表面的光輻射能收支會假設處於穩定態，因此上式推導出的黑體單位表面積自內部接收的光能功率，同時也會是黑體單位表面積向外發送的光能功率。另外，即使黑體表面有受到其他作用和鄰近系統交換能量（例如傳導或對流），只要假設黑體內近表面處處於局部熱力學平衡（因此溫度可被定義），史蒂芬－波茲曼定律仍可適用。
黑體表面輻射能量的頻率分布（普朗克定律）
普朗克定律（Planck’s Law）描述黑體單位面積、單位立體角向四周發送的光能功率隨頻率的關係，由\(\varepsilon(\nu,T)\)乘上\(\frac{c}{4 \pi}\)：
$$ B(\nu, T) = \frac{2 h}{c^2} \frac{\nu^3}{\exp\left(\frac{h\nu}{k_B T} \right) - 1} $$有時普朗克定律也會以光波波長來表示；由於\(\nu \lambda = c\)，而我們又有
$$ \int^{\nu_2}_{\nu_1} B(\nu, T) d\nu = \int^{\frac{c}{\nu_2}}_{\frac{c}{\nu_1}} B(\nu, T) \frac{d \nu}{d \lambda} d \lambda = -\int^{\frac{c}{\nu_1}}_{\frac{c}{\nu_2}} B(\lambda, T) d \lambda $$因此
$$ B(\lambda, T) = -B(\nu, T) \frac{d \nu}{d \lambda} = \frac{c}{\lambda^2} \frac{2 h}{c^2} \frac{c^3}{\lambda^3} \frac{1}{\exp\left(\frac{hc}{\lambda k_B T} \right) - 1} $$$$ B(\lambda, T) = \frac{2 h c^2}{\lambda^5} \left( \exp\left(\frac{hc}{\lambda k_B T} \right) - 1 \right)^{-1} $$黑體單位表面積接收的光子通量
在討論太陽能板的發電原理時，會需要知道黑體單位表面積接收的光子通量。這個參數的求法類似黑體表面輻射能量的求法；我們首先定義各個立體角的光子通量為
$$ dF(\Omega) = \frac{c \cos(\phi)}{4 \pi} \bar{N}(T) dA d \Omega $$定義\(\bar{F}\) 為黑體單位表面積接收到來自內部的光子通量，則可以透過下列積分式求出其：
$$ \bar{F} = \frac{c \bar{N}(T)}{4 \pi} \int^{\frac{\pi}{2}}_{0} \int^{2 \pi}_{0} \cos(\phi) \sin(\phi) d \theta d \phi = \frac{c \bar{N}(T)}{4} \int^{\frac{\pi}{2}}_{0} \sin(2 \phi) d \phi $$$$ \bar{F} = \frac{c \bar{N}(T)}{4} = \frac{4 \pi {k_B}^3 T^3}{c^2 h^3} \zeta(3) = \frac{30 \zeta(3)}{\pi^4 k_B} \sigma_{SB} T^3 $$則可以求出黑體單位表面積自內接收／向外發送的光能功率和光子通量的關係式：
$$ \bar{F} = \frac{30 \zeta(3)}{\pi^4 k_B T} \bar{W} $$太陽常數
導出史蒂芬－波茲曼定律後，我們便可以推導出太陽常數（solar constant）。太陽常數嚴格來說並不是一個真正的常數，因為如底下所要揭示的，其為太陽表面溫度以及日地距離所決定；不過一般應用上這幾個變數的變化不大，固可概略地視為定值。
假設太陽唯一黑體，則太陽表面上單位面積對外輻射功率為\(\sigma_{SB} {T_{S}}^4\)。另一方面，假設太陽半徑為\(r_S\)、而日地距離為\(L\)、而太陽輻射抵達地球過程沒有損耗，則我們可以導出在地球大氣層頂、對準太陽的方向上，單位面積接收到太陽輻射的功率為
$$ \sigma_{SB} {T_{S}}^4 \times \frac{{r_S}^2}{L^2} $$此即太陽常數\(G_{SC}\)。代入各參數數值，\(G_{SC}\) 大約為\(1360.8 Wm^{-1}\)。