Binary L2 mass-loss fraction for thin accretion disk

binary/make/makefile_base

@@ -44,6 +44,7 @@ SRCS = \
    binary_utils.f90 \
    binary_photos.f90 \
    binary_jdot.f90 \
+   binary_disk.f90 \
    binary_wind.f90 \
    binary_mdot.f90 \
    binary_tides.f90 \

binary/private/binary_disk.f90

@@ -0,0 +1,315 @@
+! ***********************************************************************
+!
+!   Copyright (C) 2025  Philip Mocz, Mathieu Renzo & The MESA Team
+!
+!   This program is free software: you can redistribute it and/or modify
+!   it under the terms of the GNU Lesser General Public License
+!   as published by the Free Software Foundation,
+!   either version 3 of the License, or (at your option) any later version.
+!
+!   This program is distributed in the hope that it will be useful,
+!   but WITHOUT ANY WARRANTY; without even the implied warranty of
+!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+!   See the GNU Lesser General Public License for more details.
+!
+!   You should have received a copy of the GNU Lesser General Public License
+!   along with this program. If not, see <https://www.gnu.org/licenses/>.
+!
+! ***********************************************************************
+
+module binary_disk
+
+   use const_def, only: dp, pi, pi2, one_third, two_thirds, standard_cgrav, Msun, Rsun, secyer, crad, boltzm, clight, mp
+   use star_lib
+   use star_def
+   use math_lib
+   use binary_def
+
+   implicit none
+
+   ! more numerical constants
+   real(dp), parameter :: one_fourth = 1.0_dp / 4.0_dp
+   real(dp), parameter :: seven_fourths = 7.0_dp / 4.0_dp
+   real(dp), parameter :: one_eigth = 1.0_dp / 8.0_dp
+   real(dp), parameter :: three_eigths = 3.0_dp / 8.0_dp
+   real(dp), parameter :: one_ninth = 1.0_dp / 9.0_dp
+   real(dp), parameter :: four_twentyseventh = 4.0_dp / 27.0_dp
+
+   ! This module includes functions for mass transfer and L2 mass loss for binary systems with a thin disk
+   ! TODO: switch Kramers opacity module to real one in eval_L2_mass_loss_fraction()
+
+contains
+
+   subroutine eval_L2_mass_loss_fraction(donor_mass, accretor_mass, mass_transfer_rate, orbital_separation, &
+                                         disk_alpha, disk_mu, &
+                                         fL2, ierr)
+      ! Calculate the (outer) L2 mass-loss fraction
+      ! according to Lu et al. (2023, MNRAS 519, 1409) "On rapid binary mass transfer -I. Physical model"
+      ! https://ui.adsabs.harvard.edu/abs/2023MNRAS.519.1409L/abstract
+      real(dp), intent(in) :: donor_mass         ! [M_sun]
+      real(dp), intent(in) :: accretor_mass      ! [M_sun]
+      real(dp), intent(in) :: mass_transfer_rate ! [M_sun/yr]
+      real(dp), intent(in) :: orbital_separation ! [R_sun]
+      real(dp), intent(in) :: disk_alpha         ! disk alpha viscosity parameter (dimensionless)
+      real(dp), intent(in) :: disk_mu            ! disk mean molecular weight (dimensionless)
+      real(dp), intent(out) :: fL2               ! L2 mass-loss fraction (dimensionless)
+      integer, intent(out) :: ierr
+
+      real(dp), parameter :: eps_small = 1d-12  ! a very small number
+      real(dp), parameter :: tol = 1d-8  ! fractional tolerance for bisection method
+      real(dp) :: M2, M1dot, a, q, log_q
+      real(dp) :: xL1, xL2, mu
+      real(dp) :: Rd_over_a, PhiL1_dimless, PhiL2_dimless, PhiRd_dimless
+      real(dp) :: GM2, Rd, Phi_units, PhiL1, PhiL2, PhiRd, omega_K
+      real(dp) :: c1, c2, c3, c4
+
+      integer :: i
+      integer, parameter :: n_the = 50  ! number of grid points for disk thickness search
+      real(dp), parameter :: the_grid_min = 0.1_dp
+      real(dp), parameter :: the_grid_max = 1.0_dp
+      real(dp) :: the_grid(n_the)  ! grid for disk thickness
+      real(dp) :: T_arr(n_the)
+      real(dp), parameter :: T_floor = 3.0d3  ! [K] the minimum value for disk temperature solution
+      real(dp) :: T, T_max, the, the_left, the_right, the_min, the_max, separation_factor, dlogthe
+      real(dp) :: T_left, T_right, f1, f1_left, f2, f2_left, f2_right, f, f_left
+      real(dp) :: logT_arr(n_the), logthe_grid(n_the)
+
+      ierr = 0
+
+      ! key parameters
+      M2 = accretor_mass * Msun
+      M1dot = mass_transfer_rate * Msun/secyer
+      a = orbital_separation * Rsun
+      q = accretor_mass / donor_mass  ! mass ratio M2/M1
+
+      log_q = log10(q)
+
+      ! positions of Lagrangian points (based on analytic fits)
+      xL1 = -0.0355_dp * log_q**2 + 0.251_dp * abs(log_q) + 0.500_dp  ! [a = SMA]
+      xL2 = 0.0756_dp * log_q**2 - 0.424_dp * abs(log_q) + 1.699_dp  ! [a]
+
+      if (log_q > 0.0_dp) then  ! m2 is more massive
+         xL1 = 1.0_dp - xL1
+         xL2 = 1.0_dp - xL2
+      end if
+      mu = q / (1.0_dp + q)
+
+      ! outer disk radius
+      Rd_over_a = pow4(1.0_dp - xL1) / mu
+      ! relavent potential energies
+      PhiL1_dimless = -((1.0_dp - mu)/abs(xL1) + mu/abs(1.0_dp - xL1) + 0.5_dp*(xL1 - mu)**2)  ! [G(M1+M2)/a]
+      PhiL2_dimless = -((1.0_dp - mu)/abs(xL2) + mu/abs(1.0_dp - xL2) + 0.5_dp*(xL2 - mu)**2)  ! [G(M1+M2)/a]
+      PhiRd_dimless = -(1.0_dp - mu + mu/Rd_over_a + 0.5_dp*(1.0_dp - mu)**2)
+
+      GM2 = standard_cgrav * M2
+
+      Rd = Rd_over_a*a
+      Phi_units = standard_cgrav * (M2 / mu) / a
+      PhiL1 = PhiL1_dimless * Phi_units
+      PhiL2 = PhiL2_dimless * Phi_units
+      PhiRd = PhiRd_dimless * Phi_units
+      ! Keplerian frequency at Rd
+      omega_K = sqrt(GM2 / pow3(Rd))
+
+
+      ! constants involved in numerical solutions
+      c1 = two_thirds * pi * crad * disk_alpha * Rd / (omega_K * M1dot)
+      c2 = boltzm * Rd / (GM2 * disk_mu * mp)
+      c3 = 8.0_dp * pi2 * crad * disk_alpha * clight * Rd**2 / (M1dot**2 * omega_K)
+      c4 = 2.0_dp * pi * disk_mu * crad * disk_alpha * omega_K * mp * pow3(Rd) / (boltzm * M1dot)
+
+      ! Create logarithmically spaced grid for grid search for disk thickness [only used at the beginning]
+      do i = 1, n_the
+         the_grid(i) = exp10(log10(the_grid_min) + (i - 1) * (log10(the_grid_max) - log10(the_grid_min)) / (n_the - 1))
+      end do
+
+      ! only T < T_max is possible to calculate
+      T_max = pow(four_twentyseventh / (c1**2 * c2), one_ninth)
+
+      do i = 1, n_the
+         T_arr(i) = 0.0_dp
+      end do
+
+      do i = 1, n_the
+         the = the_grid(i)
+         ! use bisection method
+         T_left = 0.1_dp * min(the**2 / c2, T_max)
+         f1_left = f1_the_T_fL2(the, T_left, 0.0_dp, c1, c2)
+         T_right = T_max
+         do while (abs((T_left - T_right) / T_right) > tol)
+            T = 0.5_dp * (T_left + T_right)
+            f1 = f1_the_T_fL2(the, T, 0.0_dp, c1, c2)
+            if (f1 * f1_left > 0) then
+               T_left = T
+               f1_left = f1
+            else
+               T_right = T
+            end if
+         end do
+         T_arr(i) = 0.5_dp * (T_left + T_right)
+      end do
+      ! now we have obtained numerical relation between the and T
+      do i = 1, n_the
+         logT_arr(i) = log10(T_arr(i))
+         logthe_grid(i) = log10(the_grid(i))
+      end do
+      dlogthe = logthe_grid(2) - logthe_grid(1)
+
+      ! bisection to find the numerical solution to f2(the, T, fL2=0)=0
+      the_right = 1.0_dp
+      f2_right = f2_the_T_fL2(the_right, T_the_nofL2(the_right, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
+                              c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+      separation_factor = 0.95_dp
+      the_left = separation_factor * the_right
+      f2_left = f2_the_T_fL2(the_left, T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
+                             c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+      do while (f2_left * f2_right > 0)
+         ! need to decrease the_left
+         the_right = the_left
+         f2_right = f2_left
+         the_left = the_left * separation_factor
+         f2_left = f2_the_T_fL2(the_left, T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
+                                c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+      end do
+      ! now the solution is between the_left and the_right
+      do while (abs((the_left - the_right) / the_right) > tol)
+         the = 0.5_dp * (the_left + the_right)
+         f2 = f2_the_T_fL2(the, T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
+                           c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+         if (f2 * f2_left > 0.0_dp) then
+            the_left = the
+            f2_left = f2
+         else
+            the_right = the
+         end if
+       end do
+      ! solution
+      the = 0.5_dp * (the_left + the_right)
+      T = T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
+      the_max = sqrt(three_eigths * c2 * T + one_fourth * (PhiL2 - PhiRd) / (GM2 / Rd) - one_eigth)
+
+      if (the < the_max) then
+         ! return a tiny numner
+         fL2 = eps_small
+      else if (mass_transfer_rate < eps_small) then
+         ! no mass transfer, so no L2 mass loss
+         fL2 = 0.0_dp
+      else
+         the_min = ( 0.5_dp * sqrt((PhiL2 - PhiRd) / (GM2 / Rd) - 0.5_dp) )  ! corresponding to fL2=1, T=0
+         ! need to find the maximum corresponding to fL2=0
+         ! this is given by the intersection between T_the(the), T_the_nofL2(the)
+         the_left = the_min
+         the_right = 1.0_dp
+         f_left = T_the(the_left, c2, PhiL2, PhiRd, GM2, Rd) - T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe)
+         do while (abs((the_left - the_right) / the_right) > tol)
+            the = 0.5_dp * (the_left + the_right)
+            f = T_the(the, c2, PhiL2, PhiRd, GM2, Rd) - T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
+            if (f * f_left > 0.0_dp) then
+               the_left = the
+               f_left = f
+            else
+               the_right = the
+            end if
+         end do
+         the_max = 0.5_dp * (the_left + the_right)  ! this corresponds to fL2=0
+
+         ! --- numerical solution for f2(the, T, fL2)=0 under non-zero fL2
+
+         ! -- do not use exactly the_min (corresponding to T = 0, bc. kap table breaks down)
+         ! -- define another the_min based on T_floor (kap table won't be a problem)
+         the_min = sqrt( three_eigths * c2 * T_floor &
+                       + one_fourth * (PhiL2 - PhiRd) / (GM2 / Rd) - one_eigth )
+         the_left = the_min
+         f2_left = f2_the_T_fL2(the_left, &
+                                T_the(the_left, c2, PhiL2, PhiRd, GM2, Rd), &
+                                fL2_the(the_left, c1, c2, PhiL2, PhiRd, GM2, Rd), &
+                                c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+         the_right = the_max / (1.0_dp + eps_small)
+         ! bisection again
+         do while (abs((the_left - the_right) / the_right) > tol)
+            the = 0.5_dp * (the_left + the_right)
+            f2 = f2_the_T_fL2(the, T_the(the, c2, PhiL2, PhiRd, GM2, Rd), fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd), &
+                              c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+            if (f2 * f2_left > 0.0_dp) then
+               the_left = the
+               f2_left = f2
+            else
+               the_right = the
+            end if
+         end do
+         ! solution
+         the = 0.5_dp * (the_left + the_right)
+         fL2 = fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd)
+      end if
+
+   end subroutine eval_L2_mass_loss_fraction
+
+
+   ! Helper Functions
+   real(dp) function f1_the_T_fL2(the, T, fL2, c1, c2)
+      real(dp), intent(in) :: the, T, fL2, c1, c2
+      f1_the_T_fL2 = c1 * pow4(T) * pow3(the) / (1.0_dp - fL2) - the**2 + c2 * T
+   end function f1_the_T_fL2
+
+   real(dp) function  kap(rho, T)
+      ! simplified Kramers rule (cgs; approximate)
+      real(dp), intent(in) :: rho, T
+      kap = 0.34_dp + 3.0d24 * rho * pow(T, -3.5_dp)
+   end function kap
+
+   real(dp) function f2_the_T_fL2(the, T, fL2, &
+                                  c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
+      real(dp), intent(in) :: the, T, fL2
+      real(dp), intent(in) :: c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2
+      real(dp) :: x, U_over_P, rho
+      x = c4 * pow3(T * the) / (1.0_dp - fL2)
+      U_over_P = (1.5_dp + x) / (1.0_dp + one_third * x)
+      rho = (1.0_dp - fL2) * M1dot / (2.0_dp * pi * disk_alpha * omega_K * pow3(Rd)) / pow2(the)
+      f2_the_T_fL2 = &
+         seven_fourths &
+         - (1.5_dp * U_over_P + c3 * pow4(T) / kap(rho, T) / (1.0_dp - fL2) ** 2) * the**2 &
+         - PhiRd / (GM2 / Rd) &
+         + (PhiL1 - fL2 * PhiL2) / (GM2 / Rd) / (1.0_dp - fL2)
+   end function f2_the_T_fL2
+
+   real(dp) function T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
+      real(dp), intent(in) :: the
+      integer, intent(in) :: n_the
+      real(dp), intent(in) :: logthe_grid(n_the), logT_arr(n_the)
+      real(dp), intent(in) :: dlogthe
+      real(dp) :: logthe, slope, logT
+      integer :: i_the
+      ! under the assumption fL2=0
+      logthe = log10(the)
+      if (logthe > logthe_grid(n_the-1)) then
+         ! use analytic extrapolation
+         T_the_nofL2 = exp10(logT_arr(n_the-1) - 0.25_dp * (logthe - logthe_grid(n_the-1)))
+      else if (logthe < logthe_grid(1)) then
+         ! analytic extrapolation
+         T_the_nofL2 = exp10(logT_arr(1) + 2.0_dp * (logthe - logthe_grid(1)))
+      else
+         i_the = floor((logthe - logthe_grid(1)) / dlogthe) + 1
+         slope = (logT_arr(i_the + 1) - logT_arr(i_the)) / dlogthe
+         logT = logT_arr(i_the) + (logthe - logthe_grid(i_the)) * slope
+         T_the_nofL2 = exp10(logT)
+      end if
+   end function T_the_nofL2
+
+   real(dp) function T_the(the, c2, PhiL2, PhiRd, GM2, Rd)
+      real(dp), intent(in) :: the
+      real(dp), intent(in) :: c2, PhiL2, PhiRd, GM2, Rd
+      ! only for non-zero fL2
+      T_the = (8.0_dp * the**2 + 1.0_dp - 2.0_dp * (PhiL2 - PhiRd) / (GM2 / Rd)) / (3.0_dp * c2)
+   end function T_the
+
+   real(dp) function fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd)
+      real(dp), intent(in) :: the
+      real(dp), intent(in) :: c1, c2, PhiL2, PhiRd, GM2, Rd
+      real(dp) :: T
+      ! only for non-zero fL2
+      T = T_the(the, c2, PhiL2, PhiRd, GM2, Rd)
+      fL2_the =  1.0_dp - c1 * pow4(T) * pow3(the) / (the**2 - c2 * T)
+   end function fL2_the
+
+
+end module binary_disk

binary/public/binary_lib.f90

@@ -252,5 +252,21 @@ real(dp) function binary_compute_k_div_T(binary_id, is_donor, has_convective_env
 
       end function binary_compute_k_div_T
 
+      real(dp) function binary_L2_mass_loss_fraction(donor_mass, accretor_mass, mass_transfer_rate, orbital_separation, &
+                                                     disk_alpha, disk_mu, ierr)
+         use binary_disk, only: eval_L2_mass_loss_fraction
+         real(dp), intent(in) :: donor_mass         ! [M_sun]
+         real(dp), intent(in) :: accretor_mass      ! [M_sun]
+         real(dp), intent(in) :: mass_transfer_rate ! [M_sun/yr]
+         real(dp), intent(in) :: orbital_separation ! [R_sun]
+         real(dp), intent(in) :: disk_alpha         ! disk alpha viscosity parameter (dimensionless)
+         real(dp), intent(in) :: disk_mu            ! disk mean molecular weight (dimensionless)
+         integer, intent(out) :: ierr
+
+         call eval_L2_mass_loss_fraction(donor_mass, accretor_mass, mass_transfer_rate, orbital_separation, &
+                                    disk_alpha, disk_mu, &
+                                    binary_L2_mass_loss_fraction, ierr)
+      end function binary_L2_mass_loss_fraction
+
       end module binary_lib
 

binary/test/src/test_binary.f90

@@ -14,6 +14,7 @@ module test_binary_mod
 
    subroutine do_test
       integer :: binary_id, id, ierr, res
+      real(dp) :: fL2
       type(binary_info), pointer :: b
       type(star_info), pointer :: s
 
@@ -78,6 +79,21 @@ subroutine do_test
 
       write (*, 1) 'b% max_timestep / secyer', b%max_timestep/secyer  ! should be 2.5
 
+      ! Test accretion disk functions
+      write (*, *) 'check L2_mass_loss_fraction behavior'
+
+      fL2 = binary_L2_mass_loss_fraction(20.0_dp, 10.0_dp, 0.001_dp, 2.0_dp, 0.1_dp, 1.3_dp/2.4_dp, ierr)
+
+      write (*, 1) 'fL2', fL2  ! should be ~0.9
+
+      fL2 = binary_L2_mass_loss_fraction(20.0_dp, 10.0_dp, 0.0001_dp, 2.0_dp, 0.1_dp, 1.3_dp/2.4_dp, ierr)
+
+      write (*, 1) 'fL2', fL2  ! should be ~0.2
+
+      fL2 = binary_L2_mass_loss_fraction(20.0_dp, 18.0_dp, 0.0_dp, 304.0_dp, 0.1_dp, 0.64_dp, ierr)
+
+      write (*, 1) 'fL2', fL2  ! should be 0.0
+
       write (*, '(a)') 'done'
 
    end subroutine do_test

binary/test/test_output

@@ -1,4 +1,8 @@
  check time_delta_coeff behavior
                                b% max_timestep / secyer    5.0000000000000044D+00
                                b% max_timestep / secyer    2.5000000000000067D+00
+ check L2_mass_loss_fraction behavior
+                                                    fL2    8.9701564145435098D-01
+                                                    fL2    2.1022134020418981D-01
+                                                    fL2    0.0000000000000000D+00
 done

docs/source/changelog.rst

@@ -65,6 +65,9 @@ The Fe core-collapse infall condition ``fe_core_infall_limit`` has been adjusted
 Users can switch between either choice with the new logical control ``report_max_infall_inside_fe_core``.
 See the ``&controls`` for further details.
 
+L2 mass-loss fraction according to Lu et al. (2023) is available as a public function in the bindary moduele:
+``binary_L2_mass_loss_fraction(donor_mass, accretor_mass, mass_transfer_rate, orbital_separation, disk_alpha, disk_mu, ierr)``.  
+
 .. _Bug Fixes main:
 
 Bug Fixes