Not yet reviewed by Pith; the record is open.

We introduce our approaches, in particular the modifications of the primitive recovery procedure, to handle bare quark stars in numerical relativity simulations. Reliability and convergence of our implementation are demonstrated by evolving two triaxially rotating quark star models with different mass as well as a differentially rotating quark star model which has sufficiently large kinetic energy to be dynamically unstable. These simulations allow us to verify that our method is capable of resolving the evolution of the discontinuous surface of quark stars and possible mass ejection from them. The evolution of the triaxial deformation and the properties of the gravitational-wave emission from triaxially rotating quark stars have been also studied, together with the mass ejection of the differentially rotating case. It is found that supramassive quark stars are not likely to be ideal sources of continuous gravitational wave as the star recovers axisymmetry much faster than models with smaller mass and gravitational-wave amplitude decays rapidly in a timescale of $10\,$ms, although the instantaneous amplitude from more massive models is larger. As with the differentially rotating case, our result confirms that quark stars could experience non-axisymmetric instabilities similar to the neutron star case but with quite small degree of differential rotation, which is expected according to previous initial data studies.