add cooling documentation

docs/cooling_notes.md

@@ -0,0 +1,102 @@
+---
+title: Cooling function conventions
+author: Ben Wibking
+date: February 12, 2024
+---
+
+# Summary
+
+Several different conventions exist for the definition of a collisional ionization equilibrium (CIE) cooling curve. They are explained in detail below, with references.
+
+# Conventions
+
+## $n^2 \Lambda$ convention
+
+This convention uses the *total number density* $n$ to define the volumetric cooling rate. The total particle number density is defined as the sum of the number of ions, neutrals, and electrons divided by the volume.
+
+The total number density $n$ is related to the dimensionless mean molecular weight $\mu$ by:
+
+$$n = \frac{\rho}{\mu m_H} \, ,$$ {#eq:number_density}
+
+where $\rho$ is the mass density of the fluid, $m_H$ is the mass of the hydrogen atom, and $\mu$ is the dimensionless mean molecular weight. This may be taken as the definition of $\mu$.
+
+In this convention, the volumetric cooling rate is given by
+
+$$S_{cool} = n^2 \Lambda(T) \, ,$$ {#eq:cooling_rate_nsqLambda}
+
+where $S_{cool}$ is the volumetric cooling source term, $n$ is the total number density and $\Lambda(T)$ is the CIE cooling curve in this convention.
+
+This convention is used by the PLUTO code (see the [PLUTO User Guide](http://plutocode.ph.unito.it/userguide.pdf), section 9.2, equation 9.5).
+
+## $n_i n_e \Lambda$ convention
+
+In this convention, the volumetric cooling rate is given by
+
+$$S_{cool} = n_i n_e \Lambda(T) \, ,$$ {#eq:cooling_rate_nineLambda}
+
+where $S_{cool}$ is the volumetric cooling source term, $n_i$ is the number density of ions (*usually* including neutrals by convention, but this difference should be negligible when CIE is valid), $n_e$ is the number density of electrons, and $\Lambda(T)$ is the CIE cooling curve in this convention.
+
+This convention is adopted by Sutherland and Dopita (1993) in their tabulation of CIE cooling curves (their section 5.2; see also their equations 51-58).
+
+Sutherland and Dopita additionally comment: *"The normalization factor for the cooling function $\Lambda_N$ is taken as $n_t n_e$ where $n_t$ is the total number density of ions and $n_e$ is the total number density of electrons. This definition is generally valid and independent of composition. It differs from the commonly used factor $n_H n_e$ or even $n_e^2$ commonly used in hydrogen-rich plasmas."*
+
+
+## $n_H n_e \Lambda$ convention
+
+In this convention, the volumetric cooling rate is given by
+
+$$S_{cool} = n_H n_e \Lambda(T) \, ,$$ {#eq:cooling_rate_nHneLambda}
+
+where $S_{cool}$ is the volumetric cooling source term, $n_H$ is the number density of hydrogen *nuclei* (not ions), $n_e$ is the number density of electrons, and $\Lambda(T)$ is the CIE cooling curve in this convention.
+
+This convention is adopted by Schure et al. (2009) in their definition of $\Lambda_N$. However, this is **not** the same $\Lambda_N$ as used by Sutherland and Dopita (1993)!
+
+This convention is adopted by Gnat and Sternberg (2007) in their definition of $\Lambda$.
+
+## $n_H^2 \Lambda$ convention
+
+In this convention, the volumetric cooling rate is given by
+
+$$S_{cool} = n_H^2 \Lambda(T) \, ,$$ {#eq:cooling_rate_nHsqLambda}
+
+where $S_{cool}$ is the volumetric cooling source term, $n_H$ is the number density of hydrogen *nuclei* (not ions), and $\Lambda(T)$ is the CIE cooling curve in this convention.
+
+This convention is adopted by Dalgarno and McCray (1972) (their equation 24), and Schure et al. (2009) in their definition of $\Lambda_{hd}$ (their equation 2).
+
+### AthenaPK implementation
+
+This is the convention adopted by AthenaPK, except $n_H$ is computed assuming zero metals in the mass budget:
+$$ x_H = 1 - Y \, , $$
+$$ n_H = X_H \rho / m_H  = (1 - Y) \rho / m_H \, .$$
+
+The difference in neglecting metals in the mass budget in computing $n_H$ is only a $\sim 2$ percent effect on $n_H$ (since metals are 2 percent by mass for solar metallicity) and a $\sim 4$ percent effect on the source term.
+
+The Schure table in the GitHub repository uses the using the $n_H n_e$ convention, rather than the $n_H^2$ convention. This causes an error in the cooling rate of $n_e / n_H$, or about 20 percent for temperatures above $10^{5.5}$ K.
+
+Above $\sim 10^4$ K, there is a percent-level error in the mean molecular weight due to neglecting the metal contribution to the electron fraction. *Users should be aware that at low absolute electron fraction, metals contribute a substantial proportion of the metals, so this approximation is no longer valid below $\sim 10^4$ K.* The temperature is computed as:
+$$ T = \left( \frac{\rho}{\mu m_H} \right)^{-1} \left( \frac{P}{k_B} \right) \, ,$$
+where the dimensionless mean molecular weight $\mu$ is computed assuming a zero metallicity gas:
+$$ \mu = \left( \frac{3}{4} Y + 2(1 - Y) \right)^{-1} = \left( 2 - \frac{5}{4} Y \right)^{-1} \, ,$$
+where $Y$ is the Helium mass fraction.
+
+For a more precise calculation of $\mu$ assuming a scaled-solar composition of metals, we have:
+$$ \mu = \left( 1 - x_e \right) \left( X + \frac{Y}{4} + \frac{Z}{\bar A} \right)^{-1} \, , $$
+where $\bar A$ is the mass-weighted mean atomic weight of metals:
+$$ \bar A = 16 \, .$$
+
+Additionally, for full ionization for all species, we have:
+$$ x_e = \frac{X + 2 Y (m_H/m_{He}) + (\bar A/2) Z (m_H/m_{\bar Z})}{2 X + 3 Y (m_H/m_{He}) + (1 + \bar A/2) Z (m_H/m_{\bar Z})} \approx \frac{X + Y/2 + Z/2}{2 X + 3 Y/4 + (1/{\bar A} + 1/2) Z}$$
+
+For $Y = 0.24$ and $Z = 0.02$, the mean molecular weight for complete ionization is $\mu \approx 0.60$ (if $Z = 0$, then this yields $\mu \approx 0.59$).
+
+## $n_e^2 \Lambda$ convention
+
+In this convention, the volumetric cooling rate is given by
+
+$$S_{cool} = n_e^2 \Lambda(T) \, ,$$ {#eq:cooling_rate_nHsqLambda}
+
+where $S_{cool}$ is the volumetric cooling source term, $n_e$ is the number density of electrons, and $\Lambda(T)$ is the CIE cooling curve in this convention.
+
+For the case of a hydrogen-only fully-ionized plasma *only*, this is equivalent to the $n_H^2 \Lambda$ convention, since each hydrogen atom contributes exactly one electron and there are no other sources of electrons.
+
+I am not aware of any common tables that use this convention. However, it is noted as one of the possible conventions by Sutherland and Dopita (1993).

docs/plot_lambda.py

@@ -0,0 +1,52 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+if __name__ == "__main__":
+    ## plot *solar abundance* CIE cooling curves
+
+    ## Sutherland and Dopita (1993) [NOTE: n_i n_e convention for Lambda_N, where n_t == n_i in their notation]
+    # Columns: Log(T), n_e, n_H, n_t, log(\Lambda_net), log(\Lambda_N)
+    sd_LogT, sd_ne, sd_nH, sd_nt, sd_logLambdaNet, sd_logLambdaN = np.loadtxt("sutherland_dopita_table_6.txt", unpack=True)
+    # convert to Lambda_hd:
+    sd_LambdaN = 10.**sd_logLambdaN
+    sd_LambdaHD = (sd_ne * sd_nt * sd_LambdaN) / sd_nH**2
+    # convert to nH ne Lambda(T):
+    sd_nHneLambda = (sd_ne * sd_nt * sd_LambdaN) / (sd_nH * sd_ne)
+
+    ## Schure et al. (2009) [NOTE: n_H n_e convention for Lambda_N; n_H^2 convention for Lambda_hd]
+    # Columns: log T (K), log Λ_N (erg s^−1 cm^3), log Λ_hd (erg s^−1 cm^3), n_e/n_H
+    # BEWARE: This Lambda_N is NOT the same as the Sutherland and Dopita Lambda_N!
+    s_LogT, s_logLambdaN, s_logLambdaHD, s_ne_over_nH = np.loadtxt("schure_table_2.txt", unpack=True)
+    s_LambdaHD = 10.**s_logLambdaHD
+
+    ## Gnat and Sternberg (2007) [NOTE: n_H n_e convention for Lambda]
+    # Columns: Temperature Lambda(Z=1e-3) Lambda(Z=1e-2) Lambda(Z=1e-1) Lambda(Z=1) Lambda(Z=2)
+    gs_T, gs_LambdaZem3, gs_LambdaZem2, gs_LambdaZem1, gs_LambdaZ1, gs_LambdaZ2 = np.loadtxt("gnat_sternberg_cie_table.txt", unpack=True, skiprows=23)
+    # NOTE: There is not enough information in this table alone to compute n_e!
+    #  (The electron fraction must be computed from datafile2a.txt)
+
+    # plot nH^2-normalized Lambda(T)
+    plt.figure()
+    plt.plot(sd_LogT, sd_LambdaHD, label = r"Sutherland & Dopita $\Lambda_{hd}$ for $Z_{\odot}$")
+    plt.plot(s_LogT, s_LambdaHD, label = r"Schure et al. $\Lambda_{hd}$ for $Z_{\odot}$")
+    # (Gnat & Sternberg cannot be converted to this convention without the electron fraction, which is missing.)
+    plt.yscale('log')
+    plt.ylim(1e-23, 3e-21)
+    plt.xlabel(r"$\log T$ (K)")
+    plt.ylabel(r"$\Lambda(T)$ [$n_H^2$ convention]")
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig("cooling_curves_LambdaHD.png")
+
+    # plot nH ne-normalized Lambda(T)
+    plt.figure()
+    plt.plot(sd_LogT, sd_nHneLambda, label = r"Sutherland & Dopita $\Lambda$ for $Z_{\odot}$")
+    plt.plot(s_LogT, 10.**s_logLambdaN, label = r"Schure et al. $\Lambda$ for $Z_{\odot}$")
+    plt.plot(np.log10(gs_T), gs_LambdaZ1, label = r"Gnat & Sternberg $\Lambda$ for $Z_{\odot}$")
+    plt.yscale('log')
+    plt.ylim(1e-23, 3e-21)
+    plt.xlabel(r"$\log T$ (K)")
+    plt.ylabel(r"$\Lambda(T)$ [$n_H n_e$ convention]")
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig("cooling_curves_nHneLambda.png")

docs/schure_table_2.txt

@@ -0,0 +1,112 @@
+## This is an exact reproduction of "Table 2. Cooling curve for solar metallicity" from Schure et al. (2009).
+## log T (K), log Λ_N (erg s^−1 cm^3), log Λ_hd (erg s^−1 cm^3), n_e/n_H
+3.8 -25.7331 -30.6104 1.3264e-05 
+3.84 -25.0383 -29.4107 4.2428e-05 
+3.88 -24.4059 -28.4601 8.8276e-05 
+3.92 -23.8288 -27.5743 0.00017967 
+3.96 -23.3027 -26.3766 0.00084362 
+4.0 -22.8242 -25.289 0.0034295 
+4.04 -22.3917 -24.2684 0.013283 
+4.08 -22.0067 -23.3834 0.042008 
+4.12 -21.6818 -22.5977 0.12138 
+4.16 -21.4529 -21.9689 0.30481 
+4.2 -21.3246 -21.5972 0.53386 
+4.24 -21.3459 -21.4615 0.76622 
+4.28 -21.4305 -21.4789 0.89459 
+4.32 -21.5293 -21.5497 0.95414 
+4.36 -21.6138 -21.6211 0.98342 
+4.4 -21.6615 -21.6595 1.0046 
+4.44 -21.6551 -21.6426 1.0291 
+4.48 -21.5919 -21.5688 1.0547 
+4.52 -21.5092 -21.4771 1.0767 
+4.56 -21.4124 -21.3755 1.0888 
+4.6 -21.3085 -21.2693 1.0945 
+4.64 -21.2047 -21.1644 1.0972 
+4.68 -21.1067 -21.0658 1.0988 
+4.72 -21.0194 -20.9778 1.1004 
+4.76 -20.9413 -20.8986 1.1034 
+4.8 -20.8735 -20.8281 1.1102 
+4.84 -20.8205 -20.77 1.1233 
+4.88 -20.7805 -20.7223 1.1433 
+4.92 -20.7547 -20.6888 1.1638 
+4.96 -20.7455 -20.6739 1.1791 
+5.0 -20.7565 -20.6815 1.1885 
+5.04 -20.782 -20.7051 1.1937 
+5.08 -20.8008 -20.7229 1.1966 
+5.12 -20.7994 -20.7208 1.1983 
+5.16 -20.7847 -20.7058 1.1993 
+5.2 -20.7687 -20.6896 1.1999 
+5.24 -20.759 -20.6797 1.2004 
+5.28 -20.7544 -20.6749 1.2008 
+5.32 -20.7505 -20.6709 1.2012 
+5.36 -20.7545 -20.6748 1.2015 
+5.4 -20.7888 -20.7089 1.202 
+5.44 -20.8832 -20.8031 1.2025 
+5.48 -21.045 -20.9647 1.203 
+5.52 -21.2286 -21.1482 1.2035 
+5.56 -21.3737 -21.2932 1.2037 
+5.6 -21.4573 -21.3767 1.2039 
+5.64 -21.4935 -21.4129 1.204 
+5.68 -21.5098 -21.4291 1.2041 
+5.72 -21.5345 -21.4538 1.2042 
+5.76 -21.5863 -21.5055 1.2044 
+5.8 -21.6548 -21.574 1.2045 
+5.84 -21.7108 -21.63 1.2046 
+5.88 -21.7424 -21.6615 1.2047 
+5.92 -21.7576 -21.6766 1.2049 
+5.96 -21.7696 -21.6886 1.205 
+6.0 -21.7883 -21.7073 1.2051 
+6.04 -21.8115 -21.7304 1.2053 
+6.08 -21.8303 -21.7491 1.2055 
+6.12 -21.8419 -21.7607 1.2056 
+6.16 -21.8514 -21.7701 1.2058 
+6.2 -21.869 -21.7877 1.206 
+6.24 -21.9057 -21.8243 1.2062 
+6.28 -21.969 -21.8875 1.2065 
+6.32 -22.0554 -21.9738 1.2067 
+6.36 -22.1488 -22.0671 1.207 
+6.4 -22.2355 -22.1537 1.2072 
+6.44 -22.3084 -22.2265 1.2075 
+6.48 -22.3641 -22.2821 1.2077 
+6.52 -22.4033 -22.3213 1.2078 
+6.56 -22.4282 -22.3462 1.2079 
+6.6 -22.4408 -22.3587 1.208 
+6.64 -22.4443 -22.3622 1.2081 
+6.68 -22.4411 -22.359 1.2082 
+6.72 -22.4334 -22.3512 1.2083 
+6.76 -22.4242 -22.342 1.2083 
+6.8 -22.4164 -22.3342 1.2084 
+6.84 -22.4134 -22.3312 1.2084 
+6.88 -22.4168 -22.3346 1.2085 
+6.92 -22.4267 -22.3445 1.2085 
+6.96 -22.4418 -22.3595 1.2086 
+7.0 -22.4603 -22.378 1.2086 
+7.04 -22.483 -22.4007 1.2087 
+7.08 -22.5112 -22.4289 1.2087 
+7.12 -22.5449 -22.4625 1.2088 
+7.16 -22.5819 -22.4995 1.2088 
+7.2 -22.6177 -22.5353 1.2089 
+7.24 -22.6483 -22.5659 1.2089 
+7.28 -22.6719 -22.5895 1.2089 
+7.32 -22.6883 -22.6059 1.2089 
+7.36 -22.6985 -22.6161 1.2089 
+7.4 -22.7032 -22.6208 1.209 
+7.44 -22.7037 -22.6213 1.209 
+7.48 -22.7008 -22.6184 1.209 
+7.52 -22.695 -22.6126 1.209 
+7.56 -22.6869 -22.6045 1.209 
+7.6 -22.6769 -22.5945 1.209 
+7.64 -22.6655 -22.5831 1.209 
+7.68 -22.6531 -22.5707 1.209 
+7.72 -22.6397 -22.5573 1.209 
+7.76 -22.6258 -22.5434 1.209 
+7.8 -22.6111 -22.5287 1.209 
+7.84 -22.5964 -22.514 1.209 
+7.88 -22.5816 -22.4992 1.209 
+7.92 -22.5668 -22.4844 1.209 
+7.96 -22.5519 -22.4695 1.209 
+8.0 -22.5367 -22.4543 1.209 
+8.04 -22.5216 -22.4392 1.209 
+8.08 -22.5062 -22.4237 1.2091 
+8.12 -22.4912 -22.4087 1.2091 
+8.16 -22.4753 -22.3928 1.2091 

docs/sutherland_dopita_table_6.txt

@@ -0,0 +1,97 @@
+## This is a reproduction of Table 6 of Sutherland and Dopita (1993).
+## [NOTE: Only columns 1-6 are reproduced in this file. Log(U), Log(tau_cool), and the "general plasma properties" are not included.]
+## [Headings included in the printed table:]
+##   CIE Cooling Curve : Solar Abundances
+##   Number Densities; Cooling Properties; General Plasma Properties
+## Log(T), n_e, n_H, n_t, log(\Lambda_net), log(\Lambda_N)
+4.0 0.0023 1.0 1.099 -25.65 -23.06 
+4.05 0.0142 1.0 1.099 -24.27 -22.46 
+4.1 0.0704 1.0 1.099 -23.28 -22.17 
+4.15 0.2534 1.0 1.099 -22.47 -21.92 
+4.2 0.5678 1.0 1.099 -21.99 -21.79 
+4.25 0.8177 1.0 1.099 -21.84 -21.8 
+4.3 0.9324 1.0 1.099 -21.85 -21.86 
+4.35 0.978 1.0 1.099 -21.87 -21.9 
+4.4 1.004 1.0 1.099 -21.84 -21.88 
+4.45 1.035 1.0 1.099 -21.77 -21.82 
+4.5 1.067 1.0 1.099 -21.66 -21.73 
+4.55 1.086 1.0 1.099 -21.56 -21.63 
+4.6 1.094 1.0 1.099 -21.45 -21.53 
+4.65 1.098 1.0 1.099 -21.34 -21.42 
+4.7 1.1 1.0 1.099 -21.24 -21.32 
+4.75 1.102 1.0 1.099 -21.14 -21.22 
+4.8 1.11 1.0 1.099 -21.05 -21.14 
+4.85 1.128 1.0 1.099 -20.97 -21.07 
+4.9 1.154 1.0 1.099 -20.91 -21.01 
+4.95 1.176 1.0 1.099 -20.87 -20.98 
+5.0 1.188 1.0 1.099 -20.87 -20.99 
+5.05 1.195 1.0 1.099 -20.9 -21.02 
+5.1 1.198 1.0 1.099 -20.91 -21.03 
+5.15 1.199 1.0 1.099 -20.89 -21.01 
+5.2 1.2 1.0 1.099 -20.86 -20.98 
+5.25 1.201 1.0 1.099 -20.85 -20.97 
+5.3 1.201 1.0 1.099 -20.84 -20.96 
+5.35 1.201 1.0 1.099 -20.84 -20.96 
+5.4 1.202 1.0 1.099 -20.87 -20.99 
+5.45 1.203 1.0 1.099 -21.01 -21.13 
+5.5 1.203 1.0 1.099 -21.23 -21.35 
+5.55 1.204 1.0 1.099 -21.43 -21.55 
+5.6 1.204 1.0 1.099 -21.54 -21.66 
+5.65 1.204 1.0 1.099 -21.58 -21.71 
+5.7 1.204 1.0 1.099 -21.59 -21.71 
+5.75 1.204 1.0 1.099 -21.59 -21.71 
+5.8 1.204 1.0 1.099 -21.64 -21.76 
+5.85 1.205 1.0 1.099 -21.74 -21.86 
+5.9 1.205 1.0 1.099 -21.81 -21.93 
+5.95 1.205 1.0 1.099 -21.83 -21.95 
+6.0 1.205 1.0 1.099 -21.84 -21.96 
+6.05 1.205 1.0 1.099 -21.84 -21.96 
+6.1 1.206 1.0 1.099 -21.83 -21.96 
+6.15 1.206 1.0 1.099 -21.82 -21.95 
+6.2 1.206 1.0 1.099 -21.82 -21.94 
+6.25 1.206 1.0 1.099 -21.85 -21.97 
+6.3 1.207 1.0 1.099 -21.95 -22.07 
+6.35 1.207 1.0 1.099 -22.08 -22.2 
+6.4 1.207 1.0 1.099 -22.19 -22.31 
+6.45 1.208 1.0 1.099 -22.27 -22.39 
+6.5 1.208 1.0 1.099 -22.32 -22.44 
+6.55 1.208 1.0 1.099 -22.35 -22.48 
+6.6 1.208 1.0 1.099 -22.38 -22.5 
+6.65 1.208 1.0 1.099 -22.41 -22.53
+6.7 1.208 1.0 1.099 -22.44 -22.56 
+6.75 1.208 1.0 1.099 -22.46 -22.59
+6.8 1.208 1.0 1.099 -22.48 -22.6 
+6.85 1.209 1.0 1.099 -22.48 -22.6 
+6.9 1.209 1.0 1.099 -22.46 -22.59 
+6.95 1.209 1.0 1.099 -22.45 -22.57 
+7.0 1.209 1.0 1.099 -22.45 -22.57 
+7.05 1.209 1.0 1.099 -22.46 -22.59 
+7.1 1.209 1.0 1.099 -22.49 -22.62 
+7.15 1.209 1.0 1.099 -22.53 -22.65 
+7.2 1.209 1.0 1.099 -22.56 -22.68 
+7.25 1.209 1.0 1.099 -22.58 -22.7 
+7.3 1.209 1.0 1.099 -22.6 -22.72 
+7.35 1.209 1.0 1.099 -22.61 -22.73 
+7.4 1.209 1.0 1.099 -22.61 -22.73 
+7.45 1.209 1.0 1.099 -22.61 -22.73 
+7.5 1.209 1.0 1.099 -22.6 -22.73 
+7.55 1.209 1.0 1.099 -22.6 -22.72 
+7.6 1.209 1.0 1.099 -22.59 -22.71 
+7.65 1.209 1.0 1.099 -22.37 -22.7 
+7.7 1.209 1.0 1.099 -22.36 -22.68 
+7.75 1.209 1.0 1.099 -22.34 -22.67 
+7.8 1.209 1.0 1.099 -22.33 -22.65 
+7.85 1.209 1.0 1.099 -22.51 -22.64 
+7.9 1.209 1.0 1.099 -22.49 -22.62 
+7.95 1.209 1.0 1.099 -22.48 -22.6 
+8.0 1.209 1.0 1.099 -22.46 -22.58 
+8.05 1.209 1.0 1.099 -22.44 -22.36
+8.1 1.209 1.0 1.099 -22.42 -22.34
+8.15 1.209 1.0 1.099 -22.4 -22.33
+8.2 1.209 1.0 1.099 -22.38 -22.51 
+8.25 1.209 1.0 1.099 -22.36 -22.49 
+8.3 1.209 1.0 1.099 -22.34 -22.47 
+8.35 1.209 1.0 1.099 -22.32 -22.45 
+8.4 1.209 1.0 1.099 -22.3 -22.43 
+8.45 1.209 1.0 1.099 -22.28 -22.4 
+8.5 1.209 1.0 1.099 -22.26 -22.38