# A Survey of Mathematical Models on Somitogenesis

Hanyu Song
Brandeis University

(September 09, 2019)

###### Abstract

This paper presents a comprehensive survey of various established mathematical
models pertaining to Somitogenesis, a biological process. The study begins by
revisiting and replicating the findings from prominent research papers in this
domain, subsequently offering a critical evaluation of the strengths and
weaknesses inherent in each approach. By synthesizing this knowledge, the
paper aims to contribute to a deeper understanding of Somitogenesis, and pave
the way for further advancements in the development of enhanced mathematical
models for this intricate biological process. The concluding section offers
valuable insights and directions for prospective research in this field.

###### Contents

  1. Introduction to Somitogenesis
    1. Unsolved questions
  2. Clock and Wavefront Model
    1. Summary
    2. Mathematical Equations
    3. Analysis
  3. Oscillatory-based Model
    1. Summary of the PORD Model
    2. Mathematical Equation
    3. Analysis
  4. Excitable Model
    1. Summary of the one-dimensional RD Model
    2. FhN-type system and excitability
    3. Analysis
  5. Conclusion
  6. A Code
    1. Clock and Wavefront
    2. Nagahara Discrete

## Introduction to Somitogenesis

Somites are blocks of cells that lie along the anterior-posterior (AP)
vertebrate embryonic axis of the developing embryo.

Somitogenesis is the process by which somites form by segmenting the axis into
similar morphological units such as vertebrates etc… Somitogenesis serves as
the key biological process in the embryo since it’s responsible for segmenting
the vertebrate axis and generating the prepattern that guides the formation of
the tendons, ribs, muscle, and other associated features of the body trunk.
Figure 1 illustrates the form of somites in an embryo and how segmentation
works in the AP axis.

Figure 1: Embryonic somite and the AP axis. The left picture is a human
embryo[8], where the somites are already in shape when the embryo is still
very immature. The right picture is an anterior-posterior axis where somites
segment from the PSM cell that lies on both sides of the AP axis[3]. Front the
posterior end to the anterior, the cells transform from undetermined, to
determined, then finally to somites.

Although many details about somitogenesis are still debated, there are some
scientific facts that serve as the fundamentals for further research: Somites
segment from the presomitic mesoderm (PSM): thick bands of tissue that lie on
either side of the AP axis. The segmentation begins with the establishment of
a prepattern of gene expression, and it is characterized by periodic
activations in regions where future somites will segment. Early scanning
microscope images show that the posterior PSM displays a series of cells
similar in size and structure, known as somitomeres, which seem to be the
precursors of the somites[1]. The existence of this prepattern was confirmed
by microsurgical experiments in which isolated parts of the PSM formed somites
in strict isolation[7]. Figure 2 demonstrates the wave-like gene expression in
a mouse embryo.

Figure 2: Gene expression in a mouse embryo, where the gene is marked
green[2]. The gene has a wave-like propagation from the posterior end of the
PSM to the anterior side.

Another fact about the PSM is that the PSM is not a homogeneous tissue[13].
This is supported by microsurgical experiments conducted by Dubrulle and co-
workers: AP inversions of somite-length regions of the posterior PSM resulted
in normal segmentation whilst inversions of the anterior PSM resulted in
somites with reversed polarity[10], which suggested that the anterior-most
part of the PSM is determined with regard to its segmentation program, whilst
the posterior-most part of PSM is susceptible in this respect. This proves the
PSM’s heterogeneity, which is a key feature of the models for somitogenesis.

The different regions of the PSM were found to correspond to regions of
varying FGF signaling. fgf8, which is a gradient of FGF8 (Fibroblast Growth
Factor 8). FGF8 is a gene with dynamic expression in the PSM, peaking at the
posterior end of the embryo, whilst decreasing in the direction of the
anterior end[6]. See Figure 3. The function of fgf8 is to down-regulate the
cells, meaning, higher concentration of fgf8 will prevent the segmentation of
PSM, whilst its decrease will make segmentation possible, and when fgf8
decreases past a certain threshold, the cells are then able to segment into
somites. We call that threshold ”the determination front”[10]. The uneven
distribution of fgf8 implies that the positional information of the PSM cell
is crucial. However, the role of positional information is a controversial
issue in mathematical biology, and it’s typically not possible to build robust
biological structures without additional mechanisms, such as diffusion[15].

### Unsolved questions

We know the important information that the down-regulation of fgf8 heavily
affects the somitogenesis process, and we seem to understand the logic behind
the process of somitogenesis, but it is difficult to draw any conclusion about
which specific type of model is capable of accurately recapitulating this
process. There are still many questions that must be determined: whether the
PSM cells are oscillatory or excitable with respect to fgf8 levels? Are the
cells globally controlled by the gene or do they have local interactions
between themselves as well? Does the global fgf8 down-regulation even matter?
What will happen if the fgf8 is kept constant, can a reaction-diffusion model,
that emphasizes local interactions between cells, explain this process
accurately? In the later discussions of this paper, we will look into several
kinds of mathematical models, with each of them having distinct answers to the
above questions. Admittedly, none of them are deemed to be ”perfect”, with
each of them having its own drawbacks. Although there are plenty of models out
there, understanding those models’ mechanisms is still important as it could
accelerate the process of creating a better one in the future.

Although there are plenty of models out there, to date nobody has provided any
comparison of them, and most papers on this topic don’t even reference each
other as they are in different fields: Mathematical Biology, Bio Development,
Physics, etc. Therefore, this paper’s goal is not to compare and find the
perfect model, but to see each of their distinctive advantages and try to
synthesize them if possible, while avoiding their drawbacks when we attempt to
create new models in the future.

Figure 3: PSM cell and gene gradients, FGF included[11]. The PSM elongates
posteriorly as the somites are formed, whilst the gradient of FGF8 always has
a peak at the very posterior end of the PSM cell. As it decreases to a certain
level, the prepattern arises (the blocks that are not fully green), then
somites are formed. What’s contrasting to the FGF gradient is the RA gradient,
which is another type of gene that peaks at the anterior end of the PSM, but
it’s not that relevant to the overall process compared to FGF8.

## Clock and Wavefront Model

### Summary

One of the most famous and widely studied models is the clock and wavefront (C
& W) model. As its name implies, the model proposes the existence of a
segmentation clock and a wavefront of FGF8 along the AP axis of vertebrate
embryos. This idea was first proposed in 1975 by Cooke and Zeeman, with the
gist that there is a longitudinal, global positional information, which is the
above-mentioned FGF8 gradient distribution, that interacts with a smooth
cellular oscillator, which is the so-called clock, to govern the time for the
PSM cells to segment and develop into somites. This idea was then revised by
Pourquie and co-workers, where they went into more specifics and proposed that
the clock sets the times at which new somite boundaries form whilst the
position of the determination front sets where they form[10]. For a cell at a
particular point, they assume that competence to segment is only achieved once
FGF8 signaling has decreased below a certain threshold, the position of which
is known as the determination front.

Figure 4: Representation of the vertebrate body plan during somite
formation[4]. The top part of the diagram shows the FGF8 wavefront, with a
peak in the posterior and a decrease in the direction of the anterior. When
the FGF8 decreases to a certain level, it reaches the determination front. The
middle section of the diagram shows the AP axis of the embryo with the somites
(dark grey blocks), determined region (light grey blocks), and the
undetermined region (light grey band) clearly marked. The bottom is a
visualization of a segmentation clock, which shows the time needed for cells
to gain segmentation.

Therefore, the whole somitogenesis process, according to this model, is
divided and analyzed into different parts. Before reaching a determination
front, a cell will gain the ability to segment by being able to produce a
”somitic factor”, which could be several possibilities of genes. One clock
oscillation after reaching the determination front, cells become able to
produce the ”signaling molecule”. After a cell is able to produce somitic
factor and respond to the signaling molecule, it is specified as somitic and
becomes refractory to FGF8 signaling[16]

### Mathematical Equations

C & W mathematical equations were first proposed by Collier et al. (2000) and
were developed by Mclnerny et al. in 2004, then by Baker and colleagues in
2006. One of the most important features of this model is that local
mechanisms, controlled by time points and positional information, will trigger
segmentation, which fits the C & W assumption perfectly. After segmentation,
cells adhere to each other, creating distinctive somites. When creating this
model, Collier made some further assumptions[12]: (1) The AP-axis can be seen
as fixed with respect to the cells. The PSM’s length is constant and the
segmentation pattern progresses with a constant speed. (In reality, the
posterior end is actually elongating.) (2) The signals that are emitted by
specified cells when they reach certain points are like pulses. The signaling
molecule disperses fast and diffuses rapidly. This is the key assumption since
rapid diffusion can ensure that only cells that are in certain positions will
respond to the signal, if not, all cells will segment at the same time. (3)
Somites are formed continually, and the beginning or end of this process is
not considered, which means it’s assumed that signals emitted from cells exist
all the time.

This model can be well explained by Figure 5. In this diagram, x denotes the
distance while t denotes time. There are two key components in this model:
u(x,t) and v(x,t), in which u(x,t) represents the degree of concentration of
somitic factor a cell is exposed to at a given x and t, while v(x,t)
represents the diffusive signaling molecule. A cell that has a high
concentration of u is specified as somitic, while those with a low
concentration of u are non-somitic.

Figure 5: Representation of C & W model illustrating the two time points P1
and P2 and the three key stages within the model.[12]. Cells at the posterior
end of the PSM (Region I) are less mature than those in other regions since
the somitogenesis process starts from the anterior part. As cells become more
mature in Region II, they become capable of responding to the signaling
molecule v, emitted by cells at point P2. In Region III they begin to form
somites and are no longer able to emit any signals.

The mathematical equations proposed by Collier and colleagues are also based
on these two components [12]:

$\displaystyle\partial_{t}u(x,t)=\frac{(u+\mu v)^{2}}{\gamma+\kappa
u^{2}}\chi_{u}(x,t)-\frac{u}{k}$ (1)
$\displaystyle\partial_{t}v(x,t)=\frac{\chi_{v}(x,t)}{\epsilon+u}-v-D\frac{\partial^{2}v}{\partial
x^{2}}$ (2)

$\chi_{u}$ and $\chi_{v}$ are controlled by two Heaviside step functions:

$\displaystyle\chi_{u}=H(ct-x+x_{1})$ (3)
$\displaystyle\chi_{v}=H(ct-x+x_{2})$ (4)

note that the Heaviside step functions’ rule of calculation is:

