.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/user_guide/nei_populations.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_gallery_user_guide_nei_populations.py: Non-equilibrium Ionization =========================== In this example, we'll show how to compute the ionization fractions of iron in the case where the the timescale of the temperature change is shorter than the ionization timescale. This is often referred to as *non-equilibrium ionization*. .. GENERATED FROM PYTHON SOURCE LINES 10-23 .. code-block:: Python import astropy.units as u import matplotlib.pyplot as plt import numpy as np from astropy.visualization import quantity_support from scipy.interpolate import interp1d from fiasco import Element # sphinx_gallery_thumbnail_number = 3 quantity_support() .. rst-class:: sphx-glr-script-out .. code-block:: none .MplQuantityConverter object at 0x7fc3a25aa810> .. GENERATED FROM PYTHON SOURCE LINES 24-40 The time evolution of an ion fraction :math:`f_k` for some ion :math:`k` for some element is given by, .. math:: \frac{d}{dt}f_k = n_e(\alpha_{k-1}^I f_{k-1} + \alpha_{k+1}^R f_{k+1} - \alpha_{k}^I f_k - \alpha_k^R f_k), where :math:`n_e` is the electron density, :math:`\alpha^I` is the ionization rate, and :math:`\alpha^R` is the recombination. Both the ionization and recombination rates are a function of electron temperature :math:`T_e`. Consider a very simple temperature evolution model in which the electron temperature rises from :math:`10^4` K to :math:`1.5\times10^7` K over 15 s and then decreases back to :math:`10^4` K. .. GENERATED FROM PYTHON SOURCE LINES 40-47 .. code-block:: Python t = np.arange(0, 30, 0.01) * u.s Te_min = 1e4 * u.K Te_max = 2e6 * u.K t_mid = 15 * u.s Te = Te_min + (Te_max - Te_min) / t_mid * t Te[t>t_mid] = Te_max - (Te_max - Te_min) / t_mid * (t[t>t_mid] - t_mid) .. GENERATED FROM PYTHON SOURCE LINES 48-51 Similarly, we'll model the density as increasing from :math:`10^8\,\mathrm{cm}^{-3}` to :math:`1.5^10\,\mathrm{cm}^{-3}` and then back down to :math:`10^8\,\mathrm{cm}^{-3}`. .. GENERATED FROM PYTHON SOURCE LINES 51-57 .. code-block:: Python ne_min = 1e8 * u.cm**(-3) ne_max = 1e10 * u.cm**(-3) ne = ne_min + (ne_max - ne_min) / t_mid * t ne[t>t_mid] = ne_max - (ne_max - ne_min) / t_mid * (t[t>t_mid] - t_mid) .. GENERATED FROM PYTHON SOURCE LINES 58-59 We can plot the temperature and density as a function of time. .. GENERATED FROM PYTHON SOURCE LINES 59-66 .. code-block:: Python plt.figure(figsize=(12,5)) plt.subplot(121) plt.plot(t,Te.to('MK')) plt.subplot(122) plt.plot(t,ne) plt.show() .. image-sg:: /generated/gallery/user_guide/images/sphx_glr_nei_populations_001.png :alt: nei populations :srcset: /generated/gallery/user_guide/images/sphx_glr_nei_populations_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 67-74 In the case where the ionization fraction is in equilibrium, the left hand side of the above equation is zero. As shown in previous examples, we can get the ion population fractions of any element in the CHIANTI database over some temperature array with `~fiasco.Element.equilibrium_ionization`. First, let's model the ionization fractions of C for the above temperature model, assuming ionization equilibrium. .. GENERATED FROM PYTHON SOURCE LINES 74-80 .. code-block:: Python temperature_array = np.logspace(4, 8, 1000) * u.K carbon = Element('carbon', temperature_array) func_interp = interp1d(carbon.temperature.to_value('K'), carbon.equilibrium_ionization.value, axis=0, kind='cubic', fill_value='extrapolate') carbon_ieq = u.Quantity(func_interp(Te.to_value('K'))) .. GENERATED FROM PYTHON SOURCE LINES 81-82 We can plot the population fractions as a function of time. .. GENERATED FROM PYTHON SOURCE LINES 82-89 .. code-block:: Python plt.figure(figsize=(12, 4)) for ion in carbon: plt.plot(t, carbon_ieq[:, ion.charge_state], label=ion.ion_name_roman) plt.xlim(t[[0,-1]].value) plt.legend(ncol=4, frameon=False) plt.show() .. image-sg:: /generated/gallery/user_guide/images/sphx_glr_nei_populations_002.png :alt: nei populations :srcset: /generated/gallery/user_guide/images/sphx_glr_nei_populations_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 90-106 Now, let's compute the population fractions in the case where the population fractions are out of equilibrium. We solve the above equation using the implicit method outlined in Appendix A of :cite:t:`barnes_understanding_2019`, .. math:: \mathbf{F}_{l+1} = \left(\mathbb{I} - \frac{\delta}{2}\mathbf{A}_{l+1}\right)^{-1}\left(\mathbb{I} + \frac{\delta}{2}\mathbf{A}_l\right)\mathbf{F}_l where :math:`\mathbf{F}_l` is the vector of population fractions of all states at time :math:`t_l`, :math:`\mathbf{A}_l` is the matrix of ionization and recombination rates with dimension :math:`Z+1\times Z+1` at time :math:`t_l`, and :math:`\delta` is the time step. First, we can set our initial state to the equilibrium populations. .. GENERATED FROM PYTHON SOURCE LINES 106-109 .. code-block:: Python carbon_nei = np.zeros(t.shape + (carbon.atomic_number + 1,)) carbon_nei[0, :] = carbon_ieq[0,:] .. GENERATED FROM PYTHON SOURCE LINES 110-111 Then, interpolate the rate matrix to our modeled temperatures .. GENERATED FROM PYTHON SOURCE LINES 111-115 .. code-block:: Python func_interp = interp1d(carbon.temperature.to_value('K'), carbon._rate_matrix.value, axis=0, kind='cubic', fill_value='extrapolate') fe_rate_matrix = func_interp(Te.to_value('K')) * carbon._rate_matrix.unit .. GENERATED FROM PYTHON SOURCE LINES 116-118 Finally, we can iteratively compute the non-equilibrium ion fractions using the above equation. .. GENERATED FROM PYTHON SOURCE LINES 118-129 .. code-block:: Python identity = u.Quantity(np.eye(carbon.atomic_number + 1)) for i in range(1, t.shape[0]): dt = t[i] - t[i-1] term1 = identity - ne[i] * dt/2. * fe_rate_matrix[i, ...] term2 = identity + ne[i-1] * dt/2. * fe_rate_matrix[i-1, ...] carbon_nei[i, :] = np.linalg.inv(term1) @ term2 @ carbon_nei[i-1, :] carbon_nei[i, :] = np.fabs(carbon_nei[i, :]) carbon_nei[i, :] /= carbon_nei[i, :].sum() carbon_nei = u.Quantity(carbon_nei) .. GENERATED FROM PYTHON SOURCE LINES 130-131 And plot the non-equilibrium populations as a function of time .. GENERATED FROM PYTHON SOURCE LINES 131-138 .. code-block:: Python plt.figure(figsize=(12,4)) for ion in carbon: plt.plot(t, carbon_nei[:, ion.charge_state], ls='-', label=ion.ion_name_roman,) plt.xlim(t[[0,-1]].value) plt.legend(ncol=4, frameon=False) plt.show() .. image-sg:: /generated/gallery/user_guide/images/sphx_glr_nei_populations_003.png :alt: nei populations :srcset: /generated/gallery/user_guide/images/sphx_glr_nei_populations_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 139-142 We can compare the two equilibrium and non equilbrium cases directly by overplotting the C V population fractions for the two cases. .. GENERATED FROM PYTHON SOURCE LINES 142-149 .. code-block:: Python plt.figure(figsize=(12,4)) plt.plot(t, carbon_ieq[:, 4], label='IEQ') plt.plot(t, carbon_nei[:, 4], label='NEI',) plt.xlim(t[[0,-1]].value) plt.legend(frameon=False) plt.show() .. image-sg:: /generated/gallery/user_guide/images/sphx_glr_nei_populations_004.png :alt: nei populations :srcset: /generated/gallery/user_guide/images/sphx_glr_nei_populations_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 150-159 Note that the equilibrium populations exactly track the temperature evolution while the peak of the non-equilibrium Fe XVI population lags the equilibrium case and is not symmetric. Finally let's plot the effective temperature as described in :cite:t:`bradshaw_numerical_2009`. Qualitatively, this is the temperature that one would infer if they observed the non-equilibrium populations and assumed the populations were in equilibrium. .. GENERATED FROM PYTHON SOURCE LINES 159-167 .. code-block:: Python Te_eff = carbon.temperature[[(np.fabs(carbon.equilibrium_ionization - carbon_nei[i, :])).sum(axis=1).argmin() for i in range(carbon_nei.shape[0])]] plt.plot(t, Te.to('MK'), label=r'$T$') plt.plot(t, Te_eff.to('MK'), label=r'$T_{\mathrm{eff}}$') plt.xlim(t[[0,-1]].value) plt.legend() plt.show() .. image-sg:: /generated/gallery/user_guide/images/sphx_glr_nei_populations_005.png :alt: nei populations :srcset: /generated/gallery/user_guide/images/sphx_glr_nei_populations_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 168-173 Note that the effective temperature is consistently less than the actual temperature, emphasizing the point that non-equilibrium ionization makes it difficult to detect very hot plasma when the temperature changes rapidly and/or the density is low. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.866 seconds) .. _sphx_glr_download_generated_gallery_user_guide_nei_populations.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: nei_populations.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: nei_populations.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_