Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2025 Philip Mocz, Mathieu Renzo & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module binary_disk
21 :
22 : use const_def, only: dp, pi, pi2, one_third, two_thirds, standard_cgrav, Msun, Rsun, secyer, crad, boltzm, clight, mp
23 : use star_lib
24 : use star_def
25 : use math_lib
26 : use binary_def
27 :
28 : implicit none
29 :
30 : ! more numerical constants
31 : real(dp), parameter :: one_fourth = 1.0_dp / 4.0_dp
32 : real(dp), parameter :: seven_fourths = 7.0_dp / 4.0_dp
33 : real(dp), parameter :: one_eigth = 1.0_dp / 8.0_dp
34 : real(dp), parameter :: three_eigths = 3.0_dp / 8.0_dp
35 : real(dp), parameter :: one_ninth = 1.0_dp / 9.0_dp
36 : real(dp), parameter :: four_twentyseventh = 4.0_dp / 27.0_dp
37 :
38 : ! This module includes functions for mass transfer and L2 mass loss for binary systems with a thin disk
39 : ! TODO: switch Kramers opacity module to real one in eval_L2_mass_loss_fraction()
40 :
41 : contains
42 :
43 3 : subroutine eval_L2_mass_loss_fraction(donor_mass, accretor_mass, mass_transfer_rate, orbital_separation, &
44 : disk_alpha, disk_mu, &
45 : fL2, ierr)
46 : ! Calculate the (outer) L2 mass-loss fraction
47 : ! according to Lu et al. (2023, MNRAS 519, 1409) "On rapid binary mass transfer -I. Physical model"
48 : ! https://ui.adsabs.harvard.edu/abs/2023MNRAS.519.1409L/abstract
49 : real(dp), intent(in) :: donor_mass ! [M_sun]
50 : real(dp), intent(in) :: accretor_mass ! [M_sun]
51 : real(dp), intent(in) :: mass_transfer_rate ! [M_sun/yr]
52 : real(dp), intent(in) :: orbital_separation ! [R_sun]
53 : real(dp), intent(in) :: disk_alpha ! disk alpha viscosity parameter (dimensionless)
54 : real(dp), intent(in) :: disk_mu ! disk mean molecular weight (dimensionless)
55 : real(dp), intent(out) :: fL2 ! L2 mass-loss fraction (dimensionless)
56 : integer, intent(out) :: ierr
57 :
58 : real(dp), parameter :: eps_small = 1d-12 ! a very small number
59 : real(dp), parameter :: tol = 1d-8 ! fractional tolerance for bisection method
60 6 : real(dp) :: M2, M1dot, a, q, log_q
61 6 : real(dp) :: xL1, xL2, mu
62 3 : real(dp) :: Rd_over_a, PhiL1_dimless, PhiL2_dimless, PhiRd_dimless
63 3 : real(dp) :: GM2, Rd, Phi_units, PhiL1, PhiL2, PhiRd, omega_K
64 : real(dp) :: c1, c2, c3, c4
65 :
66 : integer :: i
67 : integer, parameter :: n_the = 50 ! number of grid points for disk thickness search
68 : real(dp), parameter :: the_grid_min = 0.1_dp
69 : real(dp), parameter :: the_grid_max = 1.0_dp
70 153 : real(dp) :: the_grid(n_the) ! grid for disk thickness
71 153 : real(dp) :: T_arr(n_the)
72 : real(dp), parameter :: T_floor = 3.0d3 ! [K] the minimum value for disk temperature solution
73 6 : real(dp) :: T, T_max, the, the_left, the_right, the_min, the_max, separation_factor, dlogthe
74 3 : real(dp) :: T_left, T_right, f1, f1_left, f2, f2_left, f2_right, f, f_left
75 303 : real(dp) :: logT_arr(n_the), logthe_grid(n_the)
76 :
77 3 : ierr = 0
78 :
79 : ! key parameters
80 3 : M2 = accretor_mass * Msun
81 3 : M1dot = mass_transfer_rate * Msun/secyer
82 3 : a = orbital_separation * Rsun
83 3 : q = accretor_mass / donor_mass ! mass ratio M2/M1
84 :
85 3 : log_q = log10(q)
86 :
87 : ! positions of Lagrangian points (based on analytic fits)
88 3 : xL1 = -0.0355_dp * log_q**2 + 0.251_dp * abs(log_q) + 0.500_dp ! [a = SMA]
89 3 : xL2 = 0.0756_dp * log_q**2 - 0.424_dp * abs(log_q) + 1.699_dp ! [a]
90 :
91 3 : if (log_q > 0.0_dp) then ! m2 is more massive
92 0 : xL1 = 1.0_dp - xL1
93 0 : xL2 = 1.0_dp - xL2
94 : end if
95 3 : mu = q / (1.0_dp + q)
96 :
97 : ! outer disk radius
98 3 : Rd_over_a = pow4(1.0_dp - xL1) / mu
99 : ! relavent potential energies
100 3 : PhiL1_dimless = -((1.0_dp - mu)/abs(xL1) + mu/abs(1.0_dp - xL1) + 0.5_dp*(xL1 - mu)**2) ! [G(M1+M2)/a]
101 3 : PhiL2_dimless = -((1.0_dp - mu)/abs(xL2) + mu/abs(1.0_dp - xL2) + 0.5_dp*(xL2 - mu)**2) ! [G(M1+M2)/a]
102 3 : PhiRd_dimless = -(1.0_dp - mu + mu/Rd_over_a + 0.5_dp*(1.0_dp - mu)**2)
103 :
104 3 : GM2 = standard_cgrav * M2
105 :
106 3 : Rd = Rd_over_a*a
107 3 : Phi_units = standard_cgrav * (M2 / mu) / a
108 3 : PhiL1 = PhiL1_dimless * Phi_units
109 3 : PhiL2 = PhiL2_dimless * Phi_units
110 3 : PhiRd = PhiRd_dimless * Phi_units
111 : ! Keplerian frequency at Rd
112 3 : omega_K = sqrt(GM2 / pow3(Rd))
113 :
114 :
115 : ! constants involved in numerical solutions
116 3 : c1 = two_thirds * pi * crad * disk_alpha * Rd / (omega_K * M1dot)
117 3 : c2 = boltzm * Rd / (GM2 * disk_mu * mp)
118 3 : c3 = 8.0_dp * pi2 * crad * disk_alpha * clight * Rd**2 / (M1dot**2 * omega_K)
119 3 : c4 = 2.0_dp * pi * disk_mu * crad * disk_alpha * omega_K * mp * pow3(Rd) / (boltzm * M1dot)
120 :
121 : ! Create logarithmically spaced grid for grid search for disk thickness [only used at the beginning]
122 153 : do i = 1, n_the
123 153 : the_grid(i) = exp10(log10(the_grid_min) + (i - 1) * (log10(the_grid_max) - log10(the_grid_min)) / (n_the - 1))
124 : end do
125 :
126 : ! only T < T_max is possible to calculate
127 3 : T_max = pow(four_twentyseventh / (c1**2 * c2), one_ninth)
128 :
129 153 : do i = 1, n_the
130 153 : T_arr(i) = 0.0_dp
131 : end do
132 :
133 153 : do i = 1, n_the
134 150 : the = the_grid(i)
135 : ! use bisection method
136 150 : T_left = 0.1_dp * min(the**2 / c2, T_max)
137 150 : f1_left = f1_the_T_fL2(the, T_left, 0.0_dp, c1, c2)
138 150 : T_right = T_max
139 2855 : do while (abs((T_left - T_right) / T_right) > tol)
140 2705 : T = 0.5_dp * (T_left + T_right)
141 2705 : f1 = f1_the_T_fL2(the, T, 0.0_dp, c1, c2)
142 2855 : if (f1 * f1_left > 0) then
143 1462 : T_left = T
144 1462 : f1_left = f1
145 : else
146 : T_right = T
147 : end if
148 : end do
149 153 : T_arr(i) = 0.5_dp * (T_left + T_right)
150 : end do
151 : ! now we have obtained numerical relation between the and T
152 153 : do i = 1, n_the
153 150 : logT_arr(i) = log10(T_arr(i))
154 153 : logthe_grid(i) = log10(the_grid(i))
155 : end do
156 3 : dlogthe = logthe_grid(2) - logthe_grid(1)
157 :
158 : ! bisection to find the numerical solution to f2(the, T, fL2=0)=0
159 3 : the_right = 1.0_dp
160 : f2_right = f2_the_T_fL2(the_right, T_the_nofL2(the_right, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
161 3 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
162 3 : separation_factor = 0.95_dp
163 3 : the_left = separation_factor * the_right
164 : f2_left = f2_the_T_fL2(the_left, T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
165 3 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
166 35 : do while (f2_left * f2_right > 0)
167 : ! need to decrease the_left
168 32 : the_right = the_left
169 32 : f2_right = f2_left
170 32 : the_left = the_left * separation_factor
171 : f2_left = f2_the_T_fL2(the_left, T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
172 35 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
173 : end do
174 : ! now the solution is between the_left and the_right
175 72 : do while (abs((the_left - the_right) / the_right) > tol)
176 69 : the = 0.5_dp * (the_left + the_right)
177 : f2 = f2_the_T_fL2(the, T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
178 69 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
179 72 : if (f2 * f2_left > 0.0_dp) then
180 27 : the_left = the
181 27 : f2_left = f2
182 : else
183 42 : the_right = the
184 : end if
185 : end do
186 : ! solution
187 3 : the = 0.5_dp * (the_left + the_right)
188 3 : T = T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
189 3 : the_max = sqrt(three_eigths * c2 * T + one_fourth * (PhiL2 - PhiRd) / (GM2 / Rd) - one_eigth)
190 :
191 3 : if (the < the_max) then
192 : ! return a tiny numner
193 0 : fL2 = eps_small
194 3 : else if (mass_transfer_rate < eps_small) then
195 : ! no mass transfer, so no L2 mass loss
196 1 : fL2 = 0.0_dp
197 : else
198 2 : the_min = ( 0.5_dp * sqrt((PhiL2 - PhiRd) / (GM2 / Rd) - 0.5_dp) ) ! corresponding to fL2=1, T=0
199 : ! need to find the maximum corresponding to fL2=0
200 : ! this is given by the intersection between T_the(the), T_the_nofL2(the)
201 2 : the_left = the_min
202 2 : the_right = 1.0_dp
203 2 : f_left = T_the(the_left, c2, PhiL2, PhiRd, GM2, Rd) - T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe)
204 60 : do while (abs((the_left - the_right) / the_right) > tol)
205 58 : the = 0.5_dp * (the_left + the_right)
206 58 : f = T_the(the, c2, PhiL2, PhiRd, GM2, Rd) - T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
207 2 : if (f * f_left > 0.0_dp) then
208 19 : the_left = the
209 19 : f_left = f
210 : else
211 39 : the_right = the
212 : end if
213 : end do
214 2 : the_max = 0.5_dp * (the_left + the_right) ! this corresponds to fL2=0
215 :
216 : ! --- numerical solution for f2(the, T, fL2)=0 under non-zero fL2
217 :
218 : ! -- do not use exactly the_min (corresponding to T = 0, bc. kap table breaks down)
219 : ! -- define another the_min based on T_floor (kap table won't be a problem)
220 : the_min = sqrt( three_eigths * c2 * T_floor &
221 2 : + one_fourth * (PhiL2 - PhiRd) / (GM2 / Rd) - one_eigth )
222 2 : the_left = the_min
223 : f2_left = f2_the_T_fL2(the_left, &
224 : T_the(the_left, c2, PhiL2, PhiRd, GM2, Rd), &
225 : fL2_the(the_left, c1, c2, PhiL2, PhiRd, GM2, Rd), &
226 2 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
227 2 : the_right = the_max / (1.0_dp + eps_small)
228 : ! bisection again
229 44 : do while (abs((the_left - the_right) / the_right) > tol)
230 42 : the = 0.5_dp * (the_left + the_right)
231 : f2 = f2_the_T_fL2(the, T_the(the, c2, PhiL2, PhiRd, GM2, Rd), fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd), &
232 42 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
233 2 : if (f2 * f2_left > 0.0_dp) then
234 22 : the_left = the
235 22 : f2_left = f2
236 : else
237 20 : the_right = the
238 : end if
239 : end do
240 : ! solution
241 2 : the = 0.5_dp * (the_left + the_right)
242 2 : fL2 = fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd)
243 : end if
244 :
245 3 : end subroutine eval_L2_mass_loss_fraction
246 :
247 :
248 : ! Helper Functions
249 2855 : real(dp) function f1_the_T_fL2(the, T, fL2, c1, c2)
250 : real(dp), intent(in) :: the, T, fL2, c1, c2
251 2855 : f1_the_T_fL2 = c1 * pow4(T) * pow3(the) / (1.0_dp - fL2) - the**2 + c2 * T
252 2855 : end function f1_the_T_fL2
253 :
254 151 : real(dp) function kap(rho, T)
255 : ! simplified Kramers rule (cgs; approximate)
256 : real(dp), intent(in) :: rho, T
257 151 : kap = 0.34_dp + 3.0d24 * rho * pow(T, -3.5_dp)
258 151 : end function kap
259 :
260 151 : real(dp) function f2_the_T_fL2(the, T, fL2, &
261 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
262 : real(dp), intent(in) :: the, T, fL2
263 : real(dp), intent(in) :: c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2
264 151 : real(dp) :: x, U_over_P, rho
265 151 : x = c4 * pow3(T * the) / (1.0_dp - fL2)
266 151 : U_over_P = (1.5_dp + x) / (1.0_dp + one_third * x)
267 151 : rho = (1.0_dp - fL2) * M1dot / (2.0_dp * pi * disk_alpha * omega_K * pow3(Rd)) / pow2(the)
268 : f2_the_T_fL2 = &
269 : seven_fourths &
270 : - (1.5_dp * U_over_P + c3 * pow4(T) / kap(rho, T) / (1.0_dp - fL2) ** 2) * the**2 &
271 : - PhiRd / (GM2 / Rd) &
272 151 : + (PhiL1 - fL2 * PhiL2) / (GM2 / Rd) / (1.0_dp - fL2)
273 151 : end function f2_the_T_fL2
274 :
275 170 : real(dp) function T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
276 : real(dp), intent(in) :: the
277 : integer, intent(in) :: n_the
278 : real(dp), intent(in) :: logthe_grid(n_the), logT_arr(n_the)
279 : real(dp), intent(in) :: dlogthe
280 170 : real(dp) :: logthe, slope, logT
281 : integer :: i_the
282 : ! under the assumption fL2=0
283 170 : logthe = log10(the)
284 170 : if (logthe > logthe_grid(n_the-1)) then
285 : ! use analytic extrapolation
286 6 : T_the_nofL2 = exp10(logT_arr(n_the-1) - 0.25_dp * (logthe - logthe_grid(n_the-1)))
287 164 : else if (logthe < logthe_grid(1)) then
288 : ! analytic extrapolation
289 0 : T_the_nofL2 = exp10(logT_arr(1) + 2.0_dp * (logthe - logthe_grid(1)))
290 : else
291 164 : i_the = floor((logthe - logthe_grid(1)) / dlogthe) + 1
292 164 : slope = (logT_arr(i_the + 1) - logT_arr(i_the)) / dlogthe
293 164 : logT = logT_arr(i_the) + (logthe - logthe_grid(i_the)) * slope
294 164 : T_the_nofL2 = exp10(logT)
295 : end if
296 170 : end function T_the_nofL2
297 :
298 150 : real(dp) function T_the(the, c2, PhiL2, PhiRd, GM2, Rd)
299 : real(dp), intent(in) :: the
300 : real(dp), intent(in) :: c2, PhiL2, PhiRd, GM2, Rd
301 : ! only for non-zero fL2
302 104 : T_the = (8.0_dp * the**2 + 1.0_dp - 2.0_dp * (PhiL2 - PhiRd) / (GM2 / Rd)) / (3.0_dp * c2)
303 0 : end function T_the
304 :
305 46 : real(dp) function fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd)
306 : real(dp), intent(in) :: the
307 : real(dp), intent(in) :: c1, c2, PhiL2, PhiRd, GM2, Rd
308 46 : real(dp) :: T
309 : ! only for non-zero fL2
310 46 : T = T_the(the, c2, PhiL2, PhiRd, GM2, Rd)
311 46 : fL2_the = 1.0_dp - c1 * pow4(T) * pow3(the) / (the**2 - c2 * T)
312 46 : end function fL2_the
313 :
314 :
315 : end module binary_disk
|
