وحدات تطور النجوم الثنائية في REBOUNDx

الملخص

يقترن تطور الثنائيات المتقاربة بطفح فص روش (RLOF)، وبمقاومة الغلاف المشترك (CE)، وبالرياح النجمية، والكبح المغناطيسي، وفواقد الموجات الثقالية، بحيث تتبادل هذه العمليات الكتلة والزخم الزاوي بينما تعيد تشكيل المدارات والدورانات المغزلية. نقدّم مؤثرات قابلة للتشغيل المتبادل في امتداد REBOUNDx إلى REBOUND تدمج هذه العمليات ضمن ديناميكيات -جسمية عالية الدقة. وتشمل الحزمة: مؤثر RLOF حافظاً للزخم، بقنوات محافظة وجهازية وبفقد نوعي قابل للتهيئة لـ ؛ ونموذج مقاومة CE قائم على الاحتكاك الديناميكي المعتمد على عدد ماخ مع حدّ للدفعات؛ ورياح Reimers متناحية الخواص، ورياحاً حرارية من نمط Parker، وتدفقات محدودة بإدنغتون مشغّلة بوحدة تطور نجمي بارامترية تزوّد و معتمدين على الكتلة؛ وكبحاً مغناطيسياً عبر عزم Verbunt–Zwaan/Kawaler مع تحديث مغلق الشكل للدوران يأخذ التشبع في الحسبان؛ وتصحيحات ما بعد نيوتنية (كتل نقطية وازدواج دوران–دوران عند 2 PN؛ ورد فعل إشعاعي عند 2.5 PN). يُحفظ الزخم الخطي في النقل المحافظ، ويفرض عزم تصحيحي أصغري اتساق الزخم الزاوي، كما تثبّت الخطوات الفرعية التكيفية التطور قرب التلامس. وتنسّق رايات بين-وحدية نشاط الرياح وRLOF وCE. يتيح الإطار المستقل عن الوحدات دراسات متسقة ذاتياً ومحلولة زمنياً للثنائيات المتقاربة في أوساط معزولة أو غنية ديناميكياً. تَرِد في الملحق أمثلة متعددة ومقارنات مع شفرات أخرى. الشفرة متاحة على https://github.com/malidib/ReboundS .

1 مقدمة

تتطور النجوم الثنائية المتقاربة تحت مجموعة مترابطة بشدة من العمليات الفيزيائية التي تعيد توزيع الكتلة والزخم الزاوي، وفي بعض الأنظمة تزيل طاقة مدارية. وينطلق تبادل الكتلة عبر طفح فص روش (RLOF) عندما يملأ النجم المانح سطح روش المكافئ حجماً، والذي يُمثّل عادة بصيغة Eggleton (Eggleton, 1983)، في حين يمكن التقاط استجابة الفقد الكتلي اللحظية قرب نقطة Lagrange الداخلية بوصفات أسية معايرة على ارتفاع السلم الفوتوسفيري (Ritter, 1988). وفي أطوار أشد تطرفاً، قد يقود انتقال الكتلة غير المستقر وانكماش المدار إلى تطور الغلاف المشترك (CE)، حيث يستخلص السحب الهيدروديناميكي داخل غلاف ممتد الطاقة المدارية والزخم الزاوي (Ostriker, 1999). وتعمل قنوات مستقلة لفقد الكتلة والزخم الزاوي طوال حياة الثنائية: فالنجوم الباردة ذات الحمل الحراري تتعرض لكبح مغناطيسي عبر رياح ممغنطة (Verbunt & Zwaan, 1981; Kawaler, 1988)، وتفقد نجوم فرع العمالقة كتلتها عبر رياح نجمية مترابطة مع المعلمات الكلية (Reimers, 1975)، وتصدر الثنائيات المدمجة إشعاعاً ثقالياً يقود هبوطاً حلزونياً علمانياً على مقاييس زمنية نسبية (Peters, 1964). وتعدّل تغيرات البنية النجمية هذه المسارات عبر تغيير أنصاف الأقطار واللمعانيات وعزوم العطالة الداخلية (Hurley et al., 2002). وبمجموعها، تحدد هذه المؤثرات ما إذا كانت الثنائية المتقاربة ستخضع لتبادل كتلي مستقر، أو ستنكمش إلى حالة تماس وتندمج، أو ستنجو لتنتج بقايا مدمجة.

يبقى تمثيل هذا المشهد متعدد الفيزياء بطريقة واحدة أمراً صعباً. فالحسابات الهيدروديناميكية ثلاثية الأبعاد بالكامل تستطيع حل مجاري RLOF وتدفقات CE، لكنها مرتفعة الكلفة جداً بالنسبة للتطور طويل الأمد أو دراسات المعلمات على مستوى الجماعات. وفي المقابل، ترمّز صيغ التطور الثنائي السريعة RLOF والرياح والكبح المغناطيسي وفواقد الموجات الثقالية عبر وصفات علمانية فعالة وتنبؤية عند العزل، لكنها تفترض عموماً سياق جسمين ولا تستطيع بطبيعتها احتساب الديناميكيات الثقالية الأعلى رتبة (مثل الثلاثيات، أو اللقاءات الرنينية، أو كمونات العناقيد) التي يمكن أن تعدّل انتقال الكتلة أو حتى تطلقه. ويتمثل نهج مكمل في دمج العمليات العلمانية والتبددية الموثقة ضمن مكامل -جسمي عالي الدقة بحيث تتقدم الديناميكيات المدارية وتطور الدوران وتبادل الكتلة بصورة متسقة ذاتياً على خطوة زمنية مشتركة، مع الحفاظ على مسك الحسابات الذي تفرضه قوانين الحفظ.

تعرض هذه الورقة حزمة من وحدات تطور ثنائي جديدة نُفذت في امتداد REBOUNDx لإطار REBOUND -الجسمي (Rein & Liu, 2012; Rein & Tamayo, 2015; Tamayo et al., 2020). تستهدف الوحدات القنوات المهيمنة ذات الصلة بالثنائيات المتقاربة: (i) وحدة RLOF وCE حافظة للزخم تجمع النقل المحافظ مع فقد كتلي جهازي غير محافظ وتدعم وصفات متعددة للزخم الزاوي النوعي للغاز الهارب؛ (ii) مؤثر كبح مغناطيسي ينفذ عزم Verbunt–Zwaan/Kawaler مع فرع تشبع ملائم للأغلفة الحملية سريعة الدوران (Verbunt & Zwaan, 1981; Kawaler, 1988); (iii) فقد كتلي برياح نجمية متناحية الخواص وفق قياس Reimers (Reimers, 1975); (iv) رياحاً مدفوعة حرارياً (شبيهة بـ Parker) ممثلة بكفاءة تسخين؛ و(v) تصحيحات ما بعد نيوتنية (PN) تقدّم حدود كتل نقطية محافظة ودوران–دوران عند 2 PN، ورد فعل إشعاعي للموجات الثقالية عند 2.5 PN (Peters, 1964; Kidder, 1995)؛ وتتوفر مؤثرات 1 PN ودوران–مدار (1.5 PN) عبر وحدتي gr_full وlense_thirring. صُمم كل مؤثر ليكون مستقلاً عن الوحدات، وليكشف معلمات دنيا ذات تفسير فيزيائي يمكن معايرتها أو تغييرها منهجياً.

هدفنا هو إتاحة هذه العمليات وتوثيقها وجعلها قابلة لإعادة الإنتاج ضمن بيئة -جسمية واحدة، بحيث تستطيع دراسات تطور الثنائيات المتقاربة الانتقال بسلاسة بين الثنائيات المعزولة والسياقات الغنية ديناميكياً. وتفصّل بقية هذه الورقة النماذج الفيزيائية والتنفيذ العددي، وتعرض تطبيقات نموذجية.

2 تنفيذ المؤثرات

نعرض فيما يلي التطبيقات الفيزيائية والعددية لوحدات تطور الثنائيات الجديدة. ونلاحظ أن الملحق يتضمن أمثلة متعددة ومقارنات مع شفرات أخرى.

2.1 تطور نجمي مبسط

تعتمد بعض المؤثرات المنفذة أدناه على امتلاك أنصاف أقطار ولمعانيات نجمية متسقة ذاتياً بوصفها دوالاً في الكتلة. ولتوفير هذه الخواص، يحدّث المؤثر الاختياري stellar_evolution_sse (الجدول 1) نصف قطر كل نجم ولمعانيته من وصفات تحليلية معتمدة على الكتلة (Hurley et al., 2000). يكون التحديث جبرياً خالصاً: تُهمل الخطوة الزمنية وتُقيّم العلاقات باستخدام الكتلة الحالية للجسيم عند كل استدعاء. وتُتجاوز الجسيمات الافتراضية، وتبقى الكتل النجمية بلا تغيير.

