BCCOMICS - Baryon CDM COsMological Initial Condition generator for Small scales.

Small-scale fluctuations in the early universe, even at redshifts as high as z=200, are found to be strongly affected by large-scale density and streaming-velocity environments. BCCOMICS is an initial condition generator that allows the study of structure formation inside a simulation box of < 4 comoving Mpc, where the simulation box can have non-zero overdensity (Δ) and streaming velocity (V_cb = mean velocity of CDM - mean velocity of baryon) as its environmental condition. This allows for the study of cosmic variance of e.g., first star formation under varying large-scale environments.

Currently, it only supports enzo. We are inviting contributors to help port this to Gadget, RAMSES, and many other fabulous N-body+hydro simulation codes.

The main code is composed of two parts.

1. Realization of large-scale fluctuations (at a length resolution of 4 Mpc comoving) and transfer function calculation (bccomics_setup.m):
   Because the perturbation theory studying the dual impact from large-scale Δ and V_cb on small-scale fluctuations is relatively new (Ahn 2016), transfer functions from the usual linear Boltzmann solvers, such as CAMB, do not provide the level of accuracy of this new theory. Even when Δ variance is not considered, the suppression of high-k (around k~100/Mpc) modes due to V_cb (Tseliakhovich & Hirata 2010) is not reflected in the transfer function from CAMB. "bccomics_setup.m" makes a realization of environmental variables, and once the user chooses a specific "patch" of generically non-zero Δ and V_cb, it calculates and records the transfer function.

2. Realization of small-scale fluctuations on a selected patch (bccomics.m):
   The user is asked again to choose one from a previously-calculated set of patches, and then the corresponding transfer function is read in and used to generate 3D data of CDM and baryons. Currently, the following data are generated:
   a. CDM particle positions (cpos1, cpos2, cpos3)
   b. CDM particle velocities (vc1, vc2, vc3)
   c. baryon grid density (db)
   d. baryon grid velocities (vb1, vb2, vb3)
   e. baryon grid thermal energy (etherm)
   f. baryon grid kinetic+thermal energy (etot)
   g. baryon particle positions (bpos1, bpos2, bpos3) (optional, for SPH)
   h. baryon particle velocities (vpb1, vpb2, vpb3) (optional, for SPH)
   i. baryon particle thermal energy (eptherm) (optional, for SPH)
   Data is written in Enzo's internal unit system.
   If using MATLAB, initial conditions can be saved directly to Enzo's HDF5 format. Saving to HDF5 is not supported in OCTAVE, but a Python script for converting the native output to Enzo format is provided.
   Items g, h, and i can be written when baryonparticleflag=true and particlevelocity_accuracyflag=true in params_patch.m.

3. (OCTAVE-only) Conversion of binary data from step (2) into enzo-usable initial conditions, by "convert_enzo.py" using python+h5py.