$\displaystyle H(x)=\left\\{\begin{aligned} 1&&x\geq 0\\\
0&&x<0\end{aligned}\right.$ (5)

As mentioned above, u and v represent the concentration of ”somitic growth
factor” and ”signaling molecule” respectively, while other variables in these
equations are all positive constants. This model uses a zero flux boundary
condition. It indicates that this boundary condition prevents anything from
leaving this system, which may be an application of the third assumption made
by Collier mentioned above.

Heaviside functions,$\chi_{u}$ and $\chi_{v}$, play an important role in this
model. They can be seen as switches: the elements inside brackets, t, and x,
which are the time and location information, together determine the on and off
of the dynamics in u and v respectively. In Figure 5, the Heaviside functions
are shown along with the regions where the somitic growth factor u is,
respectively, high ($x<x_{2}+ct$) and low ($x>x_{1}+ct$). Somitic growth
factor and signaling molecules boost the somitogenesis process collectively
and they affect each other, as we can see that $\partial_{t}u$ is affected by
v and $\partial_{t}v$ is affected by u. Specifically, u inhibits v while v
activates u.

The model was further expanded by Baker and colleagues in 2006. His team made
some revisions to the two equations above and they added a third equation into
the system: $\frac{\partial w}{\partial t}$ which represents the changing
gradient of FGF8, which down-regulates the somitogenesis process[4]:

$\displaystyle\frac{\partial u}{\partial t}=\frac{(u+\mu
v)^{2}}{\gamma+u^{2}}\chi_{u}-u$ (6) $\displaystyle\frac{\partial v}{\partial
t}=k(\frac{\chi_{v}}{\epsilon+u}-v)+D_{v}\frac{\partial^{2}v}{\partial x^{2}}$
(7) $\displaystyle\frac{\partial w}{\partial t}=\chi_{w}-\eta
w+D_{w}\frac{\partial^{2}w}{\partial x^{2}}$ (8)

and $\chi_{w}=H(x-x_{n}-c_{n}t)$ where $x_{n}$ and $c_{n}$ are constants.
Based on the previous system, this system is reproducing these important
behaviors: (1) the increase of somitic factor u is activated by signaling
molecules and is self-regulating. (2) The somitic factor is an inhibiting
signaling molecule. In other words, signaling molecule is produced rapidly in
areas where somitic factor concentration is low. (3) FGF8 is produced in the
tail and regresses along the x-axis.

### Analysis

The C & W mathematical model proved to be effective in producing a
qualitatively reasonable match to reality. We recapitulated and reproduced
some of the results of the above mathematical equations shown in Mclnerney and
Baker’s papers, and they do support the gist of C & W theory. We first
analyzed the qualitative behaviors of this model, and the result can be
explained with Figure 6, which is derived from Figure 5:

Figure 6: The three stages of somitogenesis from the posterior to the anterior
end of the PSM.

In region I, since the switches $\chi_{u}=\chi_{v}=0,(u,v)\rightarrow(0,0)$,
while in Region II and III, as $\chi_{u}and\chi_{v}$ change with respect to t,
the qualitative behavior of these two regions will be different. Below are the
phase planes of u and v in Region II and III respectively, which is a
reproduction of Fig 2 and Fig 3 in Mclnerney’s paper[12] using XPP:

Figure 7: The phase planes of u and v in Region II (left) and Region III
(right).

In Region II, there are three steady states: two stable equilibriums when u is
close to 0 and 1, and a saddle in the middle. The region II carries the pulse
of the signaling molecule. In the phase plane, we can see that after the cells
pass the determination front and before they finish one clock cycle, the cells
gain the ability to respond to signals.- the somitic factor concentration u is
always above 0, while they can’t generate signaling molecules themselves, as v
remains to be 0. In Region III, after cells undergo one cycle of segmentation
clock, $\chi_{u}=\chi_{v}=1$. There is only one stable equilibrium in the
phase plane, meaning whatever u and v’s initial values are, they will all
arrive at that specific point and their values will remain stable, then cells
are identified as somitic.

We then reproduced the numerical solutions for equations (1) and (2) in
Mclnerney’s Figure 11(a) and (b)[12], shown in figure 8.

Figure 8: Numerical solution given by equation (1) and (2) for $0\leq t\leq
300$. Parameter values:
$\mu=10^{-1},\gamma=0.2,\kappa=10,c=5*10^{-3},\epsilon=10^{-3},D=100$. This
set of parameter values violated one of the conditions mentioned in the paper,
so this solution is not a perfect one since for $\mu=10^{-5}$ and
$\gamma=0.2$, a high level of v fails to activate u production[12].

We also reproduced the numerical solution of the C & W model in one spatial
dimension given by equations (5) (6) (7 )in Baker’s Fig 3[16], using the code
provided in the appendix. Figure 9 contains the numerical solutions for
$u(x,t),v(x,t)$ and $w(x,t)$ respectively:

Figure 9: Numerical solution given by equation (5) (6) and (7), showing the
spatiotemporal dynamics of the somitic factor, (a). the signaling molecule,
(b), and fgf8 (c). The regression of the FGF8 wavefront is accompanied by a
series of pulses in the signaling molecule and coherent rises in the level of
the somitic factor. Parameter values:
$\mu=10^{-4},\gamma=10^{-3},\kappa=10,\epsilon=10^{-3},\eta=1.0,D_{v}=50,D_{w}=20,x_{n}=0,c_{n}=0.5,D=100$.

However, the fact that the verification of the results of the above models, in
some ways, shows the models’ validity, can not prove the models to be
flawless. There are some issues that need to be considered before constructing
a better model. The equation set (1) and (2) is not robust because somites
depend sensitively on plenty of factors, such as mesh, the speed c, and
initial conditions. Any slight interference in those factors will prevent them
from obtaining successful results. For equations (5), (6), and (7), although
the results in Fig 9 show a clear and consistent pattern of pulses of somitic
factor and signaling molecule, they rely heavily on a very smooth gradient w.
Admittedly, in normal cases, the idea that u and v rely on a smooth gene
gradient is not, in itself, problematic. However, in this model, it’s simply
assumed that a generic FGF8 molecule makes up the gradient controlling the
position of the determination front[16], which means although the gene’s name
is FGF8, it in fact represents the aggregate influence of multiple genes that
may affect the somitogenesis process. In other words, the gradient is modeled
at a very phenomenological level. Requiring such a gradient to be perfectly
smooth becomes a drawback of this system. A stochastic FGF8 gradient or some
random fgf8 pulse will easily mess up the result. This problem is mentioned
and demonstrated in Fig 3 in Baker’s paper[16]. Also, in this paper, the
position of the determination front is prescribed, yet in reality, it will be
subject to plenty of factors such as the gradient slope, etc…

## Oscillatory-based Model

### Summary of the PORD Model

In the C & W model mentioned above, the key is that long-range molecular
gradients control the movement of the front and therefore the placement of the
stripes in the embryo. In this section, we are introducing a fundamentally
different system: the progressive oscillatory reaction-diffusion model, or
PORD model, which does not rely on a global gradient control, but is driven by
short-range interactions.

In the first section of this paper, we have introduced several ”facts” and
”unsolved questions” in the somitogenesis process. Although the oscillatory
model’s mechanism is very different from C & W, they do share a lot of
similarities - it’s just that their interpretations and understandings to
those facts are different. The PORD model admits the existence of the
posterior movement of the determination front, yet it explains that it’s not
controlled by global positional information but by interactions between cells.
In Cotterell’s paper, it’s argued that the PORD model could also explain some
other important features of somitogenesis, such as size regulation, which
previous reaction-diffusion models fail to explain. However, we did find that
controlling the FGF8 gradient such as adding random pulse in the C & W model,
will result in larger somites, which strengthens the argument that the amount
of somitic factor controls the size of somites.

Figure 10: The comparison between C & W model’s and the PORD model’s
mechanisms[5]. The left figure shows C & W model focuses on global gradient
control while the PORD model focuses on short-range interactions between
cells. The bottom figure shows the oscillation of gene expression in the PSM
cells, which are the stripes. Each stripe of gene expression will, in the
future, correspond to a subsequent somite boundary. The right figure is the
comparison of the sensitivity of the stripe positions. The positional accuracy
of the arrest front will be more sensitive to noise if defined by long-range
gradients(top) than if defined by the distance from the last-formed expression
stripe (bottom).

The PORD model argues that there is a molecular patterning process that
sequentially produces stripes of gene expression along the PSM, resulting in
the segmentations in the PSM. Figure 10 shows this mechanism and its
comparison with the C & W model. There are two dynamical systems that are
involved in this process. First, cells of the PSM exhibit oscillations of gene
expression. These oscillations are organized into traveling waves and they are
locally well synchronized: neighboring cells are in very similar phases of the
cycle [5]. Second, these oscillations are arrested in an anterior-to-posterior
progression, which means the position where the oscillations are frozen
travels posteriorly through the PSM, and that position is called the arrest
front. Note that the arrest front is similar but not equivalent to the
determination front mentioned above, and is addressed in the discussion
section in Cotterell’s paper [5].

Despite being locally self-organizing, the PORD model involves both molecular
oscillations in the PSM and a traveling wavefront. Yet, it continues to create
stripes even in the absence of a moving FGF gradient. Thus it does not rely on
positional information along the PSM. In this reaction-diffusion model, the
distance between stripes is defined by the local diffusion of a repressor
molecule, which is secreted from the stripes themselves (See the right of
Figure 10). However, the fact that the model behaves the same with and without
the gradient seems like a potential problem, since the PSM cells have been
studied without a gradient and their behaviors seem to be very different [9].

Overall, the PORD model challenges the existing clock and wavefront models by
providing a fundamentally alternative theory based on locally self-
organization. It could explain somite size scaling and have higher robustness
of somite size regulation. Some of the PORD model’s results also stand the
test in chick embryos, which shows its validity[5].

### Mathematical Equation

The exploration of mathematical equations for the PORD model is refreshing.
Cotterell and colleagues enumerated all possible topologies that are possible
for a gene regulatory network of three genes. Of the 9710 possible networks,
210 produced a multi-stripe pattern for at least one parameter set. Of all the
stalactites in the topological tree containing successful topologies, they
found two versions of the C & W model and several versions of the oscillatory
PORD model.

The simplest design of the oscillatory model is a network that contains only
two nodes ((A) of Figure 11), comprising a cell-autonomous activator (A),
which is itself activated by the FGF signal, and a diffusible repressor (R). A
and R are defined by the following equations:

$\displaystyle\frac{\partial A}{\partial
t}=\Phi(\frac{\kappa_{1}A-\kappa_{2}R+\textit{F}+\beta}{1+\kappa_{1}A+\kappa_{2}R+\textit{F
}+\beta})-\mu A$ (9) $\displaystyle\frac{\partial R}{\partial
t}=\frac{\kappa_{3}A}{1+\kappa_{3}A}-D\nabla^{2}R-\mu R$ (10)

where $\kappa_{1},\kappa_{2},\kappa_{3}$ define the strengths of regulatory
interactions between A and R. D is the diffusion constant for R, $\mu$ is a
fixed decay constant, and F is the regulatory input of the FGF gradient onto
a. $\beta$ is the background regulatory input of A. To prevent negative values
of morphogens, we use the function $\Phi(x)=xH(x)$, where $H(x)$ is the
Heaviside function.

Together they form a reaction-diffusion mechanism where R inhibition is
responsible for the spacing of adjacent stripes. Since the PORD model does not
rely on global positional information, the model does not spontaneously
generate segments everywhere but rather progresses from anterior to posterior,
which is similar to real-world biological phenomena. (B) and (C) in Figure 11
shows the wave of gene propagation and its oscillatory mechanism.

Figure 11: The PORD mechanisms [5]. (A) The minimum somite-patterning circuit
that implements the PORD mechanism. It contains an activator molecule (green)
and a diffusible repressor (red). $\kappa_{1},\kappa_{2},$ and $\kappa_{3}$
are strengths of interactions between A and R. It shows that the system relies
on the FGF to activate but does not need FGF to perpetuate, which fits the
PORD character. (B) Gene expressions are initiated at the posterior end and
they travel to the anterior end. The white stripe is the formed somite, which
is the position where the last gene expression stopped. When the next wave
propagation arrives at a certain distance from the last one, it will stop and
form the next stripe. (C) A snapshot of gene expression oscillation along the
PSM. The oval dashed arrows indicate the oscillation directions. The blue line
is the FGF gradient. It’s not directly related to stripe formation. The Buffer
Region is generated by the diffusion of the repressor from the last formed
stripe, which inhibits oscillations. Therefore, cells cannot go beyond the
Buffer Region and will exit oscillations in the Oscillating Edge Region to
form new stripes. Newly formed stripes act as the next source of repressor to
prevent oscillations. They will form new Buffer Regions and push the arrest
front posteriorly.

### Analysis

The PORD model proves to be a typical oscillatory model. Its wave propagation
theory as well as the mathematical equation both exhibit its oscillatory
nature. We used XPP and recapitulated figure E of Figure 2 in the Cotterell
paper [5], see the left of Figure 12.

Figure 12: XPP analysis of the PORD model. The left figure is the phase
portrait for the non-diffusing case of Equations (8) and (9), which means to
ignore the diffusing term D. The green and red lines are the nullclines for
the activator A and inhibitor R respectively. The right figure is the
bifurcation analysis for the same case. It shows the bifurcation of activator
A with respect to the FGF gradient.

Ignored its diffusion state and made stationary, the system reveals that
oscillations are the natural dynamic state for most cells in the PSM. The
bifurcation analysis (the right of Figure 12) also reveals its oscillatory
nature. The activator has a Hopf bifurcation which starts when the fgf8 drops
to a certain level. Meaning, when the fgf8 is high, the activator will be
stimulated. When fgf8 decreases to a certain level, the activator will
interact with the inhibitor and will start to oscillate. The part in between
the two green boundaries is where oscillation exists. And the oscillation will
stop once the cells reach the arrest front, which in this case, is when the
fgf8 decreases to 0.

However, the PORD model has received some criticisms. For one, although in the
paper, it’s been claimed multiple times that this model does not require the
moving FGF gradient, it nevertheless acts to couple the rate of embryo growth
with the integral levels of FGF8 signaling in the PSM [5]. Meaning, it can’t
ignore the fact that FGF8 plays an important role in controlling the size of
somites. Higher levels of FGF signaling will result in smaller somites. Also,
from the bifurcation analysis shown above, we can see that the character of
the PORD system is somewhat similar to C & W system in that FGF8 gradient
information could control the cell activities in both cases. The position
where the activator starts to oscillate can be seen as the determination front
in C & W model, and the new terminology ”arrest front” in the PORD system also
locates close to where the FGF8 gradient drops to a very low level. Simply
put, although the PORD system has created some new terms such as the ”buffer
region” and ”the arrest front”, its activities, similar to C & W model, could
still be explained by the FGF8 gradient control. Another problem is that, not
only in PORD model but also in some other oscillatory-based models, many of
them control the specific position of spatial stripes ”manually”, by defining
thresholds or piecewise functions [14]. Although this may help create
beautiful results, it’s in contrast to the principle of the self-organization
of biological systems.

Nonetheless, the PORD model, as one of the most famous oscillatory models for
somitogenesis, does present a very different perspective. It reveals the
possibility that cells themselves carry an oscillatory nature in the absence
of diffusion. It also produces several nice movies to show the oscillation
process clearly. We tried to reproduce that in Matlab but didn’t succeed.

## Excitable Model

### Summary of the one-dimensional RD Model

Both the C & W model and the PORD model have an important feature that has not
been mentioned in previous sections. That is, both models set spatial
continuity as a key requirement. Spatial continuity, in this case, means that
both models ignore the size of cells and see the PSM as a whole unit, or, a
spatial continuum. However, although spatial continuity is acceptable in most
chemical reaction systems, Nagahara and colleagues argue that it’s not always
the case in biological systems. It’s simply because cells in a multicellular
organism have a finite size[14]. In the initial stages of an organism’s
developmental process, when important biological structures first emerge, the
number of cells is usually small and the size of a single cell can not be
simply ignored since the size of the field where the phenomena occur is
comparable to that of a cell. As we know spatial continuity is not always met,
we may have to consider the spatial variations between cells. Yet, Nagahara
and colleagues propose that instead of thinking in that way, it’s better to
just treat cells as ”interacting discrete nodes” in a network. Since the
diffusion inside a cell is much faster than that in a membrane, treating cells
as individual nodes that compose a huge network is suitable.

Nagahara also criticized the complicity and difficulty in other models, since
models based on a continuum will have difficulty producing a narrow boundary
between different distinctive behaviors, as sharp as two or three cells. The
fact that most models assume the existence of two or more different
interactions, such as activator, inhibitor, etc., among the neighboring cells,
is also complicated. Therefore, they created a simple, one-dimensional
reaction-diffusion model that focuses on three things: (1) No diffusion among
inhibitors, (2) cells are discrete instead of continuous, and (3) spatial
inhomogeneity, where (1), (2) are new ideas, while (3) is an old one. See
Figure 13.

Figure 13: An illustration of the one-dimensional reaction-diffusion system
proposed by Nagahara and colleagues [14]. The major difference between this
one and the PORD model is that the diffusion of the inhibitor v is excluded.
There’s only one interaction between neighboring cells, which is the activator
u.

### FhN-type system and excitability

Below is a hypothetical model that describes the above features of gene
expression in somitogenesis:

$\displaystyle\frac{\partial u}{\partial t}=f(u,v)+D\mathcal{L}u$ (11)
$\displaystyle\frac{\partial v}{\partial t}=g(u,v)$ (12)

where u and v are concentrations of the activator and inhibitor respectively.
$D\mathcal{L}u$ is the diffusion term for the activator while f and g are the
reaction terms of the activator and the inhibitor. f and g are given as
follows:

$\displaystyle f(u,v)=\frac{1}{\tau_{1}}(\frac{1}{\gamma}u(u-a)(1-u)-v+\beta)$
(13) $\displaystyle g(u,v)=\frac{1}{\tau_{2}}(u-v)$ (14)

where $\tau_{1}$ and $\tau_{2}$ represent the time scales of local reaction
kinetics of u and v respectively. $\gamma$ represents the spatial gradient and
the temporal change in the concentration of a certain substance [14], which is
dependent on space x and time t. In reality,$\gamma$’s biological counterpart,
in the case of somitogenesis, could be the FGF8 gradient in the PSM. $\gamma$
thus plays an important role in this model.

We call this model the ”FhN-type” model. The reason for that is because the
reaction terms (12) and (13) resemble the ”Fitzhugh-Nagumo” model. It’s named
after Richard Fitzhugh who suggested the model in 1961 and J. Nagumo who
created the equivalent circuit the following year. The Fitzhugh-Nagumo model
is a generic model for excitable systems. Because of its simple two-variable
form and generality, it has been used widely. The Fitzhugh-Nagumo prototype
model has the following form:

$\displaystyle\dot{v}=v-\frac{v^{3}}{3}-w+I_{ext}$ (15)
$\displaystyle\tau\dot{w}=v+a-bw$ (16)

The reason why we say this model represents excitable systems, is because when
$I_{ext}$, the external stimulus, exceeds a certain threshold level, the
system will exhibit a characteristic excursion in the phase plane, before the
variables $v$ and $w$ relax back to the rest values. We say the system is
excited and will be refractory to excitability for a period of time. When
$I_{ext}$ doesn’t exceed that threshold, there will not be an excursion, and
we say the system is not yet excited and remains quiescent or excitable.
Except for the excursion, the phase plane also contains two nullclines. One is
linear and the other is a cubic nullcline, or a sigmoid. The excitability of
the system can be discovered by looking at the spatial relationship between
the two nullclines: the closer the linear nullcline to the peak of the
sigmoid, the more excitable the system is. See Figure 14 for the phase plane.

Figure 14: The phase plane for the prototype Fitzhugh-Nagumo model. $I_{ext}$
= 0.5, a = 0.8, b = 0.7. The blue line is the trajectory of the FhN model in
the phase space. The pink and yellow lines are two nullclines where the pink
line is a cubic nullcline and the yellow line is the linear nullcline.

### Analysis

In this model, $\gamma$ is defined as a linear gradient function:

$\displaystyle\gamma(x)=0.21-0.20x,\quad\quad\quad x\in[0,1]$ (17)

With the given information, we used XPP and created the bifurcation diagrams
to find the model’s qualitative features. See Figure 15. We can see that
$\gamma$ has a hopf bifurcation. When $x\in[0,1]$, $\gamma$ is invariant and
the system is stable, yet when x increases beyond the boundary, the system
tends to become oscillatory. $\alpha$’s bifurcation is a saddle-node
bifurcation. Together with the stable region of $\gamma$’s bifurcation, the
system tends to exhibit a bistable region when $\gamma$ is small. Two stable
regions will coexist and then the system will tend to be oscillatory as
$\gamma$ increases. The right figure is a cusp bifurcation. Note that the cusp
bifurcation has a normal form which resembles equation (12).

Figure 15: The bifurcation diagrams for the one-dimensional RD model. The left
figure is a hopf bifurcation of $\gamma$. The middle figure is a saddle-node
bifurcation of $\alpha$. The right figure is a co-dimension bifurcation
diagram, which is a cusp bifurcation of $\gamma$ and $\alpha$. Cusp
bifurcation is defined as two branches of saddle-node bifurcation curve that
meet tangentially, forming a semi-cubic parabola.

The Nagahara paper utilizes the idea of bistability. By changing the parameter
$\gamma$, which adjusts the amplitude of cubic function, we can vary the local
kinetics from oscillation to bistability [14]. Figure 16 shows the changes in
the nullclines $f(u,v)=0$ and $g(u,v)=0$ by manipulating $\gamma$. As we can
see, the system gains more excitability when $\gamma$ is decreased, since the
decrease in $\gamma$ will make the sigmoid, or the cubic nullcline, elongated
vertically. As the linear nullcline will not change, this elongation will
close the distance between the linear nullcline and the sigmoid, thus making
the system more easily excited. Furthermore, using plotting techniques in
Mathematica, we found that manipulating $\alpha$ and $\beta$ in equation (12)
will also change the excitability, in which the increase in $\alpha$ makes the
system less excitable, while the increase in $\beta$ makes the system more
excitable.

Figure 16: Phase planes for equation (12) and (13)[14]. They correspond to (a)
oscillatory state in big $\gamma$ case and (b) bistable state in small
$\gamma$ case.

A spatiotemporal diagram is shown in Figure 17, where (a) $\gamma$ is given as
invariant. A single pulse triggered from the left boundary propagated to the
right and generated a stationary band at a specific position. (b) If we take
into account the posterior growth of the PSM, $\gamma$ will be a function of
both space and time. The pulses will lead to a static, periodic structure. (c)
$\gamma$ sets to decrease as the wave passes, while we ignore the growth of
the PSM. The pulses can also create a static, periodic structure, but the
bandwidth will be much thicker.

Figure 17: The numerical simulations for the one-dimensional RD model,
illustrate the manner of wave propagation, depending on the geometrical
distribution of $\gamma$ [14].

However, the Nagahara model is not a typical RD model. The structures in
Figure 17 will not form unless the model is discrete, which is in opposition
to continuous. Spatial discreteness is an important feature, and also the base
for this model. Normally in the continuous case, the wave propagation
triggered from the left will not stop and generate a stationary band, rather,
it will propagate across the field without stopping[14]. The paper proposes
that, when D is small enough, the propagation of the wave will be blocked and
there exist stable steady solutions. This phenomenon is called ”wave
propagation failure”. We used Matlab and simulated what will happen to the
system as D varies, and we obtained the result shown in Figure 18, using the
code provided in the appendix.

Figure 18: The spatiotemporal diagram for the system as D varies. D decreases
from (a) to (c). In (a), $D=10^{-3}$. The system is not yet discrete and the
wave propagates through the field non-stop. In (b), $D=10^{-4}$. The system
becomes discrete, and the wave propagation leaves a steady state solution,
just like what’s shown in Fig 17 (a). In (c), $D=10^{-5}$. The system is
discrete but the diffusion is so weak that the signal fails to reach the
bistable region.

Overall, this one-dimensional RD model, which is the FhN-type model, is the
most immature model among the three. Unlike the previous two that have already
been going through in vitro experiments, this model is highly theoretical and
leans more toward physics than biology. However, the simple, fresh idea of
excitability opens a new perspective to see the whole process. Excitability in
this model, as we have found in Mathematica, can vary with respect to
$\gamma,\alpha$, and $\beta$. Depending on three variables seems a bit
capricious but the simple philosophy behind excitability - how easy it is to
cross the threshold - makes it simple to accurately manipulate the
excitability of a system or how to make it ”excited”. We look forward to
exploring more about this feature, and to implementing it in an effort to
effectively understand somitogenesis.

## Conclusion

As they represent several of the mainstream ideas in the field, it is not
surprising that all three models provide abundant insights into the hidden
mechanisms behind somitogenesis. While their core ideas are different in one
way or another, and each of them has its own flaws, each model’s existence
nevertheless greatly improves scientists’ understanding toward this field and
motivates new experiments.

The Clock and Wavefront model proposed a prescribed determination front, which
is determined by the level of FGF8 gradient, that controls the positioning of
somites. It segments the PSM cell into different regions and explains the
somitogenesis process systematically. The result is easy to recapitulate and
it shows desired characters that fit the theory. Its biggest drawback is its
hard-coded outcomes, for example, the model can not explain experiments in
which mutant embryo determination fronts change since its position is
prescribed in the model.

The PORD model proposes a local reaction-diffusion oscillatory mechanism that
can generate stripes of somites without global gradient information. The
system is simple and their video result is very impressive. However, in
contrast to its claims, the PORD model still relies on the FGF8 gradient to
control the size of somites, and many of its mechanisms can still be explained
by fgf8 which seems as if the theory is like another ”perspective” to see the
fgf8’s effects. Many of the oscillatory models control the specific position
of the stripes manually using thresholds or piecewise functions, which
violates its self-organizing nature.

The excitable model is a one-dimensional reaction-diffusion model that
resembles the ”Fitzhugh-Nagumo” model, which is an excitable, generic model
widely applied in fields of physics. Discreteness plays an important role in
the system: the model decreases the diffusion level until cells are no longer
continuous, then blocks the wave propagation generated by excitability,
creating fixed stripes of somites, which is a refreshing idea. This model is
theoretical and has not gone through in vitro experiments, and its result is
not robust as it needs a very fine-tuned diffusion level, which is hard to
achieve in reality, but the idea of an excitable system has latent potential
and much to be exploited.

The C&W and PORD models yield great results that match our desire. However,
their strict conditions often do not meet in real-life biology, and it seems
that this problem can’t simply be solved under the existing frames of
mathematical equations. In my perspective, excitability is what needs to be
studied the most rigorously in order to understand somitogenesis. Dr. Hubaud
and his colleagues’ research [9] which proposed excitability as a general
framework for oscillations in the PSM cell is a great start. While referring
to other models for inspiration, we should boldly explore the direction of
excitable models, instead of sticking to past experiences and techniques on
other models that seem more successful at the moment, since forsaking the
mindset that holds for the moment is a must to create a better one.

## References

  * [1] Martin Hrabě de Angelis Achim Gossler.  Somitogenesis.  Current Topics in Developmental Biology, 38:225 – 287, 1997.
  * [2] Alexander Aulehla, Winfried Wiegraebe, Valerie Baubet, Matthias B. Wahl, Chuxia Deng, Makoto Mark Taketo, Mark B Lewandoski, and Olivier Pourquié.  A B-catenin gradient links the clock and wavefront systems in mouse embryo segmentation.  Nature Cell Biology, 10:186–193, 2008.
  * [3] Ruth E. Baker, Santiago Schnell, and Philip K. Maini.  A clock and wavefront mechanism for somite formation.  Developmental biology, 293 1:116–26, 2006.
  * [4] Maini PK. Baker RE, Schnell S.  A mathematical investigation of a clock and wavefront model for somitogenesis.  Journal of mathematical biology, 52:458 – 482, 2006.
  * [5] James Lloyd Cotterell, Alexandre Robert-Moreno, and Jacob C. Sharpe.  A local, self-organizing reaction-diffusion model can explain somite patterning in embryos.  Cell systems, 1 4:257–69, 2015.
  * [6] Julien Dubrulle and Olivier Pourquié.  Coupling segmentation to axis formation.  Development, 131 23:5783–93, 2004.
  * [7] S.Robert Hilfer Ellen A.G. Chernoff.  Calcium dependence and contraction in somite formation.  In Tissue and Cell, volume 14, pages 435–449. IEEE, Elsevier, 1982\.
  * [8] M.A. Hill.  Hillh5 stage 16 bf04.jpg.  Embryology, 2019.
  * [9] Alexis Hubaud, Ido Regev, L. Mahadevan, and Olivier Pourquié.  Excitable dynamics and yap-dependent mechanical cues drive the segmentation clock.  Cell, 171:668–682.e11, 2017.
  * [10] Olivier Pourquié Julien Dubrulle, Michael J. McGrew.  Fgf signaling controls somite boundary position and regulates segmentation clock control of spatiotemporal hox gene activation.  Cell, 106(2):219 – 232, July 2001.
  * [11] Moises Mallo.  Revisiting the involvement of signaling gradients in somitogenesis.  The FEBS journal, 283 8:0–7, 2016.
  * [12] David McInerney, Santiago Schnell, Ruth E. Baker, and Philip K. Maini.  A mathematical formulation for the cell-cycle model in somitogenesis: analysis, parameter constraints and numerical solutions.  Mathematical medicine and biology : a journal of the IMA, 21 2:85–113, 2004.
  * [13] Stephen Meier.  Development of the chick embryo mesoblast formation of the embryonic axis and establishment of the metameric pattern.  Developmental Biology, 73(1):25 – 45, November 1979.
  * [14] Hiroki Nagahara, Yue Ma, Yoshiko Takenaka, Ryoichiro Kageyama, and Kenichi Yoshikawa.  Spatiotemporal pattern in somitogenesis: A non-turing scenario with wave propagation.  Phys. Rev. E, 80:021906, Aug 2009.
  * [15] Cristóbal Quiñinao, Alain Prochiantz, and Jonathan Touboul.  Local homeoprotein diffusion can stabilize boundaries generated by graded positional cues.  The Company of Biologists Ltd, 142(10):1860–1868, 2015.
  * [16] Philip K. Maini Ruth E. Baker, Santiago Schnell.  Mathematical models for somite formation.  Current topics in developmental biology, 81:182 – 203, 2008.

##  Appendix A Code

### Clock and Wavefront

⬇

function newclockwaves

m = 0; % a slab system (in terms of slab,cylindrical,spherical)

x = -5:0.01:15; % setting up AP axis

t = 0:0.05:25; % setting up time axis

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);

