# Machine-learning the phase diagram of a strongly-interacting Fermi gas

M. Link1, K. Gao1,2, A. Kell1, M. Breyer1, D. Eberz1, B. Rauf1, and M. Köhl1
Physikalisches Institut, University of Bonn, Wegelerstraße 8, 53115 Bonn,
Germany  Department of Physics, Renmin University of China, Beijing 100872,
China

###### Abstract

We determine the phase diagram of strongly correlated fermions in the
crossover from Bose-Einstein condensates of molecules (BEC) to Cooper pairs of
fermions (BCS) utilizing an artificial neural network. By applying advanced
image recognition techniques to the momentum distribution of the fermions, a
quantity which has been widely considered as featureless for providing
information about the condensed state, we measure the critical temperature and
show that it exhibits a maximum on the bosonic side of the crossover.
Additionally, we back-analyze the trained neural network and demonstrate that
it interprets physically relevant quantities.

When an ensemble of attractively interacting fermions is cooled to below a
critical temperature $T_{c}$ it transitions from a normal phase into a
superfluid or superconducting phase. The precise value of the phase transition
temperature is governed by the microscopic details of the system, such as the
interaction strength and interparticle correlations, and can exhibit non-
trivial dependencies. For example, in the crossover from BCS to BEC, it has
been theoretically predicted that the critical temperature depends non-
monotonically on the interaction parameter [1, 2, 3, 4, 5, 6, 7, 8], see
Figure 1a. The non-monotonic behaviour is rooted in the fundamental change of
the nature of pairing below the critical temperature. For Cooper pairing (BCS)
one expects an exponential dependence of $T_{c}$ on the interaction strength
whereas dimer pairing (BEC) implies a nearly constant $T_{c}$. The division
between the two regimes is not at unitarity but is expected to be on the BEC
side of the crossover [9, 10]. In this manuscript we study the critical
temperature across the BCS/BEC crossover using an artificial neural network to
analyze the momentum-distribution of ultracold atomic Fermi gases.

Figure 1:  A Sketch of the phase diagram across the BCS/BEC crossover
including the critical temperature for condensation (solid line). BCS ground
state is dominated by long-range Cooper pairs whereas the BEC exhibits dimer
pairing of the fermions. B Momentum distribution of an ideal Fermi gas at a
temperature of $T/T_{F}=0.15$ (dashed line) in comparison with the momentum
distribution of the BCS ground state for different interaction strengths:
$1/(k_{F}a)=-1$ (blue), $1/(k_{F}a)=0$ (yellow), and $1/(k_{F}a)=1$ (red). C
Principle of the data analysis using a neural network. Condensate fractions
are determined from time-of-flight images after a rapid-ramp. This information
is used to label time-of-flight data for equal parameters but without rapid
ramp. A neural network is trained on the labeled data to predict the
condensate fraction.

A precision determination of the critical temperature across a broad range of
interaction strengths have so far been hindered by insufficient experimental
detection capabilities. One main challenge is that upon release from the trap
in a conventional time-of-flight study, Cooper pairs break and are not
amenable for direct detection. Nevertheless, they leave a weak imprint onto
the momentum distribution of the fermions. In Figure 1b we compare the
momentum distribution of a homogeneous Fermi gas at a temperature of
$T/T_{F}=0.15$ (i.e. near the critical temperature) with the momentum
distributions of BCS ground state wave functions for different interaction
parameters. The pairing signature is by far not as pronounced as the
celebrated bimodal momentum distribution of a Bose-Einstein condensate and
therefore the detection of the condensate fraction is much more difficult.
Additionally, finite temperature, interactions and the inhomogeneity of the
harmonically trapped sample further obscure the pairing signature [11]. In
order to detect the minuscule modifications of the momentum distribution in
the time-of-flight images, we have developed and applied a neural network for
advanced image recognition. We favour neural network processing over standard
data fitting since the neural network is unbiased as compared to applying a
predetermined fitting function and therefore might detect physical signatures
beyond a model-based analysis. Recently, applications of these sophisticated
techniques have entered into the field of quantum physics for the
identification of phases of quantum matter [12, 13, 14, 15, 16, 17, 18].
However, even when being successfully trained, artificial neural networks have
acted as “black boxes” hiding their decision criteria. Specifically, whether
or not the network actually identifies physically relevant criteria for
computing its output has remained obscure. Generally, the interpretation of
neural networks and their causality is rather challenging and currently a
major topic in computer science [19]. In this work, we demonstrate that the
back-analysis of neural networks provides further details of the physics,
which are not accessible by conventional means.

Experimentally, we prepare a quantum gas of $\sim 3\cdot 10^{5}$ atoms per
spin state in the two lowest hyperfine states $\ket{1}$ and $\ket{2}$ of 6Li
in an optical dipole trap, similar to our previous work [20]. We adjust the
interaction strength of the sample by Feshbach resonance and the temperature
by changing the trap [21]. The interaction and temperature are tuned
independently of each other and the thermalized cloud is detected by
absorption imaging after ballistic expansion, see Appendix.

The neural network employed for image analysis comprises of three
convolutional and pooling layers and three fully-connected layers and is
trained through stochastic gradient descent with Adam optimizer [22], see
Appendix. In order to train and validate our neural network, we employ a
supervised learning method [23]. To this end, we measure two different density
distributions after time-of-flight, see Figure 1c: (1) The density
distribution $n_{A}$ of the atoms directly released from the optical dipole
trap. During the expansion, the Cooper pairs are broken and $n_{A}$ is related
to the momentum distribution of the fermions convolved with interaction
effects during the expansion. (2) The density distribution $n_{RR}$ after
applying the rapid ramp technique [24, 25, 26, 27], which measures the
momentum distribution of the molecules that have been created from the Cooper
pairs. Even though this technique is expected to preserve physics in many
cases, quantitatively and principally there are open questions about the
adiabaticity of the ramp and how this might affect weak signatures such as
small condensate fractions near the critical temperature. During the training
process, we label input pictures of $n_{A}$ with condensate fractions obtained
from bimodal fits to $n_{RR}$ at the same experimental parameters. We exclude
data with temperatures near the critical temperature from learning. Moreover,
in order to prevent the network from learning unwanted correlations between
directly accessible parameters (such as atom number and condensate fraction),
we use training data from different interaction values throughout the
crossover at $1/(k_{F}a)=\\{1.6,1.0,0.5,0.0,-0.5,-0.6\\}$ on a total of 7895
labeled examples. Here, $k_{F}$ denotes the Fermi wave vector calculated from
the atom number and the trap parameters and $a$ the s-wave scattering length.
We extract the critical temperature from the neural network predictions for
the direct-release time-of-flight images across the whole range of the BCS/BEC
crossover by taking a piecewise linear fit of the condensate fraction5.

Qualitatively, the behaviour of the critical temperature of the superfluid
transition across the BCS/BEC crossover can be understood by starting from the
extreme regimes: in the weakly-attractive BCS limit, the critical temperature
scales $k_{B}T_{c}\sim E_{F}\exp[-\pi/(2k_{F}|a|)]$ [1, 28]. Here, $E_{F}$
denotes the Fermi energy. In the opposite regime, far on the BEC side, we
encounter a weakly-repulsively interacting gas of bosons. The bosons have
twice the mass of the fermions $M_{B}=2m$ and half the density $n_{B}=n/2$.
The critical temperature of the ideal Bose gas is simply given by
$k_{B}T_{c}^{0}\sim\frac{\hbar^{2}n_{B}^{2/3}}{M_{B}}$. Unlike in the BCS
regime, the critical temperature of the Bose gas has a very weak dependence on
the interaction strength between the bosons
$T_{c}(a_{B})=T_{c}^{0}\left[1+cn_{B}^{1/3}a_{B}+...\right]$, where
$a_{B}=0.6\,a$ [29] denotes the s-wave scattering length between two bosons,
and $c$ is a positive constant [30, 31]. From this simple argument, we expect
an increase of the critical temperature when approaching the crossover from
the BEC side and hence a maximum critical temperature somewhere in the
crossover regime.

