-
Notifications
You must be signed in to change notification settings - Fork 11
/
Copy pathinteractions.py
1378 lines (1093 loc) · 63.5 KB
/
interactions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
from amuse.units import units, constants, quantities
import numpy as np
import sys
import scipy.integrate as integrate
from TRES_options import REPORT_BINARY_EVOLUTION, REPORT_FUNCTION_NAMES, REPORT_MASS_TRANSFER_STABILITY
#constants
numerical_error = 1.e-6
small_numerical_error = 1.e-10
minimum_eccentricity = 1.e-5
const_common_envelope_efficiency = 4.0 #1.0, 4 for now for easier testing with SeBa
const_envelope_structure_parameter = 0.5
const_common_envelope_efficiency_gamma = 1.75
stellar_types_compact_objects = [10,11,12,13,14]|units.stellar_type
stellar_types_giants = [2,3,4,5,6,8,9]|units.stellar_type
stellar_types_planetary_objects = [18,19]|units.stellar_type # planets & brown dwarfs
stellar_types_SN_remnants = [13,14,15]|units.stellar_type # remnant types created through a supernova
stellar_types_remnants = [7,8,9,10,11,12,13,14,15]|units.stellar_type
stellar_types_dr = [2,4,7,8,9,10,11,12,13,14,15]|units.stellar_type #stars which go through a instantaneous radius change at formation; hertzsprung gap stars (small envelope perturbation) + horizontal branch stars + remnants
#q_crit = 3.
#q_crit_giants_conv_env = 0.9
nucleair_efficiency = 0.007 # nuc. energy production eff, Delta E = 0.007 Mc^2
#dictionaries
bin_type = {
'unknown': 'unknown',
'merger': 'merger',
'disintegrated': 'disintegrated',
'dyn_inst': 'dynamical_instability',
'detached': 'detached',
'contact': 'contact',
'collision': 'collision',
'semisecular': 'semisecular',
'rlof': 'rlof', #only used for stopping conditions
'olof' : 'olof', #only used for stopping conditions
'stable_mass_transfer': 'stable_mass_transfer',
'common_envelope': 'common_envelope',
'common_envelope_energy_balance': 'common_envelope_energy_balance',
'common_envelope_angular_momentum_balance': 'common_envelope_angular_momentum_balance',
'double_common_envelope': 'double_common_envelope',
}
#-------------------------
#general functions
def roche_radius_dimensionless(M, m):
# Assure that the q is calculated in identical units.
unit = M.unit
# and that q itself has no unit
q = M.value_in(unit)/m.value_in(unit)
q13 = q**(1./3.)
q23 = q13**2
return 0.49*q23/(0.6*q23 + np.log(1 + q13))
def roche_radius(bin, primary, self):
if not bin.is_star and primary.is_star:
return bin.semimajor_axis * roche_radius_dimensionless(primary.mass, self.get_mass(bin)-primary.mass)
sys.exit('error in roche radius: Roche radius can only be determined in a binary')
def L2_radius_dimensionless(M,m):
# approximation for l2 overflow
# see Marchant+ 2016 equation 2
q = M/m
rl1 = roche_radius_dimensionless(M, m)
rl2_div_rl1 = 0.299 * np.arctan(1.83*q**0.397) + 1
return rl2_div_rl1 * rl1
def L2_radius(bin, primary, self):
# note: this prescription is based on the Eggleton approximation for how to adjust a circular RL to an eccentric one
# may not be consistent with Sepinsky's method for eccentric RL (L1)
if not bin.is_star and primary.is_star:
return bin.semimajor_axis * L2_radius_dimensionless(primary.mass, self.get_mass(bin)-primary.mass)*(1-bin.eccentricity)
sys.exit('Error: L2 radius can only be determined in a binary')
#for comparison with kozai timescale
def stellar_evolution_timescale(star):
if REPORT_FUNCTION_NAMES:
print("Stellar evolution timescale")
if star.stellar_type in [0,1,7]|units.stellar_type:
return (0.1 * star.mass * nucleair_efficiency * constants.c**2 / star.luminosity).in_(units.Gyr)
elif star.stellar_type in stellar_types_compact_objects:
return np.inf|units.Myr
elif star.stellar_type in stellar_types_planetary_objects:
return np.inf|units.Myr
else:
return 0.1*star.age
#for mass transfer rate
def nuclear_evolution_timescale(star):
if REPORT_FUNCTION_NAMES:
print("Nuclear evolution timescale:")
if star.stellar_type in [0,1,7]|units.stellar_type:
return (0.1 * star.mass * nucleair_efficiency * constants.c**2 / star.luminosity).in_(units.Gyr)
elif star.stellar_type in stellar_types_planetary_objects:
# print('nuclear evolution timescale for planetary objects requested')
# return np.inf|units.Myr
return dynamic_timescale(star)
else: #t_nuc ~ delta t * R/ delta R, other prescription gave long timescales in SeBa which destables the mass transfer
if star.time_derivative_of_radius <= (quantities.zero+numerical_error**2)|units.RSun/units.yr:
#when star is shrinking
# t_nuc = 0.1*main_sequence_time() # in SeBa
t_nuc = 0.1*star.age
else:
t_nuc = star.radius / star.time_derivative_of_radius #does not include the effect of mass loss on R
return t_nuc
def kelvin_helmholds_timescale(star):
if star.stellar_type in stellar_types_planetary_objects:
# print('thermal evolution timescale for planetary objects requested')
return dynamic_timescale(star)
if REPORT_FUNCTION_NAMES:
print("KH timescale:", (constants.G*star.mass**2/star.radius/star.luminosity).in_(units.Myr))
return constants.G*star.mass**2/star.radius/star.luminosity
def dynamic_timescale(star):
if REPORT_FUNCTION_NAMES:
print("Dynamic timescale:", (np.sqrt(star.radius**3/star.mass/constants.G)[0]).in_(units.yr))
return np.sqrt(star.radius**3/star.mass/constants.G)
def corotating_spin_angular_frequency_binary(semi, m1, m2):
return 1./np.sqrt(semi**3/constants.G / (m1+m2))
#Hurley, Pols en Tout 2000, eq 107-108
def lang_spin_angular_frequency(star):
v_rot = 330*star.mass.value_in(units.MSun)**3.3/(15.0+star.mass.value_in(units.MSun)**3.45)
w = 45.35 * v_rot/star.radius.value_in(units.RSun)
return w|1./units.yr
def break_up_angular_frequency(object):
return np.sqrt( constants.G * object.mass / object.radius ) / object.radius
def criticial_angular_frequency_CHE(m, Z):
#angular frequency of spin for CHE threshold
#Fitting formula for CHE from Riley+ 2021
a_coeff = np.array([5.7914 * 10 ** - 4, -1.9196 * 10 ** - 6,
-4.0602 * 10 ** - 7, 1.0150 * 10 ** - 8,
-9.1792 * 10 ** - 11, 2.9051 * 10 ** - 13])
mass_power = np.linspace(0,5,6)
omega_at_z_0d004 = np.sum(a_coeff * m.value_in(units.MSun)** mass_power / m.value_in(units.MSun) ** 0.4)
omega_crit = omega_at_z_0d004/(0.09 * np.log(Z/0.004) + 1) |1./units.s
return omega_crit
def copy_outer_orbit_to_inner_orbit(bs, self):
if REPORT_FUNCTION_NAMES:
print('Copy_outer_orbit_to_inner_orbit')
if self.is_triple():
bs.semimajor_axis = self.triple.semimajor_axis
bs.eccentricity = self.triple.eccentricity
bs.argument_of_pericenter = self.triple.argument_of_pericenter
bs.longitude_of_ascending_node = self.triple.longitude_of_ascending_node
bs.mass_transfer_rate = self.triple.mass_transfer_rate
bs.accretion_efficiency_mass_transfer, = self.triple.accretion_efficiency_mass_transfer,
bs.accretion_efficiency_wind_child1_to_child2, = self.triple.accretion_efficiency_wind_child1_to_child2,
bs.accretion_efficiency_wind_child2_to_child1, = self.triple.accretion_efficiency_wind_child2_to_child1,
bs.specific_AM_loss_mass_transfer, = self.triple.specific_AM_loss_mass_transfer,
bs.is_mt_stable = self.triple.is_mt_stable
self.triple.semimajor_axis = 1e100|units.RSun
self.triple.eccentricity = 0
def copy_outer_star_to_accretor(self):
if REPORT_FUNCTION_NAMES:
print('Copy_outer_star_to_accretor')
if self.is_triple():
if self.triple.child1.is_star:
tertiary_star = self.triple.child1
bs = self.triple.child2
else:
tertiary_star = self.triple.child2
bs = self.triple.child1
if not bs.child1.is_donor:
bs.child1 = tertiary_star
else:
bs.child2 = tertiary_star
#-------------------------
#-------------------------
# functions for mass transfer in a binary
def perform_inner_collision(self):
if self.is_triple():
if self.triple.child1.is_star:
self.triple.child2
else:
self.triple.child1
# smaller star is added to big star
if bs.child1.radius >= bs.child2.radius:
donor = bs.child1
accretor = bs.child2
else:
donor = bs.child2
accretor = bs.child1
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
accretor_in_stellar_code = accretor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
#no additional mass and Jspin loss from merged object for now
J_spin_donor_previous = self.spin_angular_momentum(donor)
J_spin_accretor_previous = self.spin_angular_momentum(accretor)
J_orbit = self.orbital_angular_momentum(bs)
J_spin_new = J_spin_donor_previous + J_spin_accretor_previous + J_orbit
#merger
donor_in_stellar_code.merge_with_other_star(accretor_in_stellar_code)
self.copy_from_stellar()
donor.moment_of_inertia_of_star = self.moment_of_inertia(donor)
#assuming conservation of total angular momentum of the inner binary
spin_angular_frequency = J_spin_new / donor.moment_of_inertia_of_star
critical_spin_angular_frequency = np.sqrt(constants.G * donor.mass/donor.radius**3)
donor.spin_angular_frequency = min(donor.spin_angular_frequency, critical_spin_angular_frequency)
self.stellar_code.particles.remove_particle(accretor)
accretor.mass = 0|units.MSun # necessary for adjust_system_after_ce_in_inner_binary
#adjust outer orbit, needs to be before the system becomes a binary
#and copy to inner orbit
# weird structure necessary for secular code -> outer orbit is redundant
adjust_system_after_ce_in_inner_binary(bs, self)
copy_outer_orbit_to_inner_orbit(bs, self)
copy_outer_star_to_accretor(self)
#functions are skipped in binaries, needs to be checked if this works well
self.secular_code.parameters.ignore_tertiary = True
self.secular_code.parameters.check_for_dynamical_stability = False
self.secular_code.parameters.check_for_outer_collision = False
self.secular_code.parameters.check_for_outer_RLOF = False
bs.bin_type = bin_type['collision']
self.save_snapshot() # just to note that it the system has merged
#use of stopping condition in this way (similar to perform inner merger) is not necessary.
#TRES.py takes care of it
# if self.check_stopping_conditions_stellar_interaction()==False:
# print('stopping conditions stellar interaction')
# return False
self.check_RLOF()
if self.has_donor():
print(self.triple.child2.child1.mass, self.triple.child2.child2.mass, self.triple.child2.semimajor_axis, self.triple.child2.eccentricity, self.triple.child2.child1.is_donor, self.triple.child2.child2.is_donor)
print(self.triple.child1.mass, self.triple.semimajor_axis, self.triple.eccentricity, self.triple.child1.is_donor)
sys.exit("error in adjusting triple after collision: RLOF")
donor.is_donor = False
bs.is_mt_stable = True
bs.bin_type = bin_type['detached']
donor.spin_angular_frequency = corotating_spin_angular_frequency_binary(bs.semimajor_axis, bs.child1.mass, bs.child2.mass)
#use of stopping condition in this way (similar to perform inner merger) is not necessary
#TRES.py takes care of it
# return True
def perform_inner_merger(bs, donor, accretor, self):
if REPORT_BINARY_EVOLUTION:
print('Merger in inner binary through common envelope phase')
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
accretor_in_stellar_code = accretor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
#no additional mass and Jspin loss from merged object for now
J_spin_donor_previous = self.spin_angular_momentum(donor)
J_spin_accretor_previous = self.spin_angular_momentum(accretor)
J_orbit = self.orbital_angular_momentum(bs)
J_spin_new = J_spin_donor_previous + J_spin_accretor_previous + J_orbit
#merger
donor_in_stellar_code.merge_with_other_star(accretor_in_stellar_code)
self.copy_from_stellar()
donor.moment_of_inertia_of_star = self.moment_of_inertia(donor)
#assuming conservation of total angular momentum of the inner binary
spin_angular_momentum = J_spin_new / donor.moment_of_inertia_of_star
critical_spin_angular_frequency = np.sqrt(constants.G * donor.mass/donor.radius**3)
donor.spin_angular_frequency = min(donor.spin_angular_frequency, critical_spin_angular_frequency)
self.stellar_code.particles.remove_particle(accretor)
accretor.mass = 0|units.MSun # necessary for adjust_system_after_ce_in_inner_binary
#adjust outer orbit, needs to be before the system becomes a binary
#and copy to inner orbit
# weird structure necessary for secular code -> outer orbit is redundant
adjust_system_after_ce_in_inner_binary(bs, self)
copy_outer_orbit_to_inner_orbit(bs, self)
copy_outer_star_to_accretor(self)
#functions are skipped in binaries, needs to be checked if this works well
self.secular_code.parameters.ignore_tertiary = True
self.secular_code.parameters.check_for_dynamical_stability = False
self.secular_code.parameters.check_for_outer_collision = False
self.secular_code.parameters.check_for_outer_RLOF = False
bs.bin_type = bin_type['merger']
self.save_snapshot() # just to note that the system has merged
if self.check_stopping_conditions_stellar_interaction()==False:
print('stopping conditions stellar interaction')
return False
else:
return True
# print(self.secular_code.give_roche_radii(self.triple),)
# print(roche_radius(self.triple.child2, self.triple.child2.child1, self), roche_radius(self.triple.child2, self.triple.child2.child2, self))
#
# print(donor.spin_angular_frequency, corotating_spin_angular_frequency_binary(bs.semimajor_axis, bs.child1.mass, bs.child2.mass), critical_spin_angular_frequency)
# donor.spin_angular_frequency = corotating_spin_angular_frequency_binary(bs.semimajor_axis, bs.child1.mass, bs.child2.mass)
def common_envelope_efficiency(donor, accretor):
return const_common_envelope_efficiency
def envelope_structure_parameter(donor):
return const_envelope_structure_parameter
def common_envelope_efficiency_gamma(donor, accretor):
return const_common_envelope_efficiency_gamma
# ang.mom balance: \Delta J = \gamma * J * \Delta M / M
# See Eq. 5 of Nelemans VYPZ 2000, 360, 1011 A&A
def common_envelope_angular_momentum_balance(bs, donor, accretor, self):
if REPORT_FUNCTION_NAMES:
print('Common envelope angular momentum balance')
if REPORT_BINARY_EVOLUTION:
if bs.eccentricity > 0.05:
print('gamma common envelope in eccentric binary')
print('Before common envelope angular momentum balance' )
self.print_binary(bs)
bs.bin_type = bin_type['common_envelope_angular_momentum_balance']
self.save_snapshot()
gamma = common_envelope_efficiency_gamma(donor, accretor)
J_init = np.sqrt(bs.semimajor_axis) * (donor.mass * accretor.mass) / np.sqrt(donor.mass + accretor.mass) * np.sqrt(1-bs.eccentricity**2)
J_f_over_sqrt_a_new = (donor.core_mass * accretor.mass) / np.sqrt(donor.core_mass + accretor.mass)
J_lost = gamma * (donor.mass-donor.core_mass) * J_init/(donor.mass + accretor.mass)
sqrt_a_new = max(0.|units.RSun**0.5, (J_init -J_lost)/J_f_over_sqrt_a_new)
a_new = pow(sqrt_a_new, 2)
Rl_donor_new = roche_radius_dimensionless(donor.core_mass, accretor.mass)*a_new
Rl_accretor_new = roche_radius_dimensionless(accretor.mass, donor.core_mass)*a_new
if REPORT_BINARY_EVOLUTION:
print('donor:', donor.radius, donor.core_radius, Rl_donor_new)
print('accretor:', accretor.radius, accretor.core_radius, Rl_accretor_new)
if (donor.core_radius > Rl_donor_new) or (accretor.radius > Rl_accretor_new):
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
else:
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
#reduce_mass not subtrac mass, want geen adjust_donor_radius
#check if star changes type
donor_in_stellar_code.change_mass(-1*(donor.mass-donor.core_mass+(small_numerical_error|units.MSun)), 0.|units.yr)
self.copy_from_stellar()
donor.moment_of_inertia_of_star = self.moment_of_inertia(donor)
accretor.moment_of_inertia_of_star = self.moment_of_inertia(accretor)
bs.semimajor_axis = a_new
bs.eccentricity = minimum_eccentricity
#set to synchronization
corotating_frequency = corotating_spin_angular_frequency_binary(a_new, donor.mass, accretor.mass)
donor.spin_angular_frequency = corotating_frequency
accretor.spin_angular_frequency = corotating_frequency
self.check_RLOF()
if self.has_donor():
print(self.triple.child2.child1.mass, self.triple.child2.child2.mass, self.triple.child2.child1.radius, self.triple.child2.child2.radius,self.triple.child2.semimajor_axis, self.triple.child2.eccentricity, self.triple.child2.child1.is_donor, self.triple.child2.child2.is_donor)
print(self.triple.child2.child1.core_mass, self.triple.child2.child1.mass-self.triple.child2.child1.core_mass, self.triple.child2.child1.stellar_type)
print(self.triple.child1.mass, self.triple.semimajor_axis, self.triple.eccentricity, self.triple.child1.is_donor)
# sys.exit("error in adjusting triple after gamma CE: RLOF")
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
# adjusting of stellar system
# in previous case of merger, the adjustment is done there as mass may be lost during the merger
adjust_system_after_ce_in_inner_binary(bs, self)
donor.is_donor = False
bs.is_mt_stable = True
bs.bin_type = bin_type['detached']
self.instantaneous_evolution = True #skip secular evolution
return True
#Following Webbink 1984
def common_envelope_energy_balance(bs, donor, accretor, self):
if REPORT_FUNCTION_NAMES:
print('Common envelope energy balance')
if REPORT_BINARY_EVOLUTION:
print('Before common envelope energy balance' )
self.print_binary(bs)
bs.bin_type = bin_type['common_envelope_energy_balance']
self.save_snapshot()
alpha = common_envelope_efficiency(donor, accretor)
lambda_donor = envelope_structure_parameter(donor)
Rl_donor = roche_radius(bs, donor, self)
donor_radius = min(donor.radius, Rl_donor)
#based on Glanz & Perets 2021 2021MNRAS.507.2659G
#eccentric CE -> end result depends on pericenter distance more than semi-major axis
pericenter_init = bs.semimajor_axis * (1-bs.eccentricity)
orb_energy_new = donor.mass * (donor.mass-donor.core_mass) / (alpha * lambda_donor * donor_radius) + donor.mass * accretor.mass/2/pericenter_init
a_new = donor.core_mass * accretor.mass / 2 / orb_energy_new
# a_new = bs.semimajor_axis * (donor.core_mass/donor.mass) / (1. + (2.*(donor.mass-donor.core_mass)*bs.semimajor_axis/(alpha_lambda*donor_radius*accretor.mass)))
Rl_donor_new = roche_radius_dimensionless(donor.core_mass, accretor.mass)*a_new
Rl_accretor_new = roche_radius_dimensionless(accretor.mass, donor.core_mass)*a_new
if REPORT_BINARY_EVOLUTION:
print('donor:', donor.radius, donor.core_radius, Rl_donor_new)
print('accretor:', accretor.radius, accretor.core_radius, Rl_accretor_new)
if (donor.core_radius > Rl_donor_new) or (accretor.radius > Rl_accretor_new):
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
else:
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
#reduce_mass not subtrac mass, want geen adjust_donor_radius
#check if star changes type
donor_in_stellar_code.change_mass(-1*(donor.mass-donor.core_mass+(small_numerical_error|units.MSun)), 0.|units.yr)
self.copy_from_stellar()
donor.moment_of_inertia_of_star = self.moment_of_inertia(donor)
accretor.moment_of_inertia_of_star = self.moment_of_inertia(accretor)
bs.semimajor_axis = a_new
bs.eccentricity = minimum_eccentricity
#set to synchronization
corotating_frequency = corotating_spin_angular_frequency_binary(a_new, donor.mass, accretor.mass)
donor.spin_angular_frequency = corotating_frequency
accretor.spin_angular_frequency = corotating_frequency
self.check_RLOF()
if self.has_donor():
print(self.triple.child2.child1.mass, self.triple.child2.child2.mass, self.triple.child2.semimajor_axis, self.triple.child2.eccentricity, self.triple.child2.child1.is_donor, self.triple.child2.child2.is_donor)
print(self.triple.child1.mass, self.triple.semimajor_axis, self.triple.eccentricity, self.triple.child1.is_donor)
# sys.exit("error in adjusting triple after alpha CE: RLOF")
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
# adjusting of stellar system
# in previous case of merger, the adjustment is done there as mass may be lost during the merger
adjust_system_after_ce_in_inner_binary(bs, self)
donor.is_donor = False
bs.is_mt_stable = True
bs.bin_type = bin_type['detached']
self.instantaneous_evolution = True #skip secular evolution
return True
# See appendix of Nelemans YPZV 2001, 365, 491 A&A
def double_common_envelope_energy_balance(bs, donor, accretor, self):
if REPORT_FUNCTION_NAMES:
print('Double common envelope energy balance')
if REPORT_BINARY_EVOLUTION:
print('Before double common envelope energy balance' )
self.print_binary(bs)
bs.bin_type = bin_type['double_common_envelope']
self.save_snapshot()
alpha = common_envelope_efficiency(donor, accretor)
lambda_donor = envelope_structure_parameter(donor)
lambda_accretor = envelope_structure_parameter(accretor)
Rl_donor = roche_radius(bs, donor, self)
donor_radius = min(donor.radius, Rl_donor)
accretor_radius = accretor.radius
#based on Glanz & Perets 2021 2021MNRAS.507.2659G
#eccentric CE -> end result depends on pericenter distance more than semi-major axis
pericenter_init = bs.semimajor_axis * (1-bs.eccentricity)
orb_energy_new = donor.mass * (donor.mass-donor.core_mass) / (alpha * lambda_donor * donor_radius) + accretor.mass * (accretor.mass-accretor.core_mass) / (alpha * lambda_accretor * accretor_radius) + donor.mass * accretor.mass/2/pericenter_init
a_new = donor.core_mass * accretor.core_mass / 2 / orb_energy_new
Rl_donor_new = roche_radius_dimensionless(donor.core_mass, accretor.core_mass)*a_new
Rl_accretor_new = roche_radius_dimensionless(accretor.core_mass, donor.core_mass)*a_new
if REPORT_BINARY_EVOLUTION:
print('donor:', donor.radius, donor.core_radius, Rl_donor_new)
print('accretor:', accretor.radius, accretor.core_radius, Rl_accretor_new)
if (donor.core_radius > Rl_donor_new) or (accretor.core_radius > Rl_accretor_new):
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
else:
#reduce_mass not subtrac mass, want geen adjust_donor_radius
#check if star changes type
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
donor_in_stellar_code.change_mass(-1*(donor.mass-donor.core_mass+(small_numerical_error|units.MSun)), 0.|units.yr)
accretor_in_stellar_code = accretor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
accretor_in_stellar_code.change_mass(-1*(accretor.mass-accretor.core_mass), 0.|units.yr)
self.copy_from_stellar()
donor.moment_of_inertia_of_star = self.moment_of_inertia(donor)
accretor.moment_of_inertia_of_star = self.moment_of_inertia(accretor)
bs.semimajor_axis = a_new
bs.eccentricity = minimum_eccentricity
#set to synchronization
corotating_frequency = corotating_spin_angular_frequency_binary(a_new, donor.mass, accretor.mass)
donor.spin_angular_frequency = corotating_frequency
accretor.spin_angular_frequency = corotating_frequency
self.check_RLOF()
if self.has_donor():
print(self.triple.child2.child1.mass, self.triple.child2.child2.mass, self.triple.child2.semimajor_axis, self.triple.child2.eccentricity, self.triple.child2.child1.is_donor, self.triple.child2.child2.is_donor)
print(self.triple.child1.mass, self.triple.semimajor_axis, self.triple.eccentricity, self.triple.child1.is_donor)
# sys.exit("error in adjusting triple after double CE: RLOF")
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
# adjusting of stellar system
# in previous case of merger, the adjustment is done there as mass may be lost during the merger
adjust_system_after_ce_in_inner_binary(bs, self)
donor.is_donor = False
bs.is_mt_stable = True
bs.bin_type = bin_type['detached']
self.instantaneous_evolution = True #skip secular evolution
return True
def common_envelope_phase(bs, donor, accretor, self):
stopping_condition = True
if REPORT_FUNCTION_NAMES:
print('Common envelope phase', self.which_common_envelope)
print('donor:', donor.stellar_type)
print('accretor:', accretor.stellar_type)
if donor.stellar_type not in stellar_types_giants and accretor.stellar_type not in stellar_types_giants:
# possible options: MS+MS, MS+remnant, remnant+remnant,
# HeMS+HeMS, HeMS+MS, HeMS+remnant
bs.bin_type = bin_type['common_envelope']
self.save_snapshot()
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
self.check_RLOF()
if self.has_donor():
print(self.triple.child2.child1.mass, self.triple.child2.child2.mass, self.triple.child2.semimajor_axis, self.triple.child2.eccentricity, self.triple.child2.child1.is_donor, self.triple.child2.child2.is_donor)
print(self.triple.child1.mass, self.triple.semimajor_axis, self.triple.eccentricity, self.triple.child1.is_donor)
print(self.triple.child2.child1.radius, self.triple.child2.child2.radius,self.triple.child1.radius)
print(self.secular_code.give_roche_radii(self.triple))
print('binary Roche lobe radii:', roche_radius(bs, bs.child1, self), roche_radius(bs, bs.child2, self))
# sys.exit("error in adjusting system after CE: RLOF")
stopping_condition = perform_inner_merger(bs, donor, accretor, self)
if not stopping_condition: #stellar interaction
return False
donor.is_donor = False
bs.is_mt_stable = True
bs.bin_type = bin_type['detached']
self.instantaneous_evolution = True #skip secular evolution
return True
if self.which_common_envelope == 0:
if donor.stellar_type in stellar_types_giants and accretor.stellar_type in stellar_types_giants:
stopping_condition = double_common_envelope_energy_balance(bs, donor, accretor, self)
else:
stopping_condition = common_envelope_energy_balance(bs, donor, accretor, self)
elif self.which_common_envelope == 1:
if donor.stellar_type in stellar_types_giants and accretor.stellar_type in stellar_types_giants:
stopping_condition = double_common_envelope_energy_balance(bs, donor, accretor, self)
else:
stopping_condition = common_envelope_angular_momentum_balance(bs, donor, accretor, self)
elif self.which_common_envelope == 2:
Js_d = self.spin_angular_momentum(donor)
Js_a = self.spin_angular_momentum(accretor)
Jb = self.orbital_angular_momentum(bs)
Js = max(Js_d, Js_a)
# print("Darwin Riemann instability? donor/accretor:", Js_d, Js_a, Jb, Jb/3.)
if donor.stellar_type in stellar_types_giants and accretor.stellar_type in stellar_types_giants:
#giant+giant
stopping_condition = double_common_envelope_energy_balance(bs, donor, accretor, self)
elif donor.stellar_type in stellar_types_compact_objects or accretor.stellar_type in stellar_types_compact_objects:
#giant+remnant
stopping_condition = common_envelope_energy_balance(bs, donor, accretor, self)
elif Js >= Jb/3. :
#darwin riemann instability
stopping_condition = common_envelope_energy_balance(bs, donor, accretor, self)
else:
#giant+normal(non-giant, non-remnant)
stopping_condition = common_envelope_angular_momentum_balance(bs, donor, accretor, self)
return stopping_condition
def adiabatic_expansion_due_to_mass_loss(a_i, Md_f, Md_i, Ma_f, Ma_i):
d_Md = Md_f - Md_i #negative mass loss rate
d_Ma = Ma_f - Ma_i #positive mass accretion rate
Mt_f = Md_f + Ma_f
Mt_i = Md_i + Ma_i
if d_Md < 0|units.MSun and d_Ma >= 0|units.MSun:
eta = d_Ma / d_Md
a_f = a_i * ((Md_f/Md_i)**eta * (Ma_f/Ma_i))**-2 * Mt_i/Mt_f
return a_f
return a_i
def adjust_system_after_ce_in_inner_binary(bs, self):
# Assumption: Unstable mass transfer (common-envelope phase) in the inner binary, affects the outer binary as a wind.
# Instanteneous effect
if REPORT_FUNCTION_NAMES:
print('Adjust system after ce in inner binary')
if self.is_triple():
M_com_after_ce = self.get_mass(bs)
M_com_before_ce = bs.previous_mass
if self.triple.child1.is_star:
tertiary_star = self.triple.child1
else:
tertiary_star = self.triple.child2
# accretion_efficiency
M_accretor_before_ce = tertiary_star.mass
M_accretor_after_ce = tertiary_star.mass
a_new = adiabatic_expansion_due_to_mass_loss(self.triple.semimajor_axis, M_com_after_ce, M_com_before_ce, M_accretor_after_ce, M_accretor_before_ce)
self.triple.semimajor_axis = a_new
# print('outer orbit', a_new)
# nice but difficult to update self.triple
# system = bs
# while True:
# try:
# system = system.parent
# if not system.child1.is_star and system.child2.is_star:
# system = adjust_triple_after_ce_in_inner_binary(system, system.child1, system.child2, self)
# elif not system.child2.is_star and system.child1.is_star:
# system = adjust_triple_after_ce_in_inner_binary(system, system.child2, system.child1, self)
# else:
# print('adjust_system_after_ce_in_inner_binary: type of system unknown')
# exit(2)
#
# except AttributeError:
# #when there is no parent
# break
#
def stable_mass_transfer(bs, donor, accretor, self):
# orbital evolution is being taken into account in secular_code
if REPORT_FUNCTION_NAMES:
print('Stable mass transfer')
if bs.bin_type != bin_type['stable_mass_transfer']:
bs.bin_type = bin_type['stable_mass_transfer']
self.save_snapshot()
else:
bs.bin_type = bin_type['stable_mass_transfer']
self.secular_code.parameters.check_for_inner_RLOF = False
self.secular_code.parameters.include_spin_radius_mass_coupling_terms_star1 = False
self.secular_code.parameters.include_spin_radius_mass_coupling_terms_star2 = False
Md = donor.mass
Ma = accretor.mass
dt = self.triple.time - self.previous_time
dm_desired = bs.mass_transfer_rate * dt
if REPORT_FUNCTION_NAMES:
print(bs.mass_transfer_rate, dt, dm_desired)
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
donor_in_stellar_code.change_mass(dm_desired+(small_numerical_error|units.MSun), dt)
# dm != dm_desired e.g. when the envelope of the star becomes empty
dm = donor_in_stellar_code.mass - Md
bs.part_dt_mt = 1.
if dm - dm_desired > numerical_error|units.MSun:
# print('WARNING:the envelope is empty, mass transfer rate should be lower or dt should be smaller... ')
bs.part_dt_mt = dm/dm_desired
# there is an implicit assumption in change_mass that the accreted mass is of solar composition (hydrogen)
accretor_in_stellar_code = accretor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
# accretor_in_stellar_code.change_mass(dm, dt)
# for now, only conservative mass transfer
accretor_in_stellar_code.change_mass(-1.*dm, -1.*dt)
#if you want seba to determine the accretion efficiency, use
#accretor_in_stellar_code.change_mass(-1.*dm, dt)
#note doesnt work perfectly, as seba is oblivious to the roche lobe radius
#to adjust radius to mass loss and increase
self.stellar_code.evolve_model(self.triple.time)
self.copy_from_stellar()
self.update_stellar_parameters()
Md_new = donor.mass
Ma_new = accretor.mass
accretion_efficiency = (Ma_new-Ma)/(Md-Md_new)
if abs(accretion_efficiency - 1.0) > numerical_error and abs(Md-Md_new - -1.*(Ma-Ma_new)) > numerical_error |units.MSun:
self.save_snapshot()
print('stable_mass_transfer: non conservative mass transfer')
print(Md, Ma, donor.previous_mass, accretor.previous_mass)
print(Md_new, Ma_new, Md-Md_new, Ma-Ma_new, accretion_efficiency)
print(donor.stellar_type, accretor.stellar_type)
sys.exit('error in stable mass transfer')
bs.accretion_efficiency_mass_transfer = accretion_efficiency
corotation_spin = corotating_spin_angular_frequency_binary(bs.semimajor_axis, donor.mass, accretor.mass)
donor.spin_angular_frequency = corotation_spin
accretor.spin_angular_frequency = corotation_spin
def semi_detached(bs, donor, accretor, self):
#only for binaries (consisting of two stars)
if REPORT_FUNCTION_NAMES:
print('Semi-detached')
print(bs.semimajor_axis, donor.mass, accretor.mass, donor.stellar_type, accretor.stellar_type, bs.is_mt_stable)
stopping_condition = True
if bs.is_mt_stable:
stable_mass_transfer(bs, donor, accretor, self)
#adjusting triple is done in secular evolution code
else:
stopping_condition = common_envelope_phase(bs, donor, accretor, self)
return stopping_condition
#possible problem if companion or tertiary accretes significantly from this
# self.update_previous_stellar_parameters() #previous_mass, previous_radius for safety check
#-------------------------
#functions for contact mass transfer in a multiple / triple
#change parameters assuming fully conservative mass transfer
def perform_mass_equalisation_for_contact(bs, donor, accretor, self):
if REPORT_FUNCTION_NAMES:
print('perform_mass_equalisation_for_contact')
if REPORT_BINARY_EVOLUTION:
print('Start of stable mass transfer of contact systems' )
if donor.mass != accretor.mass:
# if abs(donor.mass - accretor.mass)> 1e-4|units.MSun: #could be better. see if problems arise
new_mass = 0.5*(donor.mass + accretor.mass)
bs.semimajor_axis = bs.semimajor_axis * (donor.mass * accretor.mass / new_mass ** 2) ** 2
delta_mass_donor = new_mass - donor.mass
delta_mass_accretor = new_mass - accretor.mass
donor_in_stellar_code = donor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
accretor_in_stellar_code = accretor.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
donor_in_stellar_code.change_mass(delta_mass_donor, -1.0|units.yr)
accretor_in_stellar_code.change_mass(delta_mass_accretor, -1.0|units.yr)
self.stellar_code.evolve_model(minimum_time_step) #to get updates radii, not just inflation of stars due to accretion
self.copy_from_stellar()
self.update_stellar_parameters() #makes secular code use update values: mass equalisation happens instantaneously
#set to synchronization
if self.include_CHE:
corotating_frequency = corotating_spin_angular_frequency_binary(bs.semimajor_axis, donor.mass, accretor.mass)
donor.spin_angular_frequency = corotating_frequency
accretor.spin_angular_frequency = corotating_frequency
donor.rotation_period = (2*np.pi/donor.spin_angular_frequency)
accretor.rotation_period = (2*np.pi/accretor.spin_angular_frequency)
self.channel_to_stellar.copy_attributes(['rotation_period']) #only defined when include_CHE
self.secular_code.parameters.include_inner_RLOF_terms = False
self.secular_code.parameters.include_outer_RLOF_terms = False
def contact_system(bs, star1, star2, self):
#if this implementation changes, then also change the is_mt_stable
if REPORT_FUNCTION_NAMES:
print("Contact system")
bs.bin_type = bin_type['contact']
self.save_snapshot()
self.check_RLOF() #@andris: is this necessary?
#for now no W Ursae Majoris evolution
#so for now contact binaries merge in common_envelope_phase, and MS-MS contact binaries will have mass equalisation
if bs.is_mt_stable: # happens when star1 & star2 are both on MS
perform_mass_equalisation_for_contact(bs, bs.child1, bs.child2, self)
stopping_condition = True
else:
if star1.mass >= star2.mass:
stopping_condition = common_envelope_phase(bs, star1, star2, self)
else:
stopping_condition = common_envelope_phase(bs, star2, star1, self)
return stopping_condition
#-------------------------
#functions for mass transfer in a multiple / triple
def triple_stable_mass_transfer(bs, donor, accretor, self):
# mass transfer of both inner and outer orbit is not yet considered here
# orbital evolution is being taken into account in secular_code
if REPORT_FUNCTION_NAMES:
print('Triple stable mass transfer')
if bs.bin_type != bin_type['stable_mass_transfer']:
bs.bin_type = bin_type['stable_mass_transfer']
self.save_snapshot()
else:
bs.bin_type = bin_type['stable_mass_transfer']
#implementation is missing
def triple_common_envelope_phase(bs, donor, accretor, self):
# mass transfer of both inner and outer orbit is not yet considered here
# orbital evolution is being taken into account in secular_code
if REPORT_FUNCTION_NAMES:
print('Triple common envelope')
bs.bin_type = bin_type['common_envelope']
self.save_snapshot()
#implementation is missing
#when the tertiary star transfers mass to the inner binary
def outer_mass_transfer(bs, donor, accretor, self):
#only for stellar systems consisting of a star and a binary
if REPORT_FUNCTION_NAMES:
print('Triple mass transfer')
bs.semimajor_axis, donor.mass, self.get_mass(accretor), donor.stellar_type
if bs.is_mt_stable:
triple_stable_mass_transfer(bs, donor, accretor, self)
# possible the outer binary needs part_dt_mt as well.
#adjusting triple is done in secular evolution code
else:
triple_common_envelope_phase(bs, donor, accretor, self)
#stopping condition 0:False, 1:True, -1: calculate through outer mass transfer - effect on inner & outer orbit is taken care off here.
return -1
#-------------------------
#-------------------------
#Functions for detached evolution
## Calculates stellar wind velocoty.
## Steller wind velocity is 2.5 times stellar escape velocity
#def wind_velocity(star):
# v_esc2 = constants.G * star.mass / star.radius
# return 2.5*np.sqrt(v_esc2)
#}
#
#
## Bondi, H., and Hoyle, F., 1944, MNRAS 104, 273 (wind accretion.
## Livio, M., Warner, B., 1984, The Observatory 104, 152.
#def accretion_efficiency_from_stellar_wind(accretor, donor):
#velocity needs to be determined -> velocity average?
# why is BH dependent on ecc as 1/np.sqrt(1-e**2)
# alpha_wind = 0.5
# v_wind = wind_velocity(donor)
# acc_radius = (constants.G*accretor.mass)**2/v_wind**4
#
# wind_acc = alpha_wind/np.sqrt(1-bs.eccentricity**2) / bs.semimajor_axis**2
# v_factor = 1/((1+(velocity/v_wind)**2)**3./2.)
# mass_fraction = acc_radius*wind_acc*v_factor
#
# print('mass_fraction:', mass_fraction)
## mass_fraction = min(0.9, mass_fraction)
#
def detached(bs, self):
# orbital evolution is being taken into account in secular_code
if REPORT_FUNCTION_NAMES:
print('Detached')
if bs.bin_type == bin_type['detached'] or bs.bin_type == bin_type['unknown']:
bs.bin_type = bin_type['detached']
else:
bs.bin_type = bin_type['detached']
self.save_snapshot()
# wind mass loss is done by stellar_code
# wind accretion here:
# update accretion efficiency of wind mass loss
if self.is_binary(bs):
bs.accretion_efficiency_wind_child1_to_child2 = 0.0
bs.accretion_efficiency_wind_child2_to_child1 = 0.0
# child1_in_stellar_code = bs.child1.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
# child2_in_stellar_code = bs.child2.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
#
# dt = self.triple.time - self.previous_time
# dm_child1_to_child2 = -1 * child1.wind_mass_loss_rate * bs.accretion_efficiency_wind_child1_to_child2 * dt
# child2_in_stellar_code.change_mass(dm_child1_to_child2, -1*dt)
# dm_child12to_child1 = -1 * child2.wind_mass_loss_rate * bs.accretion_efficiency_wind_child2_to_child1 * dt
# child1_in_stellar_code.change_mass(dm_child2_to_child1, -1*dt)
# check if this indeed is accreted conservatively
elif bs.child1.is_star and self.is_binary(bs.child2):
#Assumption: an inner binary is not effected by wind from an outer star
bs.accretion_efficiency_wind_child1_to_child2 = 0.0
bs.accretion_efficiency_wind_child2_to_child1 = 0.0
# child1_in_stellar_code = bs.child1.as_set().get_intersecting_subset_in(self.stellar_code.particles)[0]
# dt = self.triple.time - self.previous_time
#effect of wind from bs.child2.child1 onto bs.child1
# mtr_w_in1_1 = bs.child2.child1.wind_mass_loss_rate * (1-bs.child2.accretion_efficiency_wind_child1_to_child2)
# beta_w_in1_1 = 0.0
# dm_in1_1 = -1 * mtr_w_in1_1 * beta_w_in1_1 * dt
#
#effect of wind from bs.child2.child2 onto bs.child1
# mtr_w_in2_1 = bs.child2.child2.wind_mass_loss_rate * (1-bs.child2.accretion_efficiency_wind_child2_to_child1)
# beta_w_in2_1 = 0.0
# dm_in2_1 = -1 * mtr_w_in2_1 * beta_w_in2_1 * dt
#
# dm = dm_in1_1 + dm_in2_1
# mtr = mtr_w_in1_1 + mtr_w_in2_1)
#effect of mass transfer in the binary bs.child2 onto bs.child1
# if bs.child2.child1.is_donor and bs.child2.child2.is_donor:
# print('contact binary in detached...')
# exit(1)
# elif bs.child2.child1.is_donor or bs.child2.child2.is_donor:
# #Assumption:
# #Stable mass transfer in the inner binary, affects the outer binary as a wind.
# mtr_rlof_in_1 = bs.child2.mass_transfer_rate * (1-bs.child2.accretion_efficiency_mass_transfer)
# beta_rlof_in_1 = 0.0
# dm_rlof_in_1 = -1 * mtr_rlof_in_1 * beta_rlof_in_1 * dt
# dm += dm_rlof_in_1