u = sol(:,:,1);

v = sol(:,:,2);

w = sol(:,:,3);

figure

imagesc(x,flipud(t),u);

colormap(gray)

colorbar

set(gca,’YDir’,’normal’)

title(’(a) Somitic Factor’)

xlabel(’AP axis (x)’)

ylabel(’Time (t)’)

figure

imagesc(x,flipud(t),v);

colormap(gray)

colorbar

set(gca,’YDir’,’normal’)

title(’(b) Signaling Molecule’)

xlabel(’AP axis (x)’)

ylabel(’Time (t)’)

figure

imagesc(x,flipud(t),w);

colormap(gray)

colorbar

set(gca,’YDir’,’normal’)

title(’(c) FGF8’)

xlabel(’AP axis (x)’)

ylabel(’Time (t)’)

% \--------------------------------------------------------------

function [c,f,s] = pdex4pde(x,t,u,DuDx)

x1 = 1;

x2 = 0;

xn = 0;

k = 10;

mu = 0.0001;

epsilon = 0.001;

gamma = 0.001;

Dv = 50;

Dw = 20;

n = 1;

o = 0;

cn = 0.5;

xb = 7.5;

epsilon1 = 0;

%n1=(-cn-sqrt(cn^2+4*n*Dw))/(2*Dw);