From the previous consideration it is obvious that a careful determination of
both density and temperature is very important. In the trapped gas of our
experiment, the two quantities are inversely related to each other and,
furthermore, also interparticle interactions change the density.

The calibration of density and temperature proceeds in the following way: We
take in-situ absorption images of the trapped gas along two orthogonal spatial
directions (in order to account for asymmetries of the trapped cloud) for
different interaction strengths and temperatures. On these data, we perform an
inverse Abel transform to reconstruct the density distribution inside the
trap. This serves two purposes: on the one hand, we obtain the center density
which we use for the normalisation of the data and on the other hand, the
density distribution $n_{\sigma}(r)$ feeds into the temperature calibration in
the next step. Then, we use the data from the unitary Fermi gas
[$1/(k_{F}a)=0$] and its both theoretically [2, 36, 37, 8] and experimentally
[32, 33, 34, 35] well known critical temperature of $T_{c}=0.167\,T_{F}$ to
precisely reconstruct our trapping potential. To this end, the inverse
equation of state of the unitary Fermi gas [35] is applied to the in-trap
density distribution reconstructed from in-situ high-intensity absorption
images of the cloud at $T=T_{c}$. In the final step, we use the obtained
knowledge of trap potential and measured in-situ density profiles
$n_{\sigma}(r)$ to determine the temperature by fitting a virial expansion of
the equation of state to the outermost regions of the trapped cloud where the
gas is not condensed.

In Figure 2, we show the results of the critical temperature for a homogeneous
gas in comparison with theoretical predictions as a function of the
interaction parameter $1/(k_{F}a)$. Since the condensation will initiate at
regions of highest density, i.e., at the center of the trapping potential, we
adopt a local density approximation and use the density and interaction
parameters at the center of the cloud to compare with the theory of the
uniform gas. Our results show a steady increase of the critical temperature
from the BCS side up to interaction strengths of approximately
$1/(k_{F}a)=0.5$. There, $T_{c}$ levels off and stays approximately constant
or, possibly, declines weakly for higher coupling strengths. Overall, our
results are in very good agreement with several theory predictions in
different ranges of the phase diagram. Throughout the whole crossover, the
agreement with the extended Gorkov-Melik-Bakhudarov (GMB) theory [8] is
striking and both position and value of the maximum $T_{c}$ are well
compatible with the theoretical results. On the BCS side our data are higher
than the Quantum-Monte Carlo calculations [37, 3, 4] and the extended GMB
theory, however, close to the theoretical prediction of reference [2]. The
experimental results agree with trends observed also in earlier measurements
of the phase diagram in both potassium [24] and Lithium gases [25] in which
the critical temperature on the BCS side did not fall off as rapidly with
decreasing interaction strength. Possible suggested explanations include
effects of the harmonic trap [39], the formation of dimers above the resonance
[25] and non-adiabaticities of the rapid ramp in the training data.

In order to verify our methodology, we have double-checked the performance of
the trained neural network using a different Feshbach resonance in another
spin mixture of the Lithium atoms: We prepare a condensate in the hyperfine
states $\ket{1}$ and $\ket{3}$ for which the position and width of the
Feshbach resonance are different and atom number and starting temperature as
compared to the training cases are also different. Nevertheless, the neural
network successfully predicts the condensate fraction at $1/k_{F}a=-0.41$
($k_{F}$ calculated from the atom number) with the same critical temperature.

Figure 2: Critical temperature across the BEC-BCS crossover referenced to the
homogeneous gas. Full symbols: experimental data. Error bars are calculated
from the standard error and an estimation of the systematic error caused by
non-harmonicities of the trap. The latter effect dominates and is discussed in
more detail in the appendix. Dashed line: BCS theory with GMB corrections;
solid line: extended GMB theory [8]; dash-dot line: theory from [2]; dotted
line: interacting BEC; open triangles: quantum Monte-Carlo data [3], open
circles: quantum Monte-Carlo data [4].

An important remaining question is whether the optimized neural network has
learned physically relevant quantities. In other words: does the neural
network spot hidden details in the data during its optimization and can we
extract these information to draw conclusions for the physics? In previous
applications of neural networks to analyze quantum problems this had often not
been considered, and, generally, the question of causality in machine learning
is becoming increasingly important also in computer science. We extract from
the neural network which neurons have been mostly activated. To this end, we
employ a backpropagation-based approach (DeepLIFT, [41]) that assigns
importance scores to the inputs for a given output. The importance scores can
then be identified to reveal those neurons (or, simply put, regions of the
image) that contribute most decisively to the neural network output. In Figure
3, we show the importance scores obtained for different momenta. The results
highlight that in the BCS regime [$1/(k_{F}a)<0$] the neural network output is
dominated by the momentum density near $k\simeq 0.2\,k_{F}$. In contrast, on
the bosonic side [$1/(k_{F}a)>0$], the neural network output mostly relies on
low-$k$ data. This finding is in agreement with the expectation of the effects
of pairing in the fermionic and bosonic sides of the crossover (see Figure 1b)
and indicates that the neural network optimizes indeed for physically relevant
features in the time-of-flight data.

Figure 3: Back-analysis of the neural network. For every interaction strength,
we determine which neurons activate most for determining the condensate
fraction. On the BEC side the most active neurons are at low momenta whereas
at unitarity and on the BCS side the most active are at large momenta.

In summary, we have demonstrated that a neural network can be utilized for the
detection of quantum phases of strongly-interacting systems with high
accuracy. Moreover, we show that back-analysis of the trained neural network
allows to extract insights much beyond standard fitting routines and therefore
opens a new route to the precision analysis of physical data.

This work has been supported by the Alexander-von-Humboldt Stiftung, DFG
(SFB/TR 185 project C6), Cluster of Excellence Matter and Light for Quantum
Computing (ML4Q) EXC 2004/1 – 390534769. We thank Kilian Kluge for
discussions.

