Fragmentation inside atomic cooling haloes exposed to Lyman-Werner radiation
John A. Regan Turlough P. Downes

]Fragmentation inside atomic cooling haloes exposed to Lyman-Werner radiation

J.A. Regan & T.P. Downes] 
John A. ReganE-mail:john.regan@dcu.ie, Marie Skłdowska-Curie Fellow & Turlough P. Downes 
 

Centre for Astrophysics & Relativity, School of Mathematical Sciences, Dublin City University, Glasnevin, Ireland







2016

firstpage–lastpage




[
    [
    December 30, 2023
=====================





 

 
  Supermassive stars born in pristine environments in the early Universe hold the
  promise of being the seeds for the supermassive black holes observed as high
  redshift quasars shortly after the epoch of reionisation. suppression is
  thought to be crucial in order to negate normal Population III star formation and
  allow high accretion rates to drive the formation of supermassive stars.
  Only in the cases where vigorous fragmentation is avoided
  will a monolithic collapse be successful giving rise to a single massive central object.
  We investigate the number of fragmentation sites formed in collapsing atomic cooling haloes
  subject to various levels of background Lyman-Werner flux. The background Lyman-Werner flux
  manipulates the chemical properties of the gas in the collapsing halo by destroying .
  We find that only when the collapsing gas cloud shifts from the molecular
  to the atomic cooling regime is the degree of fragmentation suppressed. In our particular case
  we find that this occurs above a critical Lyman-Werner background of J ∼ 10 J_21. The
  important criterion being the transition to the atomic cooling regime rather than the actual value of
  J, which will vary locally. Once the temperature of the gas exceeds T ≳ 10^4 K and the gas
  transitions to atomic line cooling, then vigorous fragmentation is strongly suppressed. 







Cosmology: theory – large-scale structure – first stars, methods: numerical 








§ INTRODUCTION
 




Supermassive black holes (SMBHs) populate the centres of massive galaxies
<cit.> and shine as extremely luminous quasars out to high
redshift <cit.>. Their progenitors are unknown, however,
with competing theories for whether the progenitors of SMBHs were stellar mass black holes or
supermassive stars with initial masses of M_*,init≳ 10^4 <cit.>. Stellar mass black holes formed from the
remnants of the first stars (PopIII stars) must accrete at the Eddington limit for their
entire history to reach the billion mass threshold by a redshift of 7 if they are to be the seeds
of the first quasars. This scenario appears exceedingly difficult as the first stellar
mass black holes are born “starving” <cit.>
in mini haloes which have been disrupted by both the ionising radiation from PopIII stars and the
subsequent supernova explosions <cit.>.
As a direct result, investigation of supermassive star (SMS) formation as a viable alternative has
been undertaken and appears attractive <cit.>. The
larger initial seed masses of SMSs alleviate the growth requirements somewhat <cit.>.
More fundamentally however, is the fact that the SMS is born in a much larger halo compared to
the canonical PopIII case. SMS formation is thought to be possible only in pristine
atomic cooling haloes with virial temperatures T_vir≳ 10^4 K with accretion rates exceeding
Ṁ≳ 0.1 /yr <cit.>.
The larger halo leads to a much deeper potential well relative to
the minihalo case, the SMS is also expected to collapse directly into a black hole skipping the
supernovae phase and what's more the radiation spectrum expected from a SMS peaks in the infrared
rather than the UV as is the case for PopIII stars<cit.>.
Taken together, SMSs appear to offer far more favourable conditions for being the original seeds
for SMBHs. 

The conditions for achieving SMS formation are, however, not yet fully understood. It is widely
suspected that high mass accretion rates combined with high temperature gas will lead to the
formation of a SMS. This can be achieved in a number of ways. Lyman-Werner (LW) radiation from
a nearby galaxy will disrupt abundances removing a coolant, allowing a
halo to grow sufficiently without PopIII formation <cit.>.
Haloes which accrete unusually rapidly
through a succession of mergers provide another pathway<cit.> as do streaming
velocities<cit.>
left over from the decoupling of baryons and dark matter after recombination. A combination of
one or more of these scenarios is also possible. We focus here on haloes which are exposed to a
nearby source of LW radiation as it illustrates the necessary points and is straightforward to
implement. The fragmentation of gas inside a collapsing halo could potentially disrupt the
formation of a SMS if the gas were to fragment into sufficiently distant clumps to avoid the
formation of a single (or small multiple) of very massive stars at the centre of a collapsing
atomic cooling halo. 

We note that any mechanism (with or without a LW component) which suppresses PopIII star
formation and induces the same kind of thermodynamic change to the gas as an external LW field will
also likely impact fragmentation in the same way. While the fragmentation of primordial haloes
has been thoroughly investigated for the case of PopIII
formation<cit.>
it has not been as systematically probed for the case of irradiated haloes
(but see <cit.>, <cit.> (L14b), <cit.> and
<cit.> for some preliminary work). Here we investigate this scenario
more systematically, varying the intensity of the radiation and investigating the
degree of fragmentation as a function of resolution and radiation intensity. 




The goal of this study is to identify fragmentation sites in collapsing metal-free atomic cooling
haloes which are irradiated by a uniform LW background. As noted above these pristine haloes are
expected to be ideal cradles for SMS formation if the abundances can be suppressed. If
fragmentation levels are reduced in the face of impacting LW radiation, as expected, this would be
favourable for monolithic collapse of a single massive object. We probe this scenario in this work. 

The paper is laid out as follows: in <ref> we describe the 
model setup and the numerical approach used as well as introducing our star particle formulation;
in <ref> we describe the results of our numerical simulations;
in <ref> we discuss the importance of the results and present our conclusions.  

Throughout this paper we  assume a standard ΛCDM cosmology with the following parameters 
<cit.>, Ω_Λ,0  = 0.6817, 
Ω_ m,0 = 0.3183, Ω_ b,0 = 0.0463, σ_8 = 0.8347 and h = 0.6704. 
We further assume a spectral index for the primordial density fluctuations of n=0.9616.



§ NUMERICAL FRAMEWORK
 
In this study we have used the publicly available adaptive mesh refinement code
[http://enzo-project.org/] to study the fragmentation properties of gas within
haloes irradiated by a background LW field. Into we have added a new star particle type
which we have dubbed . We now describe both components.


 §.§ Enzo

[Changeset:fedb30ff370b] <cit.> is an adaptive mesh refinement code
ideally suited for simulations of the high redshift universe. Gravity in is solved using
a fast Fourier technique <cit.> which solves the Poisson equation on the root grid
at each timestep. On subgrids, the boundary
conditions are interpolated to the subgrids and the Poisson equation is then solved at each timestep.
Dark matter is represented using particles, each particle is stored on the highest refinement grid
available to it and thus the particle has the same timestep as the gas on that grid. The
particle densities are interpolated onto the grid and solved at the same time as the gas potential.
contains several hydrodynamics schemes to solve the Euler equation. We use the piecewise
parabolic method which was originally developed by <cit.> and adapted to cosmological
flows by <cit.>. The PPM solver is an explicit, higher order accurate version of
Godunov's method for ideal gas dynamics with a spatially third accurate piecewise parabolic
monotonic interpolation scheme employed. A nonlinear Riemann solver is used for shock capturing. The
method is formally second order accurate in space and time and explicitly conserves mass, linear
momentum and energy making the scheme extremely useful for following the collapse of dense
structures. 

Chemistry is an important component is following the collapse of (ideal) gas. We use the
[https://grackle.readthedocs.org/]^,[Changeset:482876c71f73] <cit.> library to follow the evolution of ten individual species:
H,  H^+,  He,  He^+,   He^++,  e^-, 
H_2, H_2^+ H^-and HeH^+. We adopt here the 26 reaction network
determined by <cit.> as the most appropriate network for solving the chemical
equations required by gas of primordial composition with no metal pollution and exposed to an external
radiation source. The network includes the most 
up-to-date rates as described in <cit.>. The cooling mechanisms
included in the model are collisional excitation cooling, collisional ionisation cooling,
recombination cooling, bremsstrahlung and Compton cooling off the CMB.



 §.§ Simulation Setup

All simulations are run within a box of 2 (comoving), 
  the root grid size is 256^3 and we employ three levels of nested grids. The grid nesting and
  initial conditions were created using MUSIC <cit.>. Within the most refined region
  (i.e. level 3) the dark matter particle mass is ∼ 103 . In order to increase further
  the dark matter resolution of our simulations we split the dark matter particles according to the
  prescription of <cit.> and as described in <cit.>. We split particles
  centered on the position of the final collapse as found from lower resolution simulations within a
  region with a comoving side length of 43.75 h^-1 kpc. Each particle is split into 13 daughter
  particles resulting in a final high resolution region with a dark matter particle mass of
  ∼ 8 . The particle splitting is done at a redshift of 40 well before the collapse of
  the target halo. Convergence testing to study the impact of lower dark matter particle masses was
  discussed in <cit.>. 

  The baryon resolution is set by the size of the grid cells, in the highest resolution region 
  this corresponds to approximately 0.48  comoving (before adaptive refinement). We vary
  the maximum refinement level (see Table <ref>) to explore the impact of resolution on
  our results. Refinement is triggered in when the refinement criteria are exceeded. The
  refinement criteria used in this work were based on three physical measurements: (1) The dark
  matter particle over-density, (2) The baryon over-density and (3) the Jeans length. The first two
  criteria introduce additional meshes when the over-density
  (Δρρ_mean) of a grid cell with respect to the mean density exceeds 8.0
  for baryons and/or DM. Furthermore, we set the MinimumMassForRefinementExponent parameter
  to -0.1 making the simulation super-Lagrangian and therefore reducing the threshold for
  refinement as higher densities are reached. For the final criteria we set the number of cells
  per Jeans length to be 32 in these runs. 

  We use 16 levels of refinement as our fiducial refinement level. This corresponds
  to a minimum cell size of Δ x ∼ 0.01 pc/h comoving (∼ 6 × 10^-3 pc
  at z ∼ 30).  We set the effective temperature of the background radiation field to
      T_eff = 30000 K. This background temperature suitably models the spectrum of a
      population of young stars <cit.>.
      The effective temperature of the background is important as the radiation temperature
        determines the dominant photo-dissociation reaction set in the irradiated halo. This in turn
        leads to a value of J_crit - the flux above which complete isothermal collapse of the
        irradiated halo is observed due to the complete suppression of . The actual value of
        J_crit depends on the nature of the source spectrum <cit.>. 

        <cit.> proposed that the J_crit needed from a given stellar
        population modelled using realistic stellar spectra can vary widely over 2-3 orders of
        magnitude. They argue that in an external pristine atomic cooling halo, DCBH formation is
        better parameterised by using a critical curve in the H_2 and H^- photo-destruction rate
        parameter space (further confirmed by <cit.>). While our choice of
        T_eff = 30000 K falls well within the
        range advocated by these studies, including a realistic source spectrum derived from
        population synthesis models is beyond the scope of the current work and remains to be
        explored in a future study. It should also be noted that <cit.> found that
        J_crit is only very weakly dependent on the nature of the source spectrum in tension with
        the results of both <cit.> and <cit.>, however they explored a
        somewhat smaller parameter space.











 §.§ Smart Stars

As the gas density increases in high density regions hydro codes, including , require a method
to convert the high density gas into stars in many cases. This is done to
allow gas which has reached the  maximum allowed refinement level of the simulation and for which
further collapse is being artificially suppressed through
artificial pressure support to be dealt with. In this case it can often be prudent
to introduce particles into the calculation to mimic the act of real star formation. It should
however be noted that in many cases there is only a loose correspondence between the numerical
particle and resolving actual star formation. Depending on the level of resolution, primarily, the
numerical particle may represent a single star or an entire cluster of stars. In our study here
the numerical particle introduced will be a proxy for star formation (fragmentation) sites
and represent only regions in which we expect star formation to be likely to occur. Star
(sometimes also known as sink) particles were initially introduced into grid codes by
<cit.>. The implementation used here also follows this model with some
modifications. We employ the same model for gas accretion as <cit.> for our star
particles using a fixed accretion radius of four cells and accreting gas based on the Bondi-Hoyle
prescription. However, while <cit.> use only two criteria to evaluate the formation
of sink particles; that the gas is at the highest refinement level and that the cell density exceeds
the Jeans density (see equation 6 in <cit.>) we also employ additional sink particle
formation criteria. Similar to the prescriptions described by <cit.> we
only form sink particles when the following criteria are met:


  * The cell is at the highest refinement level

  * The cell exceeds the Jeans Density

  * The flow around the cell is converging along each axis

  * The cooling time of the cell is less than the freefall time

  * The cell is at the minimum of the gravitational potential

We calculate the gravitational potential in a region of twice the Jeans length around the cell.
We experimented with also including the additional conditions relating to the gas boundedness and the
Jeans instability test (see <cit.> for more details). However, we found that these
additional tests were sub-dominant compared to the criteria noted above and so in the interest of
optimisation we did not include them.

Once a is formed it can accrete gas within its accretion radius (4 cells) and it can
merge with other particles. In this study we are focused only on ascertaining the degree
of fragmentation experienced in haloes as a function of Lyman-Werner radiation. The particle type contains algorithms for different accretion modes and feedback modes as well as the
ability to differentiate between normal PopIII star formation and supermassive star formation based
on accretion rates. However, we do not employ these more sophisticated accretion or feedback
algorithms here and instead leave those for an upcoming study. Here as noted previously we focus
only on identifying fragmentation sites in collapsing atomic cooling haloes.
In Table <ref> we list the simulation parameters and some of the
fundamental parameters and results. 




















§ RESULTS
 



 §.§ The onset of fragmentation as a function of LW intensity


      In Figure <ref> we show the central 1 pc of the simulations with
      LW intensities of 0 J_21[J_21 is defined as ],
      2 J_21, 4 J_21, 10 J_21, 50 J_21, 100 J_21. Each panel has an identical
maximum refinement level of 16 levels and hence an equivalent minimum cell size in comoving units.
The projection is made 650 kyrs after the first forms. 
Each realisation shows the formation of multiple fragments within the collapsing gas with the
exception of the halo exposed to a radiation level of 100 J_21. In the final halo of Figure
<ref> (i.e. bottom right panel) only a single
fragment forms at the very centre of the halo and is surrounded by a dense accretion disk
(though mild fragmentation is also observed in this realisation at higher resolution as we will see).
The mass of the fragments show a general trend of increasing mass with increasing radiation level
as expected since the Jeans mass increases in each case due to the suppression of .

In the bottom left of each panel in Figure <ref> we give the mass of the most
massive at this time. We do not directly associate this mass with the final mass of the
star that forms since further fragmentation may well, and does (see Figure <ref>),
occur below this maximum resolution scale. That is not the goal in this study. Instead we are focused
on identifying regions within the collapsing gas cloud which are likely to experience fragmentation
and likely to foster star formation. Interestingly what we see is that, initially, as the LW
intensity increases, there is a noticeable increase in the number of fragmentation sites. 

In Figure <ref> we have plotted the number of fragments as a function of the
background LW intensity for our fiducial refinement level of 16. We choose two different times
at which to examine the number of fragments. We calculate the cooling time in a given region,
t_cool, as the average cooling time in a region. The cooling time is the time after
  which the first forms that we expect an actual star to form within the fragmentation
  site. The coarser the resolution the longer this time will be. At the other extreme with
  infinite resolution the fragmentation site would be the actual star and
  t_cool would be zero. In order to account for our limited resolution we use the
  cooling time as a estimate of the collapse of the gas within the fragmentation site into stars a
  process which occurs well below the grid scale. 

We initially select each cell and calculate
its Jeans length. We then calculate the average cooling time in a sphere surrounding that cell
with radius its Jeans length. We then loop over all cells in this manner calculating the cooling
over multiple regions to find the cooling time over a given domain. The equation for the domain
cooling time, τ_cool, is then given by:

    t_cool   =  < ∑^N_cell_i = 0 t_cool, i^J_v > 
    τ_cool   =  ∑^N_cell_i = 0 t_coolρ_i∑^N_cell_i = 0ρ_i

where N_cell is the number of cells in the domain and J_V is the Jeans volume
  (i.e. the sphere surrounding the cell i of radius the jeans length of that cell) and
  t_cool, i^J_V is the average cooling time of the gas within that cooling sphere. We then
  sum over all cooling spheres (i.e. all cells) and taking the average cooling time within each sphere
  as the cooling time for that sphere.
We further weight the contribution of each cell in the calculation by it's density with the
highest density cells (those with the shortest cooling times) having a larger contribution. Low
density cells then do not have an appreciable contribution to the calculation. 

We choose two
different domains. In the first instance we choose the Jeans length of the maximum density cell
and examine all cells within a sphere of radius the Jeans length surrounding that cell
(J_L ∼ 0.03 pc). That leads to an average cooling time of approximately 80 kyrs for the
fiducial runs (i.e. those with a maximum refinement level of 16). The second region we choose
is the entire computational domain. We found through experimentation that the cooling time of
the domain saturates after the region exceeds approximately 10 kpc as outside of this region
the densities become too low and so those cells do not contribute significantly to the computation
of the average cooling time. We find the average cooling time in the total computational domain
corresponds to approximately 650 kyrs. The cooling time will be somewhat dependent on resolution
as well. As higher densities are resolved  we are able to better resolve the actual star
formation timescale and hence our proxy for the star formation timescale (here estimated using the
cooling time) drops. We fix the computation to only calculate the cooling time for the
fiducial refinement level as a benchmark. 

The first time (80 kyrs) corresponds to roughly when we expect the first actual star to form
after the formation has been triggered. This will also correspond to the time at
which we expect radiation feedback effects to begin to impact subsequent fragmentation. The second
time is the other extreme and is the maximum time it could take the surrounding gas to cool and
initiate widespread fragmentation. These two timescales should then bracket the actual fragmentation
timescale - which is unknown and will depend on specific environmental conditions
which are either not resolved by our simulation or not included. In determining the fragmentation
  scale the Jeans length of the gas must be adequately resolved. The Jeans length
  (and cooling time/length) depend sensitively on the temperature of the
  gas. In the 0 J_21 case the Jeans length is of order J_L ∼ 0.007 pc at the centre of the collapse
  while it increases to approximately J_L ∼ 0.04 pc for the cases with radiation. For our
  fiducial runs, Δ x ∼ 0.004 pc, and so we are not well resolving fragmentation
  sites for the 0 J_21 case. The number of fragmentation sites for the 0 J_21 case can
  only be therefore taken as lower limits. For the cases with radiation, which is the goal of this study,
  we are resolving the Jeans length reasonably well. Nonetheless we include the 0 J_21 case
  for the case of completeness.

For the case of
J = 4 J_21 the number of fragmentation sites is 5 after 80 kyrs. It is only when the
background is sufficiently strong that it can effect the detailed chemistry of the gas - we then see a
reduction in the number of fragmentation sites. For the highest intensity case probed here
(J = 100 J_21) the number of fragmentation sites is one both after 80 kyrs and after 650 kyrs.
If we focus on the trend from 2 J_21 to 100 J_21 we see the number of fragmentation
  sites decrease as a function of increasing radiation. While there is also a drop from
  1  J_21 to 0 J_21 this is likely due to a lack of resolution in the fiducial case and the
  fragmentation levels for the 0 J_21 are best treated as a lower limit.

In Figure <ref> we plot the radial profiles of the gas as a function of radius
for the realisation at the same refinement level (Ref16) and varying background intensities. The
plots are made just before the first forms in each realisation. The densities reached
are approximately the same in each case since the maximum refinement level is identical and the
densities are close to those at which formation is triggered. Note that the times
(redshifts) are different in each case as the LW background delays the onset of cooling within the
halo. Note also how the enclosed mass, fraction and temperature vary according to the LW
intensity. For the highest LW backgrounds the temperature is clearly higher in the central parts of
the halo. These thermodynamic differences effect the fragmentation rates as we saw in
Figure <ref>. To further emphasise the temperature and associated thermodynamic
differences we show the phase profiles for each realisation in Figure <ref>.
As the LW intensity increases we clearly see the state of the gas change reflecting the chemical
changes. As the intensity increases up to and above 10 J_21 we see the gas transition from the
molecular to the atomic cooling regime with a significant mass fraction of the gas above
approximately 2000 K. This transition and the thermodynamic consequences of it are reflected in
the decrease in the level of fragmentation we observe. The increased temperature of a significant
fraction of the gas results in a larger Jeans mass and less fragmentation for the most highly
irradiated haloes.



 §.§ The dependence of fragmentation on resolution

In Figure <ref> we plot the fragmentation as a function of resolution for
two different values of LW intensity - 4 J_21 and 100 J_21. In the left hand column the
maximum refinement level is set to 16 (Δ x ∼ 6 × 10^-3 pc), in the middle column
the maximum refinement level is increased to 20 (Δ x ∼ 4 × 10^-4 pc) and in the
right hand column the maximum refinement level is increased to 24 (Δ x ∼ 3 × 10^-5
pc). For the left and middle column the snapshots are plotted at approximately 80 kyrs (i.e. the
cooling time with a Jeans length of the gas for the Ref16 runs) while for the right hand column the
short dynamical times make extending the runs to those times prohibitive and instead we plot them
at 5 kyrs instead (cf. <cit.>). 

The projection is made using the temperature field to illustrate the changing chemical
conditions as a function of LW intensity. Note that the spatial scale varies between columns from
left to right. We centre on the most massive in each case and as the resolution increases
we are probing fragmentation deeper and deeper inside individual fragmentation sites. As the
resolution is increased we see that the 4 J_21 realisation shows increased fragmentation, with
the colder gas and subsequently smaller Jeans masses likely facilitating the fragmentation. This is in
contrast with the bottom row where the 100 J_21 intensities show comparatively little
fragmentation for similar scales. The right hand column shows the most significant fragmentation
demonstrating the most complex structures exist at the smallest scales. Again, however, the 100
J_21 intensities show over an order of magnitude less fragmentation compared
to the 4 J_21 intensities at comparable times. The prohibitively short Courant times associated
with the highest
resolution runs (even allowing for the fact that a portion of the high density gas is
tied up in particles) prevents us from evolving the highest resolution simulations for a large number
of dynamical times. However, what is clear is that the level of fragmentation is greatly reduced at
the smallest scales as the LW intensity is increased. 



 §.§ Further evolution of fragmentation

Finally in Figure <ref> we explore the time evolution of the realisations
with LW intensity of 4 J_21 and 100 J_21 respectively. It should be noted here again that we
do not include the effects of radiation or accretion feedback in these simulations and therefore we
are neglecting some important physical processes which will effect the levels of fragmentation as
star formation proceeds. Additionally, the formation of one , may impact
  physically, on the formation of a subsequent . This would be particularly important
  if the releases ionising radiation. For the case of SMS formation with accretion rates
  around 0.1 yr^-1 the bloated stellar radii of SMSs lead to effective surface
  temperatures of approximately 6000 K. With these, relatively low, effective temperatures the
  radiation is emitted predominantly in the infrared meaning that it's effect may be quite mild and
  not impact surrounding fragmentation sites. Therefore, for the high external flux cases where
  SMS formation is expected, the formation of one is not expected to negatively
  (or positively) impact subsequent SMS formation. If, however, smaller POPIII stars form
  instead then the ionisation front may become important. However, a recent study by
  <cit.> shows that the ionisation from around POPIII stars which form in close
  proximity to SMSs have small HII regions which are closely bound to the star and do not have
  significant negative impacts at least in the first 100,000 years. 
However, since we do not take these feedback effects into account here, we simply outline the evolution of fragmentation in the absence of
these effects. In Figure <ref> we see that initially (left hand column) the
particles are quite clustered (there is only particle for the 100 J_21 case). As the
simulation proceeds up to 1 Myr (after which the first particle would be expected to go
supernova) the fragmentation sites in the 4 J_21 spread out and lie approximately
1 pc apart. For the 100 J_21 case the collapse remains completely monolithic. We can also see
that for the 4 J_21 case the surrounding gas remains comparatively cold compared to the 100
J_21 case, the temperature of the gas in the 4 J_21 case is likely to stimulate fragmentation
and hence the increased fragmentation in the lower LW intensity case is not surprising. The main
finding is that the transition to the atomic cooling regime in the 100 J_21 case leads to and
sustains complete monolithic collapse.



 §.§ Comparison to previous studies in the literature

L14b and <cit.> (L15b) performed similar studies to
  that performed here. In L14b their goal was to identify the mass scale of stars formed
  when exposed to a moderate LW background and in L15b they investigated the
  accretion rates on sink particles in the case of non-isothermal collapse. We will focus our
  comparison with the first work - L14b - as it most closely corresponds to our study.
  In L14b they investigated the flux impinging on two separate haloes with LW
  backgrounds of 10 J_21,
  100 J_21 and 500 J_21. Unlike the study performed here they did not systematically
  vary their resolution and did not explore lower LW background cases (i.e. below 10 J_21).
  Their sink particle prescription was based on that of <cit.> which is broadly similar
  to the prescription we employed although their accretion radius was significantly larger than
  ours (∼32 cells compared to our 4 cells). Nonetheless our results are broadly consistent
  for the relatively high flux cases of J =  10 J_21 and 100 J_21 where our two studies overlap.
  The main differences being that we form more sinks for the very high resolution tests
  (δ_x ∼ 1 AU) but with smaller masses. 
  L14b find that one to two sinks form in their very high resolution cases. We find,
  at comparable resolution to their study, between 3 and 7 sinks have formed. Our higher values
  are due to the fact that we have a much smaller accretion radius and hence less merging
  (since the merging radius is identical to the accretion radius in both implementations). Our
  mass scales therefore also differ somewhat, L14b have masses of 7337 and 2592 for their haloes, for the J = 100 J_21 after 10 kyr. In our work we see sink particles
  with maximum masses of 1794 (after 80 kyr) and 6 (after 5 kyr). Our
  decreased masses can be attributed to our significantly smaller accretion radius.
  Choosing the correct accretion radius is non-trivial and different prescriptions for the
  accretion/merging radius will lead to different results as noted recently by <cit.>
  in a similar context. Nonetheless, the studies
  give broadly the same result with differences attributable to slightly different subgrid
  implementations. Since their study did not investigate the lower flux cases the degree of
  fragmentation suppression is not available from their study.

  Comparing to the results of <cit.>, which is for the zero flux case,we achieve
  similar resolution and similar results to their very high resolution runs. While they investigate
  the build up of a PopIII IMF we focus on simply identifying fragmentation sites rather than
  individual star formation sites. Encouragingly their fragmentation levels are similar to ours
  (comparing our 4 J_21 case to their fiducial (0 J_21) case). After approximately 5 kyrs
  the ionisation front from the first star to form has slowed the accretion rates and subsequent
  star formation in their simulations. They find that of the order of 100 sink particles have formed
  after 5 kyr compared to our 81 sink particles. It is encouraging to see that given our
    resolution levels are comparable that we achieve similar fragmentation levels for the low flux
    cases when the dominant chemistry is still due to . 






§ SUMMARY & DISCUSSION
  
We investigated here the degree of fragmentation of atomic cooling haloes as a function of both the
strength of the LW background and as a function of resolution. The degree of fragmentation within
irradiated atomic cooling haloes is an important consideration in understanding the zero age main
sequence masses of supermassive stars <cit.>
that are candidates for the supermassive black holes that exist as high redshift
quasars <cit.>. 

Haloes that are exposed to high Lyman Werner backgrounds have delayed collapse times and
larger Jeans masses enabling them to form larger objects. However, if fragmentation of the
gas is prevalent the final masses of the (super)massive stars is likely to be greatly reduced. The
thermodynamic conditions of the gas is therefore critical in estimating the initial and
(by accounting for mass loss) the final masses of these objects.

We find that haloes exposed to a LW background of J ≳ 10 J_21 experience
progressively less fragmentation than haloes exposed to lower LW backgrounds. Haloes exposed to a
LW background of J ≳ 10 J_21 have a significant fraction of their gas pushed onto the
atomic cooling track (see Figure <ref>) and are much less prone to fragmentation.
The transition to the atomic cooling track is likely then to be the most stringent
indicator of fragmentation levels. The increased temperature of the gas in the centre
limits the degree of fragmentation significantly and for the highest irradiation levels examined
here fragmentation was eradicated almost entirely. 

Comparing to the results of <cit.> we achieve similar resolution and similar
results to their very high resolution runs. While they investigate the build up of a PopIII IMF
we focus on simply identifying fragmentation sites rather than individual star formation sites.
Encouragingly their fragmentation levels are similar to ours (comparing our 4 J_21 case to
their fiducial (0 J_21) case). After approximately 5 kyrs the ionisation front from the
first star to form has slowed the accretion rates and subsequent star formation in their
simulations. Given our resolution levels are comparable it suggests that our assertion that
once the gas chemistry switches to the atomic cooling regime fragmentation is likely to be
strongly suppressed which is something they did not investigate. 

We note also that any mechanism which suppresses formation and forces a transition
to the atomic cooling track is likely to achieve a similar outcome. For example in the simulations
conducted by <cit.>, collapsing haloes are exposed to streaming velocities of up to
18 at z = 200. The streaming velocities act in a similar way to a LW background by suppressing
production and hence impacting the gas cooling. The thermodynamic state of the gas is therefore
effected in a similar way to a LW background and hence for haloes exposed to large streaming
velocities fragmentation is also likely to be suppressed. Similarly, dynamical heating effects
<cit.> combined with mild LW radiation from nearby galaxies as recently
identified by Wise et al. (2017)(in prep) would likely not be effected by vigorous fragmentation
and would in fact be ideal candidates for a single monolithic collapse.

In summary we find that haloes exposed to LW radiation experience both delayed collapse
and suffer from decreased fragmentation. Once the gas in the core exceeds approximately
2000 K, and has transitioned chemically to predominantly atomic line cooling, then fragmentation
is almost entirely suppressed and conditions becomes ideal for monolithic collapse.
In future work we will explore the impact of feedback from the collapsing
protostar and the subsequent transition to a (massive) black hole seed.  





§ ACKNOWLEDGEMENTS

J.A.R and T.P.D thank David Hubber, Peter Johansson and Bhaskar Agarwal for useful discussion
on both the sink particle implementation and direction of the study. 
J.A.R. acknowledges the support of the EU Commission through the
Marie Skłodowska-Curie Grant - “SMARTSTARS" - grant number 699941.
Computations described in this work were performed using the 
publicly-available code (http://enzo-project.org), which is the product of a collaborative 
effort of many independent scientists from numerous institutions around the world.  Their 
commitment to open science has helped make this work possible. The freely available astrophysical 
analysis code YT <cit.> was used to construct numerous plots within this paper. The authors 
would like to extend their gratitude to Matt Turk et al. for an excellent software package.
J.A.R. would like to thank Lydia Heck and all of the support staff involved with Durham's COSMA4
and DiRAC's COSMA5 systems for their technical support. This work was supported by the Science
and Technology Facilities Council (grant numbers ST/L00075X/1 and RF040365). This work used the
DiRAC Data Centric system at Durham University,  operated  by  the  Institute  for  Computational
Cosmology on behalf of the STFC DiRAC HPC Facility  (www.dirac.ac.uk). This equipment was funded
by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1,
and STFC DiRAC Operations grant ST/K003267/1 and Durham University.  DiRAC is part of the
National E-Infrastructure. The authors also wish to acknowledge the SFI/HEA Irish Centre for
High-End Computing (ICHEC) for the provision of computational facilities and support. Finally, we
thank the anonymous referee for comments which improved both the quality and clarity of the
manuscript. 






78


[Abel et al.(2000)Abel, Bryan & Norman]Abel_2000
Abel T., Bryan G. L., Norman M. L., 2000, , 540, 39

[Agarwal & Khochfar(2015)]Agarwal_2015a
Agarwal B., Khochfar S., 2015, , 446, 160

[Agarwal et al.(2012)Agarwal, Khochfar, Johnson, Neistein,
  Dalla Vecchia & Livio]Agarwal_2012
Agarwal B., Khochfar S., Johnson J. L., Neistein E., Dalla Vecchia
  C., Livio M., 2012, , 425, 2854

[Agarwal et al.(2013)Agarwal, Davis, Khochfar, Natarajan &
  Dunlop]Agarwal_2013
Agarwal B., Davis A. J., Khochfar S., Natarajan P., Dunlop J. S.,
  2013, , 432, 3438

[Agarwal et al.(2014)Agarwal, Dalla Vecchia, Johnson,
  Khochfar & Paardekooper]Agarwal_2014b
Agarwal B., Dalla Vecchia C., Johnson J. L., Khochfar S.,
  Paardekooper J. P., 2014, , 443, 648

[Agarwal et al.(2016)Agarwal, Smith, Glover, Natarajan &
  Khochfar]Agarwal_2015b
Agarwal B., Smith B., Glover S., Natarajan P., Khochfar S., 2016,
  , 459, 4209

[Alvarez et al.(2009)Alvarez, Wise & Abel]Alvarez_2009
Alvarez M. A., Wise J. H., Abel T., 2009, , 701, L133

[Bate & Burkert(1997)]Bate_1997
Bate M. R., Burkert A., 1997, , 288, 1060

[Becerra et al.(2017)Becerra, Marinacci, Inayoshi, Bromm &
  Hernquist]Becerra_2017
Becerra F., Marinacci F., Inayoshi K., Bromm V., Hernquist L. E.,
  2017, ArXiv e-prints:1702.03941

[Begelman et al.(2006)Begelman, Volonteri &
  Rees]Begelman_2006
Begelman M. C., Volonteri M., Rees M. J., 2006, , 370, 289

[Berger & Oliger(1984)]Berger_1984
Berger M. J., Oliger J., 1984, Journal of Computational Physics, 53, 484

[Boss & Bodenheimer(1979)]Boss_1979
Boss A. P., Bodenheimer P., 1979, , 234, 289

[Bromm et al.(1999)Bromm, Coppi & Larson]Bromm_1999
Bromm V., Coppi P. S., Larson R. B., 1999, , 527, L5

[Bryan et al.(1995)Bryan, Norman, Stone, Cen &
  Ostriker]Bryan_1995
Bryan G. L., Norman M. L., Stone J. M., Cen R., Ostriker J. P., 1995,
  Computer Physics Communications, 89, 149

[Bryan et al.(2014)Bryan, Norman, O'Shea, Abel, Wise,
  Turk & The Enzo Collaboration]Enzo_2014
Bryan G. L., Norman M. L., O'Shea B. W., Abel T., Wise J. H., Turk
  M. J., The Enzo Collaboration, 2014, , 211, 19

[Burkert & Bodenheimer(1993)]Burkert_1993
Burkert A., Bodenheimer P., 1993, , 264, 798

[Burkert et al.(1997)Burkert, Bate &
  Bodenheimer]Burkert_1997
Burkert A., Bate M. R., Bodenheimer P., 1997, , 289, 497

[Chon et al.(2017)Chon, Hosokawa & Yoshida]Chon_2017b
Chon S., Hosokawa T., Yoshida N., 2017, ArXiv e-prints:1711.05262

[Clark et al.(2011)Clark, Glover, Smith, Greif, Klessen &
  Bromm]Clark_2011a
Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen
  R. S., Bromm V., 2011, Science, 331, 1040

[Coppola et al.(2011)Coppola, Longo, Capitelli, Palla &
  Galli]Coppola_2011
Coppola C. M., Longo S., Capitelli M., Palla F., Galli D., 2011,
  , 193, 7

[Coppola et al.(2012)Coppola, D'Introno, Galli, Tennyson &
  Longo]Coppola_2012
Coppola C. M., D'Introno R., Galli D., Tennyson J., Longo S., 2012,
  , 199, 16

[Fan et al.(2006)Fan, Carilli & Keating]Fan_2006
Fan X., Carilli C. L., Keating B., 2006, , 44, 415

[Federrath et al.(2010)Federrath, Banerjee, Clark &
  Klessen]Federrath_2010
Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010, , 713,
  269

[Glover(2015a)]Glover_2015a
Glover S. C. O., 2015a, , 451, 2082

[Glover(2015b)]Glover_2015b
Glover S. C. O., 2015b, , 453, 2901

[Glover & Abel(2008)]GloverAbel_2008
Glover S. C. O., Abel T., 2008, , 388, 1627

[Glover & Jappsen(2007)]GloverJappsen_2007
Glover S. C. O., Jappsen A. K., 2007, , 666, 1

[Glover & Savin(2009)]GloverSavin_2009
Glover S. C. O., Savin D. W., 2009, , 393, 911

[Haemmerlé et al.(2017)Haemmerlé, Woods, Klessen,
  Heger & Whalen]Haemmerle_2017
Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J.,
  2017, ArXiv e-prints:1705.09301

[Hahn & Abel(2011)]Hahn_2011
Hahn O., Abel T., 2011, , 415, 2101

[Haiman(2006)]Haiman_2006
Haiman Z., 2006, New Astronomy Review, 50, 672

[Hockney & Eastwood(1988)]Hockney_1988
Hockney R. W., Eastwood J. W., 1988, Computer simulation using particles.
  Bristol: Hilger, 1988

[Hosokawa et al.(2013a)Hosokawa, Omukai &
  Yorke]Hosokawa_2012
Hosokawa T., Omukai K., Yorke H. W., 2013a, , 778, 178

[Hosokawa et al.(2013b)Hosokawa, Yorke, Inayoshi,
  Omukai & Yoshida]Hosokawa_2013
Hosokawa T., Yorke H. W., Inayoshi K., Omukai K., Yoshida N.,
  2013b, , 778, 178

[Hosokawa et al.(2016)Hosokawa, Hirano, Kuiper, Yorke,
  Omukai & Yoshida]Hosokawa_2015
Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida
  N., 2016, , 824, 119

[Inayoshi & Haiman(2014)]Inayoshi_2014b
Inayoshi K., Haiman Z., 2014, , 445, 1549

[Jeon et al.(2014)Jeon, Pawlik, Bromm &
  Milosavljević]Jeon_2014
Jeon M., Pawlik A. H., Bromm V., Milosavljević M., 2014, ,
  440, 3778

[Johnson et al.(2011)Johnson, Khochfar, Greif &
  Durier]Johnson_2011
Johnson J. L., Khochfar S., Greif T. H., Durier F., 2011, , 410,
  919

[Kitsionas & Whitworth(2002)]Kitsionas_2002
Kitsionas S., Whitworth A. P., 2002, , 330, 129

[Kormendy & Ho(2013)]Kormendy_2013
Kormendy J., Ho L. C., 2013, , 51, 511

[Krumholz et al.(2004)Krumholz, McKee &
  Klein]Krumholz_2004
Krumholz M. R., McKee C. F., Klein R. I., 2004, , 611, 399

[Latif & Volonteri(2015)]Latif_2015b
Latif M. A., Volonteri M., 2015, , 452, 1026

[Latif et al.(2014)Latif, Schleicher, Bovino, Grassi &
  Spaans]Latif_2014b
Latif M. A., Schleicher D. R. G., Bovino S., Grassi T., Spaans M.,
  2014, , 792, 78

[Latif et al.(2015)Latif, Bovino, Grassi, Schleicher &
  Spaans]Latif_2015
Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M.,
  2015, , 446, 3163

[Latif et al.(2016)Latif, Schleicher & Hartwig]Latif_2016a
Latif M. A., Schleicher D. R. G., Hartwig T., 2016, , 458, 233

[Loeb & Rasio(1994)]Loeb_1994
Loeb A., Rasio F. A., 1994, , 432, 52

[Madau & Rees(2001)]Madau_2001
Madau P., Rees M. J., 2001, , 551, L27

[Milosavljević et al.(2009)Milosavljević, Couch &
  Bromm]Milosavljevic_2009
Milosavljević M., Couch S. M., Bromm V., 2009, , 696, L146

[Mortlock et al.(2011)]Mortlock_2011
Mortlock D. J. et al., 2011, , 474, 616

[Planck Collaboration et al.(2014)]Planck_2014
Planck Collaboration et al., 2014, , 571, A16

[Regan & Haehnelt(2009a)]Regan_2009b
Regan J. A., Haehnelt M. G., 2009a, , 396, 343

[Regan & Haehnelt(2009b)]Regan_2009
Regan J. A., Haehnelt M. G., 2009b, , 393, 858

[Regan et al.(2014a)Regan, Johansson &
  Haehnelt]Regan_2014a
Regan J. A., Johansson P. H., Haehnelt M. G., 2014a, ,
  439, 1160

[Regan et al.(2014b)Regan, Johansson &
  Wise]Regan_2014b
Regan J. A., Johansson P. H., Wise J. H., 2014b, , 795,
  137

[Regan et al.(2015)Regan, Johansson & Wise]Regan_2015
Regan J. A., Johansson P. H., Wise J. H., 2015, , 449, 3766

[Regan et al.(2016)Regan, Johansson & Wise]Regan_2016a
Regan J. A., Johansson P. H., Wise J. H., 2016, 

[Regan et al.(2017)Regan, Visbal, Wise, , Haiman,
  Johansson & Bryan]Regan_2017
Regan J. A., Visbal E., Wise J. H., , Haiman Z., Johansson P. H.,
  Bryan G. L., 2017, Nature Astronomy, 1, 0075

[Safranek-Shrader et al.(2012)Safranek-Shrader, Agarwal,
  Federrath, Dubey, Milosavljević & Bromm]Safranek-Shrader_2012
Safranek-Shrader C., Agarwal M., Federrath C., Dubey A.,
  Milosavljević M., Bromm V., 2012, , 426, 1159

[Safranek-Shrader et al.(2016)Safranek-Shrader, Montgomery,
  Milosavljević & Bromm]Safranek-Shrader_2016
Safranek-Shrader C., Montgomery M. H., Milosavljević M., Bromm V.,
  2016, , 455, 3288

[Schauer et al.(2017)Schauer, Regan, Glover &
  Klessen]Schauer_2017
Schauer A. T. P., Regan J., Glover S. C. O., Klessen R. S., 2017,
  , 471, 4878

[Schleicher et al.(2013)Schleicher, Palla, Ferrara, Galli &
  Latif]Schleicher_2013
Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013,
  , 558, A59

[Shang et al.(2010)Shang, Bryan & Haiman]Shang_2010
Shang C., Bryan G. L., Haiman Z., 2010, , 402, 1249

[Shu(1977)]Shu_1977
Shu F. H., 1977, , 214, 488

[Smith et al.(2017)]Grackle
Smith B. D. et al., 2017, , 466, 2217

[Stacy et al.(2016)Stacy, Bromm & Lee]Stacy_2016
Stacy A., Bromm V., Lee A. T., 2016, , 462, 1307

[Sugimura et al.(2014)Sugimura, Omukai &
  Inoue]Sugimura_2014
Sugimura K., Omukai K., Inoue A. K., 2014, , 445, 544

[Tanaka & Haiman(2009)]Tanaka_2008
Tanaka T., Haiman Z., 2009, , 696, 1798

[Tanaka & Li(2014)]Tanaka_2014
Tanaka T. L., Li M., 2014, , 439, 1092

[Turk et al.(2011)Turk, Smith, Oishi, Skory, Skillman,
  Abel & Norman]YT
Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W.,
  Abel T., Norman M. L., 2011, , 192, 9

[Volonteri(2010)]Volonteri_2010a
Volonteri M., 2010, , 18, 279

[Wang et al.(2010)Wang, Li, Abel & Nakamura]Wang_2010
Wang P., Li Z. Y., Abel T., Nakamura F., 2010, , 709, 27

[Whalen et al.(2004)Whalen, Abel & Norman]Whalen_2004
Whalen D., Abel T., Norman M. L., 2004, , 610, 14

[Wise et al.(2008)Wise, Turk & Abel]Wise_2008a
Wise J. H., Turk M. J., Abel T., 2008, , 682, 745

[Wolcott-Green & Haiman(2012)]WolcottGreen_2012
Wolcott-Green J., Haiman Z., 2012, , 425, L51

[Wolcott-Green et al.(2017)Wolcott-Green, Haiman &
  Bryan]WolcottGreen_2017
Wolcott-Green J., Haiman Z., Bryan G. L., 2017, , 469, 3329

[Woods et al.(2017)Woods, Heger, Whalen, Haemmerlé &
  Klessen]Woods_2017
Woods T. E., Heger A., Whalen D. J., Haemmerlé L., Klessen R. S.,
  2017, , 842, L6

[Wu et al.(2015)]Wu_2015
Wu X. B. et al., 2015, , 518, 512

[Yoshida et al.(2003)Yoshida, Abel, Hernquist &
  Sugiyama]Yoshida_2003a
Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003, , 592, 645








§ SMARTSTAR TESTS

As part of the integration of our particle into we rigorously tested the sink
particle technique used in this work. We tested the implementation using both the Singular
Isothermal Sphere test <cit.> and by analysing the results from a Rotating Cloud test
<cit.>.



 §.§ Collapse of a Singular Isothermal Sphere

The collapse of a singular isothermal sphere with a ρ(r) ∝ r^-2 density profile
produces a constant flux of mass through spherical shells i.e. a constant accretion rate.
We compute the collapse of a singular isothermal sphere to test whether our model of sink
particle accretion reflects this collapse behaviour. We follow the work of <cit.> in
conducting this numerical test. The sphere is set to have a truncation radius of
R = 5 × 10^16 cm, a density at this radius of ρ(R) = 3.82 × 10^-18 g cm^-3
and therefore a total mass of 3.02 . The sphere is initialized at a temperature of 10 K,
the mean molecular weight of the gas is set to 3.0 resulting in a sound speed of 0.166 . These
values lead to a large instability parameter for the sphere of A = 29.3 with
A = 4 π G ρ(R) R^2/c_s^2 <cit.>. We ran the test at a maximum refinement level of
10 leading to a minimal grid spacing of Δ x ∼ 3.05 × 10^-5 pc. No radiative cooling
or chemical reactions are calculated during the collapse. 

In the upper panel of Figure <ref> we show the mass, in solar masses,
as a function of time.
The particle mass increases monotonically as it accretes matter eventually accreting the
entire envelope (3.02 ). The slope of the line is 1.5 × 10^-4 /yr in
excellent agreement with the analytical solution predicted by <cit.>.
In the lower panel of Figure <ref> we show the mass accretion rate during this time. The test predicts
a constant accretion rate of  1.5 × 10^-4 and this is confirmed by our test. The bump
at approximately 17000 yrs corresponds to the particle accreting the last of the matter and
at this point the accretion rate fluctuates mildly before falling to zero. 

We also examine the behaviour of the gas itself and the mass that is accreted by the
particle. In Figure <ref> we see the density profile, radial velocity profile
and the mass accretion through spherical shells as a function of radius. Overlaid on each plot are
the analytical expressions for the collapse taken from <cit.>. The agreement in all three
profiles is excellent with mild deviations only appearing once the the entire envelope has been
accreted at T ∼ 17000 yrs. 











 §.§ Rotating Cloud Core Fragmentation Test

Finally, we also analyse the behaviour of the particles when the gas is subject to
rotation and collapse. This is the classic <cit.> test. It is a standard test for
examining fragmentation in hydrodynamical simulations <cit.>.
We employ a setup
very similar to that of <cit.>, <cit.> and <cit.>.
The cloud is initialised with the following parameters: radius R = 5 × 10^16 cm, constant
density  ρ_0 = 3.82 × 10^-18 g cm^-3, mass = 1 , angular velocity,
Ω = 7.2 × 10^-13 rad s^-1 (ratio of rotational to gravitational energy
β = 0.16), sound speed = 0.166 (ratio of thermal to gravitational energy α = 0.26),
global free-fall time t_ff = 1.075 × 10^12 s =  3.41 × 10^4 yr. The cloud is then
subject to a density perturbation with an m = 2 mode: ρ = ρ_0[1 + 0.1 Cos(2ϕ)], where
ϕ is the azimuthal angle. We use an isothermal equation of state throughout the test. 

The initial rotation (β = 0.16) forces the gas cloud to collapse to a disk with a
central bar due to the m = 2 initial density perturbation. particles then form at the
ends of the bar once the conditions for particle formation are met. The test was run at a maximum
refinement level of 5 with a minimum cell size of 16.32 AU (2.44 × 10^14 cm). 
The evolution of the system is shown in Figure <ref>. The figures show the
initial formation of the disk and bar structure face on. Two particles form at either
end of the bar at the same time due to the inherent symmetry of the system. A third fragment then
forms at the centre of the cloud core. Subsequently the two end fragments move along the
bar and merge with the central fragment. The results of this test are
in excellent agreement with similar tests carried out by <cit.>.  