تُبنى الكتلة النجمية عديمة الأبعاد باستخدام معلمة المؤثر sse_Msun (القيمة الافتراضية 1) لتحويل كتل الشفرة إلى كتل شمسية. وتعرّف sse_Rsun وsse_Lsun، على نحو مماثل، نصف القطر الشمسي واللمعانية الشمسية بوحدات الشفرة. يضبط حقل نصف قطر الجسيم sim.particles[index].r إلى ، وتُكتب اللمعانية في خاصية الجسيم sse_L. ويمكن أن تستخدم وحدات أخرى هذه اللمعانية المخزنة، مثل مؤثرات الرياح.

بديل أساسي بقانون قوة.

إذا وفر جسيم أيّاً من معلمات قانون القوة الأساسية sse_R_coeff أو sse_R_exp أو sse_L_coeff أو sse_L_exp، فإن المؤثر يستخدم تحويل قانون القوة

مع القيم الافتراضية و و و إذا لم تُحدَّد أي من المعلمات المناظرة على مستوى الجسيم.

مسبقات فئات الأجسام عبر sse_type.

إذا لم توجد معلمات أساسية لقانون القوة، يستخدم المؤثر العدد الصحيح لكل جسيم sse_type لاختيار واحد من خمسة مسبقات تحليلية للبنية. ويعيد مضاعفان اختياريان لكل جسيم، sse_R_mult وsse_L_mult (قيمتهما الافتراضية 1)، تحجيم قيم المسبق (والبديل) كما يلي و . والمسبقات هي:

• sse_type=1 (تسلسل رئيسي غني بالهيدروجين / شبيه بـ ZAMS). نستخدم تحجيمات مستمرة مقطعية لكل من و:

• sse_type=2 (عملاق غني بالهيدروجين؛ يجمع RGB/AGB). نعتمد تحجيماً واسع نصف القطر وضعيف الاعتماد على الكتلة مع إغلاق يقابل درجة حرارة فعالة باردة نموذجية:

  

  وتسمح المضاعفات sse_R_mult وsse_L_mult للمستخدمين بإزاحة مقياس العملاق المرجعي من غير توفير مسارات تطورية كاملة.

• sse_type=3 (نجم هيليوم منزوع الغلاف؛ شبيه بـ He-MS / WR). نستخدم تحجيمات مدمجة ومضيئة بقانون قوة:

  

• sse_type=4 (قزم أبيض). يتبع نصف القطر علاقة كتلة–نصف قطر مطبّعة بكتلة Chandrasekhar،

  

  مع (القيمة الافتراضية 1.44) و (القيمة الافتراضية 0.0112). وتُضبط اللمعانية إلى (القيمة الافتراضية ) ما لم يُعد تحجيمها بواسطة sse_L_mult.

• sse_type=5 (بقية مدمجة؛ NS/BH). يسند المؤثر نصف قطر ثابتاً لنجم نيوتروني إذا كان $m\le \text{sse\_Mns\_max}$ (القيمة الافتراضية 3)، وإلا فإنه يستخدم تحجيماً بنصف قطر Schwarzschild لثقب أسود:

  القيم الافتراضية هي sse_Rns_factor و sse_Rbh_factor.

إذا غاب sse_type أو ساوى 0 ولم توجد معلمات أساسية لقانون القوة، يعود المؤثر إلى قوانين القوة الافتراضية و.

Table 1: معلمات البنية النجمية (stellar_evolution_sse)

Name (scope)	Unit	Default	Purpose
sse_Msun (op)	mass	1	Solar mass in code units ($M_{\odot }$)
sse_Rsun (op)	length	1	Solar radius in code units ($R_{\odot }$)
sse_Lsun (op)	luminosity	1	Solar luminosity in code units ($L_{\odot }$)
sse_type (part)	—	0	Object class selector (0–5; see text)
sse_R_mult (part)	—	1	Radius multiplier for presets
sse_L_mult (part)	—	1	Luminosity multiplier for presets
sse_R_coeff (part)	—	1	basic radius prefactor (overrides sse_type)
sse_R_exp (part)	—	0.8	basic radius mass exponent
sse_L_coeff (part)	—	1	basic luminosity prefactor
sse_L_exp (part)	—	3.5	basic luminosity mass exponent
sse_Mch (op)	—	1.44	Chandrasekhar mass in $M_{\odot }$ (WD preset)
sse_Rwd_coeff (op)	—	0.0112	WD radius coefficient in $R_{\odot }$
sse_Lwd (op)	—	${10}^{-3}$	WD luminosity in $L_{\odot }$
sse_Rns_factor (op)	—	$1.724\times {10}^{-5}$	NS radius in $R_{\odot }$
sse_Rbh_factor (op)	—	$4.246\times {10}^{-6}$	Schwarzschild radius per $M_{\odot }$ in $R_{\odot }$
sse_Mns_max (op)	—	3.0	NS/BH boundary mass in $M_{\odot }$ (compact preset)
sse_L (part, out)	luminosity	—	Stored luminosity used by other operators

2.2 فقد الكتلة بالرياح النجمية وفق Reimers

في وصفة Reimers (Reimers, 1975)، تزيل الرياح متناحية الخواص كتلة من عمالقة منفردة باردة ومنخفضة إلى متوسطة الكتلة. وبالنسبة إلى نجم كتلته ، ولمعانيته ، ونصف قطره ، يكون معدل فقد الكتلة

$$\dot{M}=-4\times {10}^{-13}\eta \frac{L}{L_{\odot }}\frac{R}{R_{\odot }}\frac{M_{\odot }}{M}{M}_{\odot }{\mathrm{yr}}^{-1},$$

حيث تُحدد الكفاءة عديمة الأبعاد واللمعانية (sse_L المعرّفة أعلاه) لكل جسيم، في حين تؤخذ من حقل نصف قطر الجسيم. تُتجاوز النجوم التي تفتقر إما إلى swml_eta أو إلى sse_L. تُزال الكتلة على نحو متناحي الخواص من دون ارتداد في الزخم الخطي؛ وتُتجاهل الجسيمات الافتراضية، وبعد فقد الكتلة يُعاد تمركز النظام على مركز الكتلة. ويحد محدد أمان التغير الكسري في الكتلة عند swml_max_dlnM (القيمة الافتراضية ) لكل استدعاء. ويمكن تعديل العامل الأمامي والقيم المرجعية الشمسية وطول السنة عبر معلمات المؤثر (الجدول 2). افتراضياً، يثبط المؤثر الرياح عندما يوسم جسيم بأنه داخل غلاف مشترك أو يخضع لطفح فص روش؛ ويتحكم المفتاحان swml_disable_in_CE وswml_disable_in_RLOF في هذا السلوك عبر قراءة الرايتين لكل جسيم inside_CE وrlof_active.

Table 2: معلمات فقد الكتلة بالرياح النجمية

Name (scope)	Unit	Default	Purpose
swml_eta (part)	—	—	Wind efficiency $\eta$
swml_const (op)	$M_{\odot }$/yr	$4\times {10}^{-13}$	Reimers prefactor
swml_Msun (op)	mass	1	Solar mass in code units
swml_Rsun (op)	length	1	Solar radius in code units
swml_Lsun (op)	luminosity	1	Solar luminosity in code units
swml_year (op)	time	1	Length of Julian year in code units
swml_max_dlnM (op)	—	0.1	Max $\left|\mathrm{\Delta }M\right|/M$ per call
swml_disable_in_CE (op)	bool	1	Skip if inside_CE$>0$
swml_disable_in_RLOF (op)	bool	1	Skip if rlof_active$>0$

2.3 رياح حرارية من نمط Parker

في النجوم الباردة منخفضة الثقالة، يمكن للضغط الحراري أن يطلق رياحاً متناحية الخواص تحمل الكتلة بعيداً (Parker, 1958). وبالنسبة إلى نجم كتلته ، ولمعانيته ، ونصف قطره ، نعتمد تحجيماً شبيهاً بـ Parker

$$\dot{M}=-C_{\mathrm{th}}\eta {\left(\frac{R}{R_{\odot }}\right)}^{{\alpha }_R}{\left(\frac{L}{L_{\odot }}\right)}^{{\alpha }_L}{\left(\frac{M_{\odot }}{M}\right)}^{{\alpha }_M}{M}_{\odot }{\mathrm{yr}}^{-1},$$

حيث كفاءة تسخين عديمة الأبعاد، وللأسس القيم الافتراضية . تُقرأ اللمعانية من خاصية الجسيم sse_L، التي يجب أن يملأها إما مؤثر التطور النجمي المبسط أو المستخدم يدوياً، في حين تؤخذ و من حقلي و للجسيم. تُتجاوز النجوم التي تفتقر إما إلى tdw_eta أو إلى sse_L. تُزال الكتلة على نحو متناحي الخواص من دون ارتداد في الزخم الخطي؛ وتُتجاهل الجسيمات الافتراضية، وبعد فقد الكتلة يُعاد تمركز النظام على مركز الكتلة. ويمكن تعديل العامل الأمامي والقيم المرجعية الشمسية والأسس وأكبر تغير كسري في الكتلة وطول السنة عبر معلمات المؤثر (الجدول 3). يكون المؤثر مستقلاً عن الوحدات إذا حُددت ثوابت التحجيم هذه بصورة متسقة. افتراضياً، يوقف المؤثر الرياح للنجوم الموسومة بأنها داخل غلاف مشترك أو بأنها تخضع لطفح فص روش نشط، كما تدل الرايتان لكل جسيم inside_CE وrlof_active. ويتحكم في هذا السلوك tdw_disable_in_CE و tdw_disable_in_RLOF.

