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_eighth = 1.0_dp / 8.0_dp
34 : real(dp), parameter :: three_eighths = 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 6 : real(dp) :: GM2, Rd, Phi_units, PhiL1, PhiL2, PhiRd, omega_K
64 3 : 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 3 : if (mass_transfer_rate < eps_small) then
80 : ! no mass transfer, so no L2 mass loss
81 1 : fL2 = 0.0_dp
82 1 : return
83 : end if
84 :
85 : ! key parameters
86 2 : M2 = accretor_mass * Msun
87 2 : M1dot = mass_transfer_rate * Msun/secyer
88 2 : a = orbital_separation * Rsun
89 2 : q = accretor_mass / donor_mass ! mass ratio M2/M1
90 :
91 2 : log_q = log10(q)
92 :
93 : ! positions of Lagrangian points (based on analytic fits)
94 2 : xL1 = -0.0355_dp * log_q**2 + 0.251_dp * abs(log_q) + 0.500_dp ! [a = SMA]
95 2 : xL2 = 0.0756_dp * log_q**2 - 0.424_dp * abs(log_q) + 1.699_dp ! [a]
96 :
97 2 : if (log_q > 0.0_dp) then ! m2 is more massive
98 0 : xL1 = 1.0_dp - xL1
99 0 : xL2 = 1.0_dp - xL2
100 : end if
101 2 : mu = q / (1.0_dp + q)
102 :
103 : ! outer disk radius
104 2 : Rd_over_a = pow4(1.0_dp - xL1) / mu
105 : ! relevant potential energies
106 2 : PhiL1_dimless = -((1.0_dp - mu)/abs(xL1) + mu/abs(1.0_dp - xL1) + 0.5_dp*(xL1 - mu)**2) ! [G(M1+M2)/a]
107 2 : PhiL2_dimless = -((1.0_dp - mu)/abs(xL2) + mu/abs(1.0_dp - xL2) + 0.5_dp*(xL2 - mu)**2) ! [G(M1+M2)/a]
108 2 : PhiRd_dimless = -(1.0_dp - mu + mu/Rd_over_a + 0.5_dp*(1.0_dp - mu)**2)
109 :
110 2 : GM2 = standard_cgrav * M2
111 :
112 2 : Rd = Rd_over_a*a
113 2 : Phi_units = standard_cgrav * (M2 / mu) / a
114 2 : PhiL1 = PhiL1_dimless * Phi_units
115 2 : PhiL2 = PhiL2_dimless * Phi_units
116 2 : PhiRd = PhiRd_dimless * Phi_units
117 : ! Keplerian frequency at Rd
118 2 : omega_K = sqrt(GM2 / pow3(Rd))
119 :
120 :
121 : ! constants involved in numerical solutions
122 2 : c1 = two_thirds * pi * crad * disk_alpha * Rd / (omega_K * M1dot)
123 2 : c2 = boltzm * Rd / (GM2 * disk_mu * mp)
124 2 : c3 = 8.0_dp * pi2 * crad * disk_alpha * clight * Rd**2 / (M1dot**2 * omega_K)
125 2 : c4 = 2.0_dp * pi * disk_mu * crad * disk_alpha * omega_K * mp * pow3(Rd) / (boltzm * M1dot)
126 :
127 : ! Create logarithmically spaced grid for grid search for disk thickness [only used at the beginning]
128 102 : do i = 1, n_the
129 102 : the_grid(i) = exp10(log10(the_grid_min) + (i - 1) * (log10(the_grid_max) - log10(the_grid_min)) / (n_the - 1))
130 : end do
131 :
132 : ! only T < T_max is possible to calculate
133 2 : T_max = pow(four_twentyseventh / (c1**2 * c2), one_ninth)
134 :
135 102 : do i = 1, n_the
136 102 : T_arr(i) = 0.0_dp
137 : end do
138 :
139 102 : do i = 1, n_the
140 100 : the = the_grid(i)
141 : ! use bisection method
142 100 : T_left = 0.1_dp * min(the**2 / c2, T_max)
143 100 : f1_left = f1_the_T_fL2(the, T_left, 0.0_dp, c1, c2)
144 100 : T_right = T_max
145 2805 : do while (abs((T_left - T_right) / T_right) > tol)
146 2705 : T = 0.5_dp * (T_left + T_right)
147 2705 : f1 = f1_the_T_fL2(the, T, 0.0_dp, c1, c2)
148 2805 : if (f1 * f1_left > 0) then
149 1462 : T_left = T
150 1462 : f1_left = f1
151 : else
152 : T_right = T
153 : end if
154 : end do
155 102 : T_arr(i) = 0.5_dp * (T_left + T_right)
156 : end do
157 : ! now we have obtained numerical relation between the and T
158 102 : do i = 1, n_the
159 100 : logT_arr(i) = log10(T_arr(i))
160 102 : logthe_grid(i) = log10(the_grid(i))
161 : end do
162 2 : dlogthe = logthe_grid(2) - logthe_grid(1)
163 :
164 : ! bisection to find the numerical solution to f2(the, T, fL2=0)=0
165 2 : the_right = 1.0_dp
166 : f2_right = f2_the_T_fL2(the_right, T_the_nofL2(the_right, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
167 2 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
168 2 : separation_factor = 0.95_dp
169 2 : the_left = separation_factor * the_right
170 : f2_left = f2_the_T_fL2(the_left, T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
171 2 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
172 34 : do while (f2_left * f2_right > 0)
173 : ! need to decrease the_left
174 32 : the_right = the_left
175 32 : f2_right = f2_left
176 32 : the_left = the_left * separation_factor
177 : f2_left = f2_the_T_fL2(the_left, T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
178 34 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
179 : end do
180 : ! now the solution is between the_left and the_right
181 48 : do while (abs((the_left - the_right) / the_right) > tol)
182 46 : the = 0.5_dp * (the_left + the_right)
183 : f2 = f2_the_T_fL2(the, T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe), 0.0_dp, &
184 46 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
185 48 : if (f2 * f2_left > 0.0_dp) then
186 27 : the_left = the
187 27 : f2_left = f2
188 : else
189 19 : the_right = the
190 : end if
191 : end do
192 : ! solution
193 2 : the = 0.5_dp * (the_left + the_right)
194 2 : T = T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
195 2 : the_max = sqrt(three_eighths * c2 * T + one_fourth * (PhiL2 - PhiRd) / (GM2 / Rd) - one_eighth)
196 :
197 2 : if (the < the_max) then
198 : ! return a tiny number
199 0 : fL2 = eps_small
200 : else
201 2 : the_min = ( 0.5_dp * sqrt((PhiL2 - PhiRd) / (GM2 / Rd) - 0.5_dp) ) ! corresponding to fL2=1, T=0
202 : ! need to find the maximum corresponding to fL2=0
203 : ! this is given by the intersection between T_the(the), T_the_nofL2(the)
204 2 : the_left = the_min
205 2 : the_right = 1.0_dp
206 2 : f_left = T_the(the_left, c2, PhiL2, PhiRd, GM2, Rd) - T_the_nofL2(the_left, logthe_grid, logT_arr, n_the, dlogthe)
207 60 : do while (abs((the_left - the_right) / the_right) > tol)
208 58 : the = 0.5_dp * (the_left + the_right)
209 58 : f = T_the(the, c2, PhiL2, PhiRd, GM2, Rd) - T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
210 2 : if (f * f_left > 0.0_dp) then
211 19 : the_left = the
212 19 : f_left = f
213 : else
214 39 : the_right = the
215 : end if
216 : end do
217 2 : the_max = 0.5_dp * (the_left + the_right) ! this corresponds to fL2=0
218 :
219 : ! --- numerical solution for f2(the, T, fL2)=0 under non-zero fL2
220 :
221 : ! -- do not use exactly the_min (corresponding to T = 0, bc. kap table breaks down)
222 : ! -- define another the_min based on T_floor (kap table won't be a problem)
223 : the_min = sqrt( three_eighths * c2 * T_floor &
224 2 : + one_fourth * (PhiL2 - PhiRd) / (GM2 / Rd) - one_eighth )
225 2 : the_left = the_min
226 : f2_left = f2_the_T_fL2(the_left, &
227 : T_the(the_left, c2, PhiL2, PhiRd, GM2, Rd), &
228 : fL2_the(the_left, c1, c2, PhiL2, PhiRd, GM2, Rd), &
229 2 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
230 2 : the_right = the_max / (1.0_dp + eps_small)
231 : ! bisection again
232 44 : do while (abs((the_left - the_right) / the_right) > tol)
233 42 : the = 0.5_dp * (the_left + the_right)
234 : f2 = f2_the_T_fL2(the, T_the(the, c2, PhiL2, PhiRd, GM2, Rd), fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd), &
235 42 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
236 2 : if (f2 * f2_left > 0.0_dp) then
237 22 : the_left = the
238 22 : f2_left = f2
239 : else
240 20 : the_right = the
241 : end if
242 : end do
243 : ! solution
244 2 : the = 0.5_dp * (the_left + the_right)
245 2 : fL2 = fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd)
246 : end if
247 :
248 : end subroutine eval_L2_mass_loss_fraction
249 :
250 :
251 : ! Helper Functions
252 2805 : real(dp) function f1_the_T_fL2(the, T, fL2, c1, c2)
253 : real(dp), intent(in) :: the, T, fL2, c1, c2
254 2805 : f1_the_T_fL2 = c1 * pow4(T) * pow3(the) / (1.0_dp - fL2) - the**2 + c2 * T
255 2805 : end function f1_the_T_fL2
256 :
257 126 : real(dp) function kap(rho, T)
258 : ! simplified Kramers rule (cgs; approximate)
259 : real(dp), intent(in) :: rho, T
260 126 : kap = 0.34_dp + 3.0d24 * rho * pow(T, -3.5_dp)
261 126 : end function kap
262 :
263 126 : real(dp) function f2_the_T_fL2(the, T, fL2, &
264 : c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2)
265 : real(dp), intent(in) :: the, T, fL2
266 : real(dp), intent(in) :: c3, c4, M1dot, disk_alpha, omega_K, Rd, PhiRd, GM2, PhiL1, PhiL2
267 126 : real(dp) :: x, U_over_P, rho
268 126 : x = c4 * pow3(T * the) / (1.0_dp - fL2)
269 126 : U_over_P = (1.5_dp + x) / (1.0_dp + one_third * x)
270 126 : rho = (1.0_dp - fL2) * M1dot / (2.0_dp * pi * disk_alpha * omega_K * pow3(Rd)) / pow2(the)
271 : f2_the_T_fL2 = &
272 : seven_fourths &
273 : - (1.5_dp * U_over_P + c3 * pow4(T) / kap(rho, T) / (1.0_dp - fL2) ** 2) * the**2 &
274 : - PhiRd / (GM2 / Rd) &
275 126 : + (PhiL1 - fL2 * PhiL2) / (GM2 / Rd) / (1.0_dp - fL2)
276 126 : end function f2_the_T_fL2
277 :
278 144 : real(dp) function T_the_nofL2(the, logthe_grid, logT_arr, n_the, dlogthe)
279 : real(dp), intent(in) :: the
280 : integer, intent(in) :: n_the
281 : real(dp), intent(in) :: logthe_grid(n_the), logT_arr(n_the)
282 : real(dp), intent(in) :: dlogthe
283 144 : real(dp) :: logthe, slope, logT
284 : integer :: i_the
285 : ! under the assumption fL2=0
286 144 : logthe = log10(the)
287 144 : if (logthe > logthe_grid(n_the-1)) then
288 : ! use analytic extrapolation
289 2 : T_the_nofL2 = exp10(logT_arr(n_the-1) - 0.25_dp * (logthe - logthe_grid(n_the-1)))
290 142 : else if (logthe < logthe_grid(1)) then
291 : ! analytic extrapolation
292 0 : T_the_nofL2 = exp10(logT_arr(1) + 2.0_dp * (logthe - logthe_grid(1)))
293 : else
294 142 : i_the = floor((logthe - logthe_grid(1)) / dlogthe) + 1
295 142 : slope = (logT_arr(i_the + 1) - logT_arr(i_the)) / dlogthe
296 142 : logT = logT_arr(i_the) + (logthe - logthe_grid(i_the)) * slope
297 142 : T_the_nofL2 = exp10(logT)
298 : end if
299 144 : end function T_the_nofL2
300 :
301 150 : real(dp) function T_the(the, c2, PhiL2, PhiRd, GM2, Rd)
302 : real(dp), intent(in) :: the
303 : real(dp), intent(in) :: c2, PhiL2, PhiRd, GM2, Rd
304 : ! only for non-zero fL2
305 104 : T_the = (8.0_dp * the**2 + 1.0_dp - 2.0_dp * (PhiL2 - PhiRd) / (GM2 / Rd)) / (3.0_dp * c2)
306 0 : end function T_the
307 :
308 46 : real(dp) function fL2_the(the, c1, c2, PhiL2, PhiRd, GM2, Rd)
309 : real(dp), intent(in) :: the
310 : real(dp), intent(in) :: c1, c2, PhiL2, PhiRd, GM2, Rd
311 46 : real(dp) :: T
312 : ! only for non-zero fL2
313 46 : T = T_the(the, c2, PhiL2, PhiRd, GM2, Rd)
314 46 : fL2_the = 1.0_dp - c1 * pow4(T) * pow3(the) / (the**2 - c2 * T)
315 46 : end function fL2_the
316 :
317 :
318 : end module binary_disk
|