%n2=(-cn+sqrt(cn^2+4*n*Dw))/(2*Dw);

c = [1; 1; 1];

f = [0.00001; Dv; Dw] .* DuDx; % these are the diffusion coefficients, take D1
\approx 0 for this model

%w1 = n1/(n*(n1-n2))*exp(n2*(-xn+1));

Xu = heaviside(x1-x+cn.*t);

Xv = heaviside(x2-x+cn.*t);

Xw = heaviside(x-xn-cn.*t);

Xb = heaviside(epsilon1-xb+x)*heaviside(epsilon1+xb-x);

F1 = ((u(1)+mu*u(2)).^2)./(gamma+u(1).^2).*Xu-u(1);

F2 = k.*(Xv./(epsilon+u(1))- u(2));

F3 = Xw+o*Xb-n.*u(3);

s = [F1; F2; F3];

%% Set the initial conditions

function u0 = pdex4ic(x)

Dv = 50;

Dw = 20;

%gamma = 0.001;

k = 10;

xn = 0;

n = 1;

cn = 0.5;

epsilon = 0.001;

epsilon1 = 0;

lam = sqrt(k/Dv);

A = 1/(1+epsilon-epsilon1);

B = A*sign(x)/(2*cosh(lam*10));

%n0 = 1/2*(1+sqrt(1-4*gamma/k));

n1=(-cn-sqrt(cn^2+4*n*Dw))/(2*Dw);

n2=(-cn+sqrt(cn^2+4*n*Dw))/(2*Dw);

w0=heaviside(xn-x)*(n1/(n*(n1-n2))*exp(n2*(x-xn)))+

heaviside(x-xn)*(n2/(n*(n1-n2))*exp(n1*(x-xn))+1/n);

u0 = [heaviside(-x); A*heaviside(-x)+B*cosh(lam*(10-abs(x))); w0];

%% Set the boundary conditions

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)

epsilon = 0.001;

gamma = 0.001;

k = 10;

Dv = 50;

Dw = 20;

n = 1;

%pl = [ul(1).^2./(gamma+ul(1).^2)-ul(1); k.*(1/(epsilon+ul(1))-ul(2)); 0];

%ql = [0.00001; Dv; Dw]; % can’t be set to [0,D], need to use a small value to
approximate 0

%pr = [0; 0; 1-n.*ur(3)];

%qr = [0.00001; Dv; Dw];

pl = [0;0;0];

ql = [1;1;1];

pr = [0;0;0];

qr = [1;1;1];

### Nagahara Discrete

⬇

%% 1D Simulation of spatial savanna model with diffusion in space

tic

%% Set options for plots/movies

close all

fprintf(’\n’);

DO_MOVIE = 1; % i.e. write movie to avi file for playback, else just plots

ONE_D_PLOT = 0; % plot dynamics of a single point over time

END_STATE_PLOT = 0; % 2D plot of the final state of the system

THREE_D_PLOT = 0; % plot two spatial dimensions, time and colour

SPACE_TIME_PLOT=1;

%% Numerical method parameters

L = 1; % working on [0,L]

N = 100; % N+1 grid points

delta = L/N; % spatial discretization parameter -> 0.01 suggested

h = 0.01; % time discretisation parameter

n = 50000; % number of time steps

tau = (n-1)*h; % simulations time domain is [0,tau]

%% Function definitions

gamma_fun = @(x, a, b) a + b.*x;

f_fun = @(u,v,tau1,gamma,alpha,beta) (u.*(u-alpha).*(1-u)./gamma-v+beta)/tau1;

g_fun = @(u,v,tau2) (u-v)/tau2;

%% Model parameters (values from Nagahara et al. (2009), Phys. Rev. E)

tau1 = 0.588;

tau2 = 32.1;

alpha = 0.4;

beta = 0.33;

a = 0.21;

b = -0.2; % b = -0.2 in the paper

D1 = 1*10^(-5); % D too big or small won’t work

D2 = 0; % this is set to zero in Nagahara’s paper

%% Set up the initial distributions of the cover types on the grid

% each row is one time step of the simulation

% the initial condition is a pulse of u near x = 0

u(1,:) = 0.2*rand(1,N+1);

v(1,:) = 0.2*rand(1,N+1);

% Compute the birth and mortality matrices as a function of rainfall

X = 0:delta:L;

gamma_grad = gamma_fun(X,a,b); % compute the rainfall gradient along the
x-axis

%% The numerical scheme

for i = 2:n

% compute convolutions for this time step

progressbar(i,n);

u(i,:) = u(i-1,:) + h*( f_fun(u(i-1,:),v(i-1,:),tau1,gamma_grad,alpha,beta) +
...

D1*([0, u(i-1,1:(end-1))] - 2*u(i-1,:) + [u(i-1,2:end),0])./(delta*delta) );

% NB the zeros reflect that we are using an "open boundary"

v(i,:) = v(i-1,:) + h*( g_fun(u(i-1,:),v(i-1,:),tau2) + ...

D2*([0, v(i-1,1:(end-1))] - 2*v(i-1,:) + [v(i-1,2:end),0])./(delta*delta) );