مجال الانطباق: درجة الحرارة والجاذبية السطحية.

تهدف وصفة نمط Parker إلى تمثيل الرياح الإكليلية/الكروموسفيرية المدفوعة حرارياً. بالنسبة إلى رياح Parker متساوية الحرارة بدرجة حرارة رياح ، يكون نصف القطر الصوتي

(1)

حيث هو الوزن الجزيئي المتوسط و هي كتلة البروتون. ويتطلب حل عابر للصوت منطلق قرب سطح النجم ، أي

(2)

مع كون الجاذبية السطحية. عملياً، تقابل الرياح الشبيهة بـ Parker المجال –، أو على نحو مكافئ –. ويناظر ذلك – للأقزام الشبيهة بالشمس (–) و– للعمالقة (–). ولا يستهدف المؤثر الأجسام المدمجة (ذات عالٍ جداً) أو الأنظمة التي تكون الرياح فيها مدفوعة أساساً بالخطوط الطيفية أو بالغبار/ضغط الإشعاع.

Table 3: معلمات الرياح المدفوعة حرارياً

Name (scope)	Unit	Default	Purpose
tdw_eta (part)	—	—	Wind efficiency $\eta$
sse_L (part)	luminosity	—	Stellar luminosity $L$
tdw_const (op)	$M_{\odot }$/yr	$2\times {10}^{-14}$	Thermal-wind prefactor $C_{\mathrm{th}}$
tdw_Msun (op)	mass	1	Solar mass in code units
tdw_Rsun (op)	length	1	Solar radius in code units
tdw_Lsun (op)	luminosity	1	Reference luminosity $L_{\odot }$ in code units
tdw_year (op)	time	1	Length of Julian year in code units
tdw_alpha_R (op)	—	2	Exponent ${\alpha }_R$ on $R/R_{\odot }$
tdw_alpha_L (op)	—	$3/2$	Exponent ${\alpha }_L$ on $L/L_{\odot }$
tdw_alpha_M (op)	—	1	Exponent ${\alpha }_M$ on $M_{\odot }/M$
tdw_max_dlnM (op)	—	0.1	Max $\left|\mathrm{\Delta }M\right|/M$ per call
tdw_disable_in_CE (op)	bool	1	Skip if inside_CE$>0$
tdw_disable_in_RLOF (op)	bool	1	Skip if rlof_active$>0$

2.4 رياح محدودة بإدنغتون

تؤدي اللمعانيات التي تتجاوز حد إدنغتون لتشتت الإلكترونات في بعض الأجسام النجمية الضخمة وعالية اللمعانية إلى فقد كتلي متناحي الخواص (مثلاً Eddington, 1926; Shaviv, 2001) وفقاً لـ

$$\dot{M}=-C_{\mathrm{edd}}\max (0,\frac{L}{L_{\mathrm{Edd}}}-1)M_{\odot }{\mathrm{yr}}^{-1},$$

حيث . تُقرأ اللمعانية من خاصية الجسيم sse_L، التي يجب أن يوفرها المستخدم أو مؤثر التطور النجمي المبسط. تُزال الكتلة على نحو متناحي الخواص من دون ارتداد في الزخم الخطي؛ وتُتجاهل الجسيمات الافتراضية، وبعد فقد الكتلة يُعاد تمركز النظام على مركز الكتلة. ويحد محدد التغير الكسري في الكتلة عند edw_max_dlnM لكل استدعاء. وترد معلمات المؤثر التي تتحكم بثوابت التحجيم في الجدول 4. افتراضياً، يتجاوز المؤثر النجوم الموسومة بأنها داخل غلاف مشترك أو تخضع لطفح فص روش، ويتحكم في ذلك edw_disable_in_CE وedw_disable_in_RLOF، اللذان يقرآن inside_CE وrlof_active.

Table 4: معلمات الرياح فوق إدنغتون

Name (scope)	Unit	Default	Purpose
sse_L (part)	luminosity	—	Stellar luminosity $L$
edw_const (op)	$M_{\odot }$/yr	$1\times {10}^{-6}$	Mass-loss prefactor $C_{\mathrm{edd}}$
edw_Msun (op)	mass	1	Solar mass in code units
edw_Lsun (op)	luminosity	1	Solar luminosity in code units
edw_year (op)	time	1	Length of Julian year in code units
edw_Ledd_coeff (op)	$L_{\odot }/M_{\odot }$	$3.2\times {10}^4$	$L_{\mathrm{Edd}}$ per unit mass
edw_max_dlnM (op)	—	0.1	Max $\left|\mathrm{\Delta }M\right|/M$ per call
edw_disable_in_CE (op)	bool	1	Skip if inside_CE$>0$
edw_disable_in_RLOF (op)	bool	1	Skip if rlof_active$>0$

2.5 الكبح المغناطيسي

تفقد النجوم منخفضة الكتلة ذات الأغلفة الحملية زخماً زاوياً مغزلياً عبر رياح ممغنطة (مثلاً Skumanich, 1972; Rappaport et al., 1983). وفي الثنائيات المقترنة مدّياً (ولا سيما عندما يكون أحد المكوّنين قريباً من الدوران المتزامن)، تعوّض عزوم المد الدوران المكبوح بسحبها من الخزان المداري، مما يؤدي إلى فقد صافٍ في الزخم الزاوي المداري ونقصان علماني في نصف المحور الأكبر والدور المداري؛ أما في الأنظمة الواسعة أو الضعيفة الاقتران فيكون الرد المداري مهملاً على مقاييس زمنية مماثلة.

ننفذ قانون عزم Verbunt–Zwaan / Kawaler (Verbunt & Zwaan, 1981; Kawaler, 1988)، مطبقين عزم كبح معاكساً لاتجاه متجه الدوران (المعلمات في الجدول 5).

قانون العزم.

بالنسبة إلى نجم كتلته ، ونصف قطره ، وسرعته الزاوية $\omega =\parallel 𝛀\parallel$، يكون عزم إبطاء الدوران

(3)

حيث ثابت تطبيع (معلمة المؤثر mb_K، وقيمته الافتراضية في cgs). يعمل العزم على استقامة واحدة مع $𝛀$، ولذلك لا يغير هذا المؤثر اتجاه الدوران.

الوحدات والتحجيم.

يحدد المستخدمون mb_K في cgs. وداخلياً نحوله إلى وحدات الشفرة باستخدام معلمات المؤثر mb_Msun وmb_Rsun و mb_year، التي تعرّف و والسنة اليوليانية بوحدات الشفرة. وبتعريف الوحدات الأساسية cgs لكل وحدة شفرة ، ، ، ومع ملاحظة أن ، نحصل على

(4)

ثم نقيّم المعادلة (3) بالعوامل الصريحة عديمة الأبعاد . ولا يحتاج المستخدمون إلى إعادة تحجيم mb_K يدوياً.

التشبع.

إذا وفر جسيم زمن دوران حملي mb_tau_conv، تضبط الشفرة السرعة الزاوية الحرجة عبر عدد Rossby على النحو

(5)

مع mb_Rossby_sat (على مستوى المؤثر؛ القيمة الافتراضية 0.1). وتتجاوز القيمة التي يوفرها الجسيم mb_omega_sat هذه القيمة. ولـ يتدرج العزم خطياً مع ، بما يضمن انتقالاً مستمراً عند العتبة.

التكامل الزمني (صيغة مغلقة).

لأن العزم على استقامة واحدة مع $𝛀$، فإن المقدار يحقق معادلة تفاضلية عادية قياسية ذات حل مغلق الشكل. عرّف

$$C=\frac{K_{\mathrm{code}}}{I}{\left(\frac{R}{R_{\odot }}\right)}^{1/2}{\left(\frac{M}{M_{\odot }}\right)}^{-1/2}.$$ (6)

إذن

unsaturated: ${\displaystyle \frac{d\omega }{dt}}=-C{\omega }^3\Rightarrow \omega (t+\mathrm{\Delta }t)={\displaystyle \frac{\omega (t)}{\sqrt{1+2C\omega {(t)}^2\mathrm{\Delta }t}}},$ (7)

saturated: ${\displaystyle \frac{d\omega }{dt}}=-C{\omega }_{\mathrm{sat}}^2\omega \Rightarrow \omega (t+\mathrm{\Delta }t)=\omega (t)\exp \left[-C{\omega }_{\mathrm{sat}}^2\mathrm{\Delta }t\right].$ (8)

