ner_tags: O O O O O O O O O O O O O B-Model I-Model O O O O O B-Model I-Model I-Model I-Model I-Model O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O B-ComputingFacility I-ComputingFacility O B-ComputingFacility I-ComputingFacility O O O O O B-Software O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O B-Formula I-Formula I-Formula I-Formula O O O O B-Formula I-Formula I-Formula I-Formula I-Formula O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O B-Formula I-Formula I-Formula I-Formula O B-Citation I-Citation I-Citation I-Citation O B-Citation I-Citation I-Citation I-Citation O O O O O O O O O O O O O O O O O O O O O O O O B-Software O O O O O O O O O O O B-CelestialObject O O O O O O B-Formula I-Formula I-Formula I-Formula I-Formula O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O B-Formula I-Formula I-Formula I-Formula O O O O O O O O B-Model I-Model O O O O O O O O O O B-Model O O O O O O O O O O B-Formula I-Formula I-Formula O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O B-Citation I-Citation I-Citation O O O O O O O O O O O O O O O O O O O O B-Formula I-Formula I-Formula O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O B-Software O O O
tokens: each particle by decomposing the gravitational forces into a long-range term, computed from Particle-Mesh methods, and short scale interactions from nearest neighbours using Tree methods. The code is MPI parallel and can be used with periodic boundary conditions in the comoving frame. The scaling of the code with the number of particles is , so it can handle a large number of particles with reasonable computational resources. The simulations were run in the MareNostrum supercomputer at BSC-CNS 4 using a modified version of GADGET-2 in which the free parameters of the code were chosen such that the force computation was accurate enough for PBH binaries to form and survive without being disrupted by pure numerical errors. We took the gravitational softening length ε = 0.001 pc and a time-step of Δ t = 10 −6 in log ( a ). These parameters have a significant impact in the orbit stability, so we ensured that our choices are optimal and that a further reduction of both ε and Δ t does not yield any significant improvement in orbit stability. Previous analytic estimates suggest that the first PBH pairs begin to appear approximately at the radiation-matter energy equipartition epoch ( z ∼ 10 3 ; Nakamura et al. 1997 ; Sasaki et al. 2016 ) or, depending on the PBH fraction, even earlier – during the radiation-dominated epoch. Therefore, we had to modify the cosmological equations in the GADGET-2 code to be valid during the radiation dominated epoch of the Universe since we start our simulations at z init = 10 5 . As was discussed earlier, in this work we focus on the formation of PBH pairs for different fractions of PBH to the total DM component. So, when f PBH 1.0, we have to split the DM into two types: ‘PBH DM’ and ‘normal DM’ (see Section 1 ). As it is shown in Section 3 , in our simulations virialized structures start to form as early as z ≈ 10 3 , which could be caused by the initial Poisson distribution of PBHs. This differs significantly from the case of standard ΛCDM model, where the first virialized structures start to form at z ≃ 40. Therefore, we assume that the ‘normal DM’ behaves as an uniform medium in our simulations. Nevertheless, we will discuss later why the formation of virialized structures is important and how it affects the behaviour of PBH pairs (see Section 3.3 ). In order to quantify the range of redshifts in which this approximation is applicable, we calculate the fraction of virialized normal DM using the analytical Press–Schechter formalism for haloes more massive than 30 M ⊙ (see Press Schechter 1974 ). We find that the fraction of virialized structures made of normal DM is less than 50 per cent at z > 10. (see Fig. 1 ). Figure 1. The fraction of virialized structures is made of normal DM for haloes more massive than 30 M ⊙ , as a function of redshift z . Thus, the presence of normal DM influences the evolution of PBH DM only through its contribution to the computation of the Hubble function in the GADGET-2 code: For our