end

% The following output is useful when trying to discern whether or not

% a solution is stationary in time

fprintf(’\n’);

fprintf(’The maximum changes on the grid for each variable at the last time
step were:\n’);

fprintf([’u: ’,num2str(max(abs(u(n,:)-u(n-1,:)))),’\n’]);

fprintf([’v: ’,num2str(max(abs(v(n,:)-v(n-1,:)))),’\n’]);

toc

%% Visualise the solution...

% either in a movie...

if DO_MOVIE

subplot_vis = struct(’cdata’,[],’colormap’,[]);

w = VideoWriter(’reaction_diffusion_model.avi’);

w.FrameRate = 4;

open(w)

f = figure;

for j = 1:300:n

f.Name = [’Simulation time: t = ’, num2str((j-1)*h)];

ax1 = subplot(2,1,1);

plot(X,u(j,:));

title(’u’);

ax2 = subplot(2,1,2);

plot(X,v(j,:));

title(’v’);

xlim([ax1 ax2],[0 L])

ylim([ax1 ax2],[-0.5 1.5])

writeVideo(w,getframe(gcf))

subplot_vis(j) = getframe(gcf);

end

close(w)

end

%% Space-time plot of dynamics

if SPACE_TIME_PLOT

fST = figure;

fST.Name = ’Evolution over time’;

subplot(2,1,1)

h1 = pcolor(u);

shading interp

title(’u’);

colorbar

set(h1, ’EdgeColor’, ’none’);

ylabel(’Time’);

caxis([-0.5,1.5])

subplot(2,1,2)

h2 = pcolor(v);

shading interp

title(’v’);

colorbar

set(h2, ’EdgeColor’, ’none’);

caxis([-0.5,1.5])

end


# Letter of Intent: Muonium R&D/Physics Program at the MTA

C. Gatto,5,6 C. Izzo,2 C. J. Johnstone,2 D. M. Kaplan,11footnotemark: 1 3,4 K.
R. Lynch,2 D. C. Mancini,3 A. Mazzacane,2 B. McMorran,7 J. P. Miller,1 J. D.
Phillips,322footnotemark: 2 T. J. Phillips,3 R. D. Reasenberg,8 T. J.
Roberts,3,4 J. Terry3
1Boston U., 2Fermilab, 3Illinois Institute of Technology, 4Muons, Inc., 5INFN
Napoli,  6Northern Illinois U., 7U. Oregon, 8U. California San Diego CASS
11footnotemark: 1 Spokesperson 22footnotemark: 2 Also at Zurich Instruments

(December 5, 2022)

##  1 Introduction

There is a need for a high-efficiency source of muonium (M
$\equiv\mu^{+}e^{-}$, chemically a light isotope of hydrogen), traveling as a
beam in vacuum, for fundamental muon measurements, sensitive searches for
symmetry violation, and precision tests of theory [1]. Currently PSI in
Switzerland is the world leader for such research. With PIP-II Fermilab has
the potential to eclipse PSI and become the new world leader. It is prudent to
begin the R&D now in order to be ready when PIP-II comes online. Fermilab’s
MeV Test Area (MTA) at the 400 MeV H- Linac has a low-energy muon beamline
suitable for this R&D, with the potential to compete with PSI for this physics
in the pre-PIP-II near term as well.

Key muonium measurements include the search for M-$\overline{\rm M}$
conversion, precision measurement of the M atomic spectrum, and the study of
antimatter gravity using M. Furthermore, the J-PARC $g-2$ experiment proposes
to use a low-energy $\mu^{+}$ beam produced by photo-ionizing a slow beam of
muonium, but the needed high-intensity muonium beam has yet to be
demonstrated. The technique we propose may form a suitable muonium source for
such a $g-2$ measurement as well as for other applications of slow muon beams.

M-$\overline{\rm M}$ conversion is a double charged-lepton flavor-violating
(CLFV) reaction, allowed (albeit at an undetectably small rate) via neutrino
mixing. It may be no less likely — and in some models, more likely — than
$\mu$ to $e$ conversion [1]. Thus in a thorough CLFV research program it
should be studied as well as Mu2e. The best current limit, $P_{\rm
M\overline{M}}\leq 8.3\times 10^{-11}$ (90% C.L.) in 0.1 T field [2], was
published over 20 years ago and, given the technical progress since then, is
ripe for reexamination.

As a pure QED bound state of two point-like particles, muonium offers a more
direct test of theory than hydrogen, free of hadronic and finite-size effects.
Precise predictions and measurements have been made of its 1S–2S [3] and
hyperfine [4] splittings (to 4 and 12 ppb, respectively), again over 20 years
ago, and experiments are now under way at PSI and J-PARC to improve them.

Antimuon gravity has never been measured but its measurement now appears
feasible, thanks to the new approach [5] described below. While the weak
equivalence principle (WEP) of general relativity implies that all forms of
matter should act identically in a gravitational field, and precision
measurements supporting it have been made using torsion pendula, the
Earth–Moon–Sun system, and levitated cylinders in earth orbit [6], it has been
argued that the WEP may not hold for antimatter [7]. Moreover, in theories
that assume maximal WEP violation by antimatter (in which the gravitational
acceleration of antimatter on earth, $\overline{g}$, satisfies
$\overline{g}=-g$), major puzzles of cosmology can be resolved with no further
assumptions and no need of the as-yet-unobserved dark matter and dark energy
[8]. Or a 5th force coupling non-universally to leptons, as suggested by muon
$g-2$ and $B$-decay anomalies [9], might cause $\overline{g}$ and $g$ to
differ slightly.

##  2 Approach

Our method (proposed by PSI’s D. Taqqu) relies on the efficient conversion of
positive muons to muonium atoms in superfluid helium [10] combined with the
predicted expulsion of muonium atoms at the superfluid surface due to the
predicted large positive chemical potential of muonium in superfluid helium
[11, 12]. An electric field maintained in the superfluid helium by means of an
electron pool at the surface will drift stopped $\mu^{+}$ to the surface,
where they will combine with electrons to form muonium and be expelled into
the vacuum.

##  3 Apparatus, R&D Tasks, Budget

Figure 1: (left) Conceptual sketch of “Phase III” muonium-gravity apparatus:
degraded surface muons stop in SFHe layer, producing upward-directed M beam,
deflected into horizontal by SFHe-coated reflector [12] and entering 3-grating
interferometer; calibration soft X-ray beam enters from left; (right) muonium
detector concept: fast decay $e^{+}$ detected in barrel SciFi detector
triggered by scintillating-bar barrel hodoscope, slow $e^{-}$ in $xy$
scintillator hodoscope array with SiPM readout after electrostatic
acceleration ($xy$ hodoscope and accelerating rings not shown).

Of the three key muonium measurements described above, muonium gravity is the
one that has never been measured, uses the smallest apparatus, and will
clearly fit within the MTA; the others require more evaluation and may have to
await PIP-II. The key pieces of apparatus for our proposed program are then
(1) a small cryostat that can be operated at $\sim 0.1$ K, within which a pool
of superfluid helium (SFHe) can be created and maintained; the depth of the
pool will be $\sim$ 100 $\mu$m, sufficient to stop an appreciable fraction of
the slow $\mu^{+}$ beam incident from below after its energy has been degraded
by traversal of the cryostat wall; (2) a 3-grating atom interferometer with
$\approx$ 100 nm grating pitch and few-cm grating separation along $z$; and
(3) an efficient detector of both the M decay products and a calibration X-ray
beam; beams enter through thin spots or beryllium windows in the cryostat
walls. The apparatus is shown in Fig. 1 and will be installed within the
sample volume of a dilution refrigerator, operated at $\sim$ 0.1 K to reduce
the SFHe vapor pressure and suppress M-He scattering. Note that the actual
grating separation will be determined after a thorough optimization study;
Fig. 1(left) shows a separation corresponding to 2 muon lifetimes at the
predicted M velocity of 6.3 mm/$\mu$s. Simulations show that the sign of
$\overline{g}$ can be determined to 5$\sigma$ in a few days of running at
$10^{5}$ Mu/s incident on the interferometer; if systematics can be adequately
controlled, a 10% measurement can be made in about a month.

Figure 2: Conceptual sketch of apparatus to test and characterize charge pool
at surface of superfluid helium layer, which creates drift field in helium.
(A), (B), (C): field electrodes; (T): tungsten electron-injection tips.

As a condition to embark on the gravity measurement, we must first demonstrate
that Taqqu’s proposed muonium-beam formation method works as predicted. We are
greatly aided in this by the recent installation of a low-energy $\pi/\mu$
beamline in the Fermilab MTA. Designed to produce 40 MeV/$c$ muons for a muon-
catalyzed fusion ($\mu$CF) experiment, it can easily be tuned for efficient
transport of ($\stackrel{{\scriptstyle<}}{{{}_{\sim}}}$ 29 MeV/$c$) surface
muons, as we have shown in simulation. Our first task, which can proceed as
soon as the $\mu$CF experiment allows, is to verify these simulations by
characterizing the $\mu^{+}$ spectrum experimentally as a function of magnet
currents. Our second task (which can start in parallel with our first) is to
demonstrate sufficient control of the SFHe charge-pool technology to maintain
the needed drift field in the helium (Fig. 2). This can be done using an
existing IIT cryostat, vacuum pumps, and cryocoolers, preferably at Fermilab.
The total scale of the project is a few M$, about half of which is personnel
costs.

##  4 Draft Schedule

We anticipate that this program will require 3–5 years in total. Phase I: Use
G4beamline to simulate the muon beam, optimize its parameters for muonium
production in superfluid helium, and characterize the MTA surface-muon beam
experimentally; test and characterize the electric-field-in-helium apparatus;
design the detectors needed to characterize the produced muonium. Phase II:
Obtain a suitable dilution refrigerator to reach the $<1$ K temperature range
at which the process works best, build the detectors, install the apparatus
within the sample volume, measure muonium production, and determine the
optimal operating parameters and maximum muonium yield. Phase III: Build the
additional apparatus needed to observe and measure the effect of the earth’s
gravity on muonium and carry out a first measurement.

## References

  * [1] A. A. Petrov, R. Conlin, C. Grant, “Studying $\Delta L=2$ Lepton Flavor Violation with Muons,” Universe 8 (2022) 169, https://www.mdpi.com/2218-1997/8/3/169.
  * [2] L. Willmann et al., “New Bounds from a Search for Muonium to Antimuonium Conversion,” Phys. Rev. Lett. 82 (1998) 49.
  * [3] V. Meyer et al., “Measurement of the 1s–2s Energy Interval in Muonium,” Phys. Rev. Lett. 84 (2000) 1136.
  * [4] W. Liu et al., “High Precision Measurements of the Ground State Hyperfine Structure Interval of Muonium and of the Muon Magnetic Moment,” Phys. Rev. Lett. 82 (1999) 711.
  * [5] A. Antognini et al., “Studying Antimatter Gravity with Muonium,” Atoms 6 (2018) 17;
  * [6] T. A. Wagner, S. Schlamminger, J. H. Gundlach, E. G. Adelberger, “Torsion-balance tests of the weak equivalence principle,” Class. Quantum Gravity 29 (2012) 184002; J. G. Williams, S. G. Turyshev, D. H. Boggs, “Progress in Lunar Laser Ranging Tests of Relativistic Gravity,” Phys. Rev. Lett., 93 (2004) 261101; P. Touboul et al. [MICROSCOPE Collaboration], “MICROSCOPE Mission: Final Results of the Test of the Equivalence Principle,” Phys. Rev. Lett. 129 (2022) 121102.
  * [7] M. M. Nieto, T. Goldman, “The Arguments Against ‘Antigravity’ and the Gravitational Acceleration of Antimatter,” Phys. Rep. 205 (1991) 221–281.
  * [8] G. Chardin, “Motivations for Antigravity in General Relativity,” Hyp. Int. 109 (1997), 83; A. Benoit-Lévy, G. Chardin, “Introducing the Dirac–Milne universe,” Astron. & Astrophys. 537 (2012) A78; G. Chardin et al., “MOND-like behavior in the Dirac–Milne universe: Flat rotation curves and mass versus velocity relations in galaxies and clusters,” A&A 652 (2021) A91; G. Chardin, “Experimental and Observational Tests of Antigravity,” arXiv:2210.03445 [astro-ph.CO] (2022).
  * [9] R. Aaij et al. (LHCb Collaboration), “Test of lepton universality in beauty-quark decays,” Nat. Phys. 18 (2022) 277 and references therein.
  * [10] R. Abela et al., “Muonium in liquid helium isotopes,” JETP Lett. 57 (1993) 157.
  * [11] D. Taqqu, “Ultraslow Muonium for a Muon beam of ultra high quality,” Phys. Procedia 17 (2011) 216.
  * [12] V. G. Luppov et al., “Focusing a Beam of Ultracold Spin-Polarized Hydrogen Atoms with a Helium-Film-Coated Quasiparabolic Mirror,” Phys. Rev. Lett. 71 (1993) 2405.