## References

  * [1] C. A. R. Sá de Melo, M. Randeria, J. R. Engelbrecht, Phys. Rev. Lett. 71, 3202 (1993).
  * [2] R. Haussmann, W. Rantner, S. Cerrito, W. Zwerger, Phys. Rev. A 75, 023610 (2007).
  * [3] E. Burovski, E. Kozik, N. Prokof’ev, B. Svistunov, M. Troyer, Phys. Rev. Lett. 101, 090402 (2008).
  * [4] A. Bulgac, J. E. Drut, P. Magierski, Phys. Rev. A 78, 023625 (2008).
  * [5] S. Floerchinger, M. Scherer, S. Diehl, C. Wetterich, Phys. Rev. B 78, 174528 (2008).
  * [6] S. Floerchinger, M. M. Scherer, C. Wetterich, Phys. Rev. A 81, 063619 (2010).
  * [7] M. Pini, P. Pieri, G. C. Strinati, Phys. Rev. B 99, 094502 (2019).
  * [8] L. Pisani, A. Perali, P. Pieri, G. C. Strinati, Phys. Rev. B 97, 014528 (2018).
  * [9] D. T. Son, M. A. Stephanov, Phys. Rev. A 74, 013614 (2006).
  * [10] J. Carlson, S. Reddy, Phys. Rev. Lett. 100, 150403 (2008).
  * [11] W. Ketterle, M. W. Zwierlein, Proceedings of the International School of Physics ”Enrico Fermi” (2008), vol. 164.
  * [12] G. Carleo, M. Troyer, Science 355, 602 (2017).
  * [13] A. Tanaka, A. Tomiya, Journal of the Physical Society of Japan 86, 063001 (2017).
  * [14] E. van Nieuwenburg, Y.-H. Liu, S. Huber, Nature Physics 13, 435 (2017).
  * [15] G. Torlai, et al., Nature Physics 14, 447–450 (2018).
  * [16] L. Wang, Phys. Rev. B 94, 195105 (2016).
  * [17] P. Huembeli, A. Dauphin, P. Wittek, Phys. Rev. B 97, 134109 (2018).
  * [18] B. S. Rem, et al., Nature Physics 15, 917 (2019).
  * [19] L. H. Gilpin, et al., Explaining explanations: An overview of interpretability of machine learning (2018).
  * [20] A. Behrle, et al., Nature Physics 14, 781 (2018).
  * [21] E. D. Kuhnle, et al., Phys. Rev. Lett. 106, 170402 (2011).
  * [22] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization (2017), arXiv: 1412.6980.
  * [23] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning (MIT Press, 2016).
  * [24] C. A. Regal, M. Greiner, D. S. Jin, Phys. Rev. Lett. 92, 040403 (2004).
  * [25] M. Zwierlein, et al., Phys. Rev. Lett. 92, 120403 (2004).
  * [26] E. Altman, A. Vishwanath, Phys. Rev. Lett. 95, 110404 (2005).
  * [27] I. Tikhonenkov, E. Pazy, Y. B. Band, M. Fleischhauer, A. Vardi, Phys. Rev. A 73, 043605 (2006).
  * [28] L. P. Gorkov, T. K. Melik-Barkhudarov, Sov. Phys. JETP 13, 1018 (1961).
  * [29] D. S. Petrov, C. Salomon, G. V. Shlyapnikov, Phys. Rev. Lett. 93, 090404 (2004).
  * [30] P. Arnold, G. Moore, Phys. Rev. Lett. 87, 120401 (2001).
  * [31] V. A. Kashurnikov, N. V. Prokof’ev, B. V. Svistunov, Phys. Rev. Lett. 87, 120402 (2001).
  * [32] L. Luo, J. E. Thomas, Journal of Low Temperature Physics 154, 1 (2009).
  * [33] M. Horikoshi, S. Nakajima, M. Ueda, T. Mukaiyama, Science 327, 442 (2010).
  * [34] S. Nascimbene, N. Navon, K. Jiang, F. Chevy, C. Salomon, Nature 463, 1057 (2010).
  * [35] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, M. W. Zwierlein, Science 335, 563 (2012).
  * [36] O. Goulko, M. Wingate, Phys. Rev. A 82, 053621 (2010).
  * [37] E. Burovski, N. Prokof’ev, B. Svistunov, M. Troyer, Phys. Rev. Lett. 96, 160402 (2006).
  * [38] S. Pilati, S. Giorgini, N. Prokof’ev, Phys. Rev. Lett. 100, 140405 (2008).
  * [39] A. Perali, P. Pieri, L. Pisani, G. C. Strinati, Phys. Rev. Lett. 92, 220404 (2004).
  * [40] G. Baym, J.-P. Blaizot, M. Holzmann, F. Laloë, D. Vautherin, Phys. Rev. Lett. 83, 1703 (1999).
  * [41] A. Shrikumar, P. Greenside, A. Kundaje, Learning important features through propagating activation differences, arXiv: 1704.02685 (2017).
  * [42] J. Kinast, et al., Science 307, 1296 (2005).
  * [43] G. Reinaudi, T. Lahaye, Z. Wang, and D. Gu$\acute{\text{e}}$ry-Odelin, Strong absorption imaging of dense clouds of ultracold atoms”, Opt. Lett. 32, 3143-3145 (2007)
  * [44] M. Abadi, et al., TensorFlow: Large-scale machine learning on heterogeneous systems (2015). Software available from tensorflow.org.
  * [45] D. Lee, T. Schäfer, Phys. Rev. C 73, 015201 (2006)
  * [46] X. Leyronas, Phys. Rev. A 84, 053633 (2011)
  * [47] Y. Hou, J. E. Drut, Phys. Rev. Lett. 125, 050403 (2020)

## Appendix

### Preparation and detection the sample

The Lithium atoms are confined in a trap formed by two intersecting laser
beams of $1070\,\text{nm}$ wavelength with final trap frequencies in harmonic
approximation of $2\pi\times(168,166,238)\,\text{Hz}$. The final temperature
of the sample after evaporation close to a Feshbach resonance at
$834\,\text{G}$ is $T/T_{F}=0.08\pm 0.01$. After preparation of the condensate
in the crossover regime, we choose a desired interaction strength $1/k_{F}a$
by adiabatically ramping the magnetic field to the corresponding field value
$B_{\text{final}}$. To controllably increase the temperature, we suddenly
decompress and subsequently recompress the trap for a variable amount of time
[42, 21]. This is followed by a hold time of $50\,\text{ms}$ for
thermalisation. We perform detection of the gas by time-of-flight imaging. To
this end, the optical dipole trap is extinguished rapidly and the gas expands.
Owing to a residual inhomogeneity of our magnetic bias field, the gas expands
into a weak harmonic trapping potential with frequencies
$2\pi\times(16,16,20i)\,\text{Hz}$.

In order to image the bimodal distribution $n_{\text{RR}}$, the cloud is
subjected to a rapid projection onto Feshbach molecules by suddenly ramping
the magnetic field to the zero-crossing of the scattering length at a magnetic
field of 534 G before $15\,\text{ms}$ time-of-flight. The cloud is then imaged
close to the resonance on the repulsive side at a magnetic field strength of
760 G. To image the density distribution of the atoms after release from the
trap $n_{\text{A}}$, we perform standard time-of-flight imaging after
$5\,\text{ms}$ without changing the magnetic field.

For the reconstruction of the equation of state, we perform high-intensity
absorption imaging in-situ to resolve the very high densities in the trap
[43]. Because the size of the cloud has to be known in all spatial directions,
the high-intensity absorption imaging is done along two perpendicular lines of
sight.

### The neural network

Layer (type) | Output Shape | Parameters
---|---|---
Input | (Batch, 150, 170, 1) | 0
2D Convolutional | (Batch, 150, 170, 30) | 300
Max Pooling | (Batch, 75, 85, 30) | 0
2D Convolutional | (Batch, 75, 85, 40) | 30040
Max Pooling | (Batch, 37, 42, 40) | 0
2D Convolutional | (Batch, 37, 42, 50) | 50050
Max Pooling | (Batch, 18, 21, 50) | 0
Flatten | (Batch, 18900) | 0
Dense | (Batch, 600) | 11340600
Dropout | (Batch, 600) | 0
Dense | (Batch, 300) | 180300
Dropout | (Batch, 300) | 0
Dense | (Batch, 1) | 301
Table 1: Neural network architecture used for the phase diagrams.

The neural network architecture used to generate the phase diagram is given in
table 1 and is realised with the TensorFlow library [44]. It consists of three
convolutional layers combined with max pooling layers followed by two dense
layers with dropout regularisation and one output neuron. The total number of
tunable parameters is $11601591$. We train the network with stochastic
gradient descent using Adam optimizer with learning rate $1.5\cdot 10^{-4}$
[22] on the mean squared error loss function. We use 7895 labelled datapoints
in total from which $90\,\%$ are used for training and $10\,\%$ for
validation. We train each network for 15 epochs with a batch size of 30. The
training data is shuffled after each epoch. We tested several network
architectures and generally found the performance robust against changes of
the hyper parameters. We conclude that network architectures similar to the
one used here provide robust learners for the detection of the condensate
fraction. Moreover, we have taken data sets with different magnifications of
the optical imaging system. The images were then scaled prior to feeding the
data into the neural network and transfer learning on 2500 images was
performed, while freezing the parameters of the convolutional layers. We found
our data to be robust against this scaling operation.

### Equation of state reconstruction

We reconstruct the equation of state under local density approximation
$n_{\sigma}(\mu-V)$ for all temperatures and interaction strengths entering
the phase diagram. Here, $n_{\sigma}$ is the density distribution of the cloud
in the trap, $\mu$ the chemical potential, and $V$ the external trapping
potential. From the equation of state, we can extract the temperature of the
cloud $T$ and the density in the centre of the trap $n_{\sigma}|_{V=0}$, which
are the quantities needed to calibrate the phase diagram. The Fermi wave
vector $k_{F}$ and Fermi temperature $T_{F}$ are related to
$n_{\sigma}|_{V=0}$ via $k_{F}=(3\pi^{2}n_{\sigma}|_{V\\!=\\!0})^{1/3}$ and
$T_{F}=\frac{\hbar^{2}k_{F}^{2}}{2mk_{B}}$ where $\hbar$ denotes the reduced
Planck constant, $k_{B}$ the Boltzmann constant, and $m$ the mass of 6Li
atoms.

