  We present a systematic study of the cosmic variance that existed in the
formation of first stars and galaxies. We focus on the cosmic variance induced
by the large-scale density and velocity environment engraved at the epoch of
recombination. The density environment is predominantly determined by the
dark-matter overdensity, and the velocity environment by the dark matter-baryon
streaming velocity. Toward this end, we introduce a new cosmological initial
condition generator BCCOMICS, which solves the quasi-linear evolution of
small-scale perturbations under the large-scale density and streaming-velocity
environment and generates the initial condition for dark matter and baryons, as
either particles or grid data at a specific redshift. We also describe a scheme
to simulate the formation of first galaxies inside density peaks and voids,
where a local environment is treated as a separate universe. The resulting
cosmic variance in the minihalo number density and the amount of cooling mass
are presented as an application. Density peaks become a site for enhanced
formation of first galaxies, which compete with the negative effect from the
dark matter-baryon streaming velocity on structure formation.