Installation of BCCOMICS:
Either (1) clone this repo (git clone https://github.com/KJ-Ahn/BCCOMICS.git), or (2) download as a zip and extract its contents. Click the green "Clone or download" button and choose whichever suits you. No further installation process is required.

BCCOMICS has two main scripts that run on either MATLAB(R) or gnu OCTAVE. gnu OCTAVE is easily installed with its dependency by package managers (e.g. "apt", "rpm", "synaptic", "flatpak", ...) in usual linux distributions. Ask your system administrator for installation on a shared unix machine, unless you have a large quota.

To install a rather recent version of OCTAVE on old linux distributions (or to install it privately on a shared machine without root priviledge), you might want to try installing Anaconda, and then do conda install -c conda-forge octave for installing a recent OCTAVE with its dependencey without headache (As of 05/10/2019, the current Anaconda version of octave fails on some systems. Try instead Homebrew). Flatpak is another option, now even favored by GNU, and flatpak Octave version supercedes other binary packages like .deb, .rpm, etc. For example on Ubuntu 16.04, as of 03/16/2020, you can install Octave 5.2.0 with flatpak.

For MATLAB, in addition to the main program, it is preferred (but not essential) that following additional packages be installed.

• Image Processing Toolbox (optional; if uninstalled BCCOMICS will use padarray.m under BCCOMICS/mfiles_for_matlab/ImagProc_Tool)
• Statistics and Machine Learning Toolbox (optional; if uninstalled BCCOMICS will use functions under BCCOMICS/mfiles_for_matlab/Stat_Tool)

For gnu OCTAVE, in addition to the main program, following additional packages better be installed for the optimal performance. Use your linux distribution's package installer (e.g. sudo apt-get install octave-image for Ubuntu). It is OK not to install following "octave-*" packages if you want, BCCOMICS is shipped with necessary function files of these packages anyways.

• octave-image (optional; if uninstalled BCCOMICS will use padarray.m under BCCOMICS/mfiles_for_octave/)
• octave-statistics (optional; if uninstalled BCCOMICS will use functions under BCCOMICS/statistics-1.3.0/)
• octave-odepkg (optional; if uninstalled BCCOMICS will use octave functions or those under BCCOMICS/odepkg-0.8.5/)
• python
• h5py

You can install "octave-image" mentioned above alternatively by typing pkg install image inside Octave, but then you should remember to load the package by typing in pkg load image to use e.g. padarray fundtion. Check out octave-forge webpage for further reference.

BCCOMICS needs transfer function output files of CAMB (as provided in BCCOMICS/sample/CAMB_for_mode_finding), and also a redshift-ionization output file from RECFAST (as provided as the file BCCOMICS/sample/output_recfast). So in general, you need to install

• CAMB
• RECFAST

Best explained with an example. Let's assume that BCCOMICS is installed at /home/kjahn/BCCOMICS, whose sub-directories are src, sample, etc. (Or, assume BCCOMICS installed at c:\Documents\BCCOMICS on a Windows system)
Below $ is a linux command prompt, >> is either OCTAVE's or MATLAB's command prompt.

First, you need a work directory under which "params.m" and "params_patch.m" exist (You should stick to this naming convention!!), and also a cosmology parameter file (as provided as BCCOMICS/sample/LCDM.m). The name of the work directory can be anything. Inside OCTAVE/MATLAB, you need to go to this directory by the cd command. You also need to add the source path by the addpath commmand.
The current example "params.m" generates 151³ unigrid patches inside a (604 Mpc)³ volume with bccomics_setup.m. The example "params_patch.m" generates initial conditions of CDM and baryons with 64³ resolution with bccomics.m.
For your own setup, do a recursive copy of the sample directory to e.g. my_params, and modify "params.m" and "params_patch.m" which are both self-explanatory. The cosmology parameter file (e.g. LCDM.m) should carefully reflect parameters you used for running CAMB and RECFAST.
If you want to use the pre-generated (by bccomics_setup.m) gaussian random seed ("gaussseed.matbin"; stick to this naming convention!!), place it under setupdir specified in params.m.

$ octave

>> cd '/home/kjahn/BCCOMICS/sample'
>> addpath('/home/kjahn/BCCOMICS/src')

>> cd 'c:\Documents\BCCOMICS\sample'
>> addpath('c:\Documents\BCCOMICS\src')

You will be asked to choose a patch with your desired CDM overdensity and V_cb (absolute value).

>> bccomics_setup
... Several message outputs ...

Standard deviation of CDM overdensities (sDc) is 0.0041819
Choose CDM overdensity environment:
Input 0 for mean, 1 for overdense, 2 for underdense:

Well, if you are interested in overdense patch, enter 1

Input 0 for mean, 1 for overdense, 2 for underdense:1
What multiple of sDc away from the mean overdensity, 0? Example: for Delta_c = +1.5*sDc, Enter 1.5
Enter a floating-point number:

If you want 2sigma density environment for your patch, enter 2

Enter a floating-point number:2
CDM overdensity chosen: Delta_c = 2*sDc = 0.0083638
---------------------------------------
RMS of Vbc (rmsV) at z = 1000 is 27.528 km/s
Peak of Vbc in Maxwell-Boltzmann distribution is 22.477 km/s
Choose Vbc environment at z = 1000
Enter Vbc at z = 1000 in units of km/s: 

If you are interested in typical ones, from above V_cb=30 km/s is a good choice, so enter 30

Enter Vbc at z = 1000 in units of km/s: 30
311 patches out of total 3.44295e+06 patches satisfy your chosen condition with 1% margin.
----------------One best matching patch is being found----------------
Wanted Delta_c = 0.0083638; Selected patch's Delta_c = 0.0083689
Wanted Vcb = 30 km/s; Selected patch's Vcb = 30.011 km/s
----------------Integrating----------------
1th wavenumber out of 100 is being handled.
* 1th angle outta 21 is being handled.
* 2th angle outta 21 is being handled.
...

Then you should wait until all wavenumbers are calculated. Be patient, it takes a while.

Necessary parameters for bccomics is in "params_patch.m". It is mostly self-explanatory. Just a bit more explanation on gaussian random seed: for the final initial condition, you need to either (a) use preexisting seed or (b) generate a new seed. If (a) is the case, olseedflag = true, and you should specify diroldseed and Noldseed. The gaussian random seed is named as e.g. subgaussseed512.matbin if Noldseed=512, and the file has Noldseed*Noldseed*(Noldseed/2+1) complex values.

After setting "params_patch.m", start octave or matlab, go to the work directory, add path of the src directory, and then run bccomics:

$ octave  

>> cd '/home/kjahn/BCCOMICS/sample'
>> addpath('/home/kjahn/BCCOMICS/src')
>> bccomics

You will be asked to choose one among those patches that you have chosen in more-than-one bccomics_setup runs.

Patches ordered in calculation time, from oldest(top) to newest(bottom)
-----------------------------------------------------------------------
Patch #  ix  iy  iz  Deltac/sigma(Deltac)  V_cb(km/s)  at z=1000
  1     118 102 120      1.000e+00          3.000e+01
  2      75  42 120      1.998e+00          5.003e+01
  3      40  84  79      2.966e-04          2.299e+01
Choose a patch of your interest; default is 3 if you just hit Enter below.
Enter your choice (patch #):

Assume you want Patch # 2 above, which is 2sigma Density peak with V_cb=50.03 km/s. Then enter 2:

Enter your choice (patch #):2
Patch # 2 chosen.
----- Interpolating transfer function -----
----- Convolving transfer function with random number -----
----- Calculating CDM position x -----
... (more messages)
*********** bccomics successfully ended ************
>>

After all runs are successuful, you will have a sub-directory under ICdir, named e.g. as 1.00Mpch_64_ic75_jc42_kc120. Under this direcotyr, you will have initial conditions.

Matlab can produce enzo initial conditions directly in step (3). But Octave still lacks the functionality to do so. The initial conditions from Octave will be files named "cpos1", "cpos2", ..., "vc1", ..., "db", "vb1", ..., "ether", "etot". Go to the initial condition directory, and run the provided python script:

$ cd /home/kjahn/BCCOMICS/ICs/1.00Mpch_64_ic75_jc42_kc120
$ python /home/kjahn/BCCOMICS/src_converter/convert2enzo.py

Then you will get enzo-readable HDF5 files named "ParticlePositions", "ParticleVelocities", "GridDensity", "GridVelocities", "GasThermalSpecEnergy", and "GasTotalSpecEnergy".

Because cosmology simulation by enzo needs typically a periodic boundary condition, simulating overdense(underdense) patches requires fooling enzo or any other simulation codes (see section 3 of Ahn & Smith). For enzo, this is done by "enzo_patchcosmo.m". You first need to provide file zglobal.dat under the IC directory (e.g. /home/kjahn/BCCOMICS/ICs/1.00Mpch_64_ic75_jc42_kc120) listing "global redshifts" for enzo output dumps:

$ pwd
kjahn@machine:/home/kjahn/BCCOMICS/ICs/1.00Mpch_64_ic75_jc42_kc120
$ cat zglobal.dat
200
100
90
80
70
60
50
20
10
$

You may try zglobal.m in the source directory for easy creation of zglobal.dat to your taste.

Then, run octave(matlab) in this directory, and run "enzo_patchcosmo.m". It will then generate file enzoparam_part.enzo, which you should copy & paste to your final enzo parameter file (e.z. localcosmo.enzo). You'd better keep zglobal.dat, so that you can later match a "local redshift" in enzoparam_part.enzo to its corresponding "global redshift" in zglobal.dat. Take a look at enzoparam_part.enzo, and you will see non-zero curvature term, and also seemigly odd cosmological parameters and "local redshifts". But this is how you fool enzo (and similarly other simulation codes) to simulate overdense(underdense) patches.

$ octave
>> cd '/home/kjahn/BCCOMICS/ICs/1.00Mpch_64_ic75_jc42_kc120'
>> addpath('/home/kjahn/BCCOMICS/src')
>> enzo_patchcosmo
plotflag = 0
THflag = 0
OWRTflag =  1
mmw =  1.2195
rhocrit0 =  9.2834e-30
a_i =  0.0049751
TimeUnits =  2.3939e+14
************ Local parameters found and written ********
Contents of enzoparm_part.enzo should replace parameters in your full enzo parameter file.
*********************************************************
>>

Here is a specific example. Say you are at an IC directory named /home/kjahn/BCCOMICS/ICs/1.00Mpch_64_ic75_jc42_kc120. This patch is a 2-sigma density peak (see (2) above), as can be roughly checked out from icc_Dc_Db_Thc_Thb_Vcb1_Vcb2_Vcb3_Vcb_DT.dat as well:

$ cat icc_Dc_Db_Thc_Thb_Vcb1_Vcb2_Vcb3_Vcb_DT.dat
  75   42  120 3.191772e-02 1.762073e-02 -3.194254e-03 -2.902584e-03 2.005724e+00 -9.677555e+00 -1.796974e+00 1.004525e+01 2.651701e-03

showing Dc=3.191772e-02 at the redshift (z=200) of the initial condition, which has evolved from the value Dc= 0.0083638 at z=1000. Running enzo_patchcosmo, the output file enzoparam_part.enzo reads:

$ cat enzoparam_part.enzo
CosmologySimulationOmegaBaryonNow        = 0.603016
CosmologySimulationOmegaCDMNow           = 3.138966

CosmologyOmegaMatterNow    = 3.741982
CosmologyOmegaLambdaNow    = 0.003767
CosmologyOmegaRadiationNow = 0.015079
CosmologyHubbleConstantNow = 9.745857
CosmologyComovingBoxSize   = 1.018051   // Mpc/h
CosmologyInitialRedshift   = 13.763798
CosmologyFinalRedshift     = 0.000000

CosmologyOutputRedshift[0] = 13.763798
CosmologyOutputRedshift[1] = 10.125176
CosmologyOutputRedshift[2] = 6.489630
CosmologyOutputRedshift[3] = 5.763133
CosmologyOutputRedshift[4] = 5.036844
CosmologyOutputRedshift[5] = 4.310549
CosmologyOutputRedshift[6] = 3.584862
CosmologyOutputRedshift[7] = 2.859749
CosmologyOutputRedshift[8] = 0.695545
CosmologyOutputRedshift[9] = 0.000000

The "local" redshifts listed here correspond to those in zlobal.dat. The final redshift in zglobal.dat will always end up as "local" redshift 0 (CosmologyOutputRedshift[9] in the example above). Local redshift of 0 is the "present", and the cosmological parameters at "present", which are usually required for cosmological N-body/hydro codes, are also listed in this file; e.g. CosmologySimulationOmegaCDMNow = 3.138966 is the Omega_CDM at local present and CosmologyHubbleConstantNow = 9.745857 is h (Hubble constant in units of 100 km/s/Mpc) at local present (global redshift is 10 in the example). Omega_curvature at local present will be automatically calculated inside Enzo, which is 1-CosmologyOmegaMatterNow-CosmologyOmegaLambdaNow-CosmologyOmegaRadiationNow.

Then, copy & paste these quantities into an enzo parameter file, e.g. Unigrid.enzo you can find under sample directory, and of course these quantities should replace those already in Unigrid.enzo. Then, run enzo with this file: enzo Unigrid.enzo. Any result at local redshift e.g. 0.695545 is the result at global redshift z=20 in the above example. See Section 3 of Ahn & Smith for more details.

If you are to identify cosmological halos by e.g. the friends-of-friends algorithm in an overdense region, the linking length should differ from the "global" value and will be redshift-dependent. See Section 3 of Ahn & Smith for more details.

If enzo_patchcosmo fails to run, the redshifts in zglobal.dat is likely to contain those redshifts after the patch turns around (you know a closed universe expands, stops, and turns around to collapse, right?). In this case, try shifting the final redshift in zglobal.dat to higher redshift and try again.

• Porting to Gadget, RAMSES, etc. WE WELCOME ANYONE WHO IS INTERESTED IN THIS EFFORT. DROP US A LINE, WE'LL ADD YOU AS A CONTRIBUTOR.
• Small-wavenumber tuning: For very small wavenumbers, the number of modes are too small to correctly reproduce average power spectrum. With unlucky gaussain random seed, this might cause unwated finte-box effect. May have to restrict randomness of power spectrum amplitude (not phase though) of smallest k modes.
• Nested grid: Not yet. Again, we welcome contributors on this side.

If you use BCCOMICS in your work, please cite as:

...the BCCOMICS initial conditions generator (Ahn & Smith 2018).

BIBTEX:

@ARTICLE{0004-637X-869-1-76,
  author={Kyungjin Ahn and Britton D. Smith},
  title={Formation of First Galaxies inside Density Peaks and Voids under the Influence of Dark Matter–Baryon Streaming Velocity. I. Initial Condition and Simulation Scheme},
  journal={The Astrophysical Journal},
  volume={869},
  number={1},
  pages={76},
  url={http://stacks.iop.org/0004-637X/869/i=1/a=76},
  year={2018}
}