To determine the temperature of the cloud we utilize that, close to the
surface of the cloud, the equation of state can be approximated by a virial
expansion $n_{\sigma}\lambda^{3}=\sum_{n}{nb_{n}z^{n}}$ with
$\lambda=\sqrt{\frac{2\pi\hbar^{2}}{mk_{B}T}}$ the thermal de-Broglie
wavelength, $b_{n}$ the $n$-th virial coefficient, and
$z=\exp((\mu-V)/k_{B}T)$. By fitting the virial expansion up to order $n$=3
(or $n$=4 for clouds furthest in the BCS regime) to the dependence of
$n_{\sigma}$ on $V$ close to the surface of the cloud we get the temperature
of the cloud $T$. While $b_{1}$=1 for all scattering lengths $a$ and
temperatures $T$, higher order $b_{n}$s depend on $a$ and $T$. For $b_{2}$ an
analytic formula is known [45], $b_{3}$ has been calculated for a wide range
of $a$ and $T$ [46] and also for $a<0$ in [47] where $b_{4}$ for $a<0$ is
calculated as well.

To perform the equation of state reconstruction, intrap density and potential
have to be known. The density is determined from in-situ optical density
images by first taking the elliptic radial average - respecting the cloud’s
aspect ratio - and converting the optical density to column density
$n_{\mathrm{col}}(r)$, i.e. the density integrated along the camera’s line of
sight. The actual density $n_{\sigma}(r)$ is then reconstructed from
$n_{\mathrm{col}}(r)$ by an elliptic inverse Abel transform

$n_{\sigma}(r)=-\frac{\sigma_{r}}{\sigma_{z}}\frac{1}{\pi}\int_{r}^{\infty}\\!\mathrm{d}r^{\prime}\frac{1}{\sqrt{r^{\prime
2}-r^{2}}}\frac{\partial n_{\mathrm{col}}(r^{\prime})}{\partial r^{\prime}}$
(1)

where $\sigma_{r}\,(\sigma_{z})$ is the size of the cloud in radial direction
(the camera’s line of sight).
The external potential is only known in harmonic approximation from the trap
frequencies, but this approximation is not valid for the crossed Gaussian-beam
dipole trap in the region where we perform our thermometry. We therefore
derive the full external trapping potential $V(r)$ from the density
$n_{\sigma}(r)$ of a cloud with $1/k_{F}a=0$ and $T=T_{c}$. Since the critical
temperature of a homogeneous Fermi gas at unitarity is well known to be
$T_{c}=0.167\;T_{F}$ [35], it is sufficient to determine $T_{F}$ (resp.
$n_{\sigma}|_{V\\!=\\!0}$) instead of $T$. From reference [35] the chemical
potential of the unitary Fermi gas at $T_{c}$ can also be related to $T_{F}$
by $\mu|_{T=T_{c}}=0.416\;k_{B}T_{F}$. With known $\mu$ and $T$, combining the
local density approximation $\mu\rightarrow\mu-V(r)$ with the known equation
of state of the unitary Fermi gas [35]
$n_{\sigma}(\mu,T)=n_{\sigma}(\mu/k_{B}T)\rightarrow
n_{\sigma}(\,(\mu-V(r))/k_{B}T\,)$ yields a relation between $n_{\sigma}(r)$
and $V(r)$ without free parameters. Inversion of this relation allows the
derivation of $V(r)$ from $n_{\sigma}(r)$ of a cloud with $1/k_{F}a=0$ and
$T=T_{c}$.
Because $V(r)$ is identical for all clouds, it is possible to fit the equation
of state’s virial expansion close to the surface of the cloud, where
$\exp((\mu-V(r))/k_{B}T)\ll 1$ ensures the validity of the virial expansion,
with the free parameters $\mu$ and $T$. $T_{F}$ can again be calculated from
the density at $V=0$ to get $T/T_{F}$ for arbitrary interaction strengths and
temperatures.

It should be noted that the inverse Abel transform assumes elliptic
equipotential lines, a condition only approximately fulfilled in our crossed
Gaussian-beam dipole trap. We estimate the influence of this systematic error
by simulating the column density of a cloud with $1/k_{F}a=0$ and $T=T_{c}$ in
a trap comparable to the one used in the experiment. We then use the same
procedure as for the experimental data to reconstruct the external potential
and perform our thermometry on simulated column densities of an ideal Fermi
gas. The error in T introduced in this way stays below 5%. We therefore use
this value as an upper bound to estimate the systematic error of our
thermometry. It is the dominating contribution to the uncertainty of the
critical temperature as depicted in fig. 2.

The inverse Abel transform also relies on the differential of the column
density $\partial n_{\mathrm{col}}/\partial r$ which is very susceptible to
noise when extracting it from the experimental data. We therefore average ca.
30 identically prepared clouds and perform radial averaging. However,
experimental noise still dominates close to the center of the cloud where
radial averaging has little effect. To determine the central density more
reliably, we linearly extrapolate the measured $n_{\sigma}(V)$ data towards
$V=0$.


# The Early Evolution of Primordial Pair-Instability Supernovae

C.C. Joggerst11affiliation: Department of Astronomy and Astrophysics,
University of California at Santa Cruz, Santa Cruz, CA 95060\.
<EMAIL_ADDRESS>22affiliation: Nuclear and Particle Physics,
Theoretical Astrophysics and Cosmology (T-2), Los Alamos National Laboratory,
Los Alamos, NM 87545 and Daniel J. Whalen33affiliation: McWilliams Fellow,
Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213

###### Abstract

The observational signatures of the first cosmic explosions and their chemical
imprint on second-generation stars both crucially depend on how heavy elements
mix within the star at the earliest stages of the blast. We present numerical
simulations of the early evolution of Population III pair-instability
supernovae with the new adaptive mesh refinement code CASTRO. In stark
contrast to 15 - 40 M⊙ core-collapse primordial supernovae, we find no mixing
in most 150 - 250 M⊙ pair-instability supernovae out to times well after
breakout from the surface of the star. This may be the key to determining the
mass of the progenitor of a primeval supernova, because vigorous mixing will
cause emission lines from heavy metals such as Fe and Ni to appear much sooner
in the light curves of core-collapse supernovae than in those of pair-
instability explosions. Our results also imply that unlike low-mass Pop III
supernovae, whose collective metal yields can be directly compared to the
chemical abundances of extremely metal-poor stars, further detailed numerical
simulations will be required to determine the nucleosynthetic imprint of very
massive Pop III stars on their direct descendants.

##  1\. Introduction

The first stars in the universe form at $z\sim$ 20 and are likely very
massive, 30 - 500 M⊙ (Bromm et al., 1999; Abel et al., 2000, 2002; Bromm et
al., 2002; Nakamura & Umemura, 2001; O’Shea & Norman, 2007). The fates of
these stars depend on their masses: 15 - 50 M⊙ stars die in core collapse
supernovae (CC SNe), 140 - 260 M⊙ stars explode in far more energetic
thermonuclear pair-instability supernovae (PISNe Heger & Woosley, 2002), and
40 - 60 M⊙ stars may die as hypernovae, whose explosion mechanisms are not yet
understood but are thought to have energies intermediate to those of CC and
pair-instability SNe (Iwamoto et al., 2005). Most Pop III SNe (Bromm et al.,
2003; Kitayama & Yoshida, 2005; Greif et al., 2007) occur in low densities
(0.1 - 1 cm-1) because UV radiation from the star sweeps most of the baryons
from the dark matter halo in which it resides (Whalen et al., 2004; Kitayama
et al., 2004; Alvarez et al., 2006; Abel et al., 2007; Wise & Abel, 2008a).

