In the case of asymptotically high baryon density, the system is a color superconducor. This was proposed by Frautschi [15] and Barrois [16]. Based on the Bardeen, Cooper, and Schrieffer (BCS) theory [17], because there is a weak attractive interaction in the color antitriplet channel, the system is unstable with respect to the for- mation of particle-particle Cooper-pair condensate in the momentum space. Detailed numerical calculations of color superconducting gaps were firstly carried out by Bailin and Love [18]. They concluded that the one-gluon exchange induces gaps on the order of 1 MeV at several times of nuclear matter density. This small gap has little effect on cold dense quark matter, thus the investigation of cold quark matter lay dormant for several decades. It was only revived recently when it was found that the color superconducting gap can be of the order of 100 MeV [19], which is two orders larger than early perturbative estimates in Ref. [18]. For this reason, the topic of color superconductivity stirred a lot of interest in recent years. For review articles on the color superconductivity, see for example, Refs. [20]. roughly a constant energy density ϵ$_{cr}$ ≃ 1GeV / fm 3 [14]. In this article, I will focus on recent progress of sQGP created at RHIC and color superconducting phase struc- ture at intermediate baryon density regime. The outline of this article is as follows: I will give a brief overview on the discovery of sQGP at RHIC in Sec. II. Then introduce the status of the color superconducting phase especially the gapless color superconducting phase in Sec. III. At last, I will give a brief outlook in Sec. IV. with u µ the flow velocity, ϵ, p the energy density and pressure density, respectively. In the Navier-Stokes hydrodynamics, the energy mo- mentum tensor decomposes into ideal and dissipative parts as Studying Quantum chromodynamics (QCD) phase transition and properties of hot quark matter at high temperature has been the main target of heavy ion col- lision experiments at the Relativistic Heavy Ion collider (RHIC), the forthcoming Large Hadron Collider (LHC) and FAIR at GSI. The deconfined quark-gluon plasma, if it can be cre- ated through heavy-ion collisions, is an intermediate state and cannot be measured directly. In experiment, the detector can only measure the freeze-out hadrons. In order to extract the property of the intermediate state, hydrodynamics is often used to simulate the evolution of the fluid. The hydrodynamical equations of motion are the local conservation laws of energy-momentum and net charge In ideal hydrodynamics, the energy-momentum tensor takes the form of Where △ µν = g µν − u $^{µ}$u$^{ν}$, ∇ µ = △ $^{µν}$∂$_{ν}$ , η, ζ are the shear viscosity and bulk viscosity, respectively. It was expected that deconfined quark matter formed at high temperature should behave like a gas of weakly interacting quark-gluon plasma (wQGP). The perturba- tive QCD calculation gives a large shear viscosity in the wQGP with η/s ≃ 0 . 8 for α$_{s}$ = 0 . 3 [21]. Therefore, it turned out as a surprise that the RHIC data of el- liptic flow v$_{2}$ can be described very well by requiring a very small shear viscosity over entropy density ratio η/s [22, 23]. Lattice QCD calculation confirmed that η/s for the purely gluonic plasma is rather small and in the range of 0 . 1 − 0 . 2 [24]. It is now believed that the system created at RHIC is a strongly coupled quark-gluon plasma (sQGP) and be- haves like a nearly ”perfect” fluid [25, 26]. The AdS/CFT duality gives a lower bound η/s = 1 / 4 π [27]. Therefore, it is conjectured that the sQGP created at RHIC might be the most perfect fluid observed in nature. However, a perfect fluid should have both vanishing shear and bulk viscosities. The perturbative QCD calculation gives ζ/s = 0 . 02 α 2 s for 0 . 06 < α$_{s}$ < 0 . 3 [28]. In the hydrodynamic simula- tion used to describe the evolution of the fireball created at RHIC, the bulk viscosity ζ has often been neglected. II. STRONGLY INTERACTING QUARK-GLUON PLASMA(SQGP) A. Discovery of sQGP at RHIC FIG. 2: Elliptic flow v 2 as a function of p T for different values of Γ s /τ 0 . The figure taken from Ref.[23]. 2 (GeV) T p 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 ) T (p 2 v 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 6.8 Fm (16-24% Central) b   = 0.1 o / s  = 0.2 o / s  = 0 o / s STAR Data with T µν N S = T µν ideal + τ $^{µν}$, (3) τ µν = η ( ∇ $^{µ}$uν + ∇ $^{ν}$uµ − 2 3 △ $^{µν}$∇ $_{α}$u$^{α}$+ ζ △ $^{µν}$∇ $_{α}$u$^{α}$. (4) ∂$_{µ}$T µν = 0 , ∂$_{µ}$N µ c = 0 . (1) T µν ideal = ( ϵ + p ) u $^{µ}$uν − pg $^{µν}$, (2)