إذا بدأت خطوة في نظام التشبع وعبرت إلى النظام غير المشبع، يطبّق التحديث مقطعياً: اضمحلال أسي إلى يتبعه تطبيق الصيغة غير المشبعة لبقية الخطوة. ثم يُعاد تحجيم متجه الدوران أخيراً بواسطة $𝛀\leftarrow 𝛀\left[\omega (t+\mathrm{\Delta }t)/\omega (t)\right]$. يحفظ هذا المخطط الإيجابية ويتجنب لااستقرارات حجم الخطوة.

التفعيل والضوابط.

يطبّق الكبح فقط إذا عيّن الجسيم mb_on و mb_convective. وتتخطى المعلمات المفقودة، أو القيم غير الموجبة أو غير المنتهية أو أو ، التحديث. ويدعم المؤثر تجاوزاً لنصف القطر لكل جسيم mb_R (بوحدات طول الشفرة)؛ وإلا فيُستخدم نصف قطر الجسيم p->r. ولا يفعل المؤثر شيئاً لمتجهات الدوران المتلاشية.

Table 5: معلمات الكبح المغناطيسي

Name (scope)	Unit	Default	Purpose
mb_K (op)	cgs	$2.7\times {10}^{47}$	Braking constant $K$
mb_Msun (op)	mass	1	Solar mass in code units
mb_Rsun (op)	length	1	Solar radius in code units
mb_year (op)	time	1	Julian year in code units
mb_Rossby_sat (op)	—	0.1	Critical Rossby number
mb_on (part)	bool	0	Enable magnetic braking
mb_convective (part)	bool	0	Convective‑envelope flag
mb_omega_sat (part)	1/t	$\mathrm{\infty }$	Saturation angular velocity
mb_tau_conv (part)	time	—	Convective turnover time
mb_R (part)	length	—	Radius override (else particle radius)
I (part)	mass length2	—	Moment of inertia
Omega (part, vec)	1/t	—	Spin angular‑velocity vector (reb_vec3d)

اختبارات معقولية.

عندما يكون mb_Msun = mb_Rsun = mb_year = 1، و، و بالراديان/سنة، يتدرج العزم غير المشبع كما مع التطبيع المتوقع من الأدبيات. والتحديث مستمر عند بحكم البناء.

2.6 تصحيحات ما بعد نيوتنية (PN)

تخضع الثنائيات المدمجة لتصحيحات نسبية تقرن الدورانات المغزلية وتشع طاقة مدارية (مثلاً Blanchet, 2014). وينفذ المؤثر post_newtonian (الجدول 6) معادلات حركة الكتل النقطية بالإحداثيات التوافقية من Kidder (1995) لكل زوج ذي كتلة:

• تصحيحات محافظة للكتل النقطية (PM) والدوران–الدوران (SS) عند  PN؛

• رد فعل إشعاعي للموجات الثقالية (RR) عند  PN.

المتغيرات والكتل.

بالنسبة إلى جسمين كتلتهما و، ومتجه الفصل بينهما ، وسرعتهما النسبية ، نعرّف

$$m=m_i+m_j,\mu =\frac{m_i{m}_j}{m},\eta =\frac{\mu }{m},𝐧=\frac{𝐫}{r},\dot{r}=𝐯\cdot 𝐧.$$

الدورانات هي عزوم زاوية فيزيائية (بوحدات كتلة طول²/زمن). وإذا فضّل المستخدمون الدورانات عديمة الأبعاد ، فليحوّلوها عبر في نظام وحدات المحاكاة.

محافظ عند 2 PN.

يُكتب جزء الكتل النقطية على النحو

$${𝐚}_{2\mathrm{P}\mathrm{N}}^{\mathrm{pm}}=-\frac{Gm}{c^4{r}^2}\left[A_2𝐧+B_{v,2}𝐯\right],$$

مع بناء و من و${\dot{r}}^2$ و؛ ويتضمن التنفيذ حد $-\frac{1}{2}\eta (13-4\eta )(Gm/r)v^2$ بحيث يطابق ما في Kidder (1995). أما اقتران الدوران–الدوران فهو

$${𝐚}_{2\mathrm{P}\mathrm{N}}^{\mathrm{ss}}=-\frac{3G}{\mu c^2{r}^4}\left[\left(S_1\cdot S_2-5(𝐧\cdot S_1)(𝐧\cdot S_2)\right)𝐧+(𝐧\cdot S_2)S_1+(𝐧\cdot S_1)S_2\right],$$

ويستخدم دورانات فيزيائية والكتلة المختزلة في العامل الأمامي الكلي. ومن ثم

رد فعل إشعاعي عند 2.5 PN.

$${𝐚}_{2.5\mathrm{PN}}=\frac{8G^2{m}^2\eta }{5c^5{r}^3}\left[\dot{r}\left(18v^2+\frac{2}{3}\frac{Gm}{r}-25{\dot{r}}^2\right)𝐧-\left(6v^2-2\frac{Gm}{r}-15{\dot{r}}^2\right)𝐯\right].$$

الوحدات والمعلمات.

معلمة المؤثر c هي سرعة الضوء بوحدات الشفرة. وتفعّل المفاتيح pn_2PN وpn_25PN (القيمة الافتراضية 1) القطاعات المناظرة. ومعلمة الجسيم pn_spin هي reb_vec3d تخزن متجه الدوران الفيزيائي.

Table 6: معلمات ما بعد نيوتنية

Name (scope)	Unit	Default	Purpose
c (force)	speed	—	Speed of light (required)
pn_2PN (force)	bool	1	Enable 2 PN point-mass + spin–spin
pn_25PN (force)	bool	1	Enable 2.5 PN radiation reaction
pn_merge_dist (force)	length	0	Pre-pass merger distance: if $>0$, merge any pair with $r\le \text{pn\_merge\_dist}$ before PN; if $=0$, merge only for exact coincidence.
pn_spin (part, vec)	ang. mom.	—	Physical spin vector $S$

حفظ الزخم والنطاق.

تُبنى التسارعات للإحداثي النسبي وتُقسم بين الجسمين بنسبة كتلتيهما، بما يحفظ الزخم الخطي الكلي. ولا تُدرج معادلات بدارية الدوران؛ فإذا كانت الدورانات غير مصطفة، تضم الديناميكيات المدارية قوى دوران–دوران مع دورانات ثابتة.

الأثر في مدارات الثنائيات.

يزيل الحد 2.5 PN الطاقة المدارية والزخم الزاوي، دافعاً نقصاناً علمانياً في نصف المحور الأكبر و(عادة) في الشذوذ باتجاه الاندماج. وعند وجود الدورانات، تسبب حدود الدوران–الدوران بدارية نسبية ويمكن أن تعدّل المستوى المداري اللحظي واتجاه الحضيض، مما يغير تطور طور الموجة حتى عندما تتغير العناصر المدارية المتوسطة ببطء.

2.7 مؤثر طفح فص روش والغلاف المشترك

يمثل هذا المؤثر انتقال الكتلة في الثنائيات المدمجة عبر طفح فص روش (RLOF) إلى جانب مقاومة اختيارية للغلاف المشترك (CE) (المعلمات في الجدول 7). يستهدف التنفيذ تكاملات طويلة متينة: فالخطوات الفرعية الداخلية تقيد و، وتمنع حراس الاندماج التباعدات، وتكشف التشخيصات تغير الطاقة المتراكم من التحديثات غير الهاملتونية. بعد كل خطوة فرعية، تُنقل المحاكاة إلى مركز الكتلة المحدّث بحيث تعكس العناصر المدارية الكتلة والزخم المعاد توزيعهما.

هندسة روش وقانون فقد الكتلة.

يُقيّم نصف قطر روش المكافئ حجماً للمانح بصيغة Eggleton (Eggleton, 1983)،

(9)

حيث هو الفصل اللحظي بين المانح والمتراكم، ويتبع معدل فقد الكتلة اللحظي وصفة Ritter (Ritter, 1988)،

$${\dot{M}}_{\mathrm{d}}=-{\dot{M}}_0\exp \left[\frac{R_{\ast }-R_{\mathrm{L}}}{H_P}\right],$$ (10)

مع قص مقدار الأس ($\left|(R_{\ast }-R_{\mathrm{L}})/H_P\right|\le 80$) لتجنب الفيض العددي. ويوفر الجسيم المانح و${\dot{M}}_0$ بوصفهما rlmt_Hp و rlmt_mdot0.

القنوات المحافظة والجهازية للكتلة.

خلال خطوة فرعية حجمها ، يفقد المانح $m_{\mathrm{loss}}=-{\dot{M}}_{\mathrm{d}}\mathrm{\Delta }t>0$، وتنقسم هذه الكمية إلى جزء متراكم وجزء ريحي:

وتكون الزيادة الصافية في كتلة المتراكم خلال الخطوة الفرعية ؛ أما كتلة الريح فتغادر النظام. ولـ jloss_mode (إعادة الانبعاث متناحية الخواص)، تنفذ الشفرة ذلك كتحديث ذي مرحلتين: تُنقل أولاً كامل إلى المتراكم، ثم تُقذف من المتراكم على نحو متناحي الخواص؛ وتبقى الكتلة المحتفظ بها صافيةً .

