はじめに
空の色はなんで青いのかといった質問に対しては「レイリー散乱による光の波長依存性の影響」といった回答が多くされている。
これは実際に学研の110番や検索を入れたときトップに出てくるページなどでもあるように、わりと語りつかされていると思う。
しかし、具体的に「なにがどのくらい」レイリー散乱しているかについては意外なほど記事が少ない。というか見つけられなかったので自分で計算してみることにした。
レイリー散乱を起こす散乱体の正体
レイリー散乱とは光の波長に対して小さな物体があるときに発生する散乱の一つである。ここでおおよそ波長の1/10くらいでレイリー散乱の影響が大きくなると言われており、それより大きな散乱体はミー散乱が支配的になる。
また、空が青い理由についてはレイリー散乱係数が波長の-4乗に比例するからである。
散乱係数が大きいほど散乱されるため、その空間はその散乱された色として見える。
よって、散乱係数が大きいときは波長が短く、波長が短いということは青になるといった理屈が成り立つのである。
よって、空の色が青くなると行った説明ができる。
ここまでは前述の通り、色んなサイトや論文に書かれている事実があるが、具体的にレイリー散乱係数はどのくらいかがないので今回計算した
さて、レイリー散乱を引き起こす空の具体的な粒子というのは酸素分子や窒素分子らしい[1]。
つまり、大気中の分子組成比がレイリー散乱を起こしているようだ。
標準気圧＆温度環境下では1mol/22.4Lの密度というのは高校までの物理で明らかになっているが、実際空では高度によって圧力も温度も変化することが予想される。
高度による密度等をまとめられたサイトを参考にすると
温度と圧力を入力変数にした密度の変化は
\[
\rho = \frac{P \cdot M_0}{R^* \cdot T} = \frac{0.0034837 \cdot P}{T + 273.15}
\]
ここで、\(M_0\) はモル質量、\(R^*\) は気体定数、\(T\) は絶対温度（ケルビン）となる。
さらに温度と圧力は高度によって以下のように表される。
| ジオポテンシャル高度 H [km] | 空気の温度 T [℃] |
|---|---|
| 0 ≤ H ≤ 11 | T = 15 − 6.5 × H |
| 11 < H ≤ 20 | T = −56.5 |
| 20 < H ≤ 32 | T = −76.5 + H |
| 32 < H ≤ 47 | T = −134.1 + 2.8 × H |
| 47 < H ≤ 51 | T = −2.5 |
| 51 < H ≤ 71 | T = 140.3 − 2.8 × H |
| 71 < H ≤ 84.852 | T = 83.5 − 2 × H |
| ジオポテンシャル高度 H [km] | 空気の圧力 P [Pa] の計算式 |
|---|---|
| 0 ≤ H < 11 | P = 101325 × (288.15 / (T + 273.15))−5.256 |
| 11 ≤ H < 20 | P = 22632.064 × e−0.1577 × (H − 11) |
| 20 ≤ H < 32 | P = 5474.889 × (216.65 / (T + 273.15))34.163 |
| 32 ≤ H < 47 | P = 868.019 × (228.65 / (T + 273.15))−12.201 |
| 47 ≤ H < 51 | P = 110.906 × e−0.1262 × (H − 47) |
| 51 ≤ H < 71 | P = 66.939 × (270.65 / (T + 273.15))−12.201 |
| 71 ≤ H ≤ 84.852 | P = 3.956 × (214.65 / (T + 273.15))−17.082 |
これらの数値から窒素80%、酸素20%と仮定して実際に高度ごとの分子数を求めてみる。
ここで以下のステップで計算する。
-
空気の密度
[kg/m³] から、1m³中の空気の質量を得る
-
各成分の質量を比率で計算する（例：窒素 = 78.08%）
-
質量 ÷ モル質量（kg/mol） → モル数
その結果が以下のようになる。
ここまでのPythonコードを以下に示す。
Pythonコード
import matplotlib.pyplot as plt
import numpy as np
import japanize_matplotlib
def calculate_temperature(H):
if 0 <= H <= 11:
T = 15 - 6.5 * H
elif 11 < H <= 20:
T = -56.5
elif 20 < H <= 32:
T = -76.5 + H
elif 32 < H <= 47:
T = -134.1 + 2.8 * H
elif 47 < H <= 51:
T = -2.5
elif 51 < H <= 71:
T = 140.3 - 2.8 * H
elif 71 < H <= 84.852:
T = 83.5 - 2 * H
else:
T = None # 範囲外の場合
return T
import math
def calculate_pressure(H, T_celsius):
T_kelvin = T_celsius + 273.15 # 絶対温度に変換
if 0 <= H < 11:
P = 101325 * (288.15 / T_kelvin) ** -5.256
elif 11 <= H < 20:
P = 22632.064 * math.exp(-0.1577 * (H - 11))
elif 20 <= H < 32:
P = 5474.889 * (216.65 / T_kelvin) ** 34.163
elif 32 <= H < 47:
P = 868.019 * (228.65 / T_kelvin) ** 12.201
elif 47 <= H < 51:
P = 110.906 * math.exp(-0.1262 * (H - 47))
elif 51 <= H < 71:
P = 66.939 * (270.65 / T_kelvin) ** -12.201
elif 71 <= H <= 84.852:
P = 3.956 * (214.65 / T_kelvin) ** -17.082
else:
P = None # 範囲外
return P
def calculate_rho(T,P):
return 0.0034837*P/(T+273.15)
def air_composition_moles(rho):
# 成分の比率（モル比 ≒ 体積比）
ratio_N2 = 0.80
ratio_O2 = 0.20
# モル質量 [kg/mol]
molar_mass_N2 = 0.0280134
molar_mass_O2 = 0.0319988
# 各成分の質量 [kg]（1 m³ あたり）
mass_N2 = rho * ratio_N2
mass_O2 = rho * ratio_O2
# モル数 [mol]
mol_N2 = mass_N2 / molar_mass_N2
mol_O2 = mass_O2 / molar_mass_O2
return mol_N2, mol_O2
H = np.arange(0,84)
T = [calculate_temperature(i) for i in H]
P = [calculate_pressure(H,T[i]) for i,H in enumerate(H)]
rho = [calculate_rho(T[i],P[i]) for i in range(len(H))]
# 共通スタイル
plt.style.use("seaborn-v0_8-muted") # 落ち着いた配色
plt.rcParams['font.size'] = 12
# === 温度グラフ ===
plt.figure(figsize=(8, 5))
plt.plot(H, T, color="firebrick", linewidth=2)
plt.title("高度と気温の関係")
plt.xlabel("高度 [km]")
plt.ylabel("気温 [℃]")
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# === 圧力グラフ ===
plt.figure(figsize=(8, 5))
plt.plot(H, P, color="steelblue", linewidth=2)
plt.title("高度と気圧の関係")
plt.xlabel("高度 [km]")
plt.ylabel("圧力 [Pa]")
plt.ylim(0)
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# === 密度グラフ ===
plt.figure(figsize=(8, 5))
plt.plot(H, rho, color="darkgreen", linewidth=2)
plt.title("高度と空気密度の関係")
plt.xlabel("高度 [km]")
plt.ylabel("空気密度 [kg/m³]")
plt.ylim(0)
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
rho = np.array(rho)
mol = air_composition_moles(rho)
# === 窒素密度グラフ ===
plt.figure(figsize=(8, 5))
plt.plot(H, mol[0], color="orange", linewidth=2)
plt.title("高度と窒素mol密度の関係")
plt.xlabel("高度 [km]")
plt.ylabel("mol密度 [mol/m³]")
plt.ylim(0)
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# === 酸素密度グラフ ===
plt.figure(figsize=(8, 5))
plt.plot(H, mol[1], color="red", linewidth=2)
plt.title("高度と酸素mol密度の関係")
plt.xlabel("高度 [km]")
plt.ylabel("mol密度 [mol/m³]")
plt.ylim(0)
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()高度ごとのレイリー散乱係数を求める
レイリー散乱断面積は以下のように表される[2]。
\begin{equation}
\sigma_{Ray} = \frac{8 \pi}{3} \left( \frac{2 \pi n_{med}}{\lambda_0} \right)^4 a^6 \left( \frac{m^2-1}{m^2+2} \right)^2
\end{equation}
ここで\( N,a,\lambda_0,n_{med},m \)はそれぞれ散乱体の数、散乱体の半径、光の波長、散乱体の周囲の屈折率、散乱体の屈折率/散乱体の周囲の屈折率の比となっている。
また、ランベルト・ベールの法則(ミシガン大学の資料が割りやすいので参考にしました:transmission_lecture.ppt)より
\begin{equation}
\tau = \int_0 ^h \beta(h) dh
\end{equation}
ここで\( \beta \)は先ほどの散乱断面積を用いて
\begin{equation}
\beta_{total} = \sum_i \sigma_i N_i
\end{equation}
と表せる。ここで\(N\)は分子量(個数)を表しており、添え字の\(\sum_i\)は対象分子の総和を表している。これを用いて透過率は
\begin{equation}
T = \exp(-\tau)
\end{equation}
となる。
実際これを用いて透過率を計算してみると以下のようになる。
Pythonコード
def rayleigh_scattering(lam,H):
# all unit is cm
# should be obey: input H is 1m unit
# input wavelength unit is nm
la = lam * 10**-9 * 10**2
# O2 molecular size is 0.34nm
# a is radius
#https://klchem.co.jp/blog/2023/10/blog-12515.php
a_O2 = 0.34 / 2 * 10**-9 * 10**2
# N2 molecular size is 0.38nm
#https://pubdata.nikkan.co.jp/uploads/magazine_introduce/pdf_61d657542e1bc-2.pdf
a_N2 = 0.38 / 2 * 10**-9 * 10**2
# refractive index; セルマイヤー方程式
#https://refractiveindex.info/?shelf=main&book=O2&page=Smith
n_O2 = 1.181494*10**-4 + 9.708931*10**-3/(75.4*la**-2) + 1
n_N2 = 6.8552*10**-5 + 3.243157*10**-2/(144 - la**2) + 1
n_air = n_O2*0.2 + n_N2*0.8
m_O2 = n_O2 / n_air
m_N2 = n_N2 / n_air
sigma_O2 = (8*np.pi/3) * (2*np.pi*n_air/la) ** 4 * a_O2**6 * ((m_O2**2-1)/(m_O2**2 + 2))**2
sigma_N2 = (8*np.pi/3) * (2*np.pi*n_air/la) ** 4 * a_N2**6 * ((m_N2**2-1)/(m_N2**2 + 2))**2
T = [calculate_temperature(i) for i in H]
P = [calculate_pressure(H,T[i]) for i,H in enumerate(H)]
rho = [calculate_rho(T[i],P[i]) for i in range(len(H))]
rho = np.array(rho)
mol_O2, mol_N2 = air_composition_moles(rho)
# Avogadro constant
NA = 6.02214076e23 # mol^-1
# Height_to_scattering coefficiency β（1/m）
beta_O2 = mol_O2 * NA * sigma_O2
beta_N2 = mol_N2 * NA * sigma_N2
beta_total = beta_O2 + beta_N2
return beta_total, beta_O2, beta_N2
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import to_hex
# 高度配列（km単位 → m単位に変換）
H_km = np.arange(0, 85, 1)
H_m = H_km * 1000 # [m]
# 波長配列 [nm]
wavelengths = np.linspace(380, 2000, 200)
# 透過率を格納するリスト
transmittance = []
for lam in wavelengths:
beta_total, _, _ = rayleigh_scattering(lam, H_km)
tau = np.trapz(beta_total, x=H_m)
T_lambda = np.exp(-tau)
transmittance.append(T_lambda)
transmittance = np.array(transmittance)
# --- 波長ごとの色を設定 ---
def wavelength_to_rgb(wavelength):
"""波長[nm]をsRGBに近似変換（可視光域のみ）"""
gamma = 0.8
intensity_max = 1
if 380 <= wavelength <= 440:
r, g, b = -(wavelength - 440) / (440 - 380), 0.0, 1.0
elif 440 < wavelength <= 490:
r, g, b = 0.0, (wavelength - 440) / (490 - 440), 1.0
elif 490 < wavelength <= 510:
r, g, b = 0.0, 1.0, -(wavelength - 510) / (510 - 490)
elif 510 < wavelength <= 580:
r, g, b = (wavelength - 510) / (580 - 510), 1.0, 0.0
elif 580 < wavelength <= 645:
r, g, b = 1.0, -(wavelength - 645) / (645 - 580), 0.0
elif 645 < wavelength <= 780:
r, g, b = 1.0, 0.0, 0.0
else:
r, g, b = 1.0, 1.0, 1.0 # 非可視は白（透明にする）
# 輝度補正
def adjust(c):
return 0.0 if c == 0.0 else round(intensity_max * (c ** gamma), 3)
return (adjust(r), adjust(g), adjust(b))
# === グラフ描画 ===
plt.figure(figsize=(10, 5))
plt.title("波長と大気の透過率のグラフ")
plt.xlabel("波長 (nm)")
plt.ylabel("透過率")
plt.grid(True, linestyle='--', alpha=0.4)
plt.ylim(0.999)
for i in range(len(wavelengths) - 1):
wl = wavelengths[i]
wl_next = wavelengths[i + 1]
T1, T2 = transmittance[i], transmittance[i + 1]
color = wavelength_to_rgb(wl)
alpha = 1.0 if 380 <= wl <= 780 else 0.0 # 可視光のみ色を表示
plt.fill_between([wl, wl_next], [T1, T2], color=color, alpha=alpha)
plt.plot(wavelengths,transmittance, color="black", linewidth=2)
plt.tight_layout()
plt.show()
ここで全然レイリー散乱の影響を受けていないと思われたかもしれないが、
これは太陽光直下の太陽の色を見ているので、実際我々が太陽を見ても白く見えることからも妥当な計算となったことがわかる。
次の記事では空が青い理由をさらに追ってみようと思う。
他にも光学に関連したあほみたいな考察をしているので他の記事(虹に乗ってみたい – たこやきたべたい)を是非読んでほしい。
参考文献
[1]籔内 一博. 空が見せる多彩な色. 化学と教育 65, 28-31 (2017).
コメント