# Possible Implications of Relatively High Levels of Initial 60Fe in Iron
Meteorites for the Non-Carbonaceous – Carbonaceous Meteorite Dichotomy and
Solar Nebula Formation

Alan P. Boss  Earth & Planets Laboratory, Carnegie Institution for Science,
5241 Broad Branch Road, NW, Washington, DC 20015-1305
<EMAIL_ADDRESS>
###### Abstract

Cook et al. (2021) found that iron meteorites have an initial abundance ratio
of the short-lived isotope 60Fe to the stable isotope 56Fe of 60Fe/56Fe $\sim$
$(6.4\pm 2.0)\times 10^{-7}$. This appears to require the injection of live
60Fe from a Type II supernova (SN II) into the presolar molecular cloud core,
as the observed ratio is over a factor of ten times higher than would be
expected to be found in the ambient interstellar medium (ISM) as a result of
galactic chemical evolution. The supernova triggering and injection scenario
offers a ready explanation for an elevated initial 60Fe level, and in addition
provides a physical mechanism for explaining the non-carbonaceous –
carbonaceous (NC-CC) dichotomy of meteorites. The NC-CC scenario hypothesizes
the solar nebula first accreted material that was enriched in supernova-
derived nuclides, and then later accreted material depleted in supernova-
derived nuclides. While the NC-CC dichotomy refers to stable nuclides, not
short-lived isotopes like 60Fe, the SN II triggering hypothesis provides an
explanation for the otherwise unexplained change in nuclides being accreted by
the solar nebula. Three dimensional hydrodynamical models of SN II shock-
triggered collapse show that after triggering collapse of the presolar cloud
core, the shock front sweeps away the local ISM while accelerating the
resulting protostar/disk to a speed of several km/s, sufficient for the
protostar/disk system to encounter within $\sim$ 1 Myr the more distant
regions of a giant molecular cloud complex that might be expected to have a
depleted inventory of supernova-derived nuclides.

hydrodynamics — ISM: clouds — ISM: supernova remnants — planets and
satellites: formation — protoplanetary disks — stars: formation

##  1 Introduction

Laboratory studies of the initial abundances of the short-lived isotope 60Fe
(half-life of 2.62 Myr) have produced values differing by orders of magnitude.
Tachibana et al. (2006) determined an initial ratio of 60Fe/56Fe $\sim
5-10\times 10^{-7}$ in their studies of high Fe/Ni ferromagnesian chondrules
from the ordinary chondrites (OC) Semarkona and Bishunpur. Tang & Dauphas
(2012) performed whole rock analyses of a range of meteorites, including
unequilibrated ordinary chondrites (UOC), and found a much lower initial
abundance ratio, 60Fe/56Fe $\sim 1.15\times 10^{-8}$, a level they considered
to be representative of the general ISM. However, when Mishra & Goswami (2014)
examined seven chondrules from UOC meteorites, they found an initial ratio of
60Fe/56Fe $\sim 7\times 10^{-7}$. Mishra & Chaussidon (2014) studied three
chondrules from the OC Semarkona and the carbonaceous chondrite Efremovka,
finding initial ratios of 60Fe/56Fe in the range of $\sim 2-8\times 10^{-7}$.
Mishra et al. (2016) inferred ratios of $\sim 8-11\times 10^{-7}$ for
chondrules from an UOC, while Telus et al. (2018) found initial ratios in the
range of $\sim 0.5-3\times 10^{-7}$ in other UOC chondrules.

Telus et al. (2016) attributed the differences between whole rock analyses and
in situ studies of individual chondrules to aqueous alteration along chondrule
fracture lines, either on the parent body or prior to recovery on the Earth,
which could skew the whole rock analyses toward lower initial 60Fe/56Fe
ratios. Trappitsch et al. (2018) reanalyzed a Semarkona chondrule with a
different in situ technique, finding a low initial ratio of 60Fe/56Fe
$\sim(3.8\pm 6.9)\times 10^{-8}$ that is consistent with the low end of the
range found by some other in situ analyses and with the low values found by
the whole rock measurements of Tang & Dauphas (2012). The Trappitsch et al.
(2018) results allow a solar system initial ratio as high as $1.07\times
10^{-7}$.

The Tang & Dauphas (2012) result has been accepted by some as evidence of
meteoritic 60Fe being the result of galactic chemical evolution (e.g., Forbes
et al. 2021). Others have been more circumspect. Vescovi et al. (2018) simply
concluded that the initial 60Fe/56Fe ratio lies between $10^{-8}$ and
$10^{-6}$, while Lugaro et al. (2018) decided that the initial 60Fe/56Fe ratio
has not been determined well enough to draw any conclusions about the source
of this short-lived isotope.

Given this checkered past, the Cook et al. (2021) study represents a
potentially transformational approach to determining the initial 60Fe/56Fe
ratio. Rather than search for Fe-Ni-rich phases in the primarily silicate
mineralogy of chondritic meteorites, Cook et al. (2021) analyzed 13 samples
from two groups of magmatic iron meteorites, the common octahedrites (group
IID) and the rare ataxites (group IVB), the former with average to high Ni
contents, the latter with very high Ni. Magmatic iron meteorites offer the
advantage of samples of a well-mixed iron melt, sidestepping the concerns
about whole rock versus in situ measurements that characterize chondritic
meteorites. The study was able to rule out contamination of the Fe and Ni
isotopes by galactic cosmic-rays. Using the common assumption that the iron
meteorites formed from the chondritic reservoir, Cook et al. (2021) found that
an initial ratio of 60Fe/56Fe $\sim(6.4\pm 2.0)\times 10^{-7}$ characterized
their samples of these two groups of iron meteorites, a value considerably
higher than expected as a result of galactic chemical evolution (GCE; e.g.,
Huss et al. 2009). Cook et al. (2021) concluded that the solar system must
have been injected with additional live 60Fe, well above the level to be
expected from GCE, and noted that deriving the injected additional 60Fe from a
SN II also best explained the deficits in 60Ni and 56Fe abundances they
measured in their sample of iron meteorites.

If correct, the Cook et al. (2021) study appears to settle the brouhaha over
the initial abundance of 60Fe in the solar nebula, and thus provides
meteoritical evidence in support of the supernova triggering and injection
hypothesis that was first proposed by Cameron & Truran (1977). Lee et al.
(1976) found evidence of live initial 26Al (half-life of 0.72 Myr) in Ca, Al-
rich refractory inclusions (CAIs) from the Allende meteorite, leading Cameron
& Truran (1977) to hypothesize nucleosynthesis of 26Al in a SN II and rapid
incorporation of the 26Al into CAIs following injection into the presolar
cloud by the SN II shock wave. 26Al is also synthesized during the Wolf-Rayet
(WR) phase of massive stars, prior to a core-collapse supernova explosion, and
mixed into the interstellar medium by the WR star outflow, and so the solar
system level of 26Al has been argued to be a result of GCE (e.g., Young 2016;
Reiter 2020) or of triggering by a WR outflow (Dwarkadas et al. 2017).
However, 60Fe is not produced in significant amounts by WR stars, but it as
well as 26Al is readily produced by SN II (e.g., Tur et al. 2010), making the
initial 60Fe abundance the key meteoritical test for a SN II being the primary
source of the live 60Fe. Protosolar cosmic rays have also been suggested as a
source of 26Al (most recently by Gaches et al. 2020) but such cosmic rays fail
to produce the required 60Fe.

A long series of papers, from Foster & Boss (1996) through to Boss (2019), has
used detailed multidimensional hydrodynamics codes to show that suitable SN II
shock waves are able to trigger the self-gravitational collapse of a molecular
cloud core, while simultaneously injecting the short-lived isotopes into the
collapsing cloud through the Rayleigh-Taylor instability at the shock-cloud
boundary. These papers have proven the viability of the Cameron & Truran
(1977) hypothesis for the origin of the short-lived isotopes inferred from
meteoritical analyses such as those of Cook et al. (2021).

The purpose of this paper is to learn whether this series of shock-triggered
collapse models might have additional implications for explaining the non-
carbonaceous – carbonaceous (NC-CC) dichotomy of meteorites. Warren (2011) was
the first to show that meteoritical Cr, Ti, and O stable isotope abundances
fall into two distinct groups, the NC and CC. Recent work by Nanne et al.
(2019) and Lichtenberg et al. (2021) has developed models for creating and
preserving the NC-CC dichotomy. Both scenarios hypothesize that the solar
nebula first accreted material that was enriched in supernova-derived
nuclides, and then later accreted material depleted in supernova-derived
nuclides. No explanation is offered for a physical mechanism to explain this
significant difference in accreted matter other than isotopic heterogeneity in
the local ISM and an extended period of accretion from the ISM. Thus we seek
here to find a more dynamically-based explanation for the NC-CC dichotomy. We
extend the preferred triggered-collapse model of Boss (2019) to a much larger
calculational volume and use it to track the evolution of matter initially in
four distinct initial regions to learn the effect of shock-wave triggering on
the time evolution of matter accreted by the presolar cloud and nebula.

##  2 Numerical Hydrodynamics Code

The new model presented here was calculated in the same manner as the suite of
shock-triggered collapse and injection models of Boss (2019). The three-
dimensional hydrodynamics code used was once again FLASH 4.3, based on the
algorithms developed by Fryxell et al. (2000). The FLASH codes are adaptive-
mesh refinement (AMR) codes, ideal for following the sharp gradients in
density and temperature associated with shock fronts as they traverse more
uniformly varying regions of the ISM, such as molecular cloud cores. In
addition to the AMR feature of FLASH, the models use a sink particle (as
developed by Federrath et al. 2010) to represent the newly formed, high
density protostar as a point source of gravity able to accrete matter from its
surroundings. As in Boss (2019), the new model uses the FLASH 4.3 multigrid
Poisson solver, and as a result uses a Cartesian coordinate grid with one top
grid block with eight top grid cells in each coordinate direction. The model
begins with a maximum of six levels of refinement on the top grid block.
Because the computational volume is a rectangular cuboid with sides of length
$4\times 10^{17}$ cm in $\hat{x}$ and $\hat{z}$ and $1.12\times 10^{18}$ cm in
$\hat{y}$, the resulting cells are not cubical, being almost three times as
long in the $\hat{y}$ direction as in $\hat{x}$ and $\hat{z}$. [Note that in
Boss (2019), the rectangular cuboid had a significantly smaller length in
$\hat{y}$ of $8.2\times 10^{17}$ cm.] A single increase in the level of
refinement means that each of the three sides of the computational cells are
halved in length, leading to the increased spatial resolution that is needed
to be able to accurately follow the dynamics of shock-triggered protostar
collapse. With six levels of refinement, the smallest cell size in $\hat{x}$
and $\hat{z}$ is $4\times 10^{17}/(2^{5}\times 8)=1.6\times 10^{15}$ cm and a
size 2.8 times larger in $\hat{y}$. When the refinement level is increased to
seven, the smallest cell size is $0.8\times 10^{15}$ cm in $\hat{x}$ and
$\hat{z}$.

