COSMIC-PopSynth · TristanBruel · Jun 17, 2026 · Jun 17, 2026 · Jun 17, 2026 · Jun 17, 2026
diff --git a/examples/Params.ini b/examples/Params.ini
@@ -419,6 +419,15 @@ alpha1 = [1.0, 1.0]
 ; Default: 0.0
 lambdaf = 0.0
 
+; ce2stageflag
+; Sets whether to use the standard energy formalism or the two-stage CE formalism from Hirai & Mandel 2022
+; 
+; Options: 
+;    0 - Use the standard energy formalism
+;    1 - Use the two-stage CE formalism from Hirai & Mandel 2022
+; Default: 0
+ce2stageflag = 0
+
 ; ceflag
 ; Selects the model to set the initial orbital energy (choose whether to use total mass of the stars instead of the core masses)
 ; 

diff --git a/src/cosmic/consts.py b/src/cosmic/consts.py
@@ -42,7 +42,7 @@
         "epsnov", "eddfac", "gamma", "LBV_flag",
     ],
     "cevars": ["alpha1", "lambdaf", "qcrit_array"],
-    "ceflags": ["ceflag", "cekickflag", "cemergeflag", "cehestarflag", "ussn"],
+    "ceflags": ["ce2stageflag", "ceflag", "cekickflag", "cemergeflag", "cehestarflag", "ussn"],
     "flags": [
         "tflag", "ifflag", "wdflag", "rtmsflag", "bhflag", "remnantflag",
         "maltsev_mode", "grflag", "bhms_coll_flag", "wd_mass_lim", "aic",
@@ -63,4 +63,4 @@
     "mixvars": ["rejuv_fac"],
     "metvars": ["zsun"],
     "rand1": ["randomseed"]
-}
+}
diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json
@@ -865,6 +865,23 @@
                     }
                 ]
             }, 