Metals from Pop III SNe determine the character of second-generation stars and
the primeval galaxies they populate by enhancing cooling in the gas in which
such stars form, and therefore the mass scales on which it fragments. The
manner in which the first metals contaminate pristine gas in the primordial
IGM crucially depends on mixing processes that have only begun to be studied
numerically. Preliminary calculations indicate that metals from Pop III SNe
mix with gas in a halo on two characteristic spatial scales prior to their
emergence into cosmological flows on kpc scales (Whalen et al., 2008): 10 - 15
pc, when a reverse shock forms in the remnant, and 150 - 200 pc, when the
remnant collides with the dense shell of the relic H II region of the
progenitor. New simulations now prove that violent mixing can occur within the
star prior to shock breakout from its surface, if the star is 15 - 40 M⊙
(Joggerst et al., 2010a, hereafter JET10). Capturing mixing on the smallest
scales is prerequisite to following its cascade out to the largest ones, and
hence to determining the colors and morphologies of primitive galaxies. Early
mixing also governs which elements are imprinted on next-generation stars,
whose chemical abundances may impose indirect constraints on the masses of the
stars that enriched them. Low-mass remnants of this generation are now sought
in ongoing surveys of ancient, dim metal-poor stars in the Galactic halo
(Beers & Christlieb, 2005; Frebel et al., 2005). Since mixing in Pop III SNe
in part sets elemental abundances in second-generation stars, it is integral
to future measures of the primordial IMF.

Early mixing is also key to computing the light curves and spectra of the
first cosmic explosions, whose detection could yield the first direct measure
of the primordial IMF. Optical and UV radiation breaks free of the SN shock
when the ejecta reaches the outer envelope of the star and is exposed to the
IGM. In the frame of the shock, the photosphere from which photons escape into
space descends deeper into the ejecta as it expands outward because of the
spherical dilution of the ejecta. If mixing precedes radiation breakout, it
determines the elements that the photosphere encounters as it sinks deeper
into the ejecta, and therefore the emission lines that propagate into the IGM
over time. Accurate multidimensional models of the explosion from its earliest
stages are therefore necessary to compute lines in primordial SN light curves
and spectra.

Mixing in galactic core-collapse supernovae has been studied for over twenty
years, particularly in connection with SN 1987A (e.g. Fryer et al., 2007, and
references therein). A core objective of these studies is to understand the
premature appearance of ${}^{56}\mathrm{Ni}$ in the spectra of 1987A, whose
emission lines appear much earlier than a simple picture of segregated,
spherically-symmetric expanding mass shells would predict. Mixing in Pop III
core-collapse SNe was first examined by Joggerst et al. (2009), who found that
both vigorous mixing and fallback onto the compact remnant in 15 - 40 M⊙ Pop
III SNe govern which metals escape into the IGM at high redshifts.

PISNe may have been commonplace in the primeval universe, but their enrichment
of the early IGM is yet to be understood. To investigate the propagation of
metals into the IGM by such explosions, constrain their nucleosynthetic
imprint on second-generation stars, and to evaluate the impact of mixing on
PISN spectra, we have performed two-dimensional simulations of the explosions
of 150 - 250 M⊙ stars. In $\S\,2$ we review how blast profiles from the KEPLER
code were ported as initial conditions to CASTRO and discuss what factors
govern the presupernova structure of the star. We examine the outcomes of the
explosions in $\S\,3$ and in $\S\,4$ we conclude.

##  2\. Models

###  2.1. KEPLER

As in JET10, the simulations in our PISN survey were carried out in two
stages. First, primordial stars were evolved through all stages of stable
nuclear burning from the zero-age main sequence to initial collapse via the
pair instability in the one-dimensional Lagrangian stellar evolution code
KEPLER (Weaver et al., 1978; Woosley et al., 2002). The PISN is triggered when
this collapse incites explosive O, and some cases Si, burning. Unlike the
simulations of JET10, in which the blast is artificially launched with a
piston, the pair instability and subsequent collapse that triggers these
explosions are an emergent feature of the stellar evolution calculation. They
are genuinely spherical, barring (magneto)rotational effects, and their
energies are set by O and Si burning. The blast was followed until the end of
all nuclear burning, $\sim$ 20 s after the start of the explosion. The energy
generated was computed with a 19-isotope network up to the point of oxygen
depletion in the core of the star and with a 128-isotope quasi-equilibrium
network thereafter.

Table 1PISNe Models: Properties at time models were mapped to 2D  model |  $M_{He}$ ($\mathrm{M}_{\odot}$) |  $M_{N}$ ($\mathrm{M}_{\odot}$) |  $M_{{{}^{56}\mathrm{Ni}}}$ ($\mathrm{M}_{\odot}$) |  $M_{final}$ ($\mathrm{M}_{\odot}$) |  $R$ ($10^{13}$ cm) |  $E_{kin}$ ($10^{51}$ ergs)
---|---|---|---|---|---|---
u150 | 41 | 9.1(-5) | 0.079 | 143 | 16 | 3.7
u175 | 46 | 1.1(-4) | 0.72 | 164 | 18 | 9.5
u200 | 50 | 1.3(-4) | 5.2 | 183 | 19 | 17
u225 | 54 | 1.1(-3) | 17 | 200 | 33 | 28
u250 | 68 | 1.7(-4) | 38 | 238 | 23 | 34
z175 | 53 | 1.6(-5) | 0.24 | 175 | 4.2 | 4.7
z200 | 60 | 1.8(-5) | 2.0 | 200 | 4.5 | 12
z225 | 67 | 1.7(-5) | 8.8 | 225 | 4.9 | 17
z250 | 74 | 1.6(-5) | 23 | 250 | 6.2 | 29

###  2.2. CASTRO

The one-dimensional explosion profiles were then mapped onto two-dimensional
$RZ$ axisymmetric grids in CASTRO (Compressible ASTROphysics), a multi-
dimensional Eulerian AMR code with a high-order unsplit Godunov hydro solver
(Almgren et al., 2010). Each explosion was evolved past breakout from the
surface of the star until all mixing ceased and each element in the ejecta was
expanding homologously in mass coordinate. We smoothly join the density at the
surface of the star to a uniform circumstellar medium of 1 cm-3 with an
$r^{-3.1}$ power law, in keeping with the usual assumption of a low-density H
II region around the progenitor star with no wind-blown shell or prior mass
ejection. The medium beyond the star has no effect on the dynamics within the
star if its density falls more steeply than $r^{-3}$. In mapping the radial
profile onto the $RZ$ grid in CASTRO, care was taken to resolve the key
elements of the explosion: the shock, the shells of elements, and the high-
density core. In particular, both the ${}^{56}\mathrm{Ni}$ core and the O
shell were resolved with a minimum of 16 cells. Our mapping excludes
departures from spherical symmetry due to O burning, but such perturbations
would likely be high mode and low amplitude and therefore have minimal effect
on the evolution of instabilities. Our models thus only capture later
asymmetries of mode greater than $l=1$ or 2.

As in JET10, we adopt a monopole approximation for self-gravity. We first
compute a radial average of the density from the $RZ$ grid to create a one-
dimensional density profile. We then compute a one-dimensional gravitational
potential from this profile and map it back onto the $RZ$ grid. Since
departures from spherical symmetry in the densities are minor, this
approximation has a negligible effect on the final state of the explosion.
Because the PISN completely disperses the star, there is no compact remnant,
fallback, or thus any need to include a point potential centered at the origin
of the coordinate mesh.

We follow 15 chemical elements as individual species, each with their own
continuity equation, and calculate local energy deposition due to radioactive
decay of ${}^{56}\mathrm{Ni}$ in the same manner as JET10. However, the
formation of a nearly degenerate core at white dwarf densities in the
progenitor necessitates the use of the Helmholtz equation of state (EOS) at
early stages of the explosion. As the ejecta expands and cools, we transition
back to the ideal EOS used in JET10, which assumes that the gas is fully
ionized and includes contributions from both radiation and ideal gas pressure.

