Beta esc

src/py21cmfast/src/HaloBox.c

@@ -75,9 +75,13 @@ void set_halo_properties(double halo_mass, double M_turn_a, double M_turn_m,
     }
 
     // no rng for escape fraction yet
-    fesc = fmin(consts->fesc_10 * pow(halo_mass / 1e10, consts->alpha_esc), 1);
+    // add redshift dependence to fesc
+    fesc = fmin(consts->fesc_10 * pow(halo_mass / 1e10, consts->alpha_esc) * consts->zesc_power_law,
+                1);
     if (astro_options_global->USE_MINI_HALOS)
-        fesc_mini = fmin(consts->fesc_7 * pow(halo_mass / 1e7, consts->alpha_esc), 1);
+        fesc_mini = fmin(
+            consts->fesc_7 * pow(halo_mass / 1e7, consts->alpha_esc) * consts->zesc_power_law_mini,
+            1);
 
     n_ion_sample =
         stellar_mass * consts->pop2_ion * fesc + stellar_mass_mini * consts->pop3_ion * fesc_mini;
@@ -111,10 +115,14 @@ int get_uhmf_averages(double M_min, double M_max, double M_turn_a, double M_turn
 
     double prefactor_sfr = prefactor_stars / consts->t_star / t_h;
     double prefactor_sfr_mini = prefactor_stars_mini / consts->t_star / t_h;
-    double prefactor_nion = prefactor_stars * consts->fesc_10 * consts->pop2_ion;
-    double prefactor_nion_mini = prefactor_stars_mini * consts->fesc_7 * consts->pop3_ion;
-    double prefactor_wsfr = prefactor_sfr * consts->fesc_10 * consts->pop2_ion;
-    double prefactor_wsfr_mini = prefactor_sfr_mini * consts->fesc_7 * consts->pop3_ion;
+    double prefactor_nion =
+        prefactor_stars * consts->fesc_10 * consts->pop2_ion * consts->zesc_power_law;
+    double prefactor_nion_mini =
+        prefactor_stars_mini * consts->fesc_7 * consts->pop3_ion * consts->zesc_power_law_mini;
+    double prefactor_wsfr =
+        prefactor_sfr * consts->fesc_10 * consts->pop2_ion * consts->zesc_power_law;
+    double prefactor_wsfr_mini =
+        prefactor_sfr_mini * consts->fesc_7 * consts->pop3_ion * consts->zesc_power_law_mini;
 
     double mass_intgrl;
     double intgrl_fesc_weighted, intgrl_stars_only;

src/py21cmfast/src/IonisationBox.c

@@ -161,8 +161,9 @@ void set_ionbox_constants(double redshift, double prev_redshift, struct IonBoxCo
     consts->T_re = astro_params_global->T_RE;
 
     if (matter_options_global->SOURCE_MODEL > 0) {
-        consts->ion_eff_factor_gl = sc.pop2_ion * sc.fstar_10 * sc.fesc_10;
-        consts->ion_eff_factor_mini_gl = sc.pop3_ion * sc.fstar_7 * sc.fesc_7;
+        consts->ion_eff_factor_gl = sc.pop2_ion * sc.fstar_10 * sc.fesc_10 * sc.zesc_power_law;
+        consts->ion_eff_factor_mini_gl =
+            sc.pop3_ion * sc.fstar_7 * sc.fesc_7 * sc.zesc_power_law_mini;
     } else {
         consts->ion_eff_factor_gl = astro_params_global->HII_EFF_FACTOR;
         consts->ion_eff_factor_mini_gl = 0.;

src/py21cmfast/src/_inputparams_wrapper.h

@@ -97,6 +97,8 @@ typedef struct AstroParams {
     // Escape Fraction
     float F_ESC10;
     float ALPHA_ESC;
+    float BETA_ESC;
+    float BETA_ESC_MINI;
     float F_ESC7_MINI;
 
     float T_RE;

src/py21cmfast/src/debugging.c

@@ -126,15 +126,17 @@ void writeAstroParams(AstroParams *p) {
         "       ALPHA_STAR=%8.3f\n"
         "       F_ESC10=%8.3f\n"
         "       ALPHA_ESC=%8.3f\n"
+        "       BETA_ESC=%8.3f\n"
+        "       BETA_ESC_MINI=%8.3f\n"
         "       t_STAR=%8.3f\n"
         "       L_X=%10.3e\n"
         "       NU_X_THRESH=%8.3f\n"
         "       X_RAY_SPEC_INDEX=%8.3f,\n"
         "       UPPER_STELLAR_TURNOVER_MASS=%10.3e\n"
         "       UPPER_STELLAR_TURNOVER_INDEX=%8.3e\n",
-        p->M_TURN, p->R_BUBBLE_MAX, p->F_STAR10, p->ALPHA_STAR, p->F_ESC10, p->ALPHA_ESC, p->t_STAR,
-        p->L_X, p->NU_X_THRESH, p->X_RAY_SPEC_INDEX, p->UPPER_STELLAR_TURNOVER_MASS,
-        p->UPPER_STELLAR_TURNOVER_INDEX);
+        p->M_TURN, p->R_BUBBLE_MAX, p->F_STAR10, p->ALPHA_STAR, p->F_ESC10, p->ALPHA_ESC,
+        p->BETA_ESC, p->BETA_ESC_MINI, p->t_STAR, p->L_X, p->NU_X_THRESH, p->X_RAY_SPEC_INDEX,
+        p->UPPER_STELLAR_TURNOVER_MASS, p->UPPER_STELLAR_TURNOVER_INDEX);
     LOG_INFO(
         "\n        HaloCatalog AstroParams:\n"
         "       SIGMA_STAR=%8.3f\n"

src/py21cmfast/src/hmf.c

@@ -1038,7 +1038,6 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double lnM_co
         .delta = delta2,
         .gamma_type = -3,
     };
-
     if (lnM1 >= lnM_cond) return 0.;
     // return 1 halo at the condition mass if delta is exceeded and the condition is within the
     // integral limits

src/py21cmfast/src/scaling_relations.c

@@ -23,7 +23,8 @@ void print_sc_consts(ScalingConstants *c) {
     LOG_DEBUG("SHMR: f10 %.2e a %.2e f7 %.2e a_mini %.2e sigma %.2e", c->fstar_10, c->alpha_star,
               c->fstar_7, c->alpha_star_mini, c->sigma_star);
     LOG_DEBUG("Upper: a_upper %.2e pivot %.2e", c->alpha_upper, c->pivot_upper);
-    LOG_DEBUG("FESC: f10 %.2e a %.2e f7 %.2e", c->fesc_10, c->alpha_esc, c->fesc_7);
+    LOG_DEBUG("FESC: f10 %.2e a %.2e b %.2e f7 %.2e zesc_pl %.2e", c->fesc_10, c->alpha_esc,
+              c->beta_esc, c->beta_esc_mini, c->fesc_7, c->zesc_power_law, c->zesc_power_law_mini);
     LOG_DEBUG("SSFR: t* %.2e th %.8e sigma %.2e idx %.2e", c->t_star, c->t_h, c->sigma_sfr_lim,
               c->sigma_sfr_idx);
     LOG_DEBUG("Turnovers (nofb) ACG %.2e MCG %.2e Upper %.2e", c->mturn_a_nofb, c->mturn_m_nofb,
@@ -63,9 +64,15 @@ void set_scaling_constants(double redshift, ScalingConstants *consts, bool use_p
     consts->sigma_xray = astro_params_global->SIGMA_LX;
 
     consts->alpha_esc = astro_params_global->ALPHA_ESC;
+    consts->beta_esc = astro_params_global->BETA_ESC;
+    consts->beta_esc_mini = astro_params_global->BETA_ESC_MINI;
     consts->fesc_10 = astro_params_global->F_ESC10;
     consts->fesc_7 = astro_params_global->F_ESC7_MINI;
 
+    // Calculate the redshift-dependent power-law term for escape fraction once per redshift
+    consts->zesc_power_law = pow((1 + redshift) / 8.0, consts->beta_esc);
+    consts->zesc_power_law_mini = pow((1 + redshift) / 8.0, consts->beta_esc_mini);
+
     if (use_photoncons) {
         if (astro_options_global->PHOTON_CONS_TYPE == 2)
             consts->alpha_esc = get_fesc_fit(redshift);
@@ -108,6 +115,9 @@ ScalingConstants evolve_scaling_constants_sfr(ScalingConstants *sc) {
     sc_sfrd.fesc_10 = 1.;
     sc_sfrd.fesc_7 = 1.;
     sc_sfrd.alpha_esc = 0.;
+    sc_sfrd.beta_esc = 0.;
+    sc_sfrd.beta_esc_mini = 0.;
+    sc_sfrd.zesc_power_law = 1.;
     sc_sfrd.Mlim_Fesc = 0.;
     sc_sfrd.Mlim_Fesc_mini = 0.;
 
@@ -149,6 +159,10 @@ ScalingConstants evolve_scaling_constants_to_redshift(double redshift, ScalingCo
         sc_z.mturn_m_nofb = lyman_werner_threshold(redshift, 0., sc_z.vcb_norel);
     }
 
+    // Update the redshift-dependent power-law term for escape fraction
+    sc_z.zesc_power_law = pow((1 + redshift) / 8.0, sc_z.beta_esc);
+    sc_z.zesc_power_law_mini = pow((1 + redshift) / 8.0, sc_z.beta_esc_mini);
+
     return sc_z;
 }
 

src/py21cmfast/src/scaling_relations.h

@@ -36,10 +36,15 @@ typedef struct ScalingConstants {
 
     double fesc_10;
     double alpha_esc;
+    double beta_esc;
     double fesc_7;
+    double beta_esc_mini;
     double pop2_ion;
     double pop3_ion;
 
+    double zesc_power_law;
+    double zesc_power_law_mini;
+
     double vcb_norel;
     double acg_thresh;
     double mturn_a_nofb;

src/py21cmfast/wrapper/inputs.py

@@ -1103,11 +1103,11 @@ def _USE_MINI_HALOS_vald(self, att, val):
 
     @PHOTON_CONS_TYPE.validator
     def _PHOTON_CONS_TYPE_vld(self, att, val):
-        """Raise an error if using PHOTON_CONS_TYPE='z_photoncons' and USE_MINI_HALOS is True."""
+        """Raise an error if using PHOTON_CONS_TYPE='z-photoncons' and USE_MINI_HALOS is True."""
         if self.USE_MINI_HALOS and val == "z-photoncons":
             raise ValueError(
                 "USE_MINI_HALOS is not compatible with the redshift-based"
-                " photon conservation corrections (PHOTON_CONS_TYPE=='z_photoncons')! "
+                " photon conservation corrections (PHOTON_CONS_TYPE=='z-photoncons')! "
             )
 
     @USE_EXP_FILTER.validator
@@ -1182,6 +1182,12 @@ class AstroParams(InputStruct):
     ALPHA_ESC : float, optional
         Power-law index of escape fraction as a function of halo mass. See Sec 2.1 of
         Park+2018.
+    BETA_ESC : float, optional
+        Power-law index of escape fraction as a function of redshift. See Eq. 2 of
+        Qin+2025.
+    BETA_ESC_MINI : float, optional
+        Power-law index of escape fraction as a function of redshift for minihalos.
+        If the scaling relations are not provided explicitly, we extend the ACG ones by default.
     M_TURN : float, optional
         Turnover mass (in log10 solar mass units) for quenching of star formation in
         halos, due to SNe or photo-heating feedback, or inefficient gas accretion.
@@ -1276,6 +1282,11 @@ class AstroParams(InputStruct):
         default=-0.5,
         converter=float,
     )
+    BETA_ESC: float = field(
+        default=0.0,
+        converter=float,
+    )
+    BETA_ESC_MINI: float = field(converter=float)
     F_ESC7_MINI: float = field(
         default=-2.0,
         converter=float,
@@ -1372,6 +1383,10 @@ def _F_STAR7_MINI_default(self):
     def _ALPHA_STAR_MINI_default(self):
         return self.ALPHA_STAR
 
+    @BETA_ESC_MINI.default
+    def _BETA_ESC_MINI_default(self):
+        return self.BETA_ESC
+
     @L_X_MINI.default
     def _L_X_MINI_default(self):
         return self.L_X
@@ -1430,6 +1445,8 @@ def get_logspaced_redshifts(
 ) -> tuple[float]:
     """Compute a sequence of redshifts to evolve over that are log-spaced."""
     redshifts = [min_redshift]
+    if z_step_factor <= 1.0:
+        return (min_redshift,)
     while redshifts[-1] < max_redshift:
         redshifts.append((redshifts[-1] + 1.0) * z_step_factor - 1.0)
 
@@ -1651,6 +1668,18 @@ def _astro_params_validator(self, att, val):
                 "update of M_TURN",
                 stacklevel=2,
             )
+            if (
+                self.astro_options.USE_MASS_DEPENDENT_ZETA
+                and val.BETA_ESC != 0
+                and self.astro_options.PHOTON_CONS_TYPE
+                not in ["no-photoncons", "alpha-photoncons"]
+            ):
+                warnings.warn(
+                    f"You have set BETA_ESC != 0 but PHOTON_CONS_TYPE is {self.astro_options.PHOTON_CONS_TYPE}. "
+                    "This changes the escape fraction so it is not consistent with the manual setting of scaling."
+                    "Set PHOTON_CONS_TYPE to 'no-photoncons' or 'alpha-photoncons' if you want the scaling to be exact.",
+                    stacklevel=2,
+                )
 
         if (
             self.astro_options.HII_FILTER == "sharp-k"

tests/test_input_structs.py

@@ -368,6 +368,22 @@ def test_dim_setting_serialization(
         new = new.evolve_input_structs(**evolved)
         self.check_attributes_dim(new.simulation_options, expected)
 
+    with pytest.warns(
+        UserWarning,
+        match="You have set BETA_ESC != 0 but PHOTON_CONS_TYPE is f-photoncons.",
+    ):
+        InputParameters(
+            cosmo_params=CosmoParams(),
+            astro_params=AstroParams(BETA_ESC=0.1),
+            simulation_options=SimulationOptions(),
+            matter_options=MatterOptions(),
+            astro_options=AstroOptions(
+                PHOTON_CONS_TYPE="f-photoncons",
+                USE_MASS_DEPENDENT_ZETA=True,
+            ),
+            random_seed=1,
+        )
+
 
 class TestMatterOptions:
     """Tests of the MatterOptions class."""