{'fraction_non_alphanumeric': 0.0355564483728405, 'fraction_numerical': 0.02280945250009131, 'mean_word_length': 4.525630676084763, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rotation and Internal Structure of Population III Protostars

We analyze the cosmological simulations performed in the recent work of Greif et al. (2012), which followed the early growth and merger history of Pop III stars while resolving scales as small as 0.05 R_sol. This is the first set of cosmological simulations to self-consistently resolve the rotation and internal structure of Pop III protostars. We find that Pop III stars form under significant rotational support which is maintained for the duration of the simulations. The protostellar surfaces spin from ~50% to nearly 100% of Keplerian rotational velocity. These rotation rates persist after experiencing multiple stellar merger events. In the brief time period simulated (~ 10 yr), the protostars show little indication of convective instability, and their properties furthermore show little correlation with the properties of their host minihaloes. If Pop III protostars within this range of environments generally form with high degrees of rotational support, and if this rotational support is maintained for a sufficient amount of time, this has a number of crucial implications for Pop III evolution and nucleosynthesis, as well as the possibility for Pop III pair-instability supernovae, and the question of whether the first stars produced gamma-ray bursts.


INTRODUCTION
The formation of the first stars marked a crucial transition in the evolution of the early Universe (Bromm & Larson 2004;Loeb 2010). These stars were metal-free, and are thus also known as Population III (Pop III). They are believed to have formed at z 20 within dark matter minihaloes of mass ∼ 10 6 M⊙ (e.g. Haiman et al. 1996;Tegmark et al. 1997;Yoshida et al. 2003). Pop III stars that were sufficiently massive then formed the first Hii regions, beginning the process of reionization (e.g. Kitayama et al. 2004;Sokasian et al. 2004;Whalen et al. 2004;Alvarez et al. 2006;Johnson et al. 2007). Stars within certain mass ranges ended their lives as supernovae, thereby contributing to the initial enrichment of the intergalactic medium (IGM) with heavy elements (Madau et al. 2001;Mori et al. 2002;Bromm et al. 2003 Norman et al. 2004;Tornatore et al. 2007;Greif et al. 2007Greif et al. , 2010Wise & Abel 2008;Maio et al. 2011; recently reviewed in Karlsson et al. 2012). In particular, non-rotating primordial stars with mass 140 M⊙ < M * < 260 M⊙ are believed to have exploded as pair-instability supernovae (PISNe; Heger & Woosley 2002), releasing the entirety of their metal content into the IGM and surrounding haloes, while stars within the range 15 M⊙ < M * < 40 M⊙ ended their lives as core-collapse SNe. On the other hand, nonrotating Pop III stars with main sequence masses in the range 40 M⊙ < M * < 140 M⊙ or M * > 260 M⊙ are expected to collapse directly into black holes, thus contributing no metals to their surroundings.
The typical mass of Pop III stars is therefore a crucial factor in determining their role in the evolution of the early Universe. However, the rotation of a star can also strongly influence stellar evolution and death, (e.g. Maeder 1987; see also reviews by , 2012Langer 2012), with stronger effects for lower-metallicity stars (e.g. c 2012 RAS Brott et al. 2011). For instance, instead of a normal corecollapse SN, rapidly rotating stars may end their lives as extremely energetic hypernovae (e.g. Nomoto et al. 2003). A further consequence may be that rapid rotation lowers the minimum mass at which Pop III stars can become PISNe (Chatzopoulos & Wheeler 2012;Yoon et al. 2012), possibly to masses as low as 65 M⊙. This occurs through the process of rotationally induced mixing, which in turn leads to chemically homogeneous evolution (CHE) and a larger final oxygen core mass as compared with a non-rotating star of the same initial mass.
Theoretical studies of stellar evolution furthermore find that, compared to non-rotating stars of the same mass, the typical effective temperatures and luminosities of rotating Pop III stars are larger, and for an extended range of masses CHE shifts stellar evolutionary tracks blueward (Yoon et al. 2012). This depends sensitively, however, on how angular momentum redistribution by magnetic fields is modeled. Earlier work which does not implement the Spruit-Taylor dynamo (Spruit 2002), for instance, find that Pop III stars will not undergo CHE and will in fact end their evolution in a redder part of the Hertzsprung-Russell diagram (Ekström et al. 2008). The metal production within rotating primordial stars is generally higher, and rapidly-rotating models can lead to 14 N yields that are several orders of magnitude larger than for corresponding non-rotating or slowly-rotating models (Ekström et al. 2008;Yoon et al. 2012). However, we also note that these authors as well as Heger & Woosley (2010) somtimes find significant 14 N yields even in their nonrotating models due to chemical mixing through convection between the helium and hydrogen-burning shells.
Rotation may additionally determine a Pop III star's potential for producing a gamma-ray burst (GRB), particularly given the connection between long-duration GRBs and the death of massive stars (e.g. Woosley & Bloom 2006). This would provide a promising method of directly probing Pop III stars in their final stages, provided that Pop III GRBs occurred with sufficient frequency (Bromm & Loeb 2002Gou et al. 2004;Belczynski et al. 2007;Naoz & Bromberg 2007). The collapsar model of GRB formation requires the presence of sufficient angular momentum in the progenitor for an accretion torus to form around the remnant black hole (e.g. Woosley 1993;Lee & Ramirez-Ruiz 2006). Furthermore, to enable the escape of the accompanying relativistic jet from the star, the progenitor must also lose its hydrogen envelope (e.g. Zhang et al. 2004; but see Suwa & Ioka 2011). However, this latter requirement poses a difficulty for single-star progenitors, because removal of the hydrogen envelope will also lead to a decrease of angular momentum in the core (e.g. Spruit 2002;Heger et al. 2005;Petrovic et al. 2005). On the other hand, both conditions for collapsar GRB formation may be met in the case of close binary systems. As discussed in, e.g., MacFadyen & Woosley (1999), a binary companion can serve to remove the H envelope of the massive star, which then becomes a helium star. Subsequently, the spin-orbit tidal interactions will spin up the He star (e.g. Lee et al. 2002;Izzard et al. 2004). However, the numerical calculations of Detmers et al. (2008) question how often this scenario leads to a collapsar GRB. They find that the tidal interactions between a Wolf-Rayet (WR) star and a compact object in a binary system will more often lead to a merger event instead of a collapsing and rapidly rotating star.
Another route to GRB formation arises for Pop III stars with sufficient spin, due to the effects of rapid rotation on a star's nucleosynthesis as well as its evolution off of the main sequence (MS; e.g. Woosley & Heger 2006;Yoon et al. 2012). For instance,  and Woosley & Heger (2006) find that low-metallicity massive stars ( 20 M⊙) with rotation rates above ∼ 40-50% of their breakup speeds may undergo rotationally induced mixing and CHE, allowing the star to bypass the red giant phase and become a WR star. This may allow the star to retain enough angular momentum to become a GRB, particularly if the star also undergoes little mass loss through stellar winds, as is expected for low and zero-metallicity stars (Kudritzki 2002;Vink & de Koter 2005). Yoon et al. (2012) find that this route to GRB production through CHE applies to metal-free stars as well.
A full understanding of the nature of Pop III stars and their impact on the IGM thus requires knowledge of their typical range of spin, and recent observations lend evidence that at least some Pop III stars had significant angular momentum. For instance, it has been found that some very metal-poor stars in the Milky Way (MW) halo have an anomalous depletion of Li well below the Spite plateau (e.g. Frebel et al. 2005;Caffau et al. 2011;Bonifacio et al. 2012). This is sometimes referred to as the 'meltdown' of the Spite plateau (Aoki et al. 2009;Sbordone et al. 2010). Such Li destruction may have occurred through rotationally induced mixing (e.g. Pinsonneault et al. 2002; see also discussion and references in Asplund et al. 2006).
Further observations by Chiappini et al. (2011) also lend evidence for rapid rotation of previous generations of massive stars. They find anomalous enhancement of Ba, La, and Y within the globular cluster NGC 6522. These elements may have been produced through an enhanced s-process in rapidly rotating massive stars. Observations of N/O and C/O abundance ratios in metal-poor stars in the halo as well as damped Lyman-α systems, presented in works such as those of Spite et al. (2005) and Pettini et al. (2008), can also be more easily explained by enrichment from rapidly rotating massive stars (e.g. Chiappini et al. 2006;Hirschi 2007, see also discussion in Maeder & Meynet 2012).
Current observations of O-and B-type stars in our Galaxy reveal that massive stars can indeed be rapid rotators, with rotation rates as high as several tens of percent of break-up speed, up to over 300 km s −1 (e.g. Huang & Gies 2008;Wolff et al. 2008). Lower-metallicity massive stars such as those within the Magellanic Clouds have also been found to have faster average rotation rates than stars of higher metallicity (e.g. Hunter et al. 2008). The environments in which these stars formed, however, differs from that of Pop III stars (see, e.g., Zinnecker & Yorke 2007). While Pop III stars form in DM-dominated minihaloes, massive stars today form within non-DM-dominated giant molecular clouds. If a fraction of small Pop III stars were ejected from their minihaloes, however, their accretion would be disrupted and they would remain low-mass and long-lived (e.g. Greif et al. 2011). If these Pop III stars were later reincorporated into the Milky Way or one of its satellites, we may eventually directly observe them and their rotation rates. Identifying these as Pop III furthermore re- quires that they do not accrete too large an amount of metals from the enriched ISM (see, e.g., Frebel et al. 2009;. However, such Pop III stars have yet to be found. To examine the range of possible rotation rates for Pop III stars, and thus their potential for various spin-dependent evolutionary pathways, we therefore utilize numerical simulations, initialized in their proper cosmological context. In the work of Stacy et al. (2011), we estimated the rotation rate of Pop III stars from a cosmological simulation that resolved scales as small as 50 AU, finding these stars can attain very rapid rotation, potentially as high as their breakup speed. This was deduced from the rotational and thermal structure of the star-forming disk on the smallest resolvable scales, though the behavior of the gas as it reached the stellar surface and accreted onto the star could not be followed directly. The more recent study by Greif et al. (2012) employed simulations which reached significantly increased resolution, down to scales of 0.05 R⊙, in four different minihaloes. They found high rates of mergers between protostars, such that about half of the secondary protostars within a minihalo merge with the primary protostar. This was the first simulation to follow fragmentation and merging with this level of resolution, while previous pioneering simulations of comparable resolution were either one-dimensional (e.g. Omukai & Nishi 1998;Ripamonti et al. 2002), or threedimensional but unable to follow the gas evolution beyond the formation of the initial hydrostatic core. We analyze the Greif et al. (2012) simulations to obtain improved estimates of Pop III rotation rates and internal structure. Because these calculations resolved scales within the surface of the protostar beginning from several different cosmological realizations, we can now study the effect of mergers on the rotational structure, and the variation of spin for protostars within different host minihaloes. In Section 2 we give an overview of the numerical methodology used in the simulations. In Section 3 we present our results, assessing the main caveats in Section 4. We conclude in Section 5.  Figure 1. Note that v rad is nearly zero within R * , indicating the hydrostatic equilibrium within the newly formed protostars. Dashed lines show M turb , which reaches supersonic levels outside of Rp, but is negligible within the protostars.

Initial Setup
Our work utilizes the simulation output of Greif et al. (2012), described in detail therein. Briefly, the initial conditions for these simulations were taken from Greif et al. (2011), in which the calculations were performed with the moving-mesh code arepo (Springel 2010). The original cosmological DM simulations employed boxes of lengths 250 and 500 kpc (comoving). They were initialized with 128 3 and 256 3 particles at z = 99, assuming a Λ cold dark matter (ΛCDM) cosmology with Ωm = 0.27, Ω b = 0.046, h = 0.71, and σ8 ranging from 0.81 to 1.3. See table 1 of Greif et al. (2012) for details. The simulations are then followed until the formation of the first minihalo with virial mass greater than 5 × 10 5 M⊙.

Cut-Out and Refinement
Once the site of the first minihalo has been determined, new initial conditions at z = 99 are generated. Particles that will become part of the first minihalo, as determined by the original simulations, are replaced by 64 DM particles and 64 mesh-generating points to be used in the hydrodynamic calculations. Cells and DM particles with increasingly large distances from the high-resolution region are replaced with cells and particles of progressively higher mass. Therefore, regions farther from the target minihalo are more coarsely resolved, reducing the total number of cells and particles to be followed for these refined simulations.
After a cell in the refined simulations reaches a density of 10 9 cm −3 , the central 1 pc is extracted and used as the initial conditions for a further-refined simulation with reflective boundary conditions. These new simulations are then followed to a density of 10 19 cm −3 , after which the central 2000 AU is again extracted for further refinement, such that the final simulations resolve scales as small 0.05 R⊙.

Extraction and Identification of Protostars
Each protostar, along with its mass M * and photospheric radius Rp, was identified using the procedure described in Greif et al. (2012). In short, each new protostar was found by searching for cells with density surpassing n = 10 19 cm −3 . The densest of these cells was used as the center of the candidate protostar. If the protostellar center was located outside the radius of any previously identified protostars, it was counted as a new protostar. The protostellar boundary was defined as the edge of the photosphere, which was calcu-lated by determining the spherically averaged radius where the optical depth reached unity.
Once a candidate cell was located, the optical depth ∆τ was determined within ∼ 2 × 10 6 bins. Bins were found by dividing the area into Nang ∼ 10 4 uniformly-spaced angular regions that were further divided into N rad = 200 logarithmically-spaced radial segments between 0.01 and 10 AU, centered on the candidate cell: where j denotes the angular region, k the radial segment, ρ j,k the mass-weighted density of the bin, κ j,k the Rosseland mean opacity, and ∆r k the radial extent of the bin. Opacities were taken from Mayer & Duschl (2005). The integrated optical depth for each bin was then summed along each radius: Greif et al. (2012) next calculated the radial index kcrit at which the spherically averaged escape fraction, Comparing the solid black line with the yellow dotted lines, the protostellar density is best described by the n = 3 solution, while the temperature is better-described by the n = 1.5 solution. In panel (d), upper thin lines correspond to L eff . As expected, Lr approaches L eff near the protostellar surface Rp. Thin-lined sections of the Lr profiles represent regions where dT 4 /dr is positive, leading to negative Lr values. In these sections we thus show the modulus of Lr. Note that at these times the conditions for nuclear ignition have not yet been satisfied.
drops to βcrit = 1 − exp (−1) ≃ 0.63, which corresponds to an optical depth of unity. The photospheric radius Rp was then set equal to r k (kcrit), while M * was given by the mass enclosed within Rp.
In our current study, we furthermore define a 'hydrostatic' radius R * , defined as the distance from the protostellar center to the radial bin inside Rp which has the maximum value of κ. As will be seen in Section 3.2, the radius R * < Rp where κ reaches a peak corresponds more precisely to the location of the protostellar accretion shock.
We focus our study on four protostars, the most massive one found within each of the four minihalos simulated in Greif et al. (2012), thereby choosing the protostars most likely to later become massive stars with the greatest impact on their surroundings. This allows us to probe the early stages of rapidly accreting protostars within a range of host minihalo environments.

Rotational Profile
The significant rotational support within the protostars at the end of each simulation can be seen in Figure 1, where Figure 5. Physical and adiabatic temperature gradients with respect to distance from the center of the protostar, shown at the final snapshot time t fin of each simulation. Solid lines denote ∇ phys , and dotted lines denote ∇ ad . Dashed lines show ∇ Ω out to R * , while beyond this radius our approximation for ∇ Ω no longer applies since the gas is not undergoing solid-body rotation in these regions. Vertical dash-dotted lines show the location of the photospheric edge of the protostar at the end of the simulations, while vertical dashed lines denote the location of R * . Convection occurs where ∇ phys > ∇ ad + ∇ Ω sinθ. Thus, these protostars are generally non-convective at these early times, where the conditions for nuclear ignition have not yet been satisfied.
we show the radial profile of vrot. For each radial bin we determine its rotational velocity vrot as follows: where mi is the mass of a single hydrodynamic cell in the radial bin, vrot,x,i is the x component of cell's rotational velocity, and M is the total mass within the bin. Figure 1 furthermore displays the Keplerian velocity vKep = (GMenc/r) 1/2 , where Menc is the mass enclosed within distance r from the center of the protostar. Also shown is an example solid-body rotational velocity profile, v solid = Ωmaxr, where Ωmax = vrot,max/Rp, and vrot,max is the maximum value of vrot found within Rp. At the protostellar surface, the rotational support varies from ∼ 80% of vKep for Minihaloes 1 and 3 to > 95% for Minihaloes 2 and 4.
In panels b, c, and d, tightly bound secondary protostars are evident as extra peaks in the vrot profiles beyond the edges of the main protostars.
For comparison, the profile of radial velocity v rad can be seen in Figure 2. The large v rad outside R * signifies the radial inflow onto the protostellar surface, while v rad falls to zero inside R * due to the hydrostatic equilibrium of the gas. We additionally show the turbulent Mach number M turb over a range of radial bins in Figure 2 (see also Greif et al. 2012), defined as where cs is the sound speed of the radial bin, mi is the mass of a cell with index i contributing to the bin, and M is the total mass of the bin. Turbulence increases where there is radial inflow, and thus is approximately sonic outside of R * while becoming supersonic outside of Rp. However, it is negligible inside of R * where the gas has attained hydrostatic equilibrium.
Within the protostar, vKep roughly follows vKep ∝ r, the same as v solid , due to the nearly constant-density protostellar cores which yield Menc ∝ r 3 . The decline of vKep at radii beyond the protostar coincides with the leveling off of the enclosed mass Menc (Fig. 3), at which point vKep will scale approximately as r −1/2 (see also Clark et al. 2011b).
Along with vKep, rotation rates within the protostars also follow an approximately solid-body profile. This is apparent in Figure 3, which along with Menc also shows the radial profile of the angular velocity Ω = vrot/r, in units of km s −1 R −1 ⊙ . The protostars generally exhibit a near-constant Ω ranging between ∼ 5 × 10 −7 to 1.4 × 10 −6 s −1 . Near the protostellar surface and beyond, Ω shows the expected drop-off as roughly Ω ∝ r −3/2 . For gas within the constantdensity cores, such a profile results from vrot remaining at significant fractions of vKep throughout the protostar, because gas from the rotationally-supported disk maintains much of its rotational support as it accretes onto the protostar. Note, however, that the Ω profiles are not as flat as the inner profiles of models presented in, e.g., . Small deviation from a solid-body profile within the constant-density core (d 10 R⊙) is due to slight flattening of the protostars even at their inner radii, as well as angular momentum redistribution provided by torques from turbulent infalling gas and secondary protostars, allowing for some differential rotation such that Ω = Ω(r).

Methods of Angular Momentum Transport
Aside from these small deviations, perhaps more striking is that the protostars nevertheless maintain roughly solid body rotation even beyond the constant-density cores, nearly out to their photospheric surfaces (d ∼ 100 R⊙), where the density has declined by over an order of magnitude, or several orders of magnitude in the case of the Minihalo 1 protostar. This implies rapid angular momentum transport on scales of Rp, which is driven mainly by sonic turbulence. Figure  2 shows that between ∼ 10 and 100 R⊙ the protostellar gas transitions from non-turbulence to roughly sonic turbulence. The timescale for this hydrodynamic transport can be estimated as where Rp ∼ 100 R⊙ and cs ∼ 5 km s −1 . This yields ttrans ∼ 1 yr. Turbulent transport is therefore sufficiently rapid to account for the majority of angular momentum redistribution within the outer regions of the protostars during the time covered by the simulation. On larger scales of the protostellar disk, additional angular momentum transport occurs through gravitational torques generated by, e.g., spiral arm patterns and secondary protostars, as well as dynamical friction. This allows for protostellar mergers to occur on nearly free-fall timescales (see discussion in Greif et al. 2012). Angular momentum transport through unphysical numerical diffusion is negligible. As described in Springel (2010), the Langragian scheme used in arepo generates significantly less numerical diffusivity than corresponding Eulerian methods. For instance, when arepo was tested for its accuracy in solving the Gresho vortex problem (section 8.5 of Springel 2010), in which rotating vortex motion is applied to constant-density gas, minimal error was generated ( 1% for intermediate and high resolution cases, converging to ∼ 0.1% for the highest resolution case). These low errors were maintained even when the vortex was set in motion through the box, unlike corresponding tests with Eulerian schemes that displayed increased error due to numerical diffusivity. AREPO thus demonstrates good conservation of angular momentum and vorticity with minimal numerical diffusion effects.
Within stars that are not undergoing accretion or mass loss, redistribution of angular momentum can also occur through processes such as shear stress, meridional circulation, and convection (e.g., . Eddington-Sweet circulations occur in radiative regions of a star on the timescale tES. In a uniformly rotating star this is approximately (e.g., Maeder & Meynet 2012), where tKH is the Kelvin-Helmholtz timescale. M * , R * , and L * are the mass, radius, and luminosity of the protostar. Typical values for M * and R * at the end of this simulation are ∼ 0.5 M⊙ and 50 R⊙, while Rp is a few times larger ( 100 R⊙). The mass-weighted value of the radiative luminosity Lr within Rp, for regions where Lr is positive, ranges from ∼ 10 −3 L⊙ to 1 L⊙ (see also following section). Average Lr values within R * are smaller and range from ∼ 10 −5 L⊙ to 10 −3 L⊙ since the simulations are still in early pre-nuclear ignition stages. If we use larger photospheric values to estimate tES, 1 L⊙ and 200 R⊙, we then find typical tES values of ∼ 1000 yr, well in excess of the simulated time. For the much smaller luminosities encountered in the interior even longer times would result. Note that for stars that have contracted to the mainsequence (MS), tES will typically be several orders of magnitude longer than our estimate, e.g., tES ∼ tKH ∼ 3 × 10 7 yr for a critically rotating Sun-like star, and longer if its rotation rate is slower. Also note that our estimate for L * is several orders of magnitude lower than what might be expected in the outer regions of the star (see following section), or from the accretion luminosity Lacc: ForṀ ∼ 5 × 10 −2 M⊙ yr −1 , our typical protostar will have Lacc ∼ 4000 L⊙. These high luminosity values, however, are not representative of the radiative flux through the majority of the protostellar interior, so we use the lower values of L * quoted above.
In contrast to Eddington-Sweet circulation, shear stress and convective mixing are dynamical processes that occur on short timescales of a fraction of a year ( 0.1 yr, e.g. Zahn 1992). However, as we will describe in Section 3.3, we do not find significant convective instability in the protostars at the early times simulated. The requirement for shear instability can be expressed by the Richardson criterion: where ggrav is the gravitational acceleration, dV /dz is the change in horizontal velocity between two layers (i.e., spherical shells), and Ricrit is usually taken as equal to 1/4 (e.g. Maeder & Meynet 2012). At distances of ∼ 100 R⊙ from the protostellar center, ggrav ∼ 1 cm s −2 , n ∼ 10 17 cm −3 , and V ∼ 10 km s −1 . Taking the changes in ρ and V between the distances of 10 and 100 R⊙, we find that Ri ∼ 10 − 100. We also confirmed the accuracy of this estimate with a more precise calculation from the simulation snapshots. Thus, the density gradient strongly suppresses the shear instability. Thermal diffusivity may have an opposing destabilizing effect, but only on thermal timescales longer than followed in our simulations. In short, angular momentum is transported through turbulence and gravitational torques, and angular momentum inflow through the disk provides the protostars with substantial rotational support. These simulations confirm the prediction of Stacy et al. (2011), based on approximate modeling of the sub-grid physics, that the significant rotational support found in the large-scale ∼1000 AU starforming disk down to a few tens of AU would also persist on unresolved AU and stellar scales. Furthermore, our analysis demonstrates that the high rate of protostellar merging does not lead to significant transfer of angular momentum away from the main protostars.

Density and Temperature Profiles
In Figure 4, we show radial profiles of additional quantities for the most massive protostar in each minihalo, taken at the final simulation output. We compare the density and temperature structure of the protostar of Minihalo 1 to a range of polytropic models (yellow dotted lines in Fig. 4), where ρ(r) = ρc θ n (r), and T (r) = Tc θ(r).
θ(r) is determined through solving the Lane-Emden equation, which describes the relation between θ(r) and the dimensionless radius ǫ(r) = r/rn depending upon the polytropic index n (see, e.g., Hansen et al. 2004.) The normalization factors rn, ρc, and Tc are determined by the stellar mass and photospheric radius, M * and Rp: where ǫ1 and θ ′ 1 are the values of ǫ and dθ/dǫ at θ(ǫ1) = 0. The mass of a hydrogen atom is mH, kB is Boltzmann's constant, and µ the mean molecular weight. We take the mean molecular weight to be µ = 1.22, which is the average molecular weight of all cells with n > 5 × 10 19 cm −3 . We note that µ is not constant throughout the protostar, contributing to imperfect fits between the simulation and the polytropic solutions. However, the protostar of Minihalo 1 seems to best follow models with n ranging from 1.5 to 3. This applies to the main protostars of the other minihaloes as well, though we do not display these fits in Figure 4.

Opacity and Luminosity Profiles
The Rosseland mean opacity κ, shown in panel (c) of Fig. 4, is determined as described in Greif et al. (2012). The peak in each κ curve coincides with the protostellar accretion shock, roughly corresponding to where the temperature falls below ∼ 2 × 10 4 K and density below ∼ 10 18 cm −3 , with the opacity rapidly dropping at radii beyond the shock. The κ peaks mark where the predominant opacity contribution transitions from H bound-free (bf) absorption to H − bf absorption as the gas phase converts from ionized to neutral (e.g. Mayer & Duschl 2005). At larger radii beyond ∼ 100 R⊙, the main contribution to opacity is from H2.
Panel ( where a is the radiation constant, r is the distance from the stellar center, and dT 4 /dr is taken directly from the simulation output. The above equation describes the luminosity at a given radius using a diffusion approximation to radiative transfer and assuming that all energy is transported through the protostar by radiation. In this case the protostellar luminosity L * is equal to Lr. The simulation, however, did not include full modeling of radiative transfer, and instead assumed a local escape fraction using the Sobolev approximation. The resulting Lr profile should be interpreted only as an approximation. For comparison purposes, panel (d) of Figure 4 also shows L eff (r) = 4πr 2 σSBT 4 (r). As expected, for each protostar Lr approaches L eff (r) near r = Rp, where T approaches the effective temperature T eff . We also point out that regions of the star where dT 4 /dr is positive lead to negative values of Lr. Regions with negative Lr values correspond to the thin-lined sections of the profiles, where instead of the negative Lr values we have plotted the modulus.
The accretion flow outside of the photospheric surface Rp is subject to supersonic turbulence (see Fig. 2), which provides an additional source of heating. The values of Lr just beyond Rp are typically on the order of 1000 L⊙, similar to the value Lacc ∼ 4000 L⊙ quoted in the previous section. The main contribution to the luminosity is thus the impact of accreting gas onto the protostellar surface, though there is some contribution from the turbulence generated by the inflow as well. For M turb = v turb /cs ∼ 2, and cs ∼ 5 km s −1 , the turbulent velocity v turb is 10 km s −1 . Referring to Figure 3, the turbulent region beyond Rp (e.g., the region between 200 and 300 R⊙ for the protostar of Minihalo 2) contains a small mass of approximately 0.01 M⊙, yielding a kinetic energy of 10 43 erg. With typical radial velocities of 3 km s −1 , the timescale for gas to cross the turbulent region is ∼ 10 7 s, yielding an energy production rate of 10 36 erg s −1 , or ∼ 250 L⊙. While not negligible, this is significantly smaller than Lacc We furthermore point out that the effects of radiation pressure were not accounted for in the simulations. We expect radiation pressure from direct ionization to be a negligible effect, as the radially distended protostars have a typical effective temperature of ∼ 3000-4000 K, assuming they have surface luminosities of L * = Lacc ∼ 4000 L⊙. Note that this agrees well with the gas temperature at the photospheric surface of the protostars, ∼ 3000 -5000 K (see panel b of Fig. 4). These protostars thus do not yet ionize their surrounding gas.
For gas within the protostar, we can estimate the radiation pressure provided due to Lr using However, because the temperatures within the protostars are better-determined than our estimates for Lr, we instead approximate P rad by assuming local thermodynamic equilibrium and P rad = (1/3)aT 4 . Note that inputting this relation into Equation 17 and solving for Lr yields our Equation 16. We estimate the radiation pressure just within the photospheric surface, usually at r ∼ 100 R⊙, where typical temperatures are ∼ 4000 K, and we find P rad ∼ 1 dyne cm −2 . The radiation pressure exerted upon the outer atmosphere of the protostars is thus negligible compared to the thermal pressure, P therm = nkBT ∼ 5 × 10 4 dyne cm −2 for representative values of n = 10 17 cm −3 and T = 4000 K. Inside the protostar, typical gas temperatures and densities in the inner protostellar regions are T ∼ 3 × 10 4 K and n ∼ 3 × 10 19 cm −3 , yielding a thermal pressure of P therm ∼ 10 8 dyne cm −2 . This is again significantly greater than the radiation pressure, P rad ∼ 2000 dyne cm −2 .
We can also estimate the strength of radiation pressure due to Thomson scattering by comparing Lr within the ionized region of the protostar to the Eddington luminosity, L Edd = 4πGM * mHc/σT, where σT is the Thomson scattering cross section. Our typical protostar has L Edd ∼ 10 4 L⊙, orders of magnitude higher than Lr within the ionized gas. Radiation pressure will thus not be a significant effect until nuclear burning has commenced.

Convective Instability
In Figure 5 we estimate the regions of the protostar that are unstable to convection, determined by the relation between ∇ phys and ∇ ad . As stated by the Schwarzschild criterion, convection in non-rotating stars is expected if ∇ phys > ∇ ad , where ∇ phys = d ln T d ln P phys (18) Figure 9. Morphology of the main protostar within each minihalo, where images are oriented to show the face of the disk. Each panel has a width of 5 AU. A line of length 1 AU is also shown in each panel for comparison. Dashed white circles depict the size of the photosphere of the main protostars, while solid white circles denote R * . The color scale ranges from density of 10 12 cm −3 to 10 21 cm −3 . The diskiness and spiral structure of the accreting gas is readily apparent.
is the actual temperature gradient within the protostar, which we estimate directly from the simulation output. Furthermore, is the adiabatic variation of temperature with pressure.
To account for rotation and the restoring effect of angular momentum, we must instead use the related Solberg-Hoiland criterion for convective instability (e.g. Kippenhahn & Weigert 1990): Here ∇Ω = Hp ggravδ where HP is the pressure scale height, ggrav the gravitational acceleration, ϑ the colatitude angle from the axis of rotation, ω = rsinϑ the distance to the rotation axis, and For solid-body rotation which very roughly describes our protostars within R * , this expression for ∇Ω can be simplified to We approximate that Hp ∼ P/(ggravρ), where P is the gas pressure (see, e.g., Maeder et al. 2008). For simplicity, in our estimates we also set sinϑ equal to one. If all energy transport is through radiation, we may also write L * = Lr and ∇ phys = ∇ rad , where (24) Figure 10. Same as previous figure, except with images oriented to show the disk edge-on. Each panel has a width of 5 AU. Note the flattened structure of the protostars at these early times.
To account for the effects of rotation, we may furthermore replace Menc with (25) (von Zeipel 1924) where ρ0 is the average stellar density and is typically ∼ 10 −5 g cm −3 within R * . For our protostars the parenthetical factor in the above equation is of order unity and ranges from 0.8 to 0.9. Figure 5 shows ∇ phys , ∇ ad , and ∇Ω at the end of each simulation. Note that the estimate of ∇Ω begins to break down outside of R * , where Ω starts to deviate from a solid-body profile, so we only show ∇Ω out to this radius. Throughout most of the protostar, ∇ ad is the largest term and prevents convection within the protostars. ∇ phys and ∇Ω instead maintain values closer to zero with the exception of regions within ∼ 10 R⊙. We emphasize, however, that improved simulations, including the proper treatment of diffusive transport of radiation, will be necessary to accurately model convection within Pop III protostars. Bear-ing this caveat in mind, at the end of the simulations only one protostar, that of Minihalo 2, exhibits possible central convection. Though not shown in Figure 5, at earlier times other protostars also have short periods where the convection criterion is satisfied, giving some indication that Pop III protostars may have typically undergone periods of convection in their cores as they grew in mass. However, overall the protostars appear non-convective.
Prototars on the Hayashi track are expected to have convective zones, but if the protostar is growing then the presence of an accretion shock will return it to stability against convection (e.g. Stahler 1988b, see also Wuchterl & Klessen 2001;Wuchterl & Tscharnuter 2003). However, previous one-dimensional modeling of Pop III evolution under rapid accretion (e.g. Stahler et al. 1986;Omukai & Palla 2003;Hosokawa et al. 2010) indicate that growing protostars can still become convectively unstable after the onset of nuclear burning. In particular, Hosokawa et al. (2010) find that metal-free protostars of mass 1 M⊙ accreting through a disk may briefly have a central convective zone after deuterium burning begins.
However, this zone soon disappears while a longer-lasting surface convection zone emerges as the protostar grows and the D-burning region moves to the outer layers. For spherical accretion, however, deuterium burning and convective zones do not appear until the protostars have attained higher masses ( 10 M⊙, Omukai & Palla 2003).
Without modeling of this nuclear ignition or of radiative transfer, however, we do not expect sustained convection to appear in our simulations. We furthermore note that, despite the differences in opacity between Pop III, low-metallicity, and solar-metallicity protostars, appearance of convection zones associated with the onset of deuterium burning is expected regardless of metallicity. In all cases the main contributions to opacity are not from metals but from electron scattering and free-free absorption in stellar interiors hot enough to undergo deuterium burning, or from H − absorption and photoionization depending on the gas temperature in cooler outer regions of the star (e.g. Stahler 1988a; Durisen et al. 1989;Omukai & Palla 2003). However, the extra contribution to the electron fraction from metals does still lead to higher opacities in solar-metallicity gas as well as relatively smaller radii in Pop III stars (e.g. Alexander et al. 1983;Tumlinson & Shull 2000).

Time Evolution
Figures 6 and 7 show how various properties of the protostars evolve over time. Mass M * and radii Rp and R * of the protostars are determined as described in Greif et al. (2012) and Section 2.4. In Figure 7 we determine vrot, * similarly to the method described in Equation 4, except that the sum is performed over all hydrodynamic elements within Rp, and M is replaced with M * . When compared with vKep, * = (GM * /Rp) 1/2 , we see that each of the stars maintains a rotational velocity that is a significant fraction of its break-up speed. The protostar of Minihalo 1 (panel a of Figure 7) even appears super-Keplerian for a brief period around 4 yr. However, this is an artifact of our spherical averaging scheme, and the super-Keplerian period coincides with when the protostar is undergoing its first merger and exhibits a shape that is more flattened and bar-like than spherical. We here briefly mention that many of the secondary protostars in each minihalo also maintain very large rotation rates. The majority of secondary protostars have values of vrot, * /vKep, * ranging from ∼ 70-90% at the final simulation outputs, though some have lower values of ∼ 20-30%.
We also quantify the level of rotational support with the parameter β, which we define as the ratio of the specific angular momentum of the protostar, j * , to its Keplerian angular momentum jKep. We thus have where j * = mijx,i 2 + mijy,i 2 + mijz,i The index i refers to values for individual gas elements within the protostar, and In Figure 8, it can be seen that each of the stars remains substantially rotationally supported and maintains a relatively constant value of β ∼ 30%, with some fluctuation above and below this value. The variations in β also tend to follow the fluctuating rotational support in the surrounding disk, out to 10Rp (dotted line in Figure 8), with the exception of the Minihalo 1 protostar. In addition, the β oscillations roughly follow the fluctuations in the ratio of rotational to Keplerian velocity at the protostellar surface Rp (dashed line in Figure  8). Each protostar generally maintains vrot of 80% of vKep at its surface, demonstrating that the rapid protostellar rotation is set by the nearly Keplerian rotational speeds at the boundary between the inner disk and protostellar surface.
The rapid rotation at the inner disk boundary in turn originates from the angular momentum of the larger-scale disk. This Keplerian inner disk is likely to persist to later times (see, e.g., Stacy et al. 2011). Though these simulations only follow the initial stages of protostellar accretion, as the protostars grow by an order of magnitude their high spins may be maintained through continued accretion from a nearly Keplerian disk. This significant rotation causes the protostars to have flattened structures, as is visible in Figures 9 and 10, where we depict the density structure around the main protostars at the end of each simulation. Note that the flattening is much more significant on scales of Rp than on the smaller scales of R * . Furthermore, the main protostars of Minihalo 2 and 4 exhibit even more asymmetry within Rp than those of Minihalo 1 and 3. This is because in the former case, the protostars are undergoing mergers at the times shown, leading to greater flattening. As discussed in Lin et al. (2011), such deviations from spherical symmetry in stars undergoing disk accretion are to be expected as angular momentum is exchanged between the star and disk through gravitational torques. Similar to the stellar rotational evolution seen in their calculations, each protostar rapidly spins up over the first few Keplerian rotational periods τKep, where in our case typically τKep ∼ 1 yr. Lin et al. (2011) furthermore find a longer-term stabilization of the spin evolution, where the stellar spin may remain at ∼ 50% of the break-up value, or may undergo a slow decline, depending on the properties of the surrounding disk. Their simulations address only nonfragmenting disks. Nevertheless, even given the fragmentation and protostellar merging in our simulations, we still see a similar lack of significant evolution in β for any of the protostars.
If the very nearby secondary protostars do not rapidly merge with the main protostars, then they may also affect the longer-term evolution of the main protostars through tidal interaction and exchange of mass and angular momentum. For instance, if the binary components undergo a mass transfer through an accretion disk, the mass-gainer in the binary will be spun-up while the mass-donor spins down (see discussion and references in, e.g., Langer 2012). However, some of the rotational increase of the mass-gainer will be lost due to tidal interactions between mass-transfer events, which tend to spin down the mass-gainer and to synchronize it with the orbital motion (e.g. Zahn 1977;Petrovic et al. 2005b;Langer et al. 2003). This temporary spin down from tidal torques is visible, for instance, between 7 and 8 yr in panel b of Figure 7. However, the overall substantial rotation rate of the mass gainer is maintained through disk accretion, and this may subsequently alter its stellar evolution as previously discussed.
The close binary will additionally affect the angular momentum evolution of the circumstellar disk. In particular, the action of tidal torques between the secondary and the disk will transfer angular momentum from the disk into orbital motion, thereby enabling mass accretion onto the primary (e.g. Papaloizou & Pringle 1977). If the binary system remains wide enough, then the inflowing matter will quickly arrange into a circumbinary disk, as is already apparent in Figure 9. Spiral structure develops due to gravitational torques from the binary (e.g. Bate & Bonnell 1997), and these disk perturbations allow for further angular momentum redistribution and accretion onto the binary. Over time the binary will also tend towards an equal mass ratio between the two members. This is because angular momentum conservation causes inflowing gas to fall onto the disk at some radius beyond the more massive member, where it is more easily captured by the smaller companion.

Dependence on Minihalo Characterstics
The initial collapse of the minihalo gas from IGM to protostellar densities was discussed in detail in Greif et al. (2012). They found that the radial and rotational velocity structure of each minihalo, as well as density and temperature evolution, were consistent with previous studies (e.g. Abel et al. 2002;Yoshida et al. 2006Yoshida et al. , 2008. Furthermore, the gas collapse within each minihalo generally follows the Larson-Penston solution for isothermal and selfgravitating gas (Larson 1969;Penston 1969). The angular momentum profiles of the central few hundred solar masses of each minihalo are all the same to within a factor of a few, in addition to showing similar agreement to previous cosmological simulations (Figure 11; Abel et al. 2002;Yoshida et al. 2006;Stacy et al. 2010).
The minihaloes are not perfectly identical, however, and it is interesting to see whether variation in properties of the host minihaloes lead to any systematic variations in the properties of their protostars. Figure 12 shows the final vrot, * of each star with respect to the spin parameter, λ, and mass of the minihalo. Also denoted is the typical spin parameter of DM haloes, λ = 0.05, as measured in large-scale cosmological N-body simulations (Barnes & Efstathiou 1987;Jang-Condell & Hernquist 2001), where λ is defined as Here J, E, and M halo are the total angular momentum, energy, and mass of the halo, respectively. We have too few sample minihaloes, however, to find an obvious correlation. Furthermore, when the protostars have first formed, as well as at the end of the simulations, the rotation axes of the protostars are not aligned with that of their minihaloes. The angle between the rotation axes varies from 30 • to 160 • . This further demonstrates a lack of connection between the spins of the protostars and their host minihaloes. We additionally checked for relationships between the total number of protostellar fragments formed, N frag , and λ and vrot, * . Again, we did not simulate a sufficient number of minihaloes to find statistically significant correlations, though such an inves- Figure 11. Angular momentum profile of gas within each minihalo. Each minihalo is denoted by the same line style as in previous figures. Thin yellow lines denote angular momentum profiles found in separate cosmological simulations. Solid yellow line is taken from Stacy et al. (2010), dotted yellow line from Yoshida et al. (2006), and dashed yellow line from Abel et al. (2002). Thick black dashed line shows an approximate powerlaw fit to these profiles, j ∝ Menc. The profiles are all very similar even for a variety of cosmological realizations.
tigation with an increased number of minihaloes would be worthwhile for future work.

CAVEATS
We point out the caveat that magnetic fields, which were not included in these computations, may affect the angular momentum build-up of the protostellar disk and the subsequent protostellar rotation rate (e.g. Tan & Blackman 2004;Silk & Langer 2006;Maki & Susa 2007;Schleicher et al. 2010;Federrath et al. 2011;Sur et al. 2012;Schober et al. 2012aSchober et al. , 2012b. In addition, magnetic fields may 'disk-lock' the star to a certain rotation rate that depends upon its mass, accretion rate, magnetic field, and radius (Koenigl 1991; see also Shu et al. 1994;Matt & Pudritz 2005). This may help in estimating the longer-term angular momentum evolution of the central object. Depending on the details of such disk-locking, the large angular momentum of the disk may allow high stellar rotational velocities to be maintained beyond the phases we are modeling here, provided angular momentum is not later lost to outflows, etc. Magnetic fields may also allow for the operation of the Spruit-Taylor dynamo within the protostar (Spruit 2002). If this mechanism is active in Pop III stars, it can facilitate angular momentum transport from the core to the envelope. This allows the star to maintain solid-body rotation, spinning up its outer layers and allowing the surface to reach critical rotation speed earlier in its evolution as compared to models which do not include magnetic torques. This helps to explain the difference in the results found by, e.g., Ek-  Yoon et al. (2012). For stars that undergo SN deaths, magnetic torques may additionally lead to a rotation rate of the collapsing iron core that is over an order of magnitude smaller than those without magnetic fields . To model the formation of Pop III stars with further improved accuracy, future simulations should include the effects of magnetic fields (e.g. Turk et al. 2012).
Another important caveat to bear in mind is the short timescales followed in the simulation. The simulations followed ∼ 10 yr, which is a small fraction of the Kelvin-Helmholtz time, tKH ∼ 10 5 yr for massive stars, the typical timescale for stellar assembly. The rotational structure of a protostar only years after it has newly formed will not be the same as that when the protostar enters or later leaves the MS. Simulations that cover longer timescales will be required to see how the stellar spin as well as the inner accretion disk will evolve. Only if mass infow continues to maintain high rotational support in the inner disk, and if the inner disk continues to transfer angular momentum onto the star, can large stellar spin rates be maintained.

DISCUSSION AND CONCLUSION
We use the cosmological simulations of Greif et al. (2012), which resolve four separate minihaloes down to sub-stellar scales (0.05 R⊙), to analyze the rotation and internal structure of Pop III protostars. This is the first simulation to provide a direct view of Pop III protostellar structure from cosmological initial conditions. We find the protostars quickly develop a roughly solid-body rotation profile, while their surface rotation velocities range from ∼ 80 − 100 % of vKep at the end of each simulation. Each of the four protostars examined maintains high rotation velocities even after undergoing multiple merger events.
The protostars generally seem non-convective over the timescales simulated. However, future calculations which include a proper prescription for radiative diffusion, as op-posed to the simplified escape probablity formalism employed in Greif et al. (2012), will be necessary to properly model convection in Pop III protostars. There is also little evidence of correlation between the properties of each host minihalo and the spin of its largest protostar or the total number of protostars formed in the minihalo. More minihaloes would be necessary to derive meaningful statistics, though it would be very informative for future simulations to more thoroughly examine correlations between minihalo properties and the Pop III systems they host.
We furthermore note that future work will follow the protostellar evolution for significantly longer timescales. This will allow for a more direct determination of how the protostars evolve, and how this evolution is affected by rotation. As the protostars grow in mass and continue on to the MS, rotation could alter the protostar's life in a number of ways. Rotation rates which persist at sufficiently high speeds, for instance, may allow for mass loss through stellar winds generated at the 'ΩΓ' limit (see discussion in, e.g., Maeder & Meynet 2012), even though mass loss through line-driven winds is expected to be minimal (e.g. Kudritzki 2002). This would reduce the final mass of the star and thus may alter the stellar death it will undergo.
As described above, under the high stellar rotation rates we have inferred, the metal production during the lifetime of the star would also be generally enhanced (e.g. Ekström et al. 2008, Yoon et al. 2012. The temperature and luminosity of the star will also be altered, and possibly greatly enhanced if rotational mixing is sufficient for the star to undergo chemically homogeneous evolution (CHE; Yoon et al. 2012), though not all studies find that CHE can take place in rotating Pop III stars (e.g. Ekström et al. 2008). CHE may furthermore provide a mechanism for a Pop III star to become a WR star and eventually a GRB without being in a tight binary Woosley & Heger 2006;Yoon et al. 2012). We finally note that CHE may also lower the minimum mass at which a star will undergo a PISN death from 140 M⊙ to∼ 64 M⊙ (Chatzopoulos & Wheeler 2012;Yoon et al. 2012) The high spins seen in these protostars for a range of minihalo characteristics already demonstrate that a Pop III star of a given mass will have a number of possible evolutionary pathways depending upon its rotation rate. Computational power is progressing to the point where threedimensional simulations can begin to examine not only the mass growth but also the spin evolution of Pop III stars. In complement to numerical studies, observations of extremely metal-poor stars in the Galactic halo are key to provide constraints on the mass distribution and rotation rates of Pop III stars, their immediate progenitors (e.g. Frebel et al. 2005;Caffau et al. 2011). Much of the best evidence will come from analyzing the abundance patterns of the oldest stars found in our Milky Way or Local Group satellites, including stars within the Galactic bulge (e.g. Chiappini et al. 2011). Further important constraints on the nature of primordial stellar populations may come from abundance analysis of damped Lyman-α systems (e.g. Chiappini et al. 2006;Hirschi 2007), while any putative low-mass Pop III stars that were ejected from their minihaloes may also someday be directly detected. These continually improving numerical and obser-vational efforts will allow us to probe the early Universe in ever-greater detail.