Merge pull request #285 from SNEWS2/mori-axion-2023

Mori axion 2023

python/snewpy/models/ccsn.py

@@ -742,6 +742,7 @@ def __new__(cls, filename=None, *, progenitor_mass=None):
     # Populate Docstring with abbreviated param values
     __new__.__doc__ = __new__.__doc__.format(**_param_abbrv)
 
+
 class Fornax_2022(_RegistryModel):
     """Model based on 2D simulations of 100 progenitors from Tianshu Wang, David Vartanyan, Adam Burrows, and Matthew S.B. Coleman, MNRAS 517:543, 2022.
        Data available at https://www.astro.princeton.edu/~burrows/nu-emissions.2d.large/
@@ -804,6 +805,77 @@ def __new__(cls, *, progenitor=None, progenitor_mass=None):
     # Populate Docstring with abbreviated param values
     __new__.__doc__ = __new__.__doc__.format(**_param_abbrv)
 
+
+class Mori_2023(_RegistryModel):
+    """Model based on 2D simulations with axionlike particles, K. Mori, T.  Takiwaki, K. Kotake and S. Horiuchi, Phys. Rev. D 108:063027, 2023. All models are based on the non-rotating 20 M_sun solar metallicity progenitor model from S.E. Woolsey and A. Heger, Phys. Rep. 442:269, 2007. Data from private communication.
+        """
+    param = {'axion_mass' : [0, 100, 200]<<u.MeV,
+             'axion_coupling' : [0, 2, 4, 6, 8, 10, 12, 14, 16, 20]<<(1e-10/u.GeV) }
+
+    _param_validator = lambda p: \
+        (p['axion_mass'] == 0 and p['axion_coupling'] == 0) or \
+        (p['axion_mass'].to_value('MeV') == 100 and p['axion_coupling'].to_value('1e-10/GeV') in (2,4,10,12,14,16,20)) or \
+        (p['axion_mass'].to_value('MeV') == 200 and p['axion_coupling'].to_value('1e-10/GeV') in (2,4,6,8,10,20))
+
+    _param_abbrv = {
+        'axion_mass': '[0, 100, 200] MeV',
+        'axion_coupling' : '[0..20] 1e-10/GeV' }
+
+    def __new__(cls, axion_mass, axion_coupling):
+        """Model Initialization.
+
+        Parameters
+        ----------
+        axion_mass: int
+            Axion mass in units of MeV. Valid values are {axion_mass}.
+        axion_coupling: int
+            Axion-photon coupling, in units of 1e-10/GeV. Valid values are {axion_coupling}.
+        """
+        # Make sure axion coupling is converted to units 1e-10/GeV:
+        axion_coupling = np.round(axion_coupling.to('1e-10/GeV'))
+
+        # Load from Parameters
+        cls.check_valid_params(cls, axion_mass=axion_mass, axion_coupling=axion_coupling)
+
+        if axion_mass == 0:
+            filename = 't-prof_std.dat'
+        else:
+            filename = f't-prof_{axion_mass.to_value("MeV"):g}_{axion_coupling.to_value("1e-10/GeV"):g}.dat'
+
+        # PNS mass table, from Mori+ 2023.
+        mpns = { (0, 0):    1.78,
+                 (100, 2):  1.77,
+                 (100, 4):  1.76,
+                 (100, 10): 1.77,
+                 (100, 12): 1.77,
+                 (100, 14): 1.77,
+                 (100, 16): 1.77,
+                 (100, 20): 1.74,
+                 (200, 2):  1.77,
+                 (200, 4):  1.76,
+                 (200, 6):  1.75,
+                 (200, 8):  1.74,
+                 (200, 10): 1.73,
+                 (200, 20): 1.62 }
+
+        am = int(axion_mass.to_value('MeV'))
+        ac = int(axion_coupling.to_value('1e-10/GeV'))
+        pns_mass = mpns[(am,ac)]
+
+        # Set the metadata.
+        metadata = {
+            'Axion mass': axion_mass,
+            'Axion coupling': axion_coupling,
+            'Progenitor mass': 20*u.Msun,
+            'PNS mass' : pns_mass*u.Msun
+            }
+
+        return loaders.Mori_2023(filename, metadata)
+
+    # Populate Docstring with abbreviated param values
+    __new__.__doc__ = __new__.__doc__.format(**_param_abbrv)
+
+
 class SNOwGLoBES:
     """A model that does not inherit from SupernovaModel (yet) and imports a group of SNOwGLoBES files."""
 

python/snewpy/models/loaders.py

@@ -749,3 +749,59 @@ def __init__(self, filename, metadata={}):
             # Note factor of 0.25 in nu_x and nu_x_bar.
             factor = 1. if flavor.is_electron else 0.25
             self.luminosity[flavor] = np.sum(dLdE*dE, axis=1) * factor * 1e50 * u.erg/u.s
+
+
+class Mori_2023(PinchedModel):
+    def __init__(self, filename, metadata={}):
+        """
+        Parameters
+        ----------
+        filename : str
+            Absolute or relative path to file prefix.
+        """
+        # Set up model metadata
+        self.axion_mass = metadata['Axion mass']
+        self.axion_coupling = metadata['Axion coupling']
+        self.progenitor_mass = metadata['Progenitor mass']
+        self.pns_mass = metadata['PNS mass']
+
+        self.metadata = metadata
+
+        # Open the requested filename using the model downloader.
+        datafile = _model_downloader.get_model_data(self.__class__.__name__, filename)
+
+        # Read ASCII data.
+        simtab = Table.read(datafile, format='ascii')
+
+        # Remove the first table row, which appears to have zero input.
+        simtab = simtab[simtab['1:t_sim[s]'] > 0]
+
+        # Get grid of model times.
+        simtab['TIME'] = simtab['2:t_pb[s]'] << u.s
+        for j, (f, fkey) in enumerate(zip([Flavor.NU_E, Flavor.NU_E_BAR, Flavor.NU_X], 'ebx')):
+            simtab[f'L_{f.name}'] = simtab[f'{6+j}:Le{fkey}[e/s]'] << u.erg / u.s
+            # Compute the pinch parameter from E_rms and E_avg
+            # <E^2> / <E>^2 = (2+a)/(1+a), where
+            # E_rms^2 = <E^2> - <E>^2.
+            Eavg = simtab[f'{9+j}:Em{fkey}[MeV]']
+            Erms = simtab[f'{12+j}:Er{fkey}[MeV]']
+            E2 = Erms**2 + Eavg**2
+            x = E2 / Eavg**2
+            alpha = (2-x) / (x-1)
+
+            simtab[f'E_{f.name}'] = Eavg << u.MeV
+            simtab[f'E2_{f.name}'] = E2 << u.MeV**2
+            simtab[f'ALPHA_{f.name}'] = alpha
+
+#            simtab[f'E_{f.name}'] = simtab[f'{9+j}:Em{fkey}[MeV]'] << u.MeV
+#            Erms = simtab[f'{12+j}:Er{fkey}[MeV]'] * u.MeV
+#
+#            # Compute the pinch parameter from E_rms and E_avg
+#            simtab[f'E2_{f.name}'] = Erms**2 + simtab[f'E_{f.name}']**2
+#            x = simtab[f'E2_{f.name}'] / simtab[f'E_{f.name}']**2
+#            simtab[f'ALPHA_{f.name}'] = (2-x) / (x-1)
+
+        self.filename = os.path.basename(filename)
+
+        super().__init__(simtab, metadata)
+

python/snewpy/models/model_files.yml

@@ -50,3 +50,6 @@ models:
     Zha_2021: 
         repository: *snewpy
 
+    Mori_2023: 
+        repository: *ccsnmodels
+

python/snewpy/test/test_01_registry.py

@@ -9,7 +9,7 @@
                           Sukhbold_2015, Bollig_2016, Walk_2018, \
                           Walk_2019, Fornax_2019, Warren_2020, \
                           Kuroda_2020, Fornax_2021, Zha_2021, \
-                          Fornax_2022
+                          Fornax_2022, Mori_2023
 
 from astropy import units as u
 
@@ -33,6 +33,7 @@ def test_param_combinations(self):
                   Fornax_2021,
                   Zha_2021,
                   Fornax_2022,
+                  Mori_2023,
                   ]
 
         for model in models:
@@ -394,3 +395,36 @@ def test_Fornax_2022(self):
             self.assertEqual(len(f), len(Flavor))
             self.assertEqual(f[Flavor.NU_E].unit, 1/(u.erg * u.s))
 
+    def test_Mori_2023(self):
+        """
+        Instantiate a set of 'Mori 2023' models
+        """
+        axion_mass_coupling = [ (0,0),
+            (100,2), (100,4), (100,10), (100,12), (100,14), (100,16), (100,20),
+            (200,2), (200,4), (200,6), (200,8), (200,10), (200,20)]
+
+        pns_mass = [1.78, 1.77, 1.76, 1.77, 1.77, 1.77, 1.77, 1.74, 1.77, 1.76, 1.75, 1.74, 1.73, 1.62] * u.Msun
+
+        for ((am, ac), mpns) in zip(axion_mass_coupling, pns_mass):
+            model = Mori_2023(axion_mass=am*u.MeV, axion_coupling=ac*1e-10/u.GeV)
+
+            axion_mass = float(am) * u.MeV
+            self.assertEqual(model.metadata['Axion mass'], axion_mass)
+
+            axion_coupling = float(ac) * 1e-10/u.GeV
+            axion_coupling = np.round(axion_coupling.to('1e-10/GeV'))
+            self.assertEqual(model.metadata['Axion coupling'], axion_coupling)
+
+            self.assertEqual(model.metadata['Progenitor mass'], 20*u.Msun)
+            self.assertEqual(model.metadata['PNS mass'], mpns)
+
+            # Check that times are in proper units.
+            t = model.get_time()
+            self.assertTrue(t.unit, u.s)
+
+            # Check that we can compute flux dictionaries.
+            f = model.get_initial_spectra(0*u.s, 10*u.MeV)
+            self.assertEqual(type(f), dict)
+            self.assertEqual(len(f), len(Flavor))
+            self.assertEqual(f[Flavor.NU_E].unit, 1./(u.erg * u.s))
+