The base grid is 10242, with the initial outer boundaries set so that the
inner portion of the star is resolved as described above using no more than 6
levels of refinement. The star is centered at the lower left corner of the
mesh. We apply reflecting and outflow boundary conditions to the inner and
outer boundaries of the grid, respectively. Our refinement criteria are the
same as those in JET10. When the shock nears the edge of the grid, the
simulation is stopped, the grid is doubled, and the calculation is then
restarted, subtracting or adding levels of refinement as needed. This
procedure is repeated up to 12 times, depending on the model. We halt the
simulation when all chemical species are expanding homologously, which always
occurs by the time the ejecta has propagated a short distance into the uniform
circumstellar density.

###  2.3. Progenitor Models

JET10 found that mixing in low-mass Pop III core-collapse explosions primarily
depends on the internal structure of the progenitor. We likewise expect the
early evolution of pair-production explosions to be determined by the envelope
of the star, which is determined its mass, internal convective mixing over its
lifetime, and by its metallicity. Capturing the full range of structures for
these stars is essential to a comprehensive survey of early mixing in Pop III
supernovae.

####  2.3.1 Convective Mixing

The initial absence of metals and the large contribution of radiation to the
pressure in massive Pop III stars promotes convection within them. The CNO
cycle cannot begin in primordial stars until a threshold mass fraction of
${}^{12}\mathrm{C}$ is first created by the triple $\alpha$ process. This
trace ${}^{12}\mathrm{C}$ sets the entropy of the hydrogen layer to be just
above that in the core, without the sharp entropy gradient in the upper layer
of the helium shell that is usually present in He burning stars with metals.
This plus radiation pressure facilitates convection. In 140 - 260 M⊙ stars,
the central convection zone can approach, come in contact with, or even reach
into the lower hydrogen layers, mixing them with carbon brought up from the
core during helium burning. When these two high-temperature components mix,
they burn vigorously, elevating energy release rates in the H shell by up to
several orders of magnitude (the so-called hydrogen boost).

Convection affects the structure of the star in two ways. First, since it
raises energy production (and hence opacities) in the lower hydrogen layers,
the star can puff up by more than an order of magnitude in radius and acquire
a red supergiant structure. Second, if convection is extreme the transport of
material out of the core could reduce its size and explosion energy in
comparison to modest convection. Unfortunately, one-dimensional stellar
evolution models cannot predict these inherently three-dimensional processes
from first principles. Instead, they parametrize them with semi-convection
coefficients.

####  2.3.2 Metallicity

Gas in high-redshift halos that is enriched to metallicities below 10-3.5 Z⊙
fragments on mass scales that are essentially identical to those of pristine
gas and still forms very massive stars (e. g. Bromm et al., 2001; Mackey et
al., 2003; Smith & Sigurdsson, 2007). However, such small metal fractions are
more than enough to enhance CNO reaction and energy generation rates in the
hydrogen burning layers of the star, enlarging it in the same manner as
convective mixing (Hirschi, 2007; Ekström et al., 2008). Hence, we would
expect a strong degeneracy between the influence of metals and convection on
the envelope of the star, and that the full range of mixing in PISNe can be as
easily spanned by metallicity as by convective overshoot.