فقد (الزخم الزاوي للريح).

يُوصَف الزخم الزاوي النوعي الذي تحمله الريح بمعلمة نمط:

mode 0: ${𝐯}_{\mathrm{loss}}={𝐯}_{\mathrm{d}}(\text{wind leaves from the donor; implicit AM loss via donor mass decrement}),$

mode 1: $\text{isotropic re-emission from the accretor: }{m}_{\mathrm{loss}}\text{is transferred, then}$

$m_{\mathrm{wind}}\text{is ejected isotropically from the accretor}$

$(\text{no recoil in the accretor frame; orbital AM of the expelled mass is removed}$

$\text{implicitly by the accretor mass decrement}),$

mode 2: ${𝐯}_{\mathrm{loss}}={𝐯}_{\mathrm{CM}},(\text{wind removes centre-of-mass momentum }{m}_{\mathrm{wind}}{𝐯}_{\mathrm{CM}};$

$\text{the donor’s mass decrement still subtracts }{m}_{\mathrm{wind}}({𝐫}_{\mathrm{d}}-{𝐑}_{\mathrm{CM}})\times {𝐯}_{\mathrm{d}}),$

$|{𝐯}_{\mathrm{loss}}-{𝐯}_{\mathrm{emit}}|=f_j{j}_{\mathrm{orb}}/r\text{along }{\widehat{e}}_{⟂}⟂\widehat{𝐫},$

$j_{\mathrm{orb}}\equiv J/M=\mu |𝐫\times {𝐯}_{\mathrm{rel}}|/(M_{\mathrm{d}}+M_{\mathrm{a}}),\mu =\frac{M_{\mathrm{d}}{M}_{\mathrm{a}}}{M_{\mathrm{d}}+M_{\mathrm{a}}},$

$\mathrm{\Delta }{𝐋}_{\mathrm{wind}}=m_{\mathrm{wind}}({𝐫}_{\mathrm{emit}}-{𝐑}_{\mathrm{CM}})\times ({𝐯}_{\mathrm{loss}}-{𝐯}_{\mathrm{emit}}),$

$|\mathrm{\Delta }{𝐋}_{\mathrm{wind}}|=m_{\mathrm{wind}}{f}_j{j}_{\mathrm{orb}}{\displaystyle \frac{|{𝐫}_{\mathrm{emit}}-{𝐑}_{\mathrm{CM}}|}{r}}.$

ولا يُفرض عادة أن يصطف اتجاه العزم مع الزخم الزاوي المداري. إذا كان الاصطفاف الصارم أو المقدار الكامل لـ مطلوباً، ففسّر على أنه عامل مقياس فعال $r/|{𝐫}_{\mathrm{emit}}-{𝐑}_{\mathrm{CM}}|$ و/أو قيّد ${\widehat{e}}_{⟂}$ بالمستوى المداري.

ملاحظات. بالنسبة إلى النمط 0، لا يطبق عزم ريحي إضافي لأن ${𝐯}_{\mathrm{loss}}={𝐯}_{\mathrm{emit}}$ والزخم الزاوي المعني يُزال ضمنياً بنقصان كتلة المانح. وبالنسبة إلى النمط 1، تُزال الكتلة المقذوفة مباشرة من المتراكم (إعادة انبعاث متناحية الخواص)، ومن ثم يُزال زخمه الزاوي المداري ضمنياً بنقصان كتلة المتراكم ولا يطبّق تصحيح صريح لـ «ريح–». وبالنسبة إلى النمط 3، يعمل تصحيح صريح لـ ريح– عبر ${𝐯}_{\mathrm{loss}}-{𝐯}_{\mathrm{emit}}$. وبالنسبة إلى النمط 2، يتلاشى الحد الصريح ريح– بحكم البناء عندما تؤخذ نقطة الانبعاث عند مركز الكتلة.

الزخم الزاوي للنقل المحافظ.

في غياب العزوم الخارجية، ينبغي ألا يغير نقل الكتلة المحافظ الزخم الزاوي المداري للنظام. يحسب المؤثر فرق الزخم الزاوي بين وضع الكتلة المنقولة عند المانح ووضعها عند المتراكم (كلاهما بسرعة المانح) ويطبق تصحيح سرعة بعزم صافٍ أصغري لفرض . هنا للأنماط 0 و2 و3، أما في إعادة الانبعاث متناحية الخواص (النمط 1) فتستخدم مرحلة النقل ويُعالج القذف اللاحق على حدة عند المتراكم.

الزخم الخطي.

تُعالج المراكمة المحافظة داخلياً: تضاف الكتلة المنقولة إلى المتراكم بسرعة المانح اللحظية، في حين تُخفض كتلة المانح؛ وهذا يترك الزخم الكلي للزوج بلا تغيير في مرحلة النقل.

بالنسبة إلى الأنماط 0 و2 و3، تُزال كتلة الريح الجهازية مباشرة من المانح، ولذلك فإن نقصان كتلة المانح يزيل بالفعل $m_{\mathrm{wind}}{𝐯}_{\mathrm{d}}$. ثم يطبق المؤثر إزاحة منتظمة $-m_{\mathrm{wind}}({𝐯}_{\mathrm{loss}}-{𝐯}_{\mathrm{d}})$ على الزوج (موزعة بوصفها زيادة مشتركة في السرعة)، بحيث يكون تغير الزخم الكلي لهذه الأنماط

$$\mathrm{\Delta }𝐏=-m_{\mathrm{wind}}{𝐯}_{\mathrm{d}}-m_{\mathrm{wind}}({𝐯}_{\mathrm{loss}}-{𝐯}_{\mathrm{d}})=-m_{\mathrm{wind}}{𝐯}_{\mathrm{loss}},$$

أي إن النظام يفقد بالضبط الزخم الذي تحمله الريح الهاربة في الإطار العطالي.

بالنسبة إلى النمط 1 (إعادة الانبعاث متناحية الخواص)، تنفذ الريح بوصفها فقداً كتلياً من المتراكم: بعد نقل إلى المتراكم، تزيل الشفرة من المتراكم على نحو متناحي الخواص (من دون ارتداد في إطار المتراكم). ومن ثم يكون فقد الزخم المناظر ضمنياً عبر نقصان كتلة المتراكم، مع $\mathrm{\Delta }{𝐏}_{\mathrm{wind}}=-m_{\mathrm{wind}}{𝐯}_{\mathrm{a}}$ مقيمة عند وقت القذف.

بعد تحديثات الكتلة والسرعة في كل خطوة فرعية، تُنقل المحاكاة إلى مركز الكتلة للحفاظ على اتساق الإطار المرجعي.

رايات بين-وحدية.

بعد كل خطوة فرعية، يكتب المؤثر الخاصيتين المنطقيتين inside_CE وrlof_active على كل من المانح والمتراكم. وتشير هاتان الخاصيتان إلى ما إذا كان المتراكم يقع داخل نصف قطر المانح، وإلى ما إذا كانت قناة RLOF تزيل الكتلة بنشاط. وتفحص مؤثرات أخرى، مثل الرياح النجمية، هاتين الرايتين (بالاقتران مع مفاتيح التعطيل الخاصة بها) لتعليق عملها أثناء أطوار الغلاف المشترك أو أثناء تشغيل RLOF.

2.7.1 مقاومة الغلاف المشترك (CE)

عندما يكون المتراكم داخل المانح ()، تستطيع الشفرة تطبيق احتكاك ديناميكي وفق صيغة Ostriker (Ostriker, 1999; Ivanova et al., 2013)،

$${𝐚}_{\mathrm{DF}}=-\frac{4\pi G^2{M}_{\mathrm{a}}\rho }{v_{\mathrm{rel}}^3}I(ℳ){𝐯}_{\mathrm{rel}},ℳ=\frac{v_{\mathrm{rel}}}{c_s},$$

مع

محدوداً بـ لضمان الاستمرارية. ويمكن إضافة حد هندسي $-\pi \rho R_{\mathrm{a}}^2{Q}_d{v}_{\mathrm{rel}}{𝐯}_{\mathrm{rel}}/M_{\mathrm{a}}$. وتأتي بنية الغلاف إما من جدول يقدمه المستخدم (ce_profile_file) مع استيفاء لوغاريتمي–لوغاريتمي، أو من قانون قوة ، . وتُحد دفعة السرعة لكل خطوة فرعية بواسطة $|\mathrm{\Delta }𝐯|\le \text{ce\_kick\_cfl}{c}_s$. افتراضياً يكون الغلاف مصرفاً خارجياً: لا يطبق رد فعل معاكس على المانح؛ ويمكن تفعيل ذلك باستخدام ce_reaction_on_donor=1.

2.7.2 الجوانب العددية وسير التدفق