##  3 Initial Conditions

The model presented here is identical to one of the models (O) published in
Boss (2019), with two exceptions. First, the length of the calculation box was
extended in the direction of propagation of the shock front in order to learn
more about the effect of the shock front on the ISM gas and dust downstream
from the target molecular cloud core. Second, several new color fields were
defined, in order to better define the fate of the surrounding ISM matter. The
color fields are defined as being initially nonzero in specific regions of the
initial configuration, such as inside the target cloud or inside the shock
front, and these fields thereafter evolve and trace the location and density
of this material as the calculation proceeds. The models in Boss (2019)
focused on the material initially inside the shock front, in order to estimate
the SN II injection efficiency, and did not specify color fields that tracked
the material initially in the target cloud or the surrounding ISM. As a
result, the Boss (2019) models could not separate the evolution of the initial
target cloud material from that of the surrounding ambient ISM, which is the
key goal of this paper.

As in Boss (2019), the initial conditions consist of a stationary target
molecular cloud core that is about to be struck by a planar shock wave (Figure
1). The target cloud core and the surrounding gas are initially isothermal at
10 K, while the shock front and post-shock gas are isothermal at 1000 K. The
target cloud consists of a spherical cloud core with a radius of of 0.053 pc
and a Bonnor-Ebert radial density profile. The central density is chosen to
produce an initial cloud with a mass of 3.04 $M_{\odot}$, embedded in a
background rectangular cuboid of gas with a mass of 1.63 $M_{\odot}$ and with
random noise in the background density distribution. The cloud core is assumed
to be in solid body rotation about the direction of propagation of the shock
wave (the $-\hat{y}$ direction) at an angular frequency $\Omega_{i}=3\times
10^{-14}$ rad s-1. The initial shock wave has a speed of 40 km s-1, a width of
$3\times 10^{-4}$ pc, and a density of $7.2\times 10^{-18}$ g cm-3. These
choices are based on previous modeling, which examined a wide range of shock
parameters (Boss et al. 2010; Boss & Keiser 2010), in order to choose shock
parameters suitable for triggered collapse and injection of SN II-derived
short-lived isotopes.

##  4 Results

Figures 1 though 5 shows the initial configuration and the four different
color fields used to trace the evolution of different regions in the model.
The first color field (denoted mass scalar 1 or ms1, in the FLASH code)
depicts the matter in the target cloud, while ms2 traces the matter in the
shock front, ms3 does the same for the matter behind the shock front, and ms4
follows the ambient ISM surrounding the target cloud core.

The time evolution of this model proceeds exactly the same as the
corresponding model O in Boss (2019): the shock front smacks the top edge of
the target cloud core, leading to a Rayleigh-Taylor instability at the shock-
cloud interface. This instability allows shock front material carrying SN II-
produced short-lived isotopes such as 60Fe and 26Al to be injected into the
target cloud core, which is soon compressed sufficiently by the shock front to
initiate sustained, self-gravitational collapse. Once the collapsing region
exceeds a critical density of $10^{-15}$ g cm-3 (0.048 Myr), a sink cell is
formed at the location of the density maximum, and this sink cell thereafter
accretes the gas and dust in its vicinity, using the same sink cell parameters
as used in model O in Boss (2019). The model started with a maximum of six
levels of refinement, which was increased to seven levels after 0.050 Myr of
evolution. The portions of the shock that do not strike the target cloud exit
off the calculational grid by 0.010 Myr.

Figures 6 through 10 depict the results for the density and four color fields
at the final time calculated for the model of 0.063 Myr. This new AMR model
required a run time of three months on three 32-core nodes of the Carnegie
memex cluster. At the final time, the sink cell is located close to the center
of the density maximum evident in Figure 6 (at $y=-1.32\times 10^{17}$ cm),
has been accelerated by the shock front to a speed of 3.2 km/s in the
direction of the initial shock front, and has acquired a mass of $\sim
0.5M_{\odot}$. The protostar is accreting mass at a rate of $\sim
10^{-6}M_{\odot}$/yr, implying that it will grow to a final mass of $\sim
1M_{\odot}$ in $\sim 0.5$ Myr if that accretion rate could be sustained
indefinitely. However, at the final time, the mass accretion rate is in
decline, as the nearby gas and dust available for accretion is being depleted.

Figure 7 shows that the initial cloud core has been only partially triggered
into self-gravitational collapse, with the majority of the initial mass of
3.04 $M_{\odot}$ having been accelerated beyond the gravitational reach of the
accreting protostar formed by shock wave compression. Figure 8 indicates that,
as in the previous model O (Boss 2019), the shock front material has been
mixed into the same region occupied by the initial cloud core gas and dust,
leading to an injection efficiency into the collapsing protostar essentially
the same as was previously determined, as discussed in detail in Boss (2017).
Figure 9 shows that the low-density, hot post-shock gas follows right behind
the shock front matter with only minimal mixing into the compressed region of
high density gas and dust.

Figure 10 presents the key result of this simulation: the material initially
in the immediate vicinity of the target cloud core has been quite efficiently
swept away from the region of the protostar and of the material that it is
still accreting. Only an insignificant amount of the initial ambient gas and
dust remains close enough to the protostar to be accreted, less than 1/20 (by
mass) compared to the shock front matter at densities above $2\times 10^{-16}$
g cm-3. The model shows that SN II-derived isotopes will be injected
thereafter, with only minimal subsequent accretion of isotopes from the
initial ambient ISM in the immediate vicinity of the target molecular cloud
core, which is assumed to be initially pristine and unpolluted by SN II-shock
ejecta. This sweeping away of the local ISM allows the protostar to experience
two distinct phases of accretion, first due to SN II-shock injection, and
later during its subsequent traverse of the GMC. Figures 11 and 12 provide
close-in views of the model at the final evolution time, showing the newly-
formed protostar/protoplanet disk surrounding the location of the protostar
(i.e., sink cell). Figure 12 shows that only faint wisps of the initial
ambient ISM material still remain in the region of the protostar. Figure 13
shows that the shock wave isotopes have been injected deep within the
protostar and disk system, as a result of the Rayleigh-Taylor fingers that
pierce the target cloud core early in the evolution. Finally, Figure 14
depicts the temperature distribution at the final time, showing that the
protostar is located well behind the high temperature (1000 K) shocked region
and retains its initial temperature of 10 K as a result of molecular line
cooling in optically thin regions (Boss et al. 2008). The protostar thereafter
continues its evolution into the GMC with its ambient temperature of 10 K.

Figure 11 superficially resembles Figure 5 of Boss (2019), which depicts the
results of a model (N) with an initial cloud rotation rate 2/3 of that in the
current model, based on model O and shown in Figure 6 of Boss (2019). In Boss
(2019), model O produced a large-scale disk by 0.121 Myr, with a diameter of
order 1000 au, whereas model N did not produce a large-scale disk by the time
that the model was halted, 0.081 Myr. The present model shown in Figure 11 was
halted at 0.063 Myr, even earlier than model N in Boss (2019), as the result
of an increasingly smaller time step problem that did not afflict model O in
Boss (2019). As the present model is identical to model O, save for the
extension of the numerical grid from $y=-2\times 10^{17}$ cm to $y=-5\times
10^{17}$ cm effected in order to better study the fate of the surrounding ISM
matter, one must ascribe this small time step problem to the computational
cells being increasingly elongated in $\hat{y}$ that had to be employed in
order to use the improved multigrid Poisson solver of FLASH 4.3, as explained
by Boss (2019). Hence one can expect that the present model would produce a
large-scale disk identical to that of model O if the model could be calculated
as far as model O in Figure 6 of Boss (2019), i.e., to 0.121 Myr. A careful
examination of the present model in Figure 11 shows that a disk is beginning
to form, as there are distinct ”horns” of accreting matter both above and
below the location of the central protostar, with the horns above the
protostar being clearly bent outward compared to those of model N in Figure 5
of Boss (2019), indicative of the higher initial rotation rate of the present
model compared to model N and of the subsequent disk formation in model O.
Regardless of the small time step problem, the present model’s extension in
$\hat{y}$ allows Figure 10 to demonstrate that the shock front effectively
clears away the residual ISM matter from the vicinity of the protostar and
disk system, a fact that cannot be gleaned from model O’s Figure 6 in Boss
(2019), where the disk system is about to reach the end of the numerical grid.

The mass of the initial shock front is 0.50 $M_{\odot}$ and it is traveling at
a speed of 40 km s-1, giving it a total momentum of 20 $M_{\odot}$ km s-1. At
the final time, the protostar has a mass of $\approx 0.5M_{\odot}$ and has
been accelerated to 3.2 km s-1 in the direction of the initial shock front,
giving it a total momentum of $\approx 1.6M_{\odot}$ km s-1. Thus the
protostar has acquired less than 1/10 of the total initial momentum of the
shock front. At the rate of 3.2 km s-1, the protostar/disk system will
traverse $\sim$ 3 pc through the background giant molecular cloud complex
(GMC) in 1 Myr, accreting gas and dust from other, more distant regions of the
GMC that have not been polluted recently by a SN II explosion. Given that GMC
diameters range from $\sim$ 5 pc to $\sim$ 200 pc, there would appear to be an
adequate GMC volume for a protostar launched at $\sim$ 3 km/s to scatter off
other protostars in the GMC star-forming regions and accrete further ISM
material and non-SN II derived isotopes.

##  5 Implications for the NC-CC Dichotomy

By this time, the implications of this model should be clear: they provide a
physically reasonable explanation for the NC-CC dichotomy, as advanced by
Nanne et al. (2019, see their Figure 7) and Lichtenberg et al. (2021, see
their Figure 6). These authors hypothesized that the solar nebula initially
accreted material that was enriched in SN II-derived isotopes, and some time
later accreted material depleted in SN II-derived isotopes.

While no time scale for these two accretion phases is presented by either
Nanne et al. (2019) or Lichtenberg et al. (2021) other than the phrases ”early
infall” and ”late infall”, the present model suggests that the ”early infall”
phase lasted for $\sim 0.1$ Myr, based on the final model time reached, while
the ”late infall” phase lasted for $\sim 1$ Myr, based on the speed of the
protostar and GMC sizes. These suggestions are supported by the work on the
NC-CC dichotomy by Kruijer et al. (2020), whose Figure 5 explicitly states
that the ”early infall” occurs at ”t = 0” Myr, implying a time scale much less
than 1 Myr, along with assuming the simultaneous formation of CAIs, while
”late infall” stops by 1 Myr, when Jupiter is supposed to have formed and
prohibited mixing of the inner NC and outer CC reservoirs.

Isotopic heterogeneity throughout the GMC where the solar system was formed is
required, but the present model provides a physical reason for the change in
the SN II-derived isotopic composition of the matter being accreted by the
protosun and solar nebula. Jacquet et al. (2019) make a similar argument for
isotopic heterogeneity in the matter accreted by the solar nebula without
specifying a physical mechanism to explain the origin of the heterogeneity.
Others have argued that physicochemical processing of dust grains could be a
better explanation for certain isotopic anomalies (150Nd) than heterogeneous
infall (Saji et al. 2021).

Hopp et al. (2022) have argued that Fe isotopic abundances in iron meteorites
reflect the same NC-CC dichotomy as other meteorites. They conclude that the
Fe isotope dichotomy can be explained by nuclear statistical equilibrium in
either type Ia SN or in core-collapse SN, i.e., SN II.

