diff --git a/docs/references.bib b/docs/references.bib index 72ffdca8..fb76e507 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -157,6 +157,18 @@ @article{dere_ionization_2007 doi = {10.1051/0004-6361:20066728} } +@ARTICLE{drake_twophoton_1986, + title = "{Spontaneous two-photon decay rates in hydrogenlike and heliumlike ions}", + author = {{Drake}, G.~W.~F.}, + journal = {Phys. Rev. A}, + year = {1986}, + month = oct, + volume = {34}, + number = {4}, + pages = {2871-2880}, + doi = {10.1103/PhysRevA.34.2871} +} + @article{feldman_potential_1992, title = {The Potential for Plasma Diagnostics from Stellar Extreme-Ultraviolet Observations}, author = {Feldman, U. and Mandelbaum, P. and Seely, J. F. and Doschek, G. A. and Gursky, H.}, @@ -181,6 +193,16 @@ @article{fontes_fully_1999 doi = {10.1103/PhysRevA.59.1329} } +@ARTICLE{gronenschild_twophoton_1978, + author = {{Gronenschild}, E.~H.~B.~M. and {Mewe}, R.}, + title = "{Calculated X-radiation from optically thin plasmas. III. Abundance effects on continuum emission.}", + journal = {Astronomy and Astrophysics Supplement Series}, + year = {1978}, + month = may, + volume = {32}, + pages = {283-305} +} + @article{itoh_relativistic_2000, title = {Relativistic {{Thermal Bremsstrahlung Gaunt Factor}} for the {{Intracluster Plasma}}. {{II}}. {{Analytic Fitting Formulae}}}, author = {Itoh, Naoki and Sakamoto, Tsuyoshi and Kusano, Shugo and Nozawa, Satoshi and Kohyama, Yasuharu}, diff --git a/fiasco/collections.py b/fiasco/collections.py index 6d736486..b63485dd 100644 --- a/fiasco/collections.py +++ b/fiasco/collections.py @@ -167,6 +167,48 @@ def free_bound(self, wavelength: u.angstrom, **kwargs): free_bound += fb * abundance * ioneq[:, np.newaxis] return free_bound + @u.quantity_input + def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**-3, **kwargs): + r""" + Compute the two-photon continuum emission. + + .. note:: Both abundance and ionization equilibrium are included here. + + The combined two-photon continuum is given by + + .. math:: + + P_{2p}(\lambda,T,n_{e}) = \sum_{X,k}\mathrm{Ab}(X)f(X_{k})C_{2p, X_k}(\lambda,T,n_{e}) + + where :math:`\mathrm{Ab}(X)` is the abundance of element :math:`X`, + :math:`f(X_{k})` is the ionization equilibrium of the emitting ion :math:`X_{k}`, + and :math:`C_{fb, X_k}(\lambda,T)` is the two-photon emission of the + ion :math:`X_k` as computed by `fiasco.Ion.two_photon`. + The sum is taken over all ions in the collection. + + Parameters + ---------- + wavelength : `~astropy.units.Quantity` + electron_density: `~astropy.units.Quantity` + + See Also + -------- + fiasco.Ion.two_photon + """ + wavelength = np.atleast_1d(wavelength) + electron_density = np.atleast_1d(electron_density) + two_photon = u.Quantity(np.zeros(wavelength.shape + self.temperature.shape + electron_density.shape), + 'erg cm^3 s^-1 Angstrom^-1') + for ion in self: + try: + tp = ion.two_photon(wavelength, electron_density, **kwargs) + except MissingDatasetException as e: + self.log.warning(f'{ion.ion_name} not included in two-photon emission. {e}') + continue + else: + two_photon += tp * ion.abundance * ion.ioneq[:, np.newaxis] + return two_photon + @u.quantity_input def spectrum(self, density: u.cm**(-3), emission_measure: u.cm**(-5), wavelength_range=None, bin_width=None, kernel=None, **kwargs): diff --git a/fiasco/io/sources/continuum_sources.py b/fiasco/io/sources/continuum_sources.py index 8deb7d0f..80840c90 100644 --- a/fiasco/io/sources/continuum_sources.py +++ b/fiasco/io/sources/continuum_sources.py @@ -252,15 +252,22 @@ def to_hdf5(self, hf, df, **kwargs): class HSeqParser(GenericParser): - """ + r""" Parameters for calculating two-photon continuum for hydrogen-like ions + + Notes + ----- + * The parameter :math: `\psi_{\text{norm}}` (called :math: `A_{\text{sum}}` in CHIANTI) is a normalization factor of + the integral of the spectral distribution function :math:`\psi(y)` from 0 to 1, such + that :math: `\frac{1}{\psi_{\text{norm}}} \int_{0}^{1} \psi(y) dy = 2`. This normalization is only used For + hydrogenic ions. """ filetype = 'hseq_2photon' dtypes = [int, float, int, float, float, float] - units = [None, u.dimensionless_unscaled, None, 1/u.s, 1/u.s, u.dimensionless_unscaled] - headings = ['Z', 'y', 'Z_0', 'A', 'A_sum', 'psi'] + units = [None, u.dimensionless_unscaled, None, 1/u.s, u.dimensionless_unscaled, u.dimensionless_unscaled] + headings = ['Z', 'y', 'Z_0', 'A', 'psi_norm', 'psi'] descriptions = ['atomic number', 'fraction of energy carried by one of the two photons', - 'nominal atomic number', 'radiative decay rate', 'summed radiative decay rate', + 'nominal atomic number', 'radiative decay rate', 'normalization of the integral of psi from 0 to 1', 'spectral distribution function'] def __init__(self, filename, **kwargs): diff --git a/fiasco/ions.py b/fiasco/ions.py index d5903b74..2ed0f28f 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -6,7 +6,7 @@ import numpy as np from functools import cached_property -from scipy.interpolate import interp1d, PchipInterpolator, splev, splrep +from scipy.interpolate import CubicSpline, interp1d, PchipInterpolator, splev, splrep from scipy.ndimage import map_coordinates from fiasco import proton_electron_ratio @@ -1418,7 +1418,7 @@ def free_bound(self, C_{fb}(\lambda, T) = \frac{2}{hc^3(k_B m_e)^{3/2}\sqrt{2\pi}}\frac{E^5}{T^{3/2}}\sum_i\frac{\omega_i}{\omega_0}\sigma_i^{\mathrm{bf}}\exp{\left(-\frac{E-I_i}{k_BT}\right)} where :math:`E` is the energy of the outgoing photon, - :math:`\omega_i,\omega_0` are the statastical weights of the + :math:`\omega_i,\omega_0` are the statistical weights of the :math:`i`-th level of the recombined ion and the ground level of the recombining ion, respectively, :math:`\sigma_i^{\mathrm{bf}}` is the free-bound cross-section, and :math:`I_i` is the energy required to ionize the recombined ion from level :math:`i`. @@ -1527,3 +1527,90 @@ def _karzas_cross_section(self, photon_energy: u.erg, ionization_energy: u.erg, cross_section = prefactor * ionization_energy**2 * photon_energy**(-3) * gaunt_factor / n cross_section[np.where(photon_energy < ionization_energy)] = 0.*cross_section.unit return cross_section + + @needs_dataset('elvlc') + @u.quantity_input + def two_photon(self, + wavelength: u.angstrom, + electron_density: u.cm**(-3), + include_protons=False) -> u.Unit('erg cm3 s-1 Angstrom-1'): + r""" + Two-photon continuum emission of a hydrogenic or helium-like ion. + + .. note:: Does not include ionization equilibrium or abundance. + Unlike the equivalent IDL routine, the output here is not + expressed per steradian and as such the factor of + :math:`1/4\pi` is not included. + + In hydrogen-like ions, the transition :math:`2S_{1/2} \rightarrow 1S_{1/2} + h\nu` cannot occur + as an electric dipole transition, but only as a much slower magnetic dipole transition. + The dominant transition then becomes :math:`2S_{1/2} \rightarrow 1S_{1/2} + h\nu_{1} + h\nu_{2}`. + + In helium-like ions, the transition from :math:`1s2s ^{1}S_{0} \rightarrow 1s^{2}\ ^{1}S_{0} + h\nu` + is forbidden under quantum selection rules since :math:`\Delta J = 0`. + Similarly, the dominant transition becomes + :math:`1s2s ^{1}S_{0} \rightarrow 1s^{2}\ ^{1}S_{0} + h\nu_{1} + h\nu_{2}`. + + In both cases, the energy of the two photons emitted equals the energy difference of the two levels. + As a consequence, no photons can be emitted beneath the rest wavelength for a given transition. + See the introduction of :cite:t:`drake_1986` for a concise description of the process. + + The emission is given by + + .. math:: + + C_{2p}(\lambda, T, n_{e}) = hc \frac{n_{j}(X^{+m}) A_{ji} \lambda_{0} \psi(\frac{\lambda_{0}}{\lambda})}{\psi_{\text{norm}}\lambda^{3}} + + where :math:`\lambda_{0}` is rest wavelength of the (disallowed) transition, + :math:`A_{ji}` is the Einstein spontaneous emission coefficient, + :math:`\psi` is so-called spectral distribution function, given approximately by + :math:`\psi(y) \approx 2.623 \sqrt{\cos{\Big(\pi(y-\frac{1}{2})\Big)}}` :cite:p:`gronenschild_twophoton_1978`, + :math:`\psi_{\text{norm}}` is a normalization factor for hydrogen-like ions such + that :math:`\frac{1}{\psi_{\text{norm}}} \int_{0}^{1} \psi(y) dy = 2` (and 1 for helium-like ions), + and :math:`n_{j}(X^{+m})` is the density of ions m of element X in excited state j, given by + :math:`n_{j}(X^{+m}) = \frac{n_{j}(X^{+m})}{n(X^{+m})} \frac{n(X^{+m})}{n(X)} \frac{n(X)}{n_{H}} \frac{n_{H}}{n_{e}} n_{e}`. + + Parameters + ---------- + wavelength : `~astropy.units.Quantity` + electron_density : `~astropy.units.Quantity` + include_protons : `bool`, optional + If True, use proton excitation and de-excitation rates in the level population calculation. + """ + wavelength = np.atleast_1d(wavelength) + electron_density = np.atleast_1d(electron_density) + + final_shape = wavelength.shape + self.temperature.shape + electron_density.shape + if self.hydrogenic: + A_ji = self._hseq['A'] + psi_norm = self._hseq['psi_norm'] + cubic_spline = CubicSpline(self._hseq['y'], self._hseq['psi']) + config = '2s' # Get the index of the 2S1/2 state for H-like + J = 0.5 + elif self.helium_like: + A_ji = self._heseq['A'] + psi_norm = 1.0 * u.dimensionless_unscaled + cubic_spline = CubicSpline(self._heseq['y'], self._heseq['psi']) + config = '1s.2s' # Get the index of the 1s2s 1S0 state for He-like: + J = 0 + else: + return u.Quantity(np.zeros(final_shape), 'erg cm^3 s^-1 Angstrom^-1') + level_index = np.where((self._elvlc['config'] == config) & (np.isclose(self._elvlc['J'], J)) )[0][0] + + E_obs = self._elvlc['E_obs'][level_index] + E_th = self._elvlc['E_th'][level_index] + E_2p = E_obs if E_obs > 0.0 else E_th + rest_wavelength = 1 / E_2p + + psi_interp = cubic_spline((rest_wavelength / wavelength).decompose()) + psi_interp = np.where(wavelength < rest_wavelength, 0.0, psi_interp) + + energy_dist = (A_ji * rest_wavelength * psi_interp) / (psi_norm * wavelength**3) + + level_population = self.level_populations(electron_density, include_protons=include_protons) + level_population = level_population[..., level_index] + + level_density = level_population / electron_density + matrix = np.outer(energy_dist, level_density).reshape(*final_shape) + + return const.h * const.c * matrix diff --git a/fiasco/tests/idl/helpers.py b/fiasco/tests/idl/helpers.py index 021fc731..03d0b718 100644 --- a/fiasco/tests/idl/helpers.py +++ b/fiasco/tests/idl/helpers.py @@ -46,7 +46,7 @@ def run_idl_script(idl_env, script, input_args, save_vars, file_name, version, f file_name: `str` Name of the IDL results file version: `packaging.version.Version`, `str` - Version of IDL used to generate these test results + Version of CHIANTI used to generate these test results format_func: `dict`, optional Functions to use to format output from the IDL function. This is most useful for adding the necessary units to any of the outputs. diff --git a/fiasco/tests/test_collections.py b/fiasco/tests/test_collections.py index 8eeebb7c..d3d39742 100644 --- a/fiasco/tests/test_collections.py +++ b/fiasco/tests/test_collections.py @@ -106,6 +106,19 @@ def test_free_bound(another_collection, wavelength): index_t = 24 # This is approximately where the ioneq for Fe V peaks assert u.allclose(fb[index_t, index_w], 3.057781475607237e-36 * u.Unit('erg cm3 s-1 Angstrom-1')) +@pytest.mark.parametrize('wavelength', [wavelength, wavelength[450]]) +def test_two_photon(collection, wavelength, hdf5_dbase_root): + # add Li III to the test to include an ion that throws a MissingDatasetException + collection = collection + fiasco.Ion('Li III', collection.temperature, hdf5_dbase_root=hdf5_dbase_root) + tp = collection.two_photon(wavelength, electron_density = 1e10 * u.cm**(-3)) + if wavelength.shape: + assert tp.shape == wavelength.shape + temperature.shape + (1, ) + else: + assert tp.shape == (1, ) + temperature.shape + (1, ) + index_w = 450 if wavelength.shape else 0 + index_t = 30 + # This value has not been checked for correctness + assert u.allclose(tp[index_w, index_t, 0], 3.48580645e-27 * u.Unit('erg cm3 s-1 Angstrom-1')) def test_radiative_loss(collection): rl = collection.radiative_loss(1e9*u.cm**(-3)) diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 4b92aac1..c92ba317 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -32,6 +32,16 @@ def fe10(hdf5_dbase_root): return fiasco.Ion('Fe 10', temperature, hdf5_dbase_root=hdf5_dbase_root) +@pytest.fixture +def c4(hdf5_dbase_root): + return fiasco.Ion('C IV', temperature, hdf5_dbase_root=hdf5_dbase_root) + + +@pytest.fixture +def c5(hdf5_dbase_root): + return fiasco.Ion('C V', temperature, hdf5_dbase_root=hdf5_dbase_root) + + @pytest.fixture def c6(hdf5_dbase_root): return fiasco.Ion('C VI', temperature, hdf5_dbase_root=hdf5_dbase_root) @@ -323,6 +333,23 @@ def test_free_bound(ion): # This value has not been tested for correctness assert u.allclose(emission[0, 0], 9.7902609e-26 * u.cm**3 * u.erg / u.Angstrom / u.s) +# The two-photon test currently fails for dbase_version >= 9 because it is missing c_5.reclvl +@pytest.mark.requires_dbase_version('<=','8.0.7') +def test_two_photon(c4, c5, c6): + # test both H-like and He-like ions, and one that doesn't have two-photon emission + c4_emission = c4.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3)) + c5_emission = c5.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3)) + c6_emission = c6.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3)) + c6_emission_protons = c6.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3), include_protons=True) + assert c4_emission.shape == (1, ) + c4.temperature.shape + (1, ) + assert c5_emission.shape == (1, ) + c5.temperature.shape + (1, ) + assert c6_emission.shape == (1, ) + c6.temperature.shape + (1, ) + assert c6_emission_protons.shape == (1, ) + c6.temperature.shape + (1, ) + assert u.allclose(c4_emission[0, 30, 0], 0.0 * u.cm**3 * u.erg / u.Angstrom / u.s) + # These values have not been tested for correctness + assert u.allclose(c5_emission[0, 30, 0], 4.04634243e-25 * u.cm**3 * u.erg / u.Angstrom / u.s) + assert u.allclose(c6_emission[0, 30, 0], 8.25316887e-26 * u.cm**3 * u.erg / u.Angstrom / u.s) + assert u.allclose(c6_emission_protons[0, 30, 0], 6.79615958e-29 * u.cm**3 * u.erg / u.Angstrom / u.s) def test_free_bound_no_recombining(h1): # This is test the case where there is no data available for the recombining diff --git a/fiasco/util/data/test_file_list.json b/fiasco/util/data/test_file_list.json index 5f2c5296..2fb71a32 100644 --- a/fiasco/util/data/test_file_list.json +++ b/fiasco/util/data/test_file_list.json @@ -169,12 +169,19 @@ "c_4.diparams", "c_4.drparams", "c_4.easplups", + "c_4.elvlc", "c_4.fblvl", "c_4.rrparams", + "c_4.scups", + "c_4.wgfa", "c_5.diparams", "c_5.drparams", + "c_5.elvlc", "c_5.fblvl", + "c_5.reclvl", "c_5.rrparams", + "c_5.scups", + "c_5.wgfa", "c_6.diparams", "c_6.drparams", "c_6.elvlc", @@ -811,6 +818,8 @@ "he_2.elvlc", "he_2.fblvl", "he_2.rrparams", + "he_2.scups", + "he_2.wgfa", "he_3.rrparams", "heseq_2photon.dat", "hseq_2photon.dat",