We present theoretical approaches to high energy nuclear collisions in detail putting a special emphasis on technical aspects of numerical simulations. Models include relativistic hydrodynamics, Monte-Carlo implementation of kT -factorization formula, jet quenching in expanding fluids, a hadronic transport model and the Vlasov equation for colored particles. §1. IntroductionOne of the primary purposes of heavy ion physics at ultrarelativistic energies is to explore properties of strongly interacting matter at high temperature and densities, namely, the quark gluon plasma (QGP). 2), 3), 1) To this purpose, collider experiments of high energy pp, dA and AA collisions have been performed at both Relativistic Heavy Ion Collider (RHIC) in Brookhaven National Laboratory (BNL) and Large Hadron Collider (LHC) in European Organization for Nuclear Research (CERN). It is of particular importance to comprehensively understand space-time evolution of matter created in high energy heavy ion collisions so that one can extract information on the detailed properties of matter from experimental data.High energy nuclear collisions contain rich physics and exhibit many aspects of dynamics according to relevant energy and time scales. Collisions of two energetic nuclei can be viewed as those of highly coherent dense gluons. Universal behaviors of hadrons and nuclei at the high energy limit are called the color glass condensate (CGC). 4) Just after the collision, longitudinal color electric and magnetic fields between the two passing nuclei, which is called color flux tubes, 5) are produced from CGC initial conditions. Subsequent evolution of these flux tubes is called "glasma". 6) Long range rapidity correlations in the glasma 7), 8), 9) provide a natural explanation for the ridge structure observed at RHIC 10), 11) and LHC. 12) It was pointed out that instabilities of the gluon fields could play a significant role in the process of thermalization. 13), 14), 15)A model based on k T -factorization formula 20), 21), 16), 17), 18), 19), 22), 24), 23) or classical Yang-Mills approach 25), 26) reproduces the multiplicity distribution for charged hadrons at RHIC and LHC. However, initial transverse energy per particle is large compared with the experimental data. 27), 26), 24) Glasma evolution based on a 2+1 dimensional classical Yang-Mills simulation does not account for elliptic flow data at RHIC. 28) These facts suggest the necessity of inclusion of further evolution with much stronger interaction, e.g., hydrodynamical evolution. Indeed, a nearly perfect fluid picture 29), 30), 31), 32), 33) turns out to explain large elliptic flow 34) observed at RHIC 35), 36), 37) and LHC, 38), 39), 40) which leads to establishment of a new paradigm typeset using PTPT E X.cls Ver.0.9