Monte Carlo event generator MC-EKRT with saturated minijet production for initializing <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"><mml:mrow><mml:mo>(</mml:mo><mml:mn>3</mml:mn><mml:mo>+</mml:mo><mml:mn>1</mml:mn><mml:mo>)</mml:mo><mml:m

We present a novel Monte Carlo implementation of the EKRT model, MC-EKRT, for computing partonic initial states in high-energy nuclear collisions. Our new MC-EKRT event generator is based on collinearly factorized, dynamically fluctuating perturbative QCD (pQCD) minijet production, supplemented with a saturation conjecture that controls the low-<a:math xmlns:a="http://www.w3.org/1998/Math/MathML"><a:msub><a:mi>p</a:mi><a:mi>T</a:mi></a:msub></a:math> particle production. Previously, the EKRT model has been very successful in describing low-<b:math xmlns:b="http://www.w3.org/1998/Math/MathML"><b:msub><b:mi>p</b:mi><b:mi>T</b:mi></b:msub></b:math> observables at midrapidity in heavy-ion collisions at energies available at the CERN Large Hadron Collider (LHC) and the BNL Relativistic Heavy Ion Collider (RHIC). As novel features, our new MC implementation gives a full three-dimensional initial state event by event, includes dynamical minijet-multiplicity fluctuations in the saturation and particle production, introduces a new type of spatially dependent nuclear parton distribution functions, and accounts for the conservation of energy/momentum and valence-quark number. In this proof-of-principle study, we average a large set of event-by-event MC-EKRT initial conditions and compute the rapidity and centrality dependence of the charged hadron multiplicities and elliptic flow for LHC <c:math xmlns:c="http://www.w3.org/1998/Math/MathML"><c:mrow><c:mi>Pb</c:mi><c:mo>+</c:mo><c:mi>Pb</c:mi></c:mrow></c:math> and RHIC <d:math xmlns:d="http://www.w3.org/1998/Math/MathML"><d:mrow><d:mi>Au</d:mi><d:mo>+</d:mo><d:mi>Au</d:mi></d:mrow></d:math> collisions using <e:math xmlns:e="http://www.w3.org/1998/Math/MathML"><e:mrow><e:mo>(</e:mo><e:mn>3</e:mn><e:mo>+</e:mo><e:mn>1</e:mn><e:mo>)</e:mo><e:mi mathvariant="normal">D</e:mi></e:mrow></e:math> viscous fluid-dynamical evolution. Also event-by-event fluctuations and decorrelations of initial eccentricities are studied. The good agreement with the rapidity-dependent data suggests that the same saturation mechanism that has been very successful in explaining the midrapidity observables works well also at larger rapidities.