Ideally the present calculation could be extended indefinitely further in time
to capture both the ”early infall” and ”late infall” phases advanced by Nanne
et al. (2019) and Lichtenberg et al. (2021). This is not possible for the
present model for two reasons. First, increasingly smaller time steps are
required to push the calculation any farther in time than the 0.063 Myr shown
in Figures 6 through 12, effectively halting this particular model. Second,
even without the time step problem, the protostar and disk system seen in
Figure 12 at $y=-1.3\times 10^{17}$ cm would reach the end of the
computational volume at $y=-5\times 10^{17}$ cm within about another 0.038 My,
traveling at 3.2 km/s, too quickly to transition to the ”late infall” phase.
Note that prior to reaching the bottom of the computational volume, the dense
gas that is accreting onto the central system (orange-red colors in Figure 12)
will have been accreted, as the free fall time for collapse at a gas density
of $10^{-16}$ g cm-3 is 0.0067 Myr, effectively ending the ”early infall”
phase.

If the present scenario is deemed interesting, future work would be needed to
consider a more global model of an entire GMC struck by a SN II shock wave
that follows the progression of the resulting shock-triggered protostar(s) as
they traverse the GMC. In lieu of such an ambitious three-dimensional
hydrodynamics model, we can predict the mass accretion rate that should
characterize the ”late infall” phase. The protostar will accrete gas from the
region of the GMC it traverses at a rate that can be calculated by the Bondi-
Hoyle-Lyttleton (BHL) formula given by Ruffert & Arnett (1994):
$\dot{M}_{BHL}=\pi R_{A}^{2}\rho_{A}v_{A}$, where the accretion radius $R_{A}$
is given by $R_{A}=2GM/v_{A}^{2}$, with $G$ being the gravitational constant,
$\rho_{A}$ being the ambient GMC gas density, and $v_{A}$ is the speed of the
protostar with respect to the ambient gas. With an ambient GMC gas density of
$\rho_{A}=10^{-21}$ g cm-3, typical of GMC average densities and the same as
the background gas in the present model (see Figure 1), $v_{A}=3.2$ km/s, and
a protostar and disk system mass of $1M_{\odot}$, we obtain a mass accretion
rate of $\dot{M}_{BHL}=10^{-4}M_{\odot}$/Myr of ambient GMC gas and dust. If
the protostar system should pass through a dense molecular cloud core with a
mean density of $\sim 10^{-19}$ g cm-3, as assumed in the present model
(Figure 1), then the mass accretion rate would increase to
$\dot{M}_{BHL}\approx 10^{-2}M_{\odot}$/Myr while traversing the molecular
cloud core, which would take about 0.03 Myr at 3.2 km/s. These estimates
suggest that the mass of pre-existing GMC gas and dust that was not enriched
by the triggering SN II shock wave but was accreted by the protostar during a
”late infall” phase lasting $\sim$ 1 Myr would be in the range of
$10^{-4}M_{\odot}$ of ambient matter to $3\times 10^{-4}M_{\odot}$ of dense
cloud core gas and dust. If the ”late infall” phase lasts longer than 1 Myr,
given that a protostar moving at 3.2 km/s travels 3 pc in 1 Myr, but GMCs can
span $\sim$ 5 pc to $\sim 200$ pc in size, the protostar could continue to
accrete ambient GMC matter for longer than 1 Myr. Either way, the total
accreted mass in the ”late infall” phase is likely to be of order
$10^{-4}M_{\odot}$ to $\sim 3\times 10^{-4}M_{\odot}$.

What are the implications of these estimates for the masses of the CC and NC
component? Given that model O in Boss (2019), the model recalculated here with
a longer computational volume, produced a disk with an initial mass of
$0.05M_{\odot}$, representing the ”early infall” or CC component, whereas the
”late infall” or NC component mass accretion estimate is in the range of
$10^{-4}M_{\odot}$ to $\sim 3\times 10^{-4}M_{\odot}$, this implies that the
CC component was about 170 to 500 times as massive as the NC component. While
specific values for the masses of the ”early infall” enriched matter and ”late
infall” depleted matter are not given by Nanne et al. (2019), Lichtenberg et
al. (2021) suggest that their Reservoir I (NC) had a total mass of about a
Earth mass of planetesimals, while their Reservoir II (CC) had about a Jupiter
mass of planetesimals, implying that the ”early infall” CC matter was about
318 times as massive as the ”late infall” NC matter. Given the uncertainties
involved in both the present model and the Lichtenberg et al. (2021)
estimates, having these two estimates agree to within a factor of two is
remarkable and suggests that the scenario proposed herein is worthy of further
scrutiny.

##  6 Conclusions

Observations of star-forming regions have demonstrated that star formation can
be triggered by supernova explosions (e.g., Bialy et al. 2021). The
meteoritical evidence discussed in this paper, coupled with the results of the
new model presented here, along with the previous papers in this series,
appear to provide a reasonable argument that the solar system was formed as a
result of the interaction of a SN II shock wave with a dense molecular cloud
core residing within a GMC complex.

The computations were performed on the Carnegie Institution memex computer
cluster (hpc.carnegiescience.edu) with the support of the Carnegie Scientific
Computing Committee. I thank Floyd Fayton for his invaluable assistance with
the use of memex and Myriam Telus for discussions about the Cook et al. (2021)
paper. The referees provided numerous useful suggestions for improving the
paper. The software used in this work was in large part developed by the DOE-
supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the
University of Chicago.

## References

  * (1)
  * (2) Bialy, S., Zucker, C., Goodman, A., et al. 2021, ApJL, 919, L5
  * (3)
  * (4) Boss, A. P. 2017, ApJ, 844, 113
  * (5)
  * (6) Boss, A. P. 2019, ApJ, 870, 3
  * (7)
  * (8) Boss, A. P., Ipatov, S. I., Keiser, S. A., Myhill, E. A., & Vanhala, H. A. T. 2008, ApJL, 686, L119
  * (9)
  * (10) Boss, A. P., & Keiser, S. A. 2010, ApJL, 717, L1
  * (11)
  * (12) Boss, A. P., Keiser, S. A., Ipatov, S. I., Myhill, E. A., & Vanhala, H. A. T. 2010, ApJ, 708, 1268
  * (13)
  * (14) Cameron, A. G. W., & Truran, J. W. 1977, Icarus, 30, 447
  * (15)
  * (16) Cook, D. L., Meyer, B. S., & Schönbächler, M. 2021, ApJ, 917, 59
  * (17)
  * (18) Dwarkadas, V. V., Dauphas, N., Meyer, B., Boyajian, P., & Bojazi, M. 2017, ApJ, 851, 147
  * (19)
  * (20) Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, ApJ, 713, 269
  * (21)
  * (22) Forbes, J. C., Alves, J., & Lin, D. N. C. 2021, Nature Astron., 5, 1009
  * (23)
  * (24) Foster, P. N., & Boss, A. P. 1996, ApJ, 468, 784
  * (25)
  * (26) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
  * (27)
  * (28) Gaches, B. A. L., Walch, S., Offner, S. S. R., & Münker, C, 2020, ApJ, 898, 79
  * (29)
  * (30) Hopp, T., Dauphas, N., Spitzer, F., Burkhardt, C.. & Kleine, T. 2022, E&PSL, 577, 117245
  * (31)
  * (32) Huss, G. R., Meyer, B. S., Srinivasan, G., Goswami, J. N., & Sahijpal, S. 2009, GeCoA, 73, 4922
  * (33)
  * (34) Jacquet, E., Pignatale, F. C., Chaussidon, M., & Charnoz, S. 2019, ApJ, 884, 32
  * (35)
  * (36) Kruijer, T. S., Kleine, T., & Borg, L. E. 2022, Nature Astronomy, 4, 32
  * (37)
  * (38) Lee T., Papanastassiou, D. A., & Wasserburg, G. J. 1976, Geophys. Res. Lett., 3, 109
  * (39)
  * (40) Lichtenberg, T., Drazkowska, J., Schönbächler, M., Golabek, G. J., & Hands, T. O. 202, Science, 371, 365
  * (41)
  * (42) Lugaro, M., Ott, U., & Kereszturi, A. 2018, Prog. Part. Nuc. Phys., 102, 1
  * (43)
  * (44) Mishra, R. K., & Chaussidon, M. 2014, E&PSL, 398, 90
  * (45)
  * (46) Mishra, R. K., & Goswami, J. N. 2014, GeCoA, 132, 440
  * (47)
  * (48) Mishra, R. K., Marhas, K. K., & Sameer 2016, E&PSL, 436, 71
  * (49)
  * (50) Nanne, J. A. M., Nimmo, F., Cuzzi, J. N., & Kleine, T. 2019, E&PSL, 511, 44
  * (51)
  * (52) Reiter, M. 2020, A&A, 644, L1
  * (53)
  * (54) Ruffert, M., & Arnett, D. 1994, ApJ, 427, 351
  * (55)
  * (56) Saji, N. S., Schiller, M., Holst, J. C., & Bizzarro, M. 2021, ApJL, 919, L8
  * (57)
  * (58) Tachibana, S., Huss, G. R., Kita, N. T., Shimoda, G., & Morishita, Y. 2006, ApJ, 639, L87
  * (59)
  * (60) Tang, H., & Dauphas, N. 2012, E&PSL, 359, 248
  * (61)
  * (62) Telus, M., Huss, G. R., Ogliore, R. C., et al. 2016, GeCoA, 178, 87
  * (63)
  * (64) Telus, M., Huss, G. R., Nagashima, K., Ogliore, R. C., & Tachibana, S. 2018, GeCoA, 221, 342
  * (65)
  * (66) Trappitsch, R., Boehnke, P., Stephan, T., et al. 2018, ApJL, 857, L15
  * (67)
  * (68) Tur, C., Heger, A., & Austin, S. M. 2010, ApJ, 718, 357
  * (69)
  * (70) Vescovi, D., Busso, M., Palmerini, S., et al. 2018, ApJ, 863, 115
  * (71)
  * (72) Warren, P. H. 2011, E&PSL, 311, 93
  * (73)
  * (74) Young, E. D., 2016, ApJ, 826, 129
  * (75)

Figure 1: Initial log density cross-section ($z$ = 0) showing the entire
computational grid with a maximum of six levels of refinement.

Figure 2: Initial cross-section ($z$ = 0) of color field 1, representing the
material initially within the target cloud core, plotted as in Figure 1. The
target cloud is initially stationary.

Figure 3: Initial cross-section ($z$ = 0) of color field 2, representing the
material initially within the shock front, plotted as in Figure 1. The shock
is initially moving downwards at 40 km -1.

Figure 4: Initial cross-section ($z$ = 0) of color field 3, representing the
material initially behind the shock front, plotted as in Figure 1. The
material behind the shock is also moving downwards at 40 km -1.

Figure 5: Initial cross-section ($z$ = 0) of color field 4, representing the
material initially outside the target cloud core, plotted as in Figure 1. This
material is initially stationary.

Figure 6: Final log density cross-section ($z$ = 0) after 0.063 Myr with seven
levels of refinement. By this time, the portions of the shock front that did
not strike the target cloud have exited the bottom of the grid.

Figure 7: Final cross-section ($z$ = 0) of color field 1, representing the
material initially within the target cloud core, plotted as in Figure 6.

Figure 8: Final cross-section ($z$ = 0) of color field 2, representing the
material initially within the shock front, plotted as in Figure 6. Note the
scale change compared to Figure 3.

Figure 9: Final cross-section ($z$ = 0) of color field 3, representing the
material initially behind the shock front, plotted as in Figure 6.

Figure 10: Final cross-section ($z$ = 0) of color field 4, representing the
material initially outside the target cloud core, plotted as in Figure 6. Note
the scale change compared to Figure 5.

Figure 11: Close-up view of the final log density cross-section ($z$ = 0)
after 0.063 Myr as seen in Figure 6. The sink cell is located in the center of
the density maximum and represents the newly-formed protostar and disk system.

Figure 12: Close-up view of the final cross-section ($z$ = 0) of color field
4, representing the material initially outside the target cloud core, plotted
as in Figure 11.

Figure 13: Close-up view of the final cross-section ($z$ = 0) of color field
2, representing the material initially inside the shock front, plotted as in
Figure 11. Note the scale change compared to Figure 3.

Figure 14: Close-up view of the final cross-section ($z$ = 0) of the log
temperature distribution, plotted as in Figure 11.