يقسم كل استدعاء للمؤثر المجال المطلوب إلى ما لا يقل عن rlmt_min_substeps خطوة فرعية؛ ويتكيف الحجم لتحقيق $|\mathrm{\Delta }M|/M\le \text{rlmt\_substep\_max\_dm}$ و $|\mathrm{\Delta }r|/r\le \text{rlmt\_substep\_max\_dr}$. ويُفعّل حارس الاندماج عندما $r\le {\epsilon }_{\mathrm{merge}}$ (معلمة المستخدم merge_eps أو القيمة الافتراضية مع حد أدنى موجب صغير). وعند الاندماج يُزال المتراكم وينفصل المؤثر لتفادي انجراف الفهارس في أنظمة -جسمية. وينفصل المؤثر بنفسه إذا اختفى أي من المكوّنين (كتلة غير موجبة)، لكنه لا يحذف جسيمات التتبع عديمة الكتلة غير ذات الصلة. بعد تحديثات الكتلة أو المقاومة، يُعاد تمركز النظام على مركز الكتلة. وإذا كان مؤثر stellar_evolution_sse الاختياري موجوداً، فإنه يُستدعى بخطوة زمنية صفرية بعد تحديثات كتلة RLOF بحيث تحدّث تغيرات الكتلة أنصاف الأقطار واللمعانيات النجمية فوراً. ويُصدّر التغير المتراكم في الطاقة المدارية ثنائية الجسم بين المانح والمتراكم على مدى الاستدعاء باسم rlmt_last_dE للتشخيص؛ ولأن المخطط يتضمن مصارف فيزيائية (رياحاً ومقاومة)، فلا يُتوقع أن تتلاشى هذه القيمة.

2.7.3 الاستقرار الديناميكي والتراكم المحدود بإدنغتون

لا يفرض مؤثر roche_lobe_mass_transfer معيار استقرار تحليلياً صريحاً لـ RLOF. وبدلاً من ذلك، تنشأ طبيعة انتقال الكتلة (مستقرة، أو انفلاتية، أو مؤدية إلى التلامس) من التطور المترابط المعتمد على الزمن لنصف قطر المانح، وهندسة روش، والاستجابة المدارية تحت تحديثات الكتلة والزخم المطبقة.

تراكم اختياري محدود بإدنغتون. يمكن للمؤثر اختيارياً أن يفرض حداً أقصى لمعدل التراكم الصافي على المتراكم، بقصد تمثيل حد إدنغتون (أو أي سقف يحدده المستخدم). إذا وفر المتراكم rlmt_mdot_edd (معلمة جسيم؛ بوحدات كتلة/زمن في نظام وحدات المحاكاة)، فإن الكتلة المحتفظ بها خلال كل خطوة فرعية داخلية تُحد عند $m_{\mathrm{acc}}\le \text{rlmt\_mdot\_edd}\mathrm{\Delta }t$. وأي فائض فوق هذا السقف يُعامل تلقائياً بوصفه فقداً كتلياً غير محافظ ويُزال من المتراكم على نحو متناحي الخواص (إعادة انبعاث متناحية الخواص؛ من دون ارتداد في إطار المتراكم). يعمل هذا السقف إضافة إلى كسر الفقد الجهازي الذي يتحكم به المستخدم rlmt_loss_fraction. وإذا غاب rlmt_mdot_edd أو $\le 0$، يُعطل السقف ويعود المؤثر إلى انتقال (غير) محافظ تتحكم به بالكامل معلمات المستخدم.

2.7.4 المعلمات

Table 7: معلمات طفح فص روش وتطور الغلاف المشترك.

Name (scope)	Unit	Default	Purpose
rlmt_donor (op)	—	—	Donor index (double, cast to int)
rlmt_accretor (op)	—	—	Accretor index (double, cast to int)
rlmt_Hp (part)	length	—	Pressure scale height (donor)
rlmt_mdot0 (part)	M/t	—	Reference overflow rate (donor)
rlmt_mdot_edd (part, accretor)	M/t	0	Maximum allowed net accretion rate (Eddington cap). If $>0$, the operator caps $m_{\mathrm{acc}}\le \text{rlmt\_mdot\_edd}\mathrm{\Delta }t$ and re-emits any excess isotropically from the accretor.
rlmt_loss_fraction (op)	—	0	Wind fraction $f_{\mathrm{loss}}$
jloss_mode (op)	—	0	Wind $j$ prescription (0–3)
jloss_factor (op)	—	1	Scale $f_j$ for mode 3
rlmt_skip_in_CE (op)	bool	1	Skip RLOF if
rlmt_substep_max_dm (op)	—	${10}^{-3}$	Max $\left|\mathrm{\Delta }M\right|/M$ per sub–step
rlmt_substep_max_dr (op)	—	$5\times {10}^{-3}$	Max $\left|\mathrm{\Delta }r\right|/r$ per sub–step
rlmt_min_substeps (op)	int	3	Minimum sub–steps
ce_profile_file (op)	path	—	Table $(s,\rho ,c_s)$ for CE
ce_rho0 (op)	dens	—	CE density normalisation
ce_alpha_rho (op)	—	0	Density slope
ce_cs (op)	speed	—	CE sound-speed normalisation
ce_alpha_cs (op)	—	0	Sound-speed slope
ce_xmin (op)	—	${10}^{-4}$	Coulomb cutoff
ce_Qd (op)	—	0	Geometric drag coefficient
ce_kick_cfl (op)	—	1	Velocity–kick limiter
ce_reaction_on_donor (op)	bool	0	Apply opposite CE kick to donor
merge_eps (op)	length	$0.5\min (R_{\ast },R_{\mathrm{a}})$	Merge radius (with floor)
rlmt_last_dE (op, out)	energy	—	Energy change (diagnostic)

ملاحظات.

تُقرأ جميع المعلمات القياسية بوصفها أعداداً مضاعفة الدقة في REBOUNDx وتحوّل إلى أعداد صحيحة للفهارس/المفاتيح. وإذا كان ذلك متاحاً في البناء، يمكن توفير اسم ملف في ce_profile_file لتحميل جدول CE؛ وإلا فيُستخدم نموذج قانون القوة.

3 الملخص والاستنتاجات

لقد عرضنا حزمة من وحدات تطور الثنائيات لامتداد REBOUNDx إلى مكامل REBOUND -الجسمي، وهي تجسد الآليات المهيمنة التي تقود التطور المداري والعلماني للثنائيات المتقاربة. وتتألف الوحدات من: (i) مؤثر لطفح فص روش والغلاف المشترك يجمع انتقال الكتلة المحافظ مع فقد كتلي نوعي غير محافظ وقابل للتهيئة لـ واحتكاك ديناميكي معتمد على عدد ماخ؛ (ii) وصفة مبسطة للتطور النجمي تزود أنصاف أقطار ولمعانيات معتمدة على الكتلة؛ (iii) ثلاث قنوات رياح (Reimers، ومدفوعة حرارياً، وفوق إدنغتون) تزيل الكتلة على نحو متناحي الخواص باستخدام تلك الخواص النجمية؛ (iv) مؤثر كبح مغناطيسي ينفذ عزم Verbunt–Zwaan/Kawaler مع تحديث مغلق الشكل للدوران واعٍ بالتشبع؛ و(v) وحدة ما بعد نيوتنية تقدم تصحيحات كتل نقطية ودوران–دوران عند 2 PN مع رد فعل إشعاعي للموجات الثقالية عند 2.5 PN.

خوارزمياً، يتتبع مؤثر RLOF/CE تبادلات الزخم الخطي والزاوي بين المانح والمتراكم والكتلة الهاربة، فارضاً حفظ الزخم الخطي للنقل المحافظ ومطبقاً عزماً أدنى للحفاظ على اتساق الزخم الزاوي. وتثبّت قيود الخطوات الفرعية على التغيرات الكسرية في الكتلة والفصل، وحراس الاندماج، ومحددات الدفعات لمقاومة CE، التطور قرب التلامس. وتنقل الرايات بين-الوحدية بدء RLOF والانغمار في غلاف مشترك إلى مؤثرات الرياح، في حين يُحدّث التطور النجمي المبسط بالتزامن مع تغيرات الكتلة. وجميع المؤثرات مستقلة عن الوحدات، وتعتمد على مقاييس شمسية وزمنية يوفرها المستخدم للتحويل بين الوحدات الفيزيائية ووحدات الشفرة. وتمكّن هذه المكونات مجتمعة من إجراء بحوث متسقة ذاتياً ومحلولة زمنياً للثنائيات المتقاربة ضمن أنظمة معزولة أو ضمن بيئات -جسمية غنية ديناميكياً.

توافر البيانات

البيانات والشفرة المستخدمة في هذا العمل متاحتان للتنزيل من GitHub.

الملحق: اختبارات الشفرة وأمثلة

رد الفعل الإشعاعي ما بعد النيوتني عند 2.5 PN: REBOUNDx مقارنةً بـ phi-GPU

