""" 熱電輸送特性計算モジュール 概要: ゼーベック係数、ローレンツ因子、電気伝導度、移動度、電子熱伝導度、パワーファクター、ZT因子を計算し、 キャリア濃度依存性をプロット・Excelファイルに保存します。 詳細説明: このモジュールは、3次元放物線バンドモデルと緩和時間近似（tau(E) ∝ E^(r-1/2)）に基づき、 フェルミ・ディラック積分を用いた厳密な輸送積分によって熱電輸送特性を評価します。 コマンドライン引数を通じて、温度、有効質量、参照移動度、格子熱伝導率、キャリア種別、散乱因子rのリストなどを設定できます。 関連リンク: Seebeck_ZT_usage """ import numpy as np import matplotlib.pyplot as plt from scipy.integrate import quad from scipy.special import expit, gamma from functools import lru_cache import argparse from openpyxl import Workbook # ========================= # constants (SI) # ========================= kB = 1.380649e-23 # J/K e = 1.602176634e-19 # C h = 6.62607015e-34 # J s hbar = h / (2.0 * np.pi) m0 = 9.1093837015e-31 # kg r""" 散乱因子 r は tau(E) ∝ E^(r-1/2) で定義 """ figsize = (10, 10) # ============================================================ # Fermi–Dirac integral (NO 1/Gamma(j+1)) # Fj(j, xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx # ============================================================ @lru_cache(maxsize=40000) def Fj(j, xi): """ 概要: フェルミ・ディラック積分 Fj(j, xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx を計算します。 詳細説明: SciPyのquad関数を使用して数値積分を実行します。結果はlru_cacheでキャッシュされ、 同じ引数での再計算を回避し、パフォーマンスを向上させます。 この関数は、ガンマ関数による正規化（1/Gamma(j+1)）は含まれていません。 引数: :param j: 積分の次数。 :type j: float :param xi: 還元フェルミ準位 (reduced Fermi level)。 :type xi: float 戻り値: :returns: フェルミ・ディラック積分 Fj(j, xi) の値。 :rtype: float """ j = float(j) xi = float(np.round(xi, 12)) integrand = lambda x: (x**j) * expit(xi - x) val, _ = quad(integrand, 0.0, np.inf, epsabs=1e-10, epsrel=1e-10, limit=400) return val def Nc_3D(T, m_eff=1.0): """ 概要: 3次元放物線バンドにおける有効状態密度 Nc を計算します。 詳細説明: 以下の公式を用いて有効状態密度を計算します。 Nc = 2 * (2π m* kB T / h^2)^(3/2) ここで、m* は有効質量、kB はボルツマン定数、T は温度、h はプランク定数です。 引数: :param T: 温度 [K]。 :type T: float :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。 :type m_eff: float, optional 戻り値: :returns: 有効状態密度 Nc [m^-3]。 :rtype: float """ mstar = m_eff * m0 return 2.0 * (2.0 * np.pi * mstar * kB * T / (h**2))**1.5 def electron_density_from_xi(xi, T, m_eff=1.0): """ 概要: 還元フェルミ準位 xi から電子濃度 n を計算します。 詳細説明: 電子濃度は以下の式で計算されます。 n [m^-3] = Nc * F^{std}_{1/2}(xi) ここで、Nc は有効状態密度、F^{std}_{1/2}(xi) = Fj(1/2, xi) / Gamma(3/2) は 標準化されたフェルミ・ディラック積分です。 引数: :param xi: 還元フェルミ準位。 :type xi: float :param T: 温度 [K]。 :type T: float :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。 :type m_eff: float, optional 戻り値: :returns: 電子濃度 n [m^-3]。 :rtype: float """ Nc = Nc_3D(T, m_eff=m_eff) F12_std = Fj(0.5, xi) / gamma(1.5) return Nc * F12_std # ============================================================ # Transport integrals for tau(E) ∝ E^(r-1/2) # I0 = (r+1) F_r # I1 = (r+2) F_{r+1} - xi (r+1) F_r # I2 = (r+3) F_{r+2} - 2 xi (r+2) F_{r+1} + xi^2 (r+1) F_r # ============================================================ def transport_I0(xi, r): """ 概要: 輸送積分 I0 を計算します。 詳細説明: tau(E) ∝ E^(r-1/2) の緩和時間に対する輸送積分の一つである I0 を計算します。 I0 = (r + 1) * F_r(xi) の式で定義されます。 引数: :param xi: 還元フェルミ準位。 :type xi: float :param r: 散乱因子。緩和時間のエネルギー依存性 E^(r-1/2) を定義します。 :type r: float 戻り値: :returns: 輸送積分 I0 の値。 :rtype: float """ return (r + 1.0) * Fj(r, xi) def transport_I1(xi, r): """ 概要: 輸送積分 I1 を計算します。 詳細説明: tau(E) ∝ E^(r-1/2) の緩和時間に対する輸送積分の一つである I1 を計算します。 I1 = (r + 2) * F_{r+1}(xi) - xi * (r + 1) * F_r(xi) の式で定義されます。 引数: :param xi: 還元フェルミ準位。 :type xi: float :param r: 散乱因子。緩和時間のエネルギー依存性 E^(r-1/2) を定義します。 :type r: float 戻り値: :returns: 輸送積分 I1 の値。 :rtype: float """ return (r + 2.0) * Fj(r + 1.0, xi) - xi * (r + 1.0) * Fj(r, xi) def transport_I2(xi, r): """ 概要: 輸送積分 I2 を計算します。 詳細説明: tau(E) ∝ E^(r-1/2) の緩和時間に対する輸送積分の一つである I2 を計算します。 I2 = (r + 3) * F_{r+2}(xi) - 2 * xi * (r + 2) * F_{r+1}(xi) + xi^2 * (r + 1) * F_r(xi) の式で定義されます。 引数: :param xi: 還元フェルミ準位。 :type xi: float :param r: 散乱因子。緩和時間のエネルギー依存性 E^(r-1/2) を定義します。 :type r: float 戻り値: :returns: 輸送積分 I2 の値。 :rtype: float """ return ( (r + 3.0) * Fj(r + 2.0, xi) - 2.0 * xi * (r + 2.0) * Fj(r + 1.0, xi) + (xi**2) * (r + 1.0) * Fj(r, xi) ) def transport_prefactor_sigma(T, m_eff=1.0): """ 概要: 電気伝導度 sigma の計算における前因子 A_sigma を計算します。 詳細説明: 電気伝導度は sigma = A_sigma * tau_pref * (kB T)^(r+1) * I0 の形で表されます。 この関数は、tau(E) = tau_pref * E^(r-1/2) と定義される緩和時間において、 以下の物理定数を含む前因子 A_sigma を計算します。 A_sigma [SI] = e^2 * 2^(3/2) * sqrt(m*) / (3 π^2 ħ^3) ここで、e は素電荷、m* は有効質量、ħ は換算プランク定数です。 引数: :param T: 温度 [K]。 :type T: float :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。 :type m_eff: float, optional 戻り値: :returns: 電気伝導度の前因子 A_sigma [SI単位]。 :rtype: float """ mstar = m_eff * m0 return e**2 * (2.0**1.5) * np.sqrt(mstar) / (3.0 * np.pi**2 * hbar**3) def sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=1.0): """ 概要: 還元フェルミ準位 xi、散乱因子 r、温度 T、緩和時間の前因子 tau_pref から 電気伝導度 sigma を計算します。 詳細説明: 電気伝導度は以下の式で計算されます。 sigma = A_sigma * tau_pref * (kB T)^(r + 1.0) * I0 ここで、A_sigma は transport_prefactor_sigma で計算される前因子、 I0 は transport_I0 で計算される輸送積分です。 引数: :param xi: 還元フェルミ準位。 :type xi: float :param r: 散乱因子。 :type r: float :param T: 温度 [K]。 :type T: float :param tau_pref: 緩和時間の前因子 tau_pref [s / J^(r-1/2)]。 :type tau_pref: float :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。 :type m_eff: float, optional 戻り値: :returns: 電気伝導度 sigma [S/m]。 :rtype: float """ A = transport_prefactor_sigma(T, m_eff=m_eff) I0 = transport_I0(xi, r) return A * tau_pref * (kB * T)**(r + 1.0) * I0 def seebeck_from_xi_transport(xi, r, carrier="electron"): """ 概要: 還元フェルミ準位 xi と散乱因子 r からゼーベック係数 S を計算します。 詳細説明: ゼーベック係数は厳密な輸送積分を用いて以下の式で計算されます。 S = (kB/q) * I1/I0 q はキャリアの種類によって異なり、電子の場合は -e、正孔の場合は +e となります。 引数: :param xi: 還元フェルミ準位。 :type xi: float :param r: 散乱因子。tau(E) ∝ E^(r-1/2) で定義されます。 :type r: float :param carrier: キャリアの種類。"electron"または"hole"を指定します。デフォルトは"electron"。 :type carrier: str, optional 戻り値: :returns: ゼーベック係数 S [V/K]。 :rtype: float 例外: :raises ValueError: carrier が "electron" または "hole" 以外の場合。 """ if carrier.lower().startswith("e"): q = -e elif carrier.lower().startswith("h"): q = +e else: raise ValueError("carrier must be 'electron' or 'hole'") I0 = transport_I0(xi, r) I1 = transport_I1(xi, r) return (kB / q) * (I1 / I0) def lorenz_from_xi_transport(xi, r): """ 概要: 還元フェルミ準位 xi と散乱因子 r からローレンツ因子 L を計算します。 詳細説明: ローレンツ因子は厳密な輸送積分を用いて以下の式で計算されます。 L = (kB/e)^2 * ( I2/I0 - (I1/I0)^2 ) 引数: :param xi: 還元フェルミ準位。 :type xi: float :param r: 散乱因子。 :type r: float 戻り値: :returns: ローレンツ因子 L [V^2/K^2 (または W ohm / K^2)]。 :rtype: float """ I0 = transport_I0(xi, r) I1 = transport_I1(xi, r) I2 = transport_I2(xi, r) return (kB / e)**2 * (I2 / I0 - (I1 / I0)**2) def mobility_from_sigma_n(sigma, n): """ 概要: 電気伝導度 sigma と電子濃度 n から移動度 mu を計算します。 詳細説明: 電気伝導度 sigma、電子濃度 n、素電荷 e、移動度 mu の間に成り立つ関係式 sigma = n * e * mu を用いて移動度を計算します。 したがって、mu = sigma / (n * e) となります。 引数: :param sigma: 電気伝導度 [S/m]。 :type sigma: float :param n: 電子濃度 [m^-3]。 :type n: float 戻り値: :returns: 移動度 mu [m^2/V/s]。 :rtype: float """ return sigma / (n * e) def tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=1.0): """ 概要: 参照移動度 mu_ref と参照還元フェルミ準位 xi_ref に基づいて、 緩和時間の前因子 tau_pref を決定します。 詳細説明: 緩和時間 tau(E) が tau(E) = tau_pref * E^(r-1/2) で定義されるとき、 指定された xi_ref と T における移動度が mu_ref となるように tau_pref を逆算します。 参照移動度 mu_ref_cm2_Vs は cm^2/V/s 単位で与えられますが、計算ではSI単位に変換されます。 引数: :param mu_ref_cm2_Vs: 参照移動度 [cm^2/V/s]。 :type mu_ref_cm2_Vs: float :param xi_ref: 参照還元フェルミ準位。この xi で参照移動度が再現されます。 :type xi_ref: float :param r: 散乱因子。 :type r: float :param T: 温度 [K]。 :type T: float :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。 :type m_eff: float, optional 戻り値: :returns: 緩和時間の前因子 tau_pref [s / J^(r-1/2)]。 :rtype: float """ mu_ref_SI = mu_ref_cm2_Vs * 1e-4 # cm^2/V/s -> m^2/V/s n_ref = electron_density_from_xi(xi_ref, T=T, m_eff=m_eff) sigma_ref = n_ref * e * mu_ref_SI A = transport_prefactor_sigma(T, m_eff=m_eff) I0_ref = transport_I0(xi_ref, r) tau_pref = sigma_ref / (A * (kB * T)**(r + 1.0) * I0_ref) return tau_pref def tau_at_energy(E_J, tau_pref, r): """ 概要: 特定のエネルギー E_J における緩和時間 tau(E_J) を計算します。 詳細説明: 緩和時間は以下の式で定義されます。 tau(E) = tau_pref * E^(r-1/2) E_J が0以下の場合は np.nan を返します。 引数: :param E_J: エネルギー [J]。 :type E_J: float :param tau_pref: 緩和時間の前因子 [s / J^(r-1/2)]。 :type tau_pref: float :param r: 散乱因子。 :type r: float 戻り値: :returns: エネルギー E_J における緩和時間 tau(E_J) [s]。 :rtype: float """ if E_J <= 0.0: return np.nan return tau_pref * (E_J**(r - 0.5)) def equivalent_l0_from_tau_pref(tau_pref, r, m_eff=1.0): """ 概要: 緩和時間の前因子 tau_pref から等価な平均自由行程 l0 を計算します。 詳細説明: 一部のスライドでの慣例的な定義 tau(E) = sqrt(m*/2) * [ l0(T) * e_chg^(-r) ] * E^(r-1/2) から、l0(T) を逆算します。 したがって、l0(T) = tau_pref * sqrt(2/m*) * e_chg^r となります。 引数: :param tau_pref: 緩和時間の前因子 [s / J^(r-1/2)]。 :type tau_pref: float :param r: 散乱因子。 :type r: float :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。 :type m_eff: float, optional 戻り値: :returns: 等価な平均自由行程 l0 [m]。 :rtype: float """ mstar = m_eff * m0 return tau_pref * np.sqrt(2.0 / mstar) * (e**r) def save_results_to_excel(outfile, meta_rows, data_rows): """ 概要: 計算結果とメタデータをExcelファイルに保存します。 詳細説明: 指定された出力ファイル名で新しいExcelワークブックを作成し、 "metadata"と"transport_vs_Ne"の2つのシートにデータを書き込みます。 "metadata"シートには計算条件などのメタデータが保存され、 "transport_vs_Ne"シートにはキャリア濃度に対する輸送特性の計算結果が保存されます。 引数: :param outfile: 出力するExcelファイル名。例: "results.xlsx"。 :type outfile: str :param meta_rows: メタデータを含む行のリスト。各要素はExcelの1行に対応するリストです。 :type meta_rows: list :param data_rows: 計算結果データを含む行のリスト。各要素はExcelの1行に対応するリストです。 :type data_rows: list 戻り値: :returns: なし :rtype: None """ wb = Workbook() ws_meta = wb.active ws_meta.title = "metadata" for row in meta_rows: ws_meta.append(row) ws = wb.create_sheet("transport_vs_Ne") ws.append([ "r", "xi", "tau_pref_s_per_J^(r-1/2)", "tau_at_kBT_s", "l0_equiv_m", "Ne_m^-3", "Ne_cm^-3", "S_V_per_K", "S_uV_per_K", "L_Wohm_per_K2", "sigma_S_per_m", "mu_m2_per_Vs", "mu_cm2_per_Vs", "kappa_e_W_per_mK", "PF_W_per_mK2", "PF_mW_per_mK2", "ZT", ]) for row in data_rows: ws.append(row) wb.save(outfile) def main(): """ 概要: スクリプトの主要な実行フローを定義します。コマンドライン引数を解析し、 熱電輸送特性を計算・プロット・保存します。 詳細説明: argparse モジュールを使用してコマンドライン引数（温度、有効質量、参照移動度、 格子熱伝導率、キャリア種別、還元フェルミ準位の範囲と点数、散乱因子rのリスト、 出力ファイル名など）を解析します。 次に、指定されたパラメータに基づいて、各散乱因子 r に対して輸送特性（電子濃度、 ゼーベック係数、ローレンツ因子、電気伝導度、移動度、電子熱伝導度、パワーファクター、ZT因子）を計算します。 計算されたデータは save_results_to_excel 関数を使用してExcelファイルに保存されます。 最後に、計算結果は matplotlib を用いてキャリア濃度に対する各特性のグラフとしてプロットされ、表示されます。 戻り値: :returns: なし :rtype: None """ parser = argparse.ArgumentParser( description=( "S, L, sigma, mu, kappa_e, PF, ZT vs Ne using exact transport integrals " "for tau(E) ∝ E^(r-1/2) in a 3D parabolic band." ) ) parser.add_argument("--T", type=float, default=300.0, help="Temperature [K]") parser.add_argument("--m_eff", type=float, default=1.0, help="Effective mass m*/m0") parser.add_argument("--mu_ref", type=float, default=10.0, help="Reference mobility at xi_ref [cm^2/V/s]") parser.add_argument("--xi_ref", type=float, default=0.0, help="Reference reduced Fermi level xi where mu_ref is reproduced") parser.add_argument("--klatt", type=float, default=5.0, help="Lattice thermal conductivity [W/m/K] (default: 5)") parser.add_argument("--carrier", type=str, default="electron", choices=["electron", "hole"], help="Carrier type") parser.add_argument("--xi_min", type=float, default=-5.0, help="Minimum xi") parser.add_argument("--xi_max", type=float, default=40.0, help="Maximum xi") parser.add_argument("--nxi", type=int, default=181, help="Number of xi points") parser.add_argument("--r_list", type=str, default="2.0,1.5,1.0,0.5,0.0,-0.5", help="Comma-separated r list") parser.add_argument("--outfile", type=str, default="transport_vs_Ne.xlsx", help="Excel output filename") args = parser.parse_args() T = args.T m_eff = args.m_eff mu_ref_cm2_Vs = args.mu_ref xi_ref = args.xi_ref klatt = args.klatt carrier = args.carrier r_list = [float(x.strip()) for x in args.r_list.split(",")] xi_range = np.linspace(args.xi_min, args.xi_max, args.nxi) outfile = args.outfile print("Convention: tau(E) ∝ E^(r-1/2)") print(f"T = {T:g} K") print(f"m_eff = {m_eff:g} m0") print(f"mu_ref = {mu_ref_cm2_Vs:g} cm^2/V/s") print(f"xi_ref = {xi_ref:g}") print(f"klatt = {klatt:g} W/m/K") print(f"carrier = {carrier}") print(f"outfile = {outfile}") print() all_results = [] for r in r_list: tau_pref = tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=m_eff) tau_kBT = tau_at_energy(kB * T, tau_pref, r) l0_equiv = equivalent_l0_from_tau_pref(tau_pref, r, m_eff=m_eff) print(f"r = {r:g}") print(f" tau_pref = {tau_pref:.6e} [s / J^({r:g}-1/2)]") print(f" tau(E=kBT) = {tau_kBT:.6e} [s]") print(f" equivalent l0 = {l0_equiv:.6e} [m] (slide convention)") n_list = [] S_list = [] L_list = [] sig_list = [] mu_list = [] ke_list = [] PF_list = [] ZT_list = [] for xi in xi_range: n = electron_density_from_xi(xi, T=T, m_eff=m_eff) # [m^-3] sigma = sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff) # [S/m] mu_SI = mobility_from_sigma_n(sigma, n) # [m^2/V/s] mu_cm2_Vs = mu_SI * 1e4 S = seebeck_from_xi_transport(xi, r, carrier=carrier) # [V/K] L = lorenz_from_xi_transport(xi, r) # [V^2/K^2] kappa_e = L * sigma * T # [W/m/K] PF = S**2 * sigma # [W/m/K^2] ZT = PF * T / (kappa_e + klatt) n_list.append(n) S_list.append(S) L_list.append(L) sig_list.append(sigma) mu_list.append(mu_cm2_Vs) ke_list.append(kappa_e) PF_list.append(PF) ZT_list.append(ZT) all_results.append([ r, xi, tau_pref, tau_kBT, l0_equiv, n, n / 1e6, S, S * 1e6, L, sigma, mu_SI, mu_cm2_Vs, kappa_e, PF, PF * 1e3, ZT, ]) # 配列化して保持 # この行は、ループ内でリストに追加された要素を再度リストに代入するもので、 # 結果的にはall_resultsのリスト構造に影響を与えません。 # 論理的な意味は薄いですが、既存コードの変更禁止ルールに従い保持します。 all_results[-len(xi_range):] = all_results[-len(xi_range):] # ---------- save Excel before plotting ---------- meta_rows = [ ["Convention", "tau(E) ∝ E^(r-1/2)"], ["T_K", T], ["m_eff_over_m0", m_eff], ["mu_ref_cm2_per_Vs", mu_ref_cm2_Vs], ["xi_ref", xi_ref], ["klatt_W_per_mK", klatt], ["carrier", carrier], ["r_list", ",".join(str(r) for r in r_list)], ["xi_min", args.xi_min], ["xi_max", args.xi_max], ["nxi", args.nxi], ] save_results_to_excel(outfile, meta_rows, all_results) print(f"\nSaved calculation results to [{outfile}]\n") # ---------- plot ---------- fig, axes = plt.subplots(4, 2, figsize=figsize) axS = axes[0, 0] axL = axes[0, 1] axSig = axes[1, 0] axMu = axes[1, 1] axKe = axes[2, 0] axPF = axes[2, 1] axZT = axes[3, 0] axBlank = axes[3, 1] for r in r_list: tau_pref = tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=m_eff) n_list = [] S_list = [] L_list = [] sig_list = [] mu_list = [] ke_list = [] PF_list = [] ZT_list = [] for xi in xi_range: n = electron_density_from_xi(xi, T=T, m_eff=m_eff) sigma = sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff) mu_SI = mobility_from_sigma_n(sigma, n) S = seebeck_from_xi_transport(xi, r, carrier=carrier) L = lorenz_from_xi_transport(xi, r) kappa_e = L * sigma * T PF = S**2 * sigma ZT = PF * T / (kappa_e + klatt) n_list.append(n) S_list.append(S) L_list.append(L) sig_list.append(sigma) mu_list.append(mu_SI * 1e4) ke_list.append(kappa_e) PF_list.append(PF) ZT_list.append(ZT) n_arr_cm3 = np.array(n_list) / 1e6 S_arr_uVK = np.array(S_list) * 1e6 L_arr = np.array(L_list) sig_arr = np.array(sig_list) mu_arr = np.array(mu_list) ke_arr = np.array(ke_list) PF_arr = np.array(PF_list) * 1e3 ZT_arr = np.array(ZT_list) label = f"r = {r:g}" axS.plot(n_arr_cm3, S_arr_uVK, label=label) axL.plot(n_arr_cm3, L_arr, label=label) axSig.plot(n_arr_cm3, sig_arr, label=label) axMu.plot(n_arr_cm3, mu_arr, label=label) axKe.plot(n_arr_cm3, ke_arr, label=label) axPF.plot(n_arr_cm3, PF_arr, label=label) axZT.plot(n_arr_cm3, ZT_arr, label=label) for ax in [axS, axL, axSig, axMu, axKe, axPF, axZT]: ax.set_xscale("log") ax.grid(True, which="both", alpha=0.25) ax.set_xlabel(r"$N_e$ (cm$^{-3}$)") axS.axhline(0.0, linewidth=1.0, alpha=0.4) axS.set_ylabel(r"$S$ ($\mu$V/K)") axS.set_title("Seebeck coefficient") axL.set_ylabel(r"$L$ (V$^2$/K$^2$)") axL.set_title("Lorenz factor") axSig.set_ylabel(r"$\sigma$ (S/m)") axSig.set_title("Electrical conductivity") axMu.set_ylabel(r"$\mu$ (cm$^2$/V/s)") axMu.set_title(r"Mobility ($\tau(E)\propto E^{r-1/2}$)") axKe.set_ylabel(r"$\kappa_e$ (W/m/K)") axKe.set_title("Electronic thermal conductivity") axPF.set_ylabel(r"$PF=S^2\sigma$ (mW/m/K$^2$)") axPF.set_title("Power factor") axZT.set_ylabel(r"$ZT$") axZT.set_title(rf"$ZT$ ($\kappa_{{latt}}$={klatt:g} W/m/K)") axBlank.axis("off") handles, labels = axS.get_legend_handles_labels() fig.legend(handles, labels, loc="upper center", ncol=min(len(r_list), 6)) fig.suptitle( f"Transport vs Ne (T={T:g} K, m*/m0={m_eff:g}, " f"mu_ref={mu_ref_cm2_Vs:g} cm$^2$/V/s at xi_ref={xi_ref:g}, " r"$\tau(E)\propto E^{r-1/2}$)", y=0.995 ) plt.tight_layout(rect=[0, 0, 1, 0.96]) plt.show() if __name__ == "__main__": main()