+	    {
+                "name": "ce2stageflag",
+                "description": "Sets whether to use the standard energy formalism or the two-stage CE formalism from Hirai & Mandel 2022",
+                "type": "dropdown",
+                "options-preface": "",
+                "options": [
+                    {
+                        "name": 0,
+                        "description": "Use the standard energy formalism",
+                        "default": true
+                    },
+                    {
+                        "name": 1,
+                        "description": "Use the two-stage CE formalism from Hirai & Mandel 2022"
+                    }
+                ]
+            },
             {
                 "name": "ceflag",
                 "description": "Selects the model to set the initial orbital energy (choose whether to use total mass of the stars instead of the core masses)",

diff --git a/src/cosmic/src/SSE/SSE_zfuncs.f b/src/cosmic/src/SSE/SSE_zfuncs.f
@@ -845,6 +845,94 @@ real*8 FUNCTION lmcgbf(mc,GB)
 *
       return
       end
+***
+      real*8 FUNCTION tonset(tmin,z)
+      implicit none
+      real*8 tmin,z
+*
+* A function to evaluate Tonset (Teff when Mconv reaches 1% of its max)
+*       from Mandel, Hirai & Picker 2024
+*
+      tonset = tmin /min(0.695d0-0.057d0*log10(z),0.95d0)
+*
+      return
+      end
+***
+      real*8 FUNCTION menvmax(kw,m,z)
+      implicit none
+      integer kw
+      real*8 m,z
+      real*8 mcagbf,mifmcf
+      external mcagbf
+*
+* A function to evaluate the maximum mass of the convective envelope
+*       from Picker, Hirai & Mandel 2024
+*
+* I'm not sure how HG stars should be treated ...
+*
+      if(kw.ge.10)then
+         menvmax = 0.0d0
+      elseif(kw.lt.2)then
+         if(m.gt.1.25d0)then 
+            menvmax = 0.0d0
+         else
+            menvmax = m
+            if(m.gt.0.35d0)then
+               menvmax = 0.35d0*(1.25d0-m)*(1.25d0-m)/0.81d0
+            endif
+         endif
+      else
+         mifmcf = -0.023*log10(z)-0.0023
+         menvmax = max(m-mcagbf(m)*(1.0d0+mifmcf), 0.0d0)
+      endif
+*     
+      return
+      end
+***
+      real*8 FUNCTION mconvenv(kw,m,z,teff,tmin,tset,aj,tm)
+      implicit none
+      integer kw
+      real*8 m,z,teff,tmin,tset,aj,tm
+      real*8 menvmax
+      external menvmax
+*
+* A function to evaluate the mass of the convective envelope
+*       from Picker, Hirai & Mandel 2024
+*
+      if(kw.ge.10)then
+         mconvenv = 0.0d0
+      elseif(kw.lt.2)then
+         if(m.gt.1.25d0)then
+            mconvenv = 0.0d0
+         else
+            mconvenv = menvmax(kw,m,z)*sqrt(sqrt(1.0d0-aj/tm))
+         endif
+      else
+         mconvenv = menvmax(kw,m,z)/(1.0d0+exp(4.6d0*
+     &       (tmin+tset-2.0d0*teff)/(tmin-tset)))
+      endif
+*
+      return
+      end
+***
+      real*8 FUNCTION celamhe(m,z,mconv,mconvmax)
+      implicit none
+      real*8 m,z,mconv,mconvmax
+      real*8 m2,b1
+*
+* A function to evaluate the binding energy factor of the convective envelope
+*       from Picker, Hirai & Mandel 2024
+*
+      m2 = 0.0023d0*log10(z)**2+0.0088d0*log10(z)+0.013d0
+      b1 = m2*m-0.23d0
+      if((mconv/mconvmax).gt.0.3d0)then
+         celamhe = exp(0.42*(mconv/mconvmax)+b1)
+      else
+         celamhe = exp(0.3*0.42+b1)
+      endif
+*
+      return
+      end
 ***
       real*8 FUNCTION lHeIf(m,mhefl)
       implicit none

diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f
@@ -75,6 +75,13 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
       REAL*8 rad(2),tms(2),lumin(2),B_0(2),bacc(2),tacc(2),epoch(2)
       REAL*8 menv_bpp(2),renv_bpp(2)
       REAL*8 ALPHA_CE
+
+      REAL*8 met,reagb,tmin,tset,mconvmax
+      REAL*8 ragbf,menvmax,mconvenv,CELAMHE,tonset
+      EXTERNAL ragbf,menvmax,mconvenv,CELAMHE,tonset
+      REAL*8 mcgbtf,lmcgbf
+      EXTERNAL mcgbtf,lmcgbf
+      REAL*8 m1endstage1,m2endstage1
 *
 * Initialize
 *
@@ -120,20 +127,55 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
         RZAMS = RZAMSF(M01)
         LAMB1 = CELAMF(KW,M01,L1,R1,RZAMS,MENVD,LAMBDAF)
       ENDIF
+*
+* Two-stage energy formalism
+*
+      IF(CE2STAGEFLAG.EQ.1)THEN
+         IF (using_SSE.eq.1) THEN
+* Estimate the mass of the convective envelope 
+* from the fits in Picker, Hirai, & Mandel 2024
+           met = 10**(LOG10(zpars(14))/0.4)
+           reagb = ragbf(M1,LUMS(7),zpars(2))
+           tmin = 1000.d0*((1130.d0*LUMS(7)/
+     &                       (reagb**2.d0))**(1.d0/4.d0))
+           tset = tonset(tmin,met)
+           teff1 = 1000.d0*((1130.d0*lumin(1)/
+     &                       (rad1_bpp**2.d0))**(1.d0/4.d0))
+           MENV = mconvenv(KW,M1,met,teff1,tmin,tset,AJ1,TM1)
+           mconvmax = menvmax(KW,M1,met)
+           LAMB1 = CELAMHE(M1,met,MENV,mconvmax)
+         ENDIF
+* if > 8 Msun, the envelope participating in the first stage is only the convective one
+* if < 2 Msun, the entire envelope participates in the first stage
+* linear interpolation in between
+         IF(M1.GE.8.0d0)THEN
+            m1endstage1 = M1-MENV
+         ELSEIF(M1.LT.2.0d0)THEN
+            m1endstage1 = MC1
+         ELSE
+            m1endstage1 = MC1 + (M1-MC1 - MENV) * (M1-2.0d0) /6.0d0
+         ENDIF
+*
+* Calculate the binding energy of the convective envelope (multiplied by lambda).
+*
+         EBINDI = M1*(M1-m1endstage1)/(LAMB1*R1)
+      ELSE
+*
+* Calculate the binding energy of the giant envelope (multiplied by lambda).
+*
+      EBINDI = M1*(M1-MC1)/(LAMB1*R1)
+      ENDIF
       KW = KW2
       CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS,dtm,star2)
       CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS,
      &            R2,L2,KW2,MC2,RC2,MENV,RENV,K22,
      &            bhspin2,star2)
       OSPIN2 = JSPIN2/(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2)
 *