في الشكل 1 نعرض مثالاً لتطور نصف المحور الأكبر للثنائية بسبب رد الفعل الإشعاعي عند 2.5 PN. يُهيأ النظام بـ ، ثم يبدأ بفقد الطاقة واللولبة إلى الداخل حتى يندمج. ونُظهر كذلك مقارنة مع شفرة phi-GPU (Berczik et al., 2011; Khan et al., 2012)

Listing 1: 2.5 PN Python example (REBOUNDx)

⬇

import rebound

import reboundx

sim = rebound.Simulation()

sim.G=1

sim.add(m=1.00000001E-01,x= -2.86636483E+00,y= -3.51352660E+00,z= 1.31048588E+00,\

vx=-8.29922019E+00,vy= -2.77660036E+00,vz= -1.16069792E-01)

sim.add(m=5.99999987E-02,x= -2.86626798E+00,y= -3.51381408E+00,z= 1.31048970E+00,\

vx=1.38346340E+01,vy= 4.61131049E+00,vz= 1.91195042E-01)

sim.move_to_com()

rebx = reboundx.Extras(sim)

pn = rebx.load_force(’post_newtonian’)

rebx.add_force(pn)

pn.params[’pn_2PN’]=0

pn.params[’pn_25PN’]=1

pn.params[’c’]=457.14

sim.integrate(4.2)

Figure 1: تطور نصف المحور الأكبر للثنائية بسبب رد الفعل الإشعاعي عند 2.5 PN، في REBOUNDx (هذا العمل) وphi-GPU.

الكبح المغناطيسي مع المدود

في الشكل 2 نعرض مثالاً لتطور نظام ثنائي بسبب الكبح المغناطيسي مع المدود.

في هذا الإعداد، يبقى تشغيل “من دون كبح مغناطيسي” شبه ثابت لأن المدار يبدأ دائرياً ولأن النجمين يبدآن متزامنين تماماً، لذلك لا تكاد توجد مدود تحتاج إلى تخميد وتبقى نصف المحور الأكبر وفترات الدوران شبه ثابتة. وعندما يُفعّل الكبح المغناطيسي، فإنه يبطئ دوران النجمين فوراً، فيجعلهما دون التزامن قليلاً؛ ثم تعمل المدود على استعادة التزامن بنقل الزخم الزاوي من المدار إلى دوراني النجمين، لكن الكبح المغناطيسي يستمر في إزالة ذلك الزخم الزاوي المغزلي إلى مصرف خارجي. وتكون المحصلة استنزافاً علمانياً للزخم الزاوي المداري، ولذلك ينكمش نصف المحور الأكبر باطراد مقارنة بحالة عدم وجود كبح مغناطيسي، بينما تميل فترات دوران النجمين إلى الازدياد (إبطاء الدوران) نسبة إلى تشغيل عدم وجود الكبح المغناطيسي.

Listing 2: Magnetic braking + tides Python example

⬇

import numpy as np

import rebound

import reboundx

# ==========================================================

# Helper: build one binary simulation

# ==========================================================

def make_binary(with_mb=True):

sim = rebound.Simulation()

sim.units = (’AU’, ’yr’, ’Msun’)

sim.integrator = "ias15"

sim.ri_ias15.epsilon=0

rebx = reboundx.Extras(sim)

# ---------- Binary ICs ----------

AU_per_Rsun = 1.0/215.032

m1 = m2 = 1.0

a0 = 0.06

e0 = 0.0 # circular orbit

R1 = R2 = 1.0 * AU_per_Rsun

# Equal-mass binary in COM frame (for e0=0, this is circular)

vfac = np.sqrt(sim.G*(m1+m2)/a0) * np.sqrt((1.0 - e0)/(1.0 + e0))

sim.add(m=m1, r=R1,

x=-a0*m2/(m1+m2),

vy=-vfac*m2/(m1+m2))

sim.add(m=m2, r=R2,

x=+a0*m1/(m1+m2),

vy=+vfac*m1/(m1+m2))

sim.move_to_com()

# Smaller dt for tighter orbit

sim.dt = sim.particles[1].P / 20.0

# ---------- Spins & structure ----------

k2_moi = 0.7

I1 = k2_moi*m1*R1**2

I2 = k2_moi*m2*R2**2

# Make both stars initially synchronous with the orbit

P_orb_yr = sim.particles[1].P

omega_sync = 2.0*np.pi / P_orb_yr # rad/yr

for p, I, R in zip(sim.particles[:2], [I1, I2], [R1, R2]):

p.params["Omega"] = np.array([0.0, 0.0, omega_sync])

p.params["I"] = I

p.r = R

# Tidal parameters (used by tides_spin or tctl)

p.params["k2"] = 0.2

p.params["tau"] = 0.01/365.25 # 0.01 days; moderately strong tides

# ---------- Magnetic braking ----------

if with_mb:

mb = rebx.load_operator("magnetic_braking")

rebx.add_operator(mb)

mb.params["mb_Msun"] = 1.0

mb.params["mb_Rsun"] = AU_per_Rsun

mb.params["mb_year"] = 1.0

mb.params["mb_Rossby_sat"] = 0.1

# Stronger than physical, but good for a demo

mb.params["mb_K"] = 1e55 # cgs; default is ~2.7e47

for p in sim.particles[:2]:

p.params["mb_on"] = 1 # enable MB on this star

p.params["mb_convective"] = 1

p.params["mb_tau_conv"] = 12.5/365.25 # ~12.5 days

else:

# No MB operator at all in the "no MB" run

print("No␣magnetic_braking␣operator␣added␣for␣this␣run")

# ---------- Tides ----------

use_spin = True

try:

tides = rebx.load_force("tides_spin")

rebx.add_force(tides)

use_spin = True

print("Loaded␣tides_spin␣for", "MB" if with_mb else "no-MB", "run")

except Exception:

print("tides_spin␣not␣available;␣using␣tides_constant_time_lag␣for",

"MB" if with_mb else "no-MB", "run")

use_spin = False

try:

tides = rebx.load_force("tides_constant_time_lag")

rebx.add_force(tides)

except Exception:

tides = rebx.load_operator("tides_constant_time_lag")

rebx.add_operator(tides)

# tctl parameters

for p in sim.particles[:2]:

p.params["tctl_k2"] = 0.2

p.params["tctl_tau"] = 0.1/365.25

Om = p.params["Omega"]

p.params["OmegaMag"] = float(np.linalg.norm(Om))

return sim, use_spin, R1, R2

# ==========================================================

# Utility functions

# ==========================================================

def elems(sim):

"""Return␣(a,␣e)␣of␣secondary␣relative␣to␣primary."""

p1 = sim.particles[1]

p0 = sim.particles[0]

o = p1.orbit(primary=p0)

return o.a, o.e

def prot_days(sim, i):

"""Rotation␣period␣of␣particle␣i␣in␣days."""

w = np.linalg.norm(sim.particles[i].params["Omega"])

return (2.0*np.pi/w)*365.25

# ==========================================================

# Build two simulations: with MB and without MB

# ==========================================================

sim_mb, use_spin_mb, R1_mb, R2_mb = make_binary(with_mb=True)

sim_nom, use_spin_nom, R1_nom, R2_nom = make_binary(with_mb=False)

# ==========================================================

# Integrate both and compare

# ==========================================================

t_end = 500 # years

N = 200

times = np.linspace(0.0, t_end, N)

for t in times:

# --- With magnetic braking ---

sim_mb.integrate(t)

# --- Without magnetic braking ---

sim_nom.integrate(t)

Figure 2: اللوحة اليسرى: تطور نصف المحور الأكبر للثنائية بسبب مزيج من المدود والكبح المغناطيسي (_mb) مقارنة بحالة اسمية من دون كبح مغناطيسي (_nom). اللوحة اليمنى: تطور فترة دوران النجم الأولي في وجود الكبح المغناطيسي وغيابه.

الكبح المغناطيسي مع المدود: REBOUNDx مقارنةً بـ MESA

في الشكل 3 نعرض مثالاً لتطور نظام ثنائي بسبب الكبح المغناطيسي مع المدود، مقارنةً بـ MESA (Paxton et al., 2011, 2013, 2015, 2019).

لإنشاء مسار MESA المعروض في الشكل 3، نفذنا وحدة “extras” ثنائية مخصصة (run_binary_extras) تتجاوز فقد الزخم الزاوي بالكبح المغناطيسي في MESA بربط نداء راجع يحدده المستخدم بسائق الثنائية. في extras_binary_controls نربط other_jdot_mb، الذي يحسب في كل خطوة عزماً للكبح المغناطيسي يعامل كحد لفقد الزخم الزاوي المداري وفق تقريب الاقتران المتزامن ويسنده مباشرة إلى b%jdot_mb. وتقيّم الروتين compute_mb القيمة من الدور المداري اللحظي، وتعرّف عتبة تشبع إما من تجاوز صريح أو من وصفة Rossby ، مع تمرير التطبيع و و عبر x_ctrl (وتُعتمد قيم افتراضية معقولة إن لم تُضبط). ثم تُشكّل القيمة الكلية بجمع مساهمات النجمين، متدرجة مثل ومتحولة من نظام غير مشبع إلى نظام مشبع عندما . ومن أجل الشفافية والتصحيح، نسجل بالإضافة إلى ذلك ثلاثة أعمدة إضافية في binary_history هي (jdot_mb_custom وomega_orb وomega_sat) لتمكين مقارنة تشخيصية مباشرة بين قانون الكبح المفروض والتطور الثنائي الناتج.

