-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnotes-lecture-neutrinos.html
1292 lines (1142 loc) · 55.7 KB
/
notes-lecture-neutrinos.html
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
<!DOCTYPE html>
<html lang="en">
<head>
<!-- 2025-02-06 Thu 09:04 -->
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>400A - Neutrinos & computational stellar evolution</title>
<meta name="author" content="Mathieu Renzo" />
<meta name="generator" content="Org Mode" />
<link rel="stylesheet" href="./css/style.css" />
<link rel="stylesheet" href="./fontawesome-free-6.7.2-web/css/all.min.css">
<meta name="keywords" content="Mathieu, Renzo, Mathieu Renzo,
stellar evolution, 400A, University of
Arizona, Steward Observatory, stars,
theoretical astrophysics">
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
displayAlign: "center",
displayIndent: "0em",
"HTML-CSS": { scale: 100,
linebreaks: { automatic: "false" },
webFont: "TeX"
},
SVG: {scale: 100,
linebreaks: { automatic: "false" },
font: "TeX"},
NativeMML: {scale: 100},
TeX: { equationNumbers: {autoNumber: "AMS"},
MultLineWidth: "85%",
TagSide: "right",
TagIndent: ".8em"
}
});
</script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS_HTML"></script>
</head>
<body>
<div id="preamble" class="status">
<!-- Preamble -->
<!-- The header -->
<div class="header">
<!-- The site name -->
<div class="site-name">
<a id="top" href="./index.html">Stellar Evolution</a>
</div>
<!-- The hamburger -->
<div class="hamburger">
<div id="myLinks" class="menu">
<a href="./index.html">Home</a>
<a href="./syllabus.html">Syllabus</a>
<a href="./lectures.html">Lectures</a>
<a href="./projects.html">Projects</a>
</div>
<a href="javascript:void(0);" class="icon" onclick="HamburgerMenuFunction()">
<i class="fa fa-bars"></i>
</a>
</div>
<div class="navbar">
<a href="./syllabus.html">Syllabus</a>
<a href="./lectures.html">Lectures</a>
<a href="./projects.html">Projects</a>
</div>
</div>
<!-- scripts -->
<script>
function HamburgerMenuFunction() {
var x = document.getElementById("myLinks");
if (x.style.display === "block") {
x.style.display = "none";
} else {
x.style.display = "block";
}
}
</script>
<script>
window.onload = () => {
const toggleButton = document.getElementById("light-dark-toggle");
const body = document.body;
const html = document.documentElement;
// Check localStorage for user preference
if (localStorage.getItem("theme") === "dark") {
body.classList.add("dark-mode");
html.classList.toggle("dark-mode");
}
// Toggle theme on click
toggleButton.addEventListener("click", () => {
body.classList.toggle("dark-mode");
html.classList.toggle("dark-mode");
// Save user preference
if (body.classList.contains("dark-mode")) {
localStorage.setItem("theme", "dark");
} else {
localStorage.setItem("theme", "light");
}
});
};
</script>
<!-- end scripts-->
</div>
<div id="content" class="content">
<p>
<b>Materials:</b> Onno Pols' notes Sec. 6.5 and 12.3.1, Chapter 7, Kippenhahn chapter
9, 10, 11
</p>
<div id="outline-container-org7bf990c" class="outline-2">
<h2 id="org7bf990c"><a href="#org7bf990c">Neutrino physics in stars</a></h2>
<div class="outline-text-2" id="text-org7bf990c">
<p>
We have almost seen all the pieces needed to make a one-dimensional
hydrostatic model of a star, including the input physics in the
quantities κ (opacity) and ε<sub>nuc</sub>. The aim of this lecture is to unpack
the last input quantities entering in the energy conservation which
connects stellar physics with neutrino physics: ε<sub>ν.</sub>
</p>
</div>
<div id="outline-container-orgaf3db0c" class="outline-3">
<h3 id="orgaf3db0c"><a href="#orgaf3db0c">Summary of equations we have derived</a></h3>
<div class="outline-text-3" id="text-orgaf3db0c">
</div>
<div id="outline-container-org560ca5d" class="outline-4">
<h4 id="org560ca5d"><a href="#org560ca5d">Mass conservation</a></h4>
<div class="outline-text-4" id="text-org560ca5d">
<div class="latex" id="org05a0a90">
\begin{equation}\label{eq:mass_cont}
\frac{dr}{dm} = \frac{1}{4\pi r^{2}\rho}\ \ .
\end{equation}
</div>
</div>
</div>
<div id="outline-container-org8933510" class="outline-4">
<h4 id="org8933510"><a href="#org8933510">Hydrostatic equilibrium</a></h4>
<div class="outline-text-4" id="text-org8933510">
<div class="latex" id="org96ca659">
\begin{equation}\label{eq:HSE}
\frac{dP}{dm} = -\frac{Gm}{4\pi r^{4}} \ \ ,
\end{equation}
</div>
<p>
which follows from the momentum conservation equation assuming the
acceleration to be negligible <b>N.B.:</b> this last assumption could be
relaxed to add an extra term in this equation analogous to ma in F=ma.
</p>
</div>
</div>
<div id="outline-container-orgf5ff9d9" class="outline-4">
<h4 id="orgf5ff9d9"><a href="#orgf5ff9d9">Equation of state</a></h4>
<div class="outline-text-4" id="text-orgf5ff9d9">
<div class="latex" id="org8d0a09c">
\begin{equation}
P_\mathrm{tot} = P_\mathrm{gas} + P_\mathrm{rad} = \frac{\rho}{\mu m_{u}}k_{B}T + P_{QM} + \frac{1}{3}aT^{4} \ \ .
\end{equation}
</div>
</div>
</div>
<div id="outline-container-org9b66410" class="outline-4">
<h4 id="org9b66410"><a href="#org9b66410">Energy transport</a></h4>
<div class="outline-text-4" id="text-org9b66410">
<div class="latex" id="org7f1274d">
\begin{equation}
\frac{dT}{dm} = \frac{T}{P}\frac{dP}{dm}\nabla
\end{equation}
</div>
<p>
where ∇ =∂ log(T)/∂ log(P) is the local temperature gradient, equal to
the radiative gradient in stably stratified regions:
</p>
<div class="latex" id="org1970364">
\begin{equation}
\nabla \equiv \nabla_\mathrm{rad} = \frac{3 P}{14\pi acGm T^{4}}\kappa L
\end{equation}
</div>
<p>
with κ = (1/κ<sub>rad</sub> + 1/κ<sub>cond</sub>)<sup>-1</sup> the combination "in parallel" of the
radiative and conductive opacity (assumed to be known from atomic
physics), and ∇ ≡ ∇<sub>ad</sub> the adiabatic gradient (within ∼10<sup>-7-8</sup>
precision) for convective regions. We also have a criterion
(Schwarzschild criterion) to determine which region is which.
</p>
</div>
</div>
<div id="outline-container-org7e03897" class="outline-4">
<h4 id="org7e03897"><a href="#org7e03897">Energy conservation</a></h4>
<div class="outline-text-4" id="text-org7e03897">
<div class="latex" id="org8c2ff9c">
\begin{equation}
\frac{dL}{dm} = \varepsilon_\mathrm{nuc} -\varepsilon_{\nu} + \varepsilon_\mathrm{grav} \ \ .
\end{equation}
</div>
<p>
with ε<sub>grav</sub> = T∂ s/∂ t the change in internal energy and ε<sub>ν</sub>>0
(⇒ it is always a <i>loss</i> term during stellar evolution).
</p>
</div>
</div>
<div id="outline-container-orgd542d5f" class="outline-4">
<h4 id="orgd542d5f"><a href="#orgd542d5f">Composition evolution</a></h4>
<div class="outline-text-4" id="text-orgd542d5f">
<div class="latex" id="orga36c32a">
\begin{equation}
\frac{dX_{i}}{dt} = A_{i}\frac{m_{u}}{\rho}\left(\sum_{k,l} r_{k,l} - \sum_{i,j} (1+\delta_{ij})r_{ij}\right) + D_\mathrm{mix}\frac{dX_{i}}{dr} \equiv \frac{dX_{i}}{dt}(T,\rho, {X_{j}}) \ \ ,
\end{equation}
</div>
<p>
where the first two terms represent production and destruction of the
i-th isotope by nuclear processes and the third term represent mixing
by various instabilities treated in diffusion approximation.
</p>
</div>
</div>
</div>
<div id="outline-container-org9b9ab60" class="outline-3">
<h3 id="org9b9ab60"><a href="#org9b9ab60">Two types of neutrinos</a></h3>
<div class="outline-text-3" id="text-org9b9ab60">
<p>
For the density typical in stellar interior (think of ⟨ρ<sub>☉</sub>⟩ and the other
average densities you estimated in the various homeworks), neutrinos
can safely be considered as non-interacting. This means we can
consider neutrinos as ultra-relativistic mass-less particles
propagating at v ≅ c. The mass less approximation also implies that we
neglect neutrino oscillations: even if a neutrino changes leptonic
flavor, it's still going to leave the star carrying away its energy!
</p>
<p>
<b>N.B.:</b> recall that the leptonic number is +1 for the leptons electron
e<sup>-</sup>, muon μ<sup>-</sup>, tau τ<sup>-</sup> and the corresponding neutrinos ν<sub>e</sub>, ν<sub>μ</sub>,
ν<sub>τ</sub> and -1 for their antiparticles positron e<sup>+</sup>, positive muon μ<sup>+</sup>,
and positive τ<sup>+</sup> and the corresponding antineutrinos.
</p>
<p>
Therefore, any neutrino produced in the stellar interior leaves the
star in the light-crossing time R/c ≤ 1000 R<sub>☉</sub>/c ≅ 1000 ×
2 sec = 2000 sec which is much shorter than evolutionary timescales
(and even during shorter evolutionary phases such as silicon shell
burning, the energy carried by the neutrinos is on its way out and
nothing is going to stop it from leaving the star).
</p>
<p>
These arguments apply to both kinds of neutrinos that we encounter in
stellar evolution (and are about to introduce), and stop working only
when densities get extremely high (ρ≥10<sup>10</sup> g cm<sup>-3</sup>): at this point
neutrinos can and do start interacting with matter. Such densities are
only encountered in a collapsing core that has exhausted nuclear fuel
and in neutron stars: neutrino interactions are crucial to the physics
of the supernova explosions, but can safely be neglected before.
</p>
<p>
During the evolution of the star that sets the stage for the
explosion, we encounter two "classes" of neutrinos.
</p>
</div>
<div id="outline-container-org2af9c50" class="outline-4">
<h4 id="org2af9c50"><a href="#org2af9c50">Nuclear neutrinos</a></h4>
<div class="outline-text-4" id="text-org2af9c50">
<p>
We have already encountered processes producing neutrinos occurring in
stars: <i>weak nuclear reactions</i>! These are needed since the main
sequence (both in the pp chain and CNO) to convert some protons into
neutrons for the helium nuclei. On the scale of r<sub>nuc</sub>, the neutrino
crossing time (that is the light crossing time) is so short that we
can neglect neutrino oscillations and assume conservation of leptonic
number to write down the equations.
</p>
<p>
These neutrinos do carry away energy as described above, <i>but</i> that is
just some of the energy that the exo-energetic thermonuclear reaction
would have released. Effectively, they decrease the energy gain from
releasing nuclear binding energy but they are not a net energy loss.
</p>
<p>
This allows to incorporate them in models by re-writing ε<sub>nuc</sub>
→ ε<sub>nuc</sub> - |ε<sub>ν, nuc</sub>|, where the neutrino term is always negative
(hence we can write it as minus an absolute value), and needs to be
calculated from nuclear physics, determining the spectrum of the
neutrino produced by a given reaction (a non-trivial task which
stellar physics leaves to neutrinos and nuclear physicists).
</p>
<p>
Particularly important for the late evolution of massive stars are
neutrinos from the so-called <i>URCA processes</i> (named after a casino in
Rio de Janeiro because when these processes kick in the energy of the
nuclear reaction goes the same way as the money in the casino!):
</p>
<div class="latex" id="orgd352bc6">
\begin{equation}
^{A}Z+e^{-}\rightarrow^{A}(Z-1) + \nu_{e }\\
^{A}(Z-1)\rightarrow ^{A}Z+ e^{+} +\bar{\nu_{e}}
\end{equation}
</div>
<p>
which produce one neutrino and one anti-neutrino without changing the
composition of the star. This requires that the nucleus \(^{A}(Z-1)\) is
unstable to β<sup>-</sup>-decay and the cross section for electron capture on
\(^AZ\) is non-negligible, which can happen during Si core burning and
the subsequent gravitational collapse of the core once the nuclear
fuel runs out.
</p>
<p>
The nuclear neutrinos are mostly sensitive to the core temperature for
the activation of certain thermonuclear reaction chains (except for
pycno-nuclear reaction at extremely high densities where the electron
screening effectively makes the Coulomb barrier negligible).
</p>
</div>
</div>
<div id="outline-container-org1c43654" class="outline-4">
<h4 id="org1c43654"><a href="#org1c43654">Thermal neutrinos</a></h4>
<div class="outline-text-4" id="text-org1c43654">
<p>
After helium core burning, the density in the stellar cores become
sufficiently high (because of the gravothermal collapse) that
non-nuclear processes producing neutrinos start occurring. After
carbon depletion, the neutrinos produced by these processes can take
away more energy than is locally lost to photons by each stellar
layer: <i>evolved massive stars are neutrino stars</i> L<sub>ν</sub>
≫ L<sub>rad</sub> (<a href="https://ui.adsabs.harvard.edu/abs/1968Ap%26SS...2...96F/abstract">Fraley 1968</a>).
</p>
<figure id="org6317da4">
<img src="./images/neutrino_HRD.png" alt="neutrino_HRD.png" width="100%">
<figcaption><span class="figure-number">Figure 1: </span>Left: evolutionary tracks on a photon HR diagram. Right: corresponding evolutionary tracks on a "neutrino" HR diagram. This is Fig. 2 from <a href="https://ui.adsabs.harvard.edu/abs/2020ApJ...893..133F/abstract">Farag et al. 2020</a>. Note the y-axis scale in both panels.</figcaption>
</figure>
<p>
This also effectively means that the stellar core of evolved (∼ during
and after carbon core burning) massive stars is <i>decoupled</i> from the
stellar envelope: the <i>gravothermal collapse of the core occurs to
compensate the neutrino losses from the core itself</i>! The thermal
timescale of the core becomes τ<sub>KH,ν</sub>≅ GM<sub>core</sub><sup>2</sup>/(2R<sub>core</sub> L<sub>ν</sub>) and
the nuclear timescale becomes τ<sub>nuc,ν</sub> = φ f<sub>burn</sub> Mc<sup>2</sup>/L<sub>ν</sub> both of which
are much shorter than the timescales in the low density, photon-cooled
envelope: in the late stages of stellar evolution the envelope should
be <i>frozen</i> and the core evolves driven by neutrino losses.
</p>
<p>
<b>N.B.:</b> it is still the energy losses driving the gravothermal collapse
because of the virial theorem that govern the evolution, but the
envelope does not have time to keep up with the core.
</p>
<p>
<b>N.B.:</b> recently, observations of early signals of stellar explosion
have questioned this picture of <i>frozen</i> envelope: there <i>may</i> be some
presently unknown phenomena happening on a <i>dynamical</i> timescale of the
envelope in the final years/months of a massive star evolution that
affect the envelope. The fact that they <i>need to be dynamical</i> to do
anything is related to the fact that the evolution of the core is sped
up by the thermal neutrino losses.
</p>
<p>
The neutrinos that do <i>not</i> come from nuclear reactions are a real
energy <i>loss</i> term for the star that enter in the local energy
conservation equation ε<sub>ν</sub> (which is also always negative!): dL/dm =
ε<sub>nuc</sub> - |ε<sub>ν</sub>| +ε<sub>grav</sub>.
</p>
<p>
The <i>thermal</i> processes producing these neutrinos typically will produce
a neutrino-antineutrino <i>pair</i> to conserve the leptonic number.
Typically only electron neutrinos will be relevant: leptons other than
e<sup>±</sup> are unstable and are not commonly found in stars. They fall into
several categories (here x={e,μ, or τ} but in stellar conditions it's
typically e):
</p>
<ul class="org-ul">
<li>plasmon processes: \(\gamma + \tilde{\gamma} \rightarrow \nu_x + \bar{\nu}_x\),</li>
<li>bremsstrahlung: \(e^{-} + ^{A}Z \rightarrow e^{-} + ^{A}Z + \nu_x + \bar{\nu}_x\),</li>
<li>pair-production: \(\gamma + \gamma \rightarrow \nu_x + \bar{\nu}_x\),</li>
<li>pair annihilation: \(e^{+}+e^{-} \rightarrow \nu_x + \bar{\nu}_x\) ,</li>
<li>photo-processes: \(\gamma +e^{-} \rightarrow e^{-} + \nu_x + \bar{\nu}_x\),</li>
</ul>
<figure id="orgf673fa6">
<img src="./images/feynman_diagram_neutrinos.png" alt="feynman_diagram_neutrinos.png" width="100%">
<figcaption><span class="figure-number">Figure 2: </span>Feynman diagrams for the dominant neutrino cooling processes. The top row shows the photo-emission processes, the middle row shows the e± annihilation processes, the bottom row shows the plasmon processes. The neutral (charged) current ractions are in the left (right) column. This is Figure 1 from <a href="https://ui.adsabs.harvard.edu/abs/1993ApJ...411..813A/abstract">Aufderheide 1993</a>, and Z and W represent the boson that mediate weak interactions: the left column shows interactions mediated by the neutral boson Z, while the right column shows interactions mediated by the charged boson W<sup>±</sup>.</figcaption>
</figure>
<p>
<b>N.B.:</b> "plasmons" are collective excitations of the stellar plasma that
propagate (the analogy is with "solitons" in fluid dynamics), and can
be quantized and coupled via quantum electrodynamics to e<sup>±</sup>.
</p>
<p>
<b>N.B.:</b> The Feynman diagrams of some of these processes above are
illustrating the processes (and the various pieces that contribute in
the quantum field theory calculation of the cross section), the
particles that are not "free legs" are mediators and not real
particles, do not over-interpret these diagrams as physical pictures!
</p>
<p>
Neutrino cooling processes are mostly sensitive to the core density ρ.
The typical energy carried away per neutrino-antineutrino pair is of
the order of the thermal energy of the electrons (i.e., their Fermi
energy if the region from where the neutrinos are emitted is partially
degenerate).
</p>
<p>
Lke κ and ε<sub>nuc</sub>, the energy losses to thermal neutrinos are usually
tabulated in stellar physics codes (see especially the widely used
<a href="https://ui.adsabs.harvard.edu/abs/1996ApJS..102..411I/abstract">Itoh et al. 1996</a>), and the figure below shows the |ε<sub>ν</sub>|
≡ |ε<sub>ν</sub>|(T,ρ) resulting from these tables:
</p>
<figure id="orgca84fd3">
<img src="./images/Trho_neutrinos.png" alt="Trho_neutrinos.png" width="100%">
<figcaption><span class="figure-number">Figure 3: </span>Neutrino energy losses on the T(ρ) plane. White lines mark the separation between regions where different neutrino emission processes dominate, the colored lines mark the T(ρ) tracks of helium cores of the labeled masses (in M<sub>☉</sub> units) computed with MESA. Credits: R. Farmer.</figcaption>
</figure>
</div>
</div>
</div>
</div>
<div id="outline-container-org78dbf0b" class="outline-2">
<h2 id="org78dbf0b"><a href="#org78dbf0b">Principles of computational stellar evolution</a></h2>
<div class="outline-text-2" id="text-org78dbf0b">
<p>
"<i>Traditional scientific knowledge has generally taken the form of</i>
<i>either theory or experimental data. However, where theory and</i>
<i>experiment stumble, simulations may offer a third way.</i>" - Simulation,
Johannes Lenhard et al.
</p>
<p>
We now have finally derived/discussed all the pieces of physics
necessary to compute a stellar <i>structure</i> model, treat its nuclear
energy generation and thus driving its <i>evolution</i>.
</p>
<p>
The description we have obtained for a spherically symmetric star is
made of <b>four non-linear, coupled, ordinary differential equations</b>
(ODE) <b>plus the equation of state that acts as closure condition</b> for
the system. Solving these equations we can:
</p>
<ol class="org-ol">
<li>study the interior structure of modelled stars and try to learn about
the parts of the star that are not accessible to direct
observations (hidden inside the photosphere)</li>
<li>study the time evolution of modelled stars, which we cannot
observe for real stars since most evolve way too slow for us to
follow within our lifetimes.</li>
</ol>
<p>
However, this system of coupled, non-linear ODEs is not easily solved
by hand. Since the early days of the availability of computers in the
1960s, people have been designing and leveraging <i>computational
techniques</i> to <i>numerically</i> solve this system of equations (see e.g.,
<a href="https://ui.adsabs.harvard.edu/abs/1959ApJ...129..628H/abstract">Henyey 1959</a>, <a href="https://ui.adsabs.harvard.edu/abs/1962ApJ...135..770I/abstract">Iben & Erhman 1962</a>, <a href="https://ui.adsabs.harvard.edu/abs/1970ApJ...159..619S/abstract">Sujimoto 1970</a> <a href="https://ui.adsabs.harvard.edu/abs/1971MNRAS.151..351E/abstract">Eggleton 1971</a>, <a href="https://ui.adsabs.harvard.edu/abs/1978ApJ...225.1021W/abstract">Weaver
et al. 1978</a>). For the rest of this lecture, we are going to discuss
some general principles behind these computational techniques.
</p>
</div>
<div id="outline-container-org9626394" class="outline-3">
<h3 id="org9626394"><a href="#org9626394">The most important thing: <i>Computer simulations are not empirical evidence</i></a></h3>
<div class="outline-text-3" id="text-org9626394">
<p>
Computational techniques take a system of equations describing a
physical model (for example the equations we have derived), which
already rely on a whole variety of physical approximations (e.g.,
spherical symmetry, LTE, free-streaming neutrinos, etc.), and apply a
whole new set of <i>numerical</i> approximations to obtain a numerical
solution.
</p>
<p>
Presumably, nature does not do <i>any</i> of this: without opening the
philosophical question of whether nature writes down equations to
solve, it certainly does not need to make physical approximations
(while we need to do it to reason on a problem and keep it
manageable), and even less make <i>numerical</i> approximations needed to
solve equations. At ontological level, numerical simulations are <i>not</i>
equivalent to empirical evidence! Never trust numerical results as if
they were ground truth: a computer will only do what it is told, and
we tell it our physically and numerically approximated best guess for
what we are trying to simulate which is <i>not</i> what nature empirically
provides. I emphasize this as a person who spends most of his time
making numerical simulations! Note also how this is more general than
<i>just</i> stellar physics: this applies to <i>any</i> computational physics field.
</p>
<p>
Because of these, it is always crucially important to <i>do resolution
tests</i>: when performing a simulation of a physical phenomenon, you
should always test that the scientific results you obtain do not
depend on how you discretize your equations in time and space and on how
you numerically solve them. This is often a painful task, but very
important to not fool ourselves!
</p>
<p>
"<i>All models are wrong, but some are useful!</i>" - G. Box
</p>
<p>
In stellar physics, there is another problem: the <i>high non-linearity</i>
of the coupled system of ODEs means that there can be chaotic
behavior! A small perturbation (maybe because of a numerical error or
a small inconsistency in the tabulated input physics!) can cascade
into dramatic consequences, similar to the famous "butterfly effect"
(which was also theorized from numerical computations, but of
atmospheric physics, by <a href="https://en.wikipedia.org/wiki/Edward_Norton_Lorenz">E. Lorenz</a>).
</p>
</div>
</div>
<div id="outline-container-orgc2b86dd" class="outline-3">
<h3 id="orgc2b86dd"><a href="#orgc2b86dd">Why did we limit ourselves to spherical symmetry?</a></h3>
<div class="outline-text-3" id="text-orgc2b86dd">
<p>
Throughout the course so far we have explicitly assumed spherical
symmetry, although since the beginning we have discussed some
phenomena that can break the global spherical symmetry (e.g.,
rotation, magnetic fields, or the presence of gravity and/or
irradiation from companion stars), and we have seen phenomena (e.g.,
convection) that break the spherical symmetry locally. Without the
assumption of spherical symmetry, the equations describing our system
would become partial differential equations (PDE), changing the
mathematical and thus computational approach.
</p>
<p>
In some cases, it is possible to treat some of the non-spherical
effects by casting them in a form that still allows a 1D formalism
with ODEs (e.g., "shellular rotation" assuming that all quantities are
constant along isobars and uses P as the independent coordinate, e.g.,
<a href="https://ui.adsabs.harvard.edu/abs/2000A%26A...361..159M/abstract">Maeder & Meynet 2000</a>, or using volume-weighted equivalent potentials
accounting for the gravity of a companion star, e.g., <a href="https://ui.adsabs.harvard.edu/abs/2015ApJS..220...15P/abstract">Paxton et al.
2015</a>). But ultimately the physics is non-spherical, so why do we
insist with this limited approximation?
</p>
<p>
The main reason is <i>necessity</i> (see for example <a href="https://ui.adsabs.harvard.edu/abs/2022ApJS..262...19J/abstract">Jermyn et al. 2022
section 5.4</a>) <b>the contrast of scale in both space and time</b> required to
follow the <i>evolution</i> of a star from birth to death (or to the current
age of the Universe, whichever comes first!) <b>is just too large</b> to be
manageable with present day <i>and</i> with foreseeable computational
capabilities.
</p>
</div>
<div id="outline-container-org95a4de0" class="outline-4">
<h4 id="org95a4de0"><a href="#org95a4de0">Spatial scale limitations</a></h4>
<div class="outline-text-4" id="text-org95a4de0">
<p>
Consider the convective envelope in the Sun. It has a Reynolds number
Re= v<sub>conv</sub> ℓ/ν ∼ 10<sup>12</sup>, which, using Kolmogorov's model of
turbulence implies that the turbulent cascade spans a range of scales
from L to ℓ with a contrast L/ℓ ∼ Re<sup>3/4</sup> ∼ 10<sup>9</sup>. To resolve this
contrast in a 3D simulation we would need this number of cells <i>in
each direction</i>, meaning ∼10<sup>27</sup> points.
</p>
</div>
</div>
<div id="outline-container-orgf3f2af4" class="outline-4">
<h4 id="orgf3f2af4"><a href="#orgf3f2af4">Temporal scale limitations</a></h4>
<div class="outline-text-4" id="text-orgf3f2af4">
<p>
Consider again the Sun. The free fall timescale is generously
τ<sub>ff</sub>∼1h, and its thermal timescale is τ<sub>KH</sub>∼1.5×10<sup>7</sup> years,
and current age is ∼ 5×10<sup>9</sup> years ≅ τ<sub>nuc</sub>. In 3D we would need to
resolve the dynamics of the gas, with timesteps Δ t ≤ τ<sub>ff</sub>. Even
considering an equal sign Δ t=τ<sub>ff</sub> (certainly insufficient to <i>resolve</i>
the dynamics of the stellar plasma), it would take ∼ τ<sub>Kh</sub> / τ<sub>ff</sub> = 10<sup>10</sup>
timesteps to calculate a thermal timescale and ∼ τ<sub>nuc</sub> / τ<sub>ff</sub> = 10<sup>12</sup>
timesteps for a nuclear timescale.
</p>
</div>
</div>
<div id="outline-container-orgcd80636" class="outline-4">
<h4 id="orgcd80636"><a href="#orgcd80636">Present day and future prospects</a></h4>
<div class="outline-text-4" id="text-orgcd80636">
<p>
The largest present day multi-dimensional stellar hydrodynamics
simulations reach maybe up to 10<sup>12</sup> resolution points, and about 10<sup>8</sup>
timesteps and using lower-than-realistic Reynolds number and/or
boosting the luminosity to drive dynamical effects faster than in
reality (see <a href="https://ui.adsabs.harvard.edu/abs/2022ApJ...928L..10A/abstract">Anders et al. 2022</a> or <a href="https://ui.adsabs.harvard.edu/abs/2024MNRAS.531.1316T/abstract">Thompson et al. 2024</a>).
</p>
<p>
Assuming that Moore's law (the number of transistors in a CPU roughly
doubles every year) - an assumptions that may clash with the
engineering reality, it will take several decades to have a
numerically resolved 3D simulation of convection in a star for a
convective turnover timescale, and longer for thermal and nuclear
timescale: 1D will remain necessary, and we need to keep representing
the complex multi-scale, multi-physics, and multi-dimensional problem
of stellar evolution considering only variations along the radial
direction!
</p>
<p>
Nevertheless, 3D stellar hydrodynamics of restricted problems ("box in
a star" approach - e.g., <a href="https://ui.adsabs.harvard.edu/abs/2024MNRAS.533..687R/abstract">Rizzuti et al. 2024</a>, short timescales such as
stellar explosions - e.g., <a href="https://ui.adsabs.harvard.edu/abs/2021ApJ...921...28F/abstract">Fields & Couch 2021</a>) are possible, and
stellar <i>structure</i> simulations in 2D are also at the forefront of
possibilities (see e.g., <a href="https://ui.adsabs.harvard.edu/abs/2023A%26A...677L...5M/abstract">Mombarag et al. 2023</a>).
</p>
</div>
</div>
</div>
<div id="outline-container-org05728e9" class="outline-3">
<h3 id="org05728e9"><a href="#org05728e9">Boundary conditions</a></h3>
<div class="outline-text-3" id="text-org05728e9">
<p>
So, for the foreseeable future we will need 1D stellar structure and
evolution calculation using the set of equations we derived. To solve
them, we need boundary conditions. These are crucial in determining
the solution of the structure at each timestep!
</p>
<p>
For the stellar problem, we need to specify boundary conditions at two
locations.
</p>
</div>
<div id="outline-container-org3944e35" class="outline-4">
<h4 id="org3944e35"><a href="#org3944e35">Center</a></h4>
<div class="outline-text-4" id="text-org3944e35">
<p>
These are fairly intuitive:
</p>
<ul class="org-ul">
<li>r(m=0) = 0. This is necessary to keep the local density ρ ∼ m/r<sup>3</sup> finite.</li>
<li>L(m=0) = 0. This is necessary to keep the energy generation per
unit volume finite (as the volume goes to zero)</li>
</ul>
<p>
However, we have 4 ODEs and these are only two boundary conditions, so
the system is still undetermined. The central pressure and density are
not a priori known (we can estimate them, but that is not precise
enough!).
</p>
</div>
</div>
<div id="outline-container-org96dec03" class="outline-4">
<h4 id="org96dec03"><a href="#org96dec03">Surface</a></h4>
<div class="outline-text-4" id="text-org96dec03">
<p>
We need to turn to the surface to get observationally informed
boundary conditions. Remember that from the point of view of a
detailed stellar evolution code, which assumes LTE to solve for the
internal structure, the "surface" (meaning: the outer boundary) is
usually defined at the photosphere. This is by definition the
"idealized" surface where T=T<sub>eff</sub> which corresponds to the location
outside of which LTE is not a good assumption anymore, the radiation
field is <i>not</i> isotropic, and the problem becomes the calculation of a
stellar <i>atmosphere</i>.
</p>
<p>
Accepting to externalize the problem of stellar atmosphere (which we
will treat in more detail <a href="./notes-lecture-radTrans.html">in a future lecture</a>), we already have
written one of the missing outer boundary conditions at mass
coordinate equal to the total mass of the star (m=M):
</p>
<ul class="org-ul">
<li>T(m=M) = T<sub>eff</sub>, where T<sub>eff</sub> is <i>defined</i> as the temperature of a black
body producing the same luminosity as the star: L=4π R<sup>2</sup>σ T<sub>eff</sub>
with R=r(M).</li>
</ul>
<p>
For the final missing boundary condition, we need to specify the outer
pressure P(r=R) at R=r(M). This typically comes from imposing a <i>smooth</i>
transition from the stellar interior (inside, T ≥ T<sub>eff</sub>) and the
stellar atmosphere (outside, T ≤ T<sub>edd</sub>), which requires calculating the
pressure in the stellar atmosphere where the assumptions we made so
far, and consequently the equations we wrote do not hold.
</p>
<ul class="org-ul">
<li><b>Q</b>: Why should P be smooth? (<b>Hint:</b> think of dP/dr!)</li>
</ul>
</div>
</div>
</div>
<div id="outline-container-orgee55e5d" class="outline-3">
<h3 id="orgee55e5d"><a href="#orgee55e5d">Solving strategies</a></h3>
<div class="outline-text-3" id="text-orgee55e5d">
<p>
Let's assume we have a way to specify P(r=R), and postpone the
discussion of stellar atmosphere to a <a href="./notes-lecture-radTrans.html">future lecture</a>.
We then have 4 coupled non-linear ODEs with 4 boundary conditions (2
at the center and 2 at the surface), and the EOS as a closure
condition. How can we solve them numerically?
</p>
</div>
<div id="outline-container-orga3c8797" class="outline-4">
<h4 id="orga3c8797"><a href="#orga3c8797">Discretization</a></h4>
<div class="outline-text-4" id="text-orga3c8797">
<p>
First, to represent the system of equations in a computer we want to
discretize them, that is convert every derivative into a <i>finite
difference</i>. This can be done in various ways, each with specific
advantages and disadvantages. The simplest is a first-order forward
discretization of the form:
</p>
<div class="latex" id="org4d4b005">
\begin{equation}
\frac{df}{dx} \rightarrow \frac{f(x_{k+1})-f(x_{k})}{x_{k+1}-x_{k}} \ \ ,
\end{equation}
</div>
<p>
where f and x are a generic function and variable, and the index k
labels the discretized points.
</p>
<p>
For example, if x=m, the index k will label which cell of the
mass-coordinate "mesh" we are considering and m<sub>k+1</sub> - m<sub>k</sub> = Δ m<sub>k</sub> is the
local resolution in mass at location k, while if x=t, then k will
label the "timestep" we are considering, and t<sub>k+1</sub>-t<sub>k</sub> = Δ t<sub>k+1</sub> is the
timestep size.
</p>
<p>
<b>N.B.:</b> typically both the spatial resolution Δ m<sub>k</sub> and the temporal
resolution Δ t<sub>k</sub> are <i>adaptive</i>, meaning the stellar evolution code will
put more mesh points where / take more timesteps when quantities vary
more rapidly. This is necessary to deal with the large dynamic range
of multiple quantities.
</p>
<p>
<b>N.B.:</b> Nature ≫ numerical models, because nature does not need to do
this!
</p>
<figure id="org458547d">
<img src="./images/mesh.png" alt="mesh.png" width="100%">
<figcaption><span class="figure-number">Figure 4: </span>Schematic representation of the spatial mesh in MESA. Intensive quantities (T, ρ, P) that do not depend on the amount matter are defined at the cell "center", while extensive quantities (m, L) are defined at the cell "boundaries". This is Figure 9 in <a href="https://ui.adsabs.harvard.edu/abs/2011ApJS..192....3P/abstract">Paxton et al. 2011</a>.</figcaption>
</figure>
<p>
By discretizing the ODEs we can rewrite them as algebraic equations:
</p>
<div class="latex" id="orgb419e81">
\begin{equation}
\frac{dm}{dr}=4\pi r^2\rho \Leftrightarrow \ln(r_k) = \frac{1}{3}\ln\left( r_{k+1}^{3} +\frac{3}{4\pi}\frac{dm_k}{\rho_k}\right)
\end{equation}
\begin{equation}
\frac{dP}{dr}=-\frac{Gm(r)\rho}{r^2} \Leftrightarrow \frac{P_{k-1} - P_k}{0.5(dm_{k-1} - dm_k)} = - \frac{G m_k}{4\pi r_k^4} \\
\end{equation}
\begin{equation}
\frac{dT}{dr} = - \frac{3}{16\pi a c}\frac{\kappa\rho L}{r^2 T^3} \,\mathrm{or}\, \frac{dT}{dr}\bigg|_\mathrm{ad} \Leftrightarrow \frac{T_{k-1} - T_k}{(dm_{k-1} - dm_k)/2} = -\nabla_{T,k}\left(\frac{dP}{dm}\bigg|_k\right)\frac{\tilde{T_k}}{\tilde{P_k}} \\
\end{equation}
\begin{equation}
\frac{dL}{dr}=4\pi r^{2}\rho (\varepsilon_\mathrm{nuc} - \varepsilon_{\nu} +\varepsilon_\mathrm{grav}) \Leftrightarrow L_k-L_{k+1} = dm_k (\varepsilon_\mathrm{nuc}-\varepsilon_\nu + \varepsilon_\mathrm{grav})
\end{equation}
\begin{equation}
P\equiv P(\rho,\mu,T) \Leftrightarrow P\equiv P(\rho,\mu,T)
\end{equation}
\begin{equation}
\frac{dX_i}{dt}\bigg|_{r} = \left[ \sum\limits_j \mathcal{P}_{j,i}(T,\rho)
-\sum\limits_k \mathcal{D}_{i,k}(T,\rho)\right] + \left( D_i
\nabla^2X_i\phantom{\bigg|} \right)\\
\Updownarrow \\
X_{i,k}(t_n+\Delta t_{n+1}) = X_{i,k}(t_n) + \Delta
t_{n+1}\left(\frac{dX_{i,k}}{dt}\right)_\mathrm{nuc} +
\frac{\left(X_{i,k}-X_{i,k-1} \right)D_{k}\Delta t_{n+1}}{0.5(dm_{k-1} - dm_k)}
\end{equation}
</div>
<p>
Where, when going from the physical model to the
numerical implementation:
</p>
<ul class="org-ul">
<li>we use m rather than r as independent coordinate <a href="./notes-lecture-HSE.html#org3e45872">as expected</a>;</li>
<li>the mass continuity if re-formulated in terms of natural logarithm
or r instead or r itself (to keep numbers smaller);</li>
<li>the EOS, κ, ε<sub>nuc,</sub>, and ε<sub>ν</sub> represent input physics that we borrow
from other fields, and typically are <i>tabulated</i> as a function of T,
ρ or equivalent pairs of variables.</li>
</ul>
<p>
Now that we have transformed a set of ODEs into a set of algebraic
equations, the stellar structure and evolution problem is reduced to
solving a matrix equation:
</p>
<figure id="orgcbf6569">
<img src="./images/matrix_mesa.png" alt="matrix_mesa.png" width="100%">
<figcaption><span class="figure-number">Figure 5: </span>Section of an example of the matrix to solve. Black dots are non-zero entries, meaning the variables are coupled by an equation. Vertical dashed lines denote blocks of the matrix for cell <code>k-1</code>, <code>k</code>, and <code>k+1</code> respectively, red dotted lines separate structural variables and compositional variables. This is Fig. 47 of <a href="https://ui.adsabs.harvard.edu/abs/2013ApJS..208....4P/abstract">Paxton et al. 2013</a>.</figcaption>
</figure>
<p>
Writing this in a symbolic compact form we an array of quantities X(t)
of length equal to the number of mesh points for our spatial
discretization times the number of variables, and a matrix A which
represents how all the entries of X are coupled to each other in the
algebraic form of the equations: AX(t<sub>k</sub>) = X(t<sub>k+1</sub>).
</p>
<p>
To solve this we could consider two approaches:
</p>
<ul class="org-ul">
<li><b>Explicit methods</b>: in this case X(t<sub>k+1</sub>) is expressed as a function of
X(t<sub>k</sub>). This however will require resolving the <i>local</i> dynamics (and
be limited by the sound crossing time), making it impossible to perform a
calculation for the <i>evolution</i> of the star</li>
<li><b>Implicit methods:</b> in this case we write a function F(X(t<sub>k</sub>), X(t<sub>k+1</sub>))
= 0 (e.g., F(X(t<sub>k</sub>), X(t<sub>k+1</sub>)) ≡ AX(t<sub>k</sub>) - X(t<sub>k+1</sub>) = 0), and by solving the homogeneous
equation we have obtained we derive X(t<sub>k+1</sub>). Implicit methods are
preferred in stellar evolution as they are not limited by the sound
crossing time.</li>
</ul>
</div>
</div>
<div id="outline-container-org667cef7" class="outline-4">
<h4 id="org667cef7"><a href="#org667cef7">Solver</a></h4>
<div class="outline-text-4" id="text-org667cef7">
<p>
For each of these methods, we still need to specify a numerical
solver. Typically stellar evolution codes rely on generalized first
order Newton Raphson solvers, which find the zeros of a function
starting from an initial guess and then computing its derivative to
create a better second guess:
</p>
<figure id="org10af7ef">
<img src="https://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif" alt="NewtonIteration_Ani.gif" width="100%">
<figcaption><span class="figure-number">Figure 6: </span>Animation of the Newton-Raphson method to find the zero of a function y=f(x) in one dimension through tangent approximations (from <a href="https://en.wikipedia.org/wiki/Newton's_method">Wikipedia</a>).</figcaption>
</figure>
</div>
</div>
<div id="outline-container-org22a395b" class="outline-4">
<h4 id="org22a395b"><a href="#org22a395b">Shooting method vs. Henyey method</a></h4>
<div class="outline-text-4" id="text-org22a395b">
<p>
There is a peculiarity in the way we built the system of equations to
be solved: we have two boundary conditions at one edge of the (one dimensional) domain
(the center), and one set at the other edge (the surface).
</p>
<p>
Historically, the first computations of stellar models would use a
"shooting method": integrate numerically from the center outward, and
verify how far from the outer boundary conditions the solution lands,
based on the distance, calculate a correction on the initial guess,
and iterate. This was computationally expensive, and would require
many iterations before the solution build "inside out" would be
sufficiently good.
</p>
<p>
In the late 1960s, <a href="https://en.wikipedia.org/wiki/Louis_G._Henyey">L. Henyey</a> designed an alternative solving strategy,
which was rapidly adopted by others and is still widely used today.
The key idea was to not start computing the variables inside out, but
instead, make an initial guess for the entire domain (typically the
initial guess can be the solution for the previous timestep!), and
calculate a correction not based on the distance between the solution
and the outer boundary condition only, but based on both boundaries at
the same time!
</p>
<p>
<b>N.B.:</b> for a biased overview of the evolution of the field of
computational stellar physics by one of the "founding fathers" of the
field, see this <a href="https://www.aip.org/history-programs/niels-bohr-library/oral-histories/5091">1978 interview to R. Kippenhahn</a>.
</p>
</div>
</div>
<div id="outline-container-org86a225d" class="outline-4">
<h4 id="org86a225d"><a href="#org86a225d">Evolution splitting</a></h4>
<div class="outline-text-4" id="text-org86a225d">
<p>
One final peculiarity of the stellar structure and evolution problem
is that the evolution is <i>slow</i> and driven by nuclear physics, which
provides some of the hardest equations to solve, because they are
extremely <i>stiff</i>. Mathematically, this means the eigenvalues of the
nuclear physics part of the problem have very different norms, making
it challenging to design numerical algorithm that can find all the
eigenvalues efficiently. Physically, this is because many reactions
have extreme temperature sensitivities (because of the Coulomb
barriers and the tunneling probability).
</p>
<p>
The clear separation of timescales for the structure and evolution
(τ<sub>ff</sub> ≪ τ<sub>KH </sub>≪ τ<sub>nuc</sub>) allows for one more trick, so called
"evolution splitting": one can compute the <i>structure</i> of the star
assuming fixed composition, then <i>evolve</i> the composition with that
fixed structure (meaning T, and ρ), and iterate computing a new
structure with the updated composition. This allows to separate the
problem and can make it more tractable.
</p>
</div>
</div>
</div>
<div id="outline-container-org625c91c" class="outline-3">
<h3 id="org625c91c"><a href="#org625c91c">Open science: <a href="https://docs.mesastar.org/en/latest/">MESA</a> and <code>MESA-web</code></a></h3>
<div class="outline-text-3" id="text-org625c91c">
<p>
"<i>An algorithm must be seen to be believed</i>" - D. Knuth
</p>
<p>
The discussion above is an attempt to highlight the general principles
and technical aspects specific of computational stellar structure and
evolution without narrowing down to a specific code. Historically,
many groups independently developed computer codes to solve the
stellar structure and evolution problem, and only few universities
would have the computing power to run them. In the 1960s-1980s, even
if the group leaders may have been willing to share their codes upon
request, all the algorithmic details would typically be hidden from
the scientific community, and considered as technicalities. Many of
the people actually <i>writing</i> the code and <i>implementing the algorithms</i>
were women, and their contribution has been under-appreciated because
of this attitude of considering algorithms as technicalities.
</p>
<p>
Still today, the algorithmic details of many codes are only known to
the inner circle of the group owning the code. In 2011, an