Split ThreeFlavorDecoherence

python/snewpy/flavor_transformation/__init__.py

@@ -16,9 +16,6 @@
 .. autoclass:: snewpy.flavor_transformation.CompleteExchange
     :members:
 
-.. autoclass:: snewpy.flavor_transformation.ThreeFlavorDecoherence
-    :members:
-
 .. autoclass:: snewpy.flavor_transformation.TransformationChain
     :members:
 
@@ -71,6 +68,7 @@ def construct(mixing_params:ThreeFlavorMixingParameters|FourFlavorMixingParamete
 AdiabaticMSWes = construct_chain(in_sn.AdiabaticMSWes)
 NonAdiabaticMSWes = construct_chain(in_sn.NonAdiabaticMSWes)
 TwoFlavorDecoherence = construct_chain(in_sn.TwoFlavorDecoherence)
+ThreeFlavorDecoherence = construct_chain(in_sn.ThreeFlavorDecoherence)
 NeutrinoDecay = construct_chain(in_sn.AdiabaticMSW, 
                                 in_vacuum.NeutrinoDecay)
 QuantumDecoherence = construct_chain(in_sn.AdiabaticMSW, 
@@ -114,12 +112,3 @@ def P(f1, f2):
         return P
 
 
-class ThreeFlavorDecoherence(FlavorTransformation):
-    """Equal mixing of all threen eutrino matter states and antineutrino matter states"""
-
-    def P_ff(self, t, E):
-        """Equal mixing so Earth matter has no effect"""
-        @FlavorMatrix.from_function(ThreeFlavor)
-        def P(f1, f2):
-            return (f1.is_neutrino == f2.is_neutrino)*1/3.
-        return P

python/snewpy/flavor_transformation/in_sn.py

@@ -8,7 +8,7 @@
 import numpy as np
 
 from snewpy.flavor import FlavorMatrix
-from snewpy.neutrino import MassHierarchy
+from snewpy.neutrino import MassHierarchy, ThreeFlavorMixingParameters, FourFlavorMixingParameters
 from dataclasses import dataclass
 
 from .base import ThreeFlavorTransformation, FourFlavorTransformation
@@ -22,16 +22,15 @@ class SNTransformation(ABC):
     """Base class for all transformations in SN"""
     @abstractmethod
     def P_mf(self, t, E)->FlavorMatrix:
-        pass
-
-
+        pass        
 
 ###############################################################################
 class AdiabaticMSW(SNTransformation, ThreeFlavorTransformation):
     """Adiabatic MSW effect."""
 
     def P_mf(self, t, E):
         return self.mixing_params.Pmf_HighDensityLimit()
+
 ###############################################################################
 class NonAdiabaticMSWH(SNTransformation, ThreeFlavorTransformation):
     """Nonadiabatic MSW H resonance. 
@@ -49,6 +48,7 @@ def P_mf(self, t, E):
         else:
             Pmf[['NU_1_BAR','NU_3_BAR'],:] = Pmf[['NU_3_BAR','NU_1_BAR'],:]
         return Pmf
+
 ###############################################################################
 class TwoFlavorDecoherence(SNTransformation, ThreeFlavorTransformation):
     """Equal mixing of whatever two matter states form the MSW H resonance.
@@ -72,6 +72,20 @@ def P_mf(self, t, E):
 
         return Pmf
 
+###############################################################################
+class ThreeFlavorDecoherence(SNTransformation, ThreeFlavorTransformation):
+    """Equal mixing of all three neutrino states, and all three antineutrinos
+
+    ThreeFlavorDecoherence is relevant when the size of the density fluctuations is >= 10% for densities 
+    around the H resonance density —see Kneller (2010); Kneller & Mauney (2013). 
+    """
+
+    def P_mf(self, t, E):
+        @FlavorMatrix.from_function(ThreeFlavorMixingParameters.basis_mass, ThreeFlavorMixingParameters.basis_flavor)
+        def P(f1, f2):
+            return (f1.is_neutrino == f2.is_neutrino)*1/3.
+        return P
+
 ###############################################################################       
 class MSWEffect(SNTransformation, ThreeFlavorTransformation):
     """The MSW effect using a density profile and electron 
@@ -187,4 +201,15 @@ def P_mf(self, t, E):
             Pmf[['NU_2','NU_4'],:] = Pmf[['NU_4','NU_2'],:]
             Pmf[['NU_1_BAR','NU_2_BAR','NU_4_BAR'],:] = Pmf[['NU_2_BAR','NU_4_BAR','NU_1_BAR'],:]
 
-        return Pmf
+        return Pmf
+
+###############################################################################
+class FourFlavorDecoherence(SNTransformation, FourFlavorTransformation):
+    """Equal mixing of all four neutrino states, and all four antineutrinos
+    """
+
+    def P_mf(self, t, E):
+        @FlavorMatrix.from_function(FourFlavorMixingParameters.basis_mass, FourFlavorMixingParameters.basis_flavor)
+        def P(f1, f2):
+            return (f1.is_neutrino == f2.is_neutrino)*1/4.
+        return P  

python/snewpy/snowglobes.py

@@ -71,7 +71,17 @@ def generate_time_series(model_path, model_type, transformation_type, d, output_
     model_class = getattr(snewpy.models.ccsn, model_type)
 
     # Choose flavor transformation. Use dict to associate the transformation name with its class.
-    flavor_transformation_dict = {'NoTransformation': NoTransformation(), 'AdiabaticMSW_NMO': AdiabaticMSW(mh=MassHierarchy.NORMAL), 'AdiabaticMSW_IMO': AdiabaticMSW(mh=MassHierarchy.INVERTED), 'NonAdiabaticMSWH_NMO': NonAdiabaticMSWH(mh=MassHierarchy.NORMAL), 'NonAdiabaticMSWH_IMO': NonAdiabaticMSWH(mh=MassHierarchy.INVERTED), 'TwoFlavorDecoherence': TwoFlavorDecoherence(), 'ThreeFlavorDecoherence': ThreeFlavorDecoherence(), 'NeutrinoDecay_NMO': NeutrinoDecay(mh=MassHierarchy.NORMAL), 'NeutrinoDecay_IMO': NeutrinoDecay(mh=MassHierarchy.INVERTED), 'QuantumDecoherence_NMO': QuantumDecoherence(mh=MassHierarchy.NORMAL), 'QuantumDecoherence_IMO': QuantumDecoherence(mh=MassHierarchy.INVERTED)}
+    flavor_transformation_dict = {'NoTransformation': NoTransformation(), 
+                                  'AdiabaticMSW_NMO': AdiabaticMSW(NMO), 
+                                  'AdiabaticMSW_IMO': AdiabaticMSW(IMO), 
+                                  'NonAdiabaticMSWH_NMO': NonAdiabaticMSWH(NMO), 
+                                  'NonAdiabaticMSWH_IMO': NonAdiabaticMSWH(IMO), 
+                                  'TwoFlavorDecoherence': TwoFlavorDecoherence(NMO), 
+                                  'ThreeFlavorDecoherence': ThreeFlavorDecoherence(NMO), 
+                                  'NeutrinoDecay_NMO': NeutrinoDecay(NMO), 
+                                  'NeutrinoDecay_IMO': NeutrinoDecay(IMO), 
+                                  'QuantumDecoherence_NMO': QuantumDecoherence(NMO), 
+                                  'QuantumDecoherence_IMO': QuantumDecoherence(IMO)}
     flavor_transformation = flavor_transformation_dict[transformation_type]
 
     model_dir, model_file = os.path.split(os.path.abspath(model_path))
@@ -140,7 +150,7 @@ def generate_fluence(model_path, model_type, transformation_type, d, output_file
                                   'NonAdiabaticMSWH_NMO': NonAdiabaticMSWH(NMO), 
                                   'NonAdiabaticMSWH_IMO': NonAdiabaticMSWH(IMO), 
                                   'TwoFlavorDecoherence': TwoFlavorDecoherence(NMO), 
-                                  'ThreeFlavorDecoherence': ThreeFlavorDecoherence(), 
+                                  'ThreeFlavorDecoherence': ThreeFlavorDecoherence(NMO), 
                                   'NeutrinoDecay_NMO': NeutrinoDecay(NMO), 
                                   'NeutrinoDecay_IMO': NeutrinoDecay(IMO), 
                                   'QuantumDecoherence_NMO': QuantumDecoherence(NMO), 

python/snewpy/test/test_04_xforms.py

@@ -424,20 +424,16 @@ def test_ThreeFlavorDecoherence(self):
         """
         Three flavor decoherence
         """
-        P = xforms.ThreeFlavorDecoherence().P_ff(t,E)
+        P = xforms.ThreeFlavorDecoherence(MixingParameters(MassHierarchy.NORMAL)).P_ff(t,E)
         #convert to TwoFlavor case
-        P = (TwoFlavor<<P<<TwoFlavor)
-
-        assert np.isclose(P['e','e'], 1./3)
-        assert np.isclose(P['e','x'] , 2./3)
-        assert np.isclose(P['x','x'] , 2./3)
-        assert np.isclose(P['x','e'] , 1./3)
-
-        assert np.isclose(P['e_bar','e_bar'] , 1./3)
-        assert np.isclose(P['e_bar','x_bar'] , 2./3)
-        assert np.isclose(P['x_bar','x_bar'] , 2./3)
-        assert np.isclose(P['x_bar','e_bar'] , 1./3)
-
+        P = (ThreeFlavor<<P<<ThreeFlavor)
+        for orig_flv in ThreeFlavor:
+            for dest_flv in ThreeFlavor:
+                if orig_flv.is_neutrino == dest_flv.is_neutrino:
+                    assert np.isclose(P[orig_flv,dest_flv], 1/3)
+                else:
+                    assert np.isclose(P[orig_flv,dest_flv], 0)
+
     def test_NeutrinoDecay_NMO(self, t, E):
         """
         Neutrino decay with NMO and new mixing angles