Listing 3: Magnetic braking + tides Python example (REBOUNDx)

⬇

import numpy as np

import rebound

import reboundx

# --- Parameters ---

K2 = 0.07

K2_MOI = 0.25

TAU_CONV_YR = 12.5 / 365.25

TAU_TIDE_YR = 0.1 / 365.0

AU_PER_RSUN = 1.0 / 215.032

def build_binary_mb() -> rebound.Simulation:

"""Equal-mass␣circular␣binary␣with␣tides_spin␣+␣magnetic␣braking."""

sim = rebound.Simulation()

sim.units = ("AU", "yr", "Msun")

sim.integrator = "whfast"

rebx = reboundx.Extras(sim)

# --- Binary ICs ---

m1 = m2 = 1.0

a0 = 0.06

e0 = 0.0

R1 = R2 = 1.0 * AU_PER_RSUN

v = np.sqrt(sim.G * (m1 + m2) / a0) * np.sqrt((1.0 - e0) / (1.0 + e0))

sim.add(m=m1, r=R1, x=-a0 * m2 / (m1 + m2), vy=-v * m2 / (m1 + m2))

sim.add(m=m2, r=R2, x=+a0 * m1 / (m1 + m2), vy=+v * m1 / (m1 + m2))

sim.move_to_com()

orb0 = sim.particles[1].orbit(primary=sim.particles[0])

sim.dt = orb0.P / 10.0

# --- Spins / structure (start synchronous) ---

omega_sync = 2.0 * np.pi / orb0.P

I1 = K2_MOI * m1 * R1**2

I2 = K2_MOI * m2 * R2**2

for p, I in zip(sim.particles[:2], (I1, I2)):

p.params["Omega"] = np.array([0.0, 0.0, omega_sync])

p.params["OmegaMag"] = float(np.linalg.norm(p.params["Omega"]))

p.params["I"] = I

p.params["k2"] = K2

p.params["tau"] = TAU_TIDE_YR

p.params["tctl_k2"] = K2

p.params["tctl_tau"] = TAU_TIDE_YR

# --- Magnetic braking ---

mb = rebx.load_operator("magnetic_braking")

rebx.add_operator(mb)

mb.params["mb_Msun"] = 1.0

mb.params["mb_Rsun"] = AU_PER_RSUN

mb.params["mb_year"] = 1.0

mb.params["mb_Rossby_sat"] = 0.1

mb.params["mb_K"] = 1e55 # demo-strength

for p in sim.particles[:2]:

p.params["mb_on"] = 1

p.params["mb_convective"] = 1

p.params["mb_tau_conv"] = TAU_CONV_YR

# --- Tides ---

tides = rebx.load_force("tides_spin")

rebx.add_force(tides)

return sim

if __name__ == "__main__":

sim = build_binary_mb()

t_end_yr = 500.0

n_steps = 200

for t in np.linspace(0.0, t_end_yr, n_steps):

sim.integrate(t)

Figure 3: تطور الدور المداري للثنائية بسبب مزيج من المدود والكبح المغناطيسي، في شفرتنا وفي MESA

طفح فص روش

في الشكل 4 نعرض مثالاً لتطور نظام ثنائي بسبب طفح فص روش. في هذه التجربة نستعيد الاستجابة المميزة ذات الطورين للشذوذ المتوقعة لانتقال كتلة معتمد بقوة على الطور في ثنائية شاذة. ولأن الطفح (هنا، وهو فعلياً قناة “تسرب” أسية بالنظر إلى ارتفاع السلم المختار) متمركز قرب الحضيض، فإن التطور المداري تقوده سلسلة من النبضات المتمركزة عند الحضيض بدلاً من عزم علماني متوسط على الطور، وتعتمد إشارة التغير الصافي في الشذوذ بحساسية على نسبة الكتلتين اللحظية. وبينما يبقى المانح هو المكوّن الأكثر كتلة، يميل النقل المتركز عند الحضيض إلى تدوير المدار: إذ إن إعادة توزيع الكتلة والزخم (مع تصحيحنا الأصغري لاتساق الزخم الزاوي) تخفض تفضيلياً الشذوذ التماسكي، منتجة الانخفاض الأولي. وحالما يعبر النظام انقلاب نسبة الكتلتين، ينتقل النقل نفسه الموضع عند الحضيض إلى نظام مثير للشذوذ: فتعمل التحديثات الآن على زيادة التباين بين حركة الحضيض وحركة الأوج، ولأن مؤثرنا غير هاملتوني ولا يفرض حفظ الطاقة المدارية حتى عندما يكون النقل محافظاً شكلياً، ولعدم وجود تدوير مدّي يقاوم الإثارة، يمكن للشذوذ التماسكي أن ينمو علمانياً. عملياً، يقود ذلك إلى انفلات نحو شذوذات كبيرة جداً، مع اقتراب العناصر المدارية من الحد شبه المكافئ بينما تتراكم تغيرات صغيرة في توازن الطاقة–الزخم الزاوي تودع مراراً عند الطور المداري نفسه.

Listing 4: RLOF Python example

⬇

import math

import rebound

import reboundx

import numpy as np

# ---- Build a simple 2-body setup ----

sim = rebound.Simulation()

sim.units=[’AU’,’yr’,’Msun’]

sim.integrator = "ias15"

sim.dt= 0.1

sim.add(m=1.0, r=0.00465) # donor (index 0)

sim.add(m=0.5, a=2., e=0.4, r=0.00465) # accretor (index 1)

sim.move_to_com()

rx = reboundx.Extras(sim)

rlmt = rx.load_operator("roche_lobe_mass_transfer")

rx.add_operator(rlmt)

# Required operator params

rlmt.params["rlmt_donor"] = 0 # donor is particle 0

rlmt.params["rlmt_accretor"] = 1 # accretor is particle 1

rlmt.params["rlmt_skip_in_CE"] = 1

# RLOF controls

rlmt.params["rlmt_loss_fraction"] = 0.0

rlmt.params["jloss_mode"] = 0

sim.particles[0].params["rlmt_Hp"] = 0.1

sim.particles[0].params["rlmt_mdot0"] = 1e-2

for t in np.logspace(1,6,100):

sim.integrate(t)

Figure 4: اللوحة اليسرى: تطور كتل مكوّنات الثنائية بسبب طفح فص روش. اللوحتان الوسطى واليمنى: تطور شذوذ الثانوي ونصف المحور الأكبر.

طفح فص روش: Rebound مقارنةً بـ MESA

في الشكل 5 نعرض مثالاً لتطور نظام ثنائي بسبب طفح فص روش، مقارنةً بـ MESA.

Listing 5: RLOF Python example for REBOUNDx–MESA comparison

⬇

import rebound

import reboundx

import numpy as np

# ---- Initiate the simulation ----

sim = rebound.Simulation()

sim.units=[’AU’,’yr’,’Msun’]

sim.integrator = "ias15"

# Add donor and accretor

sim.add(m=1.1, r=0.00465) # donor (index 0)

sim.add(m=0.4, P=0.2 / 365.0, e=1e-4, r=0.00465) # accretor (index 1)

# move the simulation to the center of mass

sim.move_to_com()

# -----------------------

# Roche lobe mass transfer

# -----------------------

rx = reboundx.Extras(sim)

rlmt = rx.load_operator("roche_lobe_mass_transfer")

rx.add_operator(rlmt)

# RLOF operator params

rlmt.params["rlmt_donor"] = 0

rlmt.params["rlmt_accretor"] = 1

rlmt.params["rlmt_skip_in_CE"] = 1 # Skip common envelope evolution

rlmt.params["rlmt_loss_fraction"] = 0.0 # No mass loss to winds

rlmt.params["jloss_mode"] = 0

rlmt.params["rlmt_substep_max_dm"] = 0.1 # To speed up the simulation

rlmt.params["rlmt_substep_max_dr"] = 1000.0 # To speed up the simulation

rlmt.params["rlmt_min_substeps"] = 3 # default value

# Donor particle parameters (index 0)

# Hp is set to a large number to keep mdot = mdot0

sim.particles[0].params["rlmt_Hp"] = 10_000.0

sim.particles[0].params["rlmt_mdot0"] = 1e-7

# Integrate

for time in np.logspace(1, np.log10(7.39e6), 1000):

sim.integrate(time)

Figure 5: اللوحة اليسرى: التطور الزمني لكتلتي المانح والمتراكم في نماذج Rebound وMESA قابلة للمقارنة، بسبب RLOF. اللوحة اليمنى: التطور الزمني للدور المداري للثنائية في الحالة نفسها.