-* Calculate the binding energy of the giant envelope (multiplied by lambda).
-*
-      EBINDI = M1*(M1-MC1)/(LAMB1*R1)
-*
 * If the secondary star is also giant-like add its envelopes energy.
 * Determine EORBI based on CEFLAG (CEFLAG=1 for de Kool prescription)
 *
+      m2endstage1 = M2
       IF(KW2.GE.2.AND.KW2.LE.9.AND.KW2.NE.7)THEN
          MENVD = MENV/(M2-MC2)
          IF (using_METISSE.eq.1) THEN
@@ -142,7 +184,44 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
             RZAMS = RZAMSF(M02)
             LAMB2 = CELAMF(KW,M02,L2,R2,RZAMS,MENVD,LAMBDAF)
          ENDIF
-         EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2)
+*
+* Two-stage energy formalism
+*
+         IF(CE2STAGEFLAG.EQ.1)THEN
+            IF (using_SSE.eq.1) THEN
+* Estimate the mass of the convective envelope 
+* from the fits in Picker, Hirai, & Mandel 2024
+              met = 10**(LOG10(zpars(14))/0.4)
+              reagb = ragbf(M2,LUMS(7),zpars(2))
+              tmin = 1000.d0*((1130.d0*LUMS(7)/
+     &                          (reagb**2.d0))**(1.d0/4.d0))
+              tset = tonset(tmin,met)
+              teff2 = 1000.d0*((1130.d0*lumin(2)/
+     &                          (rad2_bpp**2.d0))**(1.d0/4.d0))
+              MENV = mconvenv(KW,M2,met,teff2,tmin,tset,AJ2,TM2)
+              mconvmax = menvmax(KW,M2,met)
+              LAMB2 = CELAMHE(M2,met,MENV,mconvmax)
+            ENDIF
+* if > 8 Msun, the envelope participating in the first stage is only the convective one
+* if < 2 Msun, the entire envelope participates in the first stage
+* linear interpolation in between
+            IF(M2.GE.8.0d0)THEN
+               m2endstage1 = M2-MENV
+            ELSEIF(M1.LT.2.0d0)THEN
+               m2endstage1 = MC2
+            ELSE
+               m2endstage1 = MC2 + (M2-MC2 - MENV) * (M2-2.0d0) /6.0d0
+            ENDIF
+*
+* Calculate the binding energy of the convective envelope (multiplied by lambda).
+*
+            EBINDI = EBINDI + M2*(M2-m2endstage1)/(LAMB2*R2)
+         ELSE
+*
+* Calculate the binding energy of the giant envelope (multiplied by lambda).
+*
+            EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2)
+         ENDIF
 *
 * Calculate the initial orbital energy
 *
@@ -167,7 +246,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
       ENDIF
       EORBF = EORBI + EBINDI/ALPHA_CE
 *
-* If the secondary is on the main sequence see if it fills its Roche lobe.
+* If the secondary is on the main sequence.
 *
       IF(KW2.LE.1.OR.KW2.EQ.7)THEN
          SEPF = MC1*M2/(2.D0*EORBF)
@@ -176,6 +255,63 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
          RL1 = RL(Q1)
          RL2 = RL(Q2)
 *