We considered 150, 175, 200, 225 and 250 M⊙ non-rotating progenitors at $Z=0$
(the z-series) and $Z=10^{-4}Z_{\odot}$ (the u-series), which we summarize in
Table 1. The $Z=0$ 150 M⊙ star collapses to a black hole without exploding, so
we exclude it from our CASTRO models. We employ metallicity rather than
convective overshoot to cover the range of plausible progenitor structures in
our study, because there is uncertainty about how much semi-convection there
is in a given star but none about the range of metallicities over which it can
form ($Z$ = 0 and 10-4. Comparison of Table 1 with Table 1 of Scannapieco et
al. (2005) confirms that for a given progenitor mass these two metallicities
do yield upper and lower limits to stellar radius similar to those for all
reasonable values of semi-convection coefficients.

##  3\. Results

As expected, the u-series models die as red giants, with radii more than an
order of magnitude larger than those of the z-series, which die as blue
giants. As shown in Table 1, the u-series models in general have more
energetic explosions than the z-series models for a given mass. As we show in
Figure 1, convective mixing has completely disrupted the helium layer in model
u225 and mixed it with the hydrogen envelope. The helium layer has also been
mixed with the hydrogen envelope, although to a lesser extent, in model u200
and to an even slighter degree in model u250. These u-series models also
experienced some mass loss due to pulsations.

The $\rho r^{3}$ values through which the shock passes are an effective
predictor of post-explosion hydrodynamics. In regions where $\rho r^{3}$
increases, the shock must decelerate, and a reverse shock forms that inverts
the pressure gradient and induces Rayleigh-Taylor (RT) instabilities. Models
with steeper values of $\rho r^{3}$ will experience more mixing because the
forward shock slows down more abruptly, which leads to a stronger reverse
shock. We show $\rho r^{3}$ values for all models in Figure 1, scaled to the
maximum value in model u225. It is clear from this figure that models u225 and
u200 exhibit the largest increases in $\rho r^{3}$ near the outer edge of the
star. This is due to slight bumps in density that are connected to the
pulsations that ejected mass from these stars earlier.

Figure 1.— The structure of helium and oxygen shells (delineated by dashed
lines for clarity) in all the progenitors, with $\rho r^{3}$ superimposed on
them. The $\rho r^{3}$ profiles have all been scaled to the maximum $\rho
r^{3}$ in model u225. $\rho r^{3}$ increases dramatically near the edge of
some models that die as red giants, while none of the blue giants show a
similar structure. The plateaus in helium abundance above cosmological values,
which are most apparent in model u225, denote convective regions in which the
hydrogen envelope mixed with a portion of the helium layer. This process
completely disrupted the helium layer in model u225, and mixed a smaller
fraction of the helium layer with the hydrogen envelope in model u200. A small
amount of the helium shell convectively mixed with the hydrogen envelope in
model u250. Figure 2.— Images of density and oxygen abundance for all models
after the forward shock exits the star or the reverse shock (if one forms)
traverses the star, whichever is later. Density is scaled to the minimum and
maximum values within a given simulation. A dense shell formed in all models
with masses greater than 200 $\mathrm{M}_{\odot}$, but the RT instability only
grows appreciably in the u200 and u225 models. These were the largest
progenitors in radius and had the steepest $\rho r^{3}$ curves. Only slight RT
growth is evident in model u250.

We show final states of mixing for the PISNe in Figure 2. The models are shown
at different times, but always after the shock emerged from the surface of the
star and the reverse shock, if one formed, dissipated. Density is shown on the
left in each panel, with values scaled to the maximum and minimum densities in
the simulation, while oxygen is shown on the right, with the values indicated
in the color bar. Some similarities between the models exist across all
stellar structures. In the higher-mass compact models (z225 and z250) and in
the red u-series models above 150 solar masses, the bulk of the oxygen layer
was swept up into a shell of higher density than the material on either side.

In general, however, the z-series explosions evolved quite differently than
the u-series SNe. In the z-series, no reverse shock formed. A reverse shock
formed in all the u-series models, but dissipated by the time it reached the
density contrast at the oxygen layer (u150 and u175) or was too weak by the
time it reached this dense shell to induce rapid RT growth (u250). In two
models, a reverse shock formed that was strong enough to drive rapid growth of
RT instabilities (u200 and u225). The reverse shock was strongest in model
u225, which had the steepest increase in $\rho r^{3}$ and the second highest
explosion energy, after model u250. The RT instability is clearly visible in
the u200 and u225 panels in Figure 2. In these models the instability became
nonlinear because individual RT fingers strongly interacted. In model u250,
the RT instability only reached the early linear phase before its growth was
halted.

The extent to which the structure of the PISN at the time of explosion was
preserved or disrupted by mixing is shown in Figure 3. We plot the initial
shell structure of the PISN at the time of mapping into CASTRO as a dashed
line, and overlay a solid line indicating the abundances of these elements
after the shock exited the star or the reverse shock traversed the star,
whichever was later. The slight smoothing in the final profiles in comparison
to the initial ones is due to numerical diffusion over the course of the
simulation. The u200 and u225 panels demonstrate the extent of mixing in these
explosions. The oxygen layer has been completely disrupted, and mixing has
penetrated just to the top of the silicon layer (in u200) or through the
silicon layer (in u225). The most vigorous mixing occured in model u225, where
hydrogen shell boost led to convection that completely mixed the helium shell
and the hydrogen envelope. The only other models to experience RT mixing,
models u200 and u250, are also the only other models that show evidence of
convective mixing in the outer region of the star prior to explosion. Model
u200 manifests more RT mixing and a larger region (in radius and mass) in
which convection mixed part of the helium layer with the hydrogen envelope
prior to explosion. Model u250, which has the smallest amount of convection
extending into the helium layer, also exhibits the least RT growth. Models
with no convective mixing between the helium layer and the hydrogen envelope
manifest no RT mixing during the explosion, so mixing seems tied to the depth
to which this convective envelope has penetrated the helium layer of the star.
The deeper the convective envelope, the more RT mixing occurs during the
supernova shock.

Figure 3.— Initial and final abundances of He, O, Si, and Fe for all models.
The RT instability is only present in models u200, u225, and to a small degree
in u250. The slight disagreement between initial and final profiles for models
which experienced no mixing is due to numerical diffusion over the evolution
times of these models. In none of the mixed models did mixing break through
the Si layer of the star, so ${}^{56}\mathrm{Ni}$ never reaches the surface of
the explosion as it does in core-collapse supernova explosions.

RT-induced mixing in PISNe, if it occurs at all, is unlikely to leave the
remnant in a state that resembles a core-collapse supernova remnant. In
particular, it is unlikely to dredge up significant amounts of oxygen, let
alone silicon or ${}^{56}\mathrm{Ni}$, into the outer regions of the star, or
draw much hydrogen back towards the center of the blast. The most vigorous
mixing occurs in the model with the helium shell that has been completely
disrupted by convective mixing with the hydrogen envelope, and likely
represents the most mixing that can occur in a PISN. Even in this model,
${}^{56}\mathrm{Ni}$ cannot reach the upper layers of lighter elements.

As noted earlier, our CASTRO simulations exclude the first stages of the
explosion, in which oxygen burning drives convective mixing that may perturb
the star (Chen et al., 2010). Since these perturbations are expected to be of
low mode and high amplitude, it is unlikely they would alter the essential
conclusions of our survey. Reprising our two-dimensional calculations in three
dimensions is unlikely to yield significantly different results. The RT
instability initially grows about $30\%$ faster in three dimensions than in
two because of artificial drag forces arising from two-dimensional geometries
(Hammer et al., 2010). However, once the fingers of the instability begin to
interact with one another, they mix more efficiently in three dimensions than
in two, reducing the Atwood number and hence their growth rate (Joggerst et
al., 2010b). These two effects cancel each other, and the width of the mixed
region in two and three dimensions is the same. For spherical simulations like
ours in which the RT fingers significantly interact, two-dimensional and
three-dimensional simulations would exhibit comparable degrees of mixing.
Also, the RT instability will not appear in three dimensions when it is not
manifest in two.

##  4\. Discussion and Conclusions

Unlike 15 - 40 M⊙ core-collapse Pop III SNe, 150 - 250 M⊙ Pop III PISNe
experience either no internal mixing prior to shock breakout from surface of
the star or only modest mixing between the O and He shells. Minor mixing
occurs in only two of the explosions and is due to the formation of a reverse
shock that is strong enough to trigger the RT instability at the dense shell
created by the forward shock at the top of the oxygen layer. The degree of
mixing is principally a function of how well hydrogen shell boost mixes the
helium shell with the hydrogen envelope. The model in which the helium shell
is completely mixed with the hydrogen envelope exhibits the most mixing. The
general lack of internal mixing in PISNe has several consequences for early
chemical enrichment of the IGM, the formation of second generation stars, and
the observational signatures of such explosions.

First, the elements that are expelled by low-mass Pop III SNe (which are
governed by both mixing and fallback onto the compact remnant) are later
imprinted on new stars in essentially the proportions in which they are
created by the explosions, regardless of how and where the stars form. This is
because the metals are already highly mixed by the time the shock exits the
star and are merely diluted upon further expansion into the halo and IGM,
where new stars might form. This implies that IMF averages of nucleosynthetic
yields of primordial core-collapse SNe can be directly compared to chemical
abundances in ancient metal-poor stars, without regard for any intervening
hydrodynamical processes. Indeed, the fact that elemental yields from Salpeter
IMF averages of 15 - 40 M⊙ progenitors in JET10 match those found in two
observational surveys of extremely metal-poor (EMP) stars suggests that the
bulk of early chemical enrichment may have been due to low-mass Pop III stars.
This is at odds with the current state of the art in Pop III star formation
simulations, which suggest that the first stars were predominantly 100 - 500
M⊙.

Determining the nucleosynthetic imprint of PISNe on second-generation stars is
much more problematic because their metals may have differentially
contaminated new stars. This is because metals in PISNe mixed with each other
and with the surrounding IGM on much larger spatial scales, at radii where new
stars may have formed. Except in models u200 and u225 at the interface of the
O and He shells, and in model u225 at the interface of the O and Si shells,
the shells of elements in the ejecta of other PISNe expand homologously until
they sweep up their own mass in the ambient H II region, at radii of 10 - 15
pc. At this point, a reverse shock detaches from the forward shock and a
contact discontinuity forms between them (the Chevalier phase– Whalen et al.,
2008). Both the reverse shock and contact discontinuity are prone to dynamical
instabilities that will mix elements from the interior of the remnant with the
surrounding medium. Later, on scales of 100 \- 200 pc, further mixing will
occur upon collision of the remnant with the dense shell of the relic H II
region.

In this case, the elements that are taken up into new stars depends on how and
where the stars form. If metals migrate out into cosmological flows and then
fall back into the the halo via accretion and mergers on timescales of 50 -
100 Myr, they likely will be well-mixed, with all elements formed in the
explosion appearing in the new star (Greif et al., 2007; Wise & Abel, 2008b;
Greif et al., 2010). However, new stars may form at much earlier times in the
SN remnant at radii where mixing takes place. One way this could happen is if
metals at the interface between the ejecta and swept-up shell suffuse into and
cool the shell, causing it to fragment into clumps that are unstable to
gravitational collapse (e.g. Mackey et al., 2003). In this scenario, such
clumps would be enriched only by the elements residing in the relatively
narrow zone in which the stars form. This still does not explain why some
hyper metal-poor stars have such high [C, N, and O/Fe] ratios. While C and O
would be preferentially deposited on the clumps in which such stars form
because they are predominant in the outer layers of the ejecta, PISNe do not
produce enough N to account for measured abundances in these stars. However,
little iron from deep in the interior of the remnant would reach the clumps
because it is not mixed, which is consistent with the abundances measured in
EMP stars to date. The failure to uncover the characteristic odd-even
nucleosynthetic signature of PISNe predicted by Heger & Woosley (2002) metal-
poor star surveys thus far has led some to suppose that primordial stars may
not have been very massive, but it is possible that such signatures have been
masked by observational selection effects (Karlsson et al., 2008). Unlike for
low-mass stars, detailed numerical simulations that follow differential
enrichment are required to determine the true chemical imprint of PISNe on
second-generation stars.

Second, our results imply that, unlike in SN 1987A or low-mass Pop III SNe, Ni
and Fe emission lines will not appear immediately after radiation breakout
from the shock because there is not enough mixing to transport these elements
out to the photosphere of the fireball. This may be key to distinguishing
between Pop III core-collapse explosions and PISNe. The light curves of these
supernovae are characterized by a sharp intense initial transient that decays
over several hours into a dimmer extended plateau that persists for 2 - 3
months in core-collapse events and 2 - 3 years in PISNe (Fryer et al., 2010;
Whalen & Fryer, 2010). The inital pulse is powered by the thermal energy of
the shock and the plateau is energized by radioactive decay of
${}^{56}\mathrm{Ni}$ (the long life of the plateau is due to the longer
radiation diffusion timescales through the massive ejecta). Preliminary
radiation hydrodynamical calculations of PISN light curves and spectra
indicate that their peak bolometric luminosities are similar to those of Type
Ia and core-collapse SNe, making determination of the mass of the progenitor
from the magnitude of the initial transient problematic. If, however, Ni and
Fe are detected just after radiation breakout, one can be confident that the
explosion is due to a low-mass primordial star, and a Pop III IMF could begin
to be built up by samples of such detections.

Our models also demonstrate that one-dimensional radiation hydrodynamical
models are sufficient to capture most features of PISN light curves and
spectra because the mixing such calculations exclude, which would alter the
order in which emission lines appear over time, is minor. The picture for
core-collapse Pop III SNe is quite different because vigorous mixing prior to
the eruption of the shock from the surface of the star mandates its inclusion
in light curve models. Two-dimensional multigroup radiation hydrodynamical
calculations of such spectra lie within the realm of current petascale
platforms, but a less costly approach can incorporate mixing in one-
dimensional models. If most mixing in low-mass Pop III explosions occurs
before radiation breaks free from the shock, two-dimensional models such as
those in JET10 can be used to compute the distribution of elements in the
ejecta just before breakout. These explosions can then be azimuthally averaged
onto the one-dimensional grid of the light curve calculation and evolved to
compute spectra. On average, along any given line of sight out of the SN, this
method will produce emission lines in the likely order they would be observed.

These simulations, together with our previous survey of mixing and fallback in
low-mass Pop III SNe, are the first of a numerical campaign to model the
chemical enrichment of the early cosmos from its smallest relevant spatial
scales. The eventual goal of this campaign is to understand the contribution
of the first SNe to the formation of new stars and the assembly of primeval
galaxies, which will soon be probed by the James Webb Space Telescope (JWST)
and the Atacama Large Millimeter Array (ALMA). The next stage of these
numerical simulations will incorporate mixing on subparsec scales to determine
if new stars directly form in the remnants of the first supernova explosions
and follow their congregation into the first galaxies.

The authors thank Stan Woosley for helpful discussions and the use of his
KEPLER progenitor models. CCJ was supported in part by the SciDAC Program
under contract DE-FC02-06ER41438. DJW acknowledges support from the Bruce and
Astrid McWilliams Center for Cosmology at Carnegie Mellon University. Work at
LANL was done under the auspices of the National Nuclear Security
Administration of the U.S. Department of Energy at Los Alamos National
Laboratory under Contract No. DE-AC52-06NA25396. All simulations were
performed on the open cluster Coyote at Los Alamos National Laboratory.

## References

  * Abel et al. (2000) Abel, T., Bryan, G. L., & Norman, M. L. 2000, ApJ, 540, 39
  * Abel et al. (2002) —. 2002, Science, 295, 93
  * Abel et al. (2007) Abel, T., Wise, J. H., & Bryan, G. L. 2007, ApJ, 659, L87
  * Almgren et al. (2010) Almgren, A. S., Beckner, V. E., Bell, J. B., Day, M. S., Howell, L. H., Joggerst, C. C., Lijewski, M. J., Nonaka, A., Singer, M., & Zingale, M. 2010, ApJ, 715, 1221
  * Alvarez et al. (2006) Alvarez, M. A., Bromm, V., & Shapiro, P. R. 2006, ApJ, 639, 621
  * Beers & Christlieb (2005) Beers, T. C. & Christlieb, N. 2005, ARA&A, 43, 531
  * Bromm et al. (1999) Bromm, V., Coppi, P. S., & Larson, R. B. 1999, ApJ, 527, L5
  * Bromm et al. (2002) —. 2002, ApJ, 564, 23
  * Bromm et al. (2001) Bromm, V., Ferrara, A., Coppi, P. S., & Larson, R. B. 2001, MNRAS, 328, 969
  * Bromm et al. (2003) Bromm, V., Yoshida, N., & Hernquist, L. 2003, ApJ, 596, L135
  * Chen et al. (2010) Chen, K., Heger, A., & Almgren, A. 2010, ArXiv e-prints
  * Ekström et al. (2008) Ekström, S., Meynet, G., Chiappini, C., Hirschi, R., & Maeder, A. 2008, A&A, 489, 685
  * Frebel et al. (2005) Frebel, A., Aoki, W., Christlieb, N., Ando, H., Asplund, M., Barklem, P. S., Beers, T. C., Eriksson, K., Fechner, C., Fujimoto, M. Y., Honda, S., Kajino, T., Minezaki, T., Nomoto, K., Norris, J. E., Ryan, S. G., Takada-Hidai, M., Tsangarides, S., & Yoshii, Y. 2005, Nature, 434, 871
  * Fryer et al. (2007) Fryer, C. L., Hungerford, A. L., & Rockefeller, G. 2007, International Journal of Modern Physics D, 16, 941
  * Fryer et al. (2010) Fryer, C. L., Whalen, D. J., & Frey, L. 2010, ArXiv e-prints
  * Greif et al. (2010) Greif, T. H., Glover, S. C. O., Bromm, V., & Klessen, R. S. 2010, ApJ, 716, 510
  * Greif et al. (2007) Greif, T. H., Johnson, J. L., Bromm, V., & Klessen, R. S. 2007, ApJ, 670, 1
  * Hammer et al. (2010) Hammer, N. J., Janka, H., & Müller, E. 2010, ApJ, 714, 1371
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * Hirschi (2007) Hirschi, R. 2007, A&A, 461, 571
  * Iwamoto et al. (2005) Iwamoto, N., Umeda, H., Tominaga, N., Nomoto, K., & Maeda, K. 2005, Science, 309, 451
  * Joggerst et al. (2010a) Joggerst, C. C., Almgren, A., Bell, J., Heger, A., Whalen, D., & Woosley, S. E. 2010a, ApJ, 709, 11
  * Joggerst et al. (2010b) Joggerst, C. C., Almgren, A., & Woosley, S. E. 2010b, ArXiv e-prints
  * Joggerst et al. (2009) Joggerst, C. C., Woosley, S. E., & Heger, A. 2009, ApJ, 693, 1780
  * Karlsson et al. (2008) Karlsson, T., Johnson, J. L., & Bromm, V. 2008, ApJ, 679, 6
  * Kitayama & Yoshida (2005) Kitayama, T. & Yoshida, N. 2005, ApJ, 630, 675
  * Kitayama et al. (2004) Kitayama, T., Yoshida, N., Susa, H., & Umemura, M. 2004, ApJ, 613, 631
  * Mackey et al. (2003) Mackey, J., Bromm, V., & Hernquist, L. 2003, ApJ, 586, 1
  * Nakamura & Umemura (2001) Nakamura, F. & Umemura, M. 2001, ApJ, 548, 19
  * O’Shea & Norman (2007) O’Shea, B. W. & Norman, M. L. 2007, ApJ, 654, 66
  * Scannapieco et al. (2005) Scannapieco, E., Madau, P., Woosley, S., Heger, A., & Ferrara, A. 2005, ApJ, 633, 1031
  * Smith & Sigurdsson (2007) Smith, B. D. & Sigurdsson, S. 2007, ApJ, 661, L5
  * Weaver et al. (1978) Weaver, T. A., Zimmerman, G. B., & Woosley, S. E. 1978, ApJ, 225, 1021
  * Whalen et al. (2004) Whalen, D., Abel, T., & Norman, M. L. 2004, ApJ, 610, 14
  * Whalen et al. (2008) Whalen, D., van Veelen, B., O’Shea, B. W., & Norman, M. L. 2008, ApJ, 682, 49
  * Whalen & Fryer (2010) Whalen, D. J. & Fryer, C. 2010, ArXiv e-prints
  * Wise & Abel (2008a) Wise, J. H. & Abel, T. 2008a, ApJ, 684, 1
  * Wise & Abel (2008b) —. 2008b, ApJ, 685, 40
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015