+* If degenerate or giant secondary.
+*
+      ELSE
+         SEPF = MC1*MC2/(2.D0*EORBF)
+         Q1 = MC1/MC2
+         Q2 = 1.D0/Q1
+         RL1 = RL(Q1)
+         RL2 = RL(Q2)
+      ENDIF
+*
+* Two-stage energy formalism
+*
+      IF(CE2STAGEFLAG.EQ.1)THEN
+         SEPF = m1endstage1*m2endstage1/(2.D0*EORBF)
+* log end of the first stage with evolve_type=3
+         q1_bpp = m1endstage1/m2endstage1
+         q2_bpp = 1.d0/q1_bpp
+         rrl1_bpp = RC1/(RL(q1_bpp)*SEP_postCE)
+         rrl2_bpp = R2/(RL(q2_bpp)*SEP_postCE)
+         TB = (SEPF/AURSUN)*
+     &        SQRT(SEPF/(AURSUN*(m1endstage1+m2endstage1)))
+         CALL writetab(jp,tphys,3.d0,
+     &        m1endstage1,m2endstage1,KW1,KW2,
+     &        SEPF,TB,ECC,
+     &        rrl1_bpp,rrl2_bpp,
+     &        AJ1,AJ2,tms1_bpp,tms2_bpp,
+     &        mc_he(1),mc_he(2),mc_co(1),mc_co(2),
+     &        rad1_bpp,rad2_bpp,
+     &        M02,M01,lumin(1),lumin(2),
+     &        teff1,teff2,
+     &        RC2,RC1,menv_bpp(1),menv_bpp(2),
+     &        renv_bpp(1),
+     &        renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2),
+     &        bacc(1),bacc(2),tacc(1),tacc(2),epoch(1),
+     &        epoch(2),bhspin2,bhspin1,
+     &        deltam_2,deltam_1,formation2,formation1,
+     &        binstate,mergertype,zpars(14)**2.d5,'bpp')
+*
+* Second stage: stable mass transfer of the radiative intershell
+* Binary hardening formula from Picker, Hirai & Mandel 2024
+*
+         IF(m1endstage1.GT.MC1)THEN
+            SEPF = SEPF*((m1endstage1+m2endstage1)/(MC1+m2endstage1))
+     &          *(m1endstage1/MC1)**2
+     &          *EXP(-2*(m1endstage1-MC1)/m2endstage1)
+            EORBF = MC1*m2endstage1/(2.D0*SEPF)
+         ENDIF
+         Q1 = MC1/m2endstage1
+         Q2 = 1.D0/Q1
+         RL1 = RL(Q1)
+         RL2 = RL(Q2)
+      ENDIF
+*
+* If the secondary is on the main sequence see if it fills its Roche lobe.
+*
+      IF(KW2.LE.1.OR.KW2.EQ.7)THEN
+*
 * If cemergeflag is set, cause kstars without clear core-envelope
 * structure to merge automatically if they enter a CE
 *

diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h
@@ -18,7 +18,9 @@
 
       REAL*8 don_lim,acc_lim(2),Mbh_initial,smt_periastron_check
       COMMON /MTVARS/ don_lim,acc_lim,Mbh_initial,smt_periastron_check
+      INTEGER ce2stageflag
       INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn
+      COMMON /CEFLAGS/ ce2stageflag
       COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn
       INTEGER pisn_track(2)
       COMMON /TRACKERS/ pisn_track

diff --git a/src/cosmic/src/test_bse.f b/src/cosmic/src/test_bse.f
@@ -123,6 +123,7 @@ PROGRAM bse
 
       WRITE(*,*)mass0(1),mass0(2),tphysf,tb,kstar(1),kstar(2),z,ecc
       WRITE(*,*)neta,bwind,hewind,alpha1(1),alpha1(2),lambdaf,windflag
+      WRITE(*,*)ce2stageflag
       WRITE(*,*)rtmsflag,ceflag,tflag,ifflag,wdflag,bhflag,remnantflag
       WRITE(*,*)mxns,idum,pts1,pts2,pts3
       WRITE(*,*)sigma,beta,xi,acc2,epsnov,eddfac,gamma

diff --git a/src/cosmic/tests/data/Params.ini b/src/cosmic/tests/data/Params.ini
@@ -194,6 +194,11 @@ alpha1 = [1.0, 1.0]
 ; default=0.5
 lambdaf = 0.5
 
+; ce2stageflag=0 uses the standard energy formalism for CEE
+; ce2stageflag=1 uses the two-stage formalism CEE from Hirai & Mandel 2022
+; default=0 
+ce2stageflag=1
+
 ; ceflag=1 used the method from de Kool 1990 for setting the initial orbital energy
 ; ceflag=0 does not use this method (uses the core mass to calculate initial orbital energy)
 ; default=0

diff --git a/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 b/src/cosmic/tests/data/initial_conditions_for_testing.hdf5
diff --git a/src/cosmic/tests/data/kick_initial_conditions.h5 b/src/cosmic/tests/data/kick_initial_conditions.h5