represents a stationary state of the full GR+MHD equa- tions. These evolutions are computed with a general rel- ativity code employing the generalized harmonic scheme (allowing for black hole excision). Further, no symme- try assumptions are made. High resolutions are achieved using a distributed adaptive mesh refinement infrastruc- ture. These evolutions demonstrate that our code can evolve a stable rotating star for many periods accurately. Similarly, unstable stars evolve to black holes with no evidence of any significant disk forming. Finally, we give evidence that unstable, rotating, magnetized stars rep- resent minimally unstable solutions which could serve as Type I critical solutions. In Section II we provide details about the formulation of the equations. In Section III we discuss aspects of our numerical implementation, and describe the diagnostic quantities evaluated in Section IV. In Section V we dis- cuss the initial data we use. We present our results in Section VI, and conclude in Section VII. Neutron stars can be modeled by relativistic fluids (possibly with the inclusion of magnetic fields) under the action of strong gravitational fields [26]. These systems are governed both by the Einstein equations for the ge- ometry and by the relativistic equations of magnetohy- drodynamics for the matter. We write both systems as first order hyperbolic equations. This form of the equa- tions is convenient in order to take advantage of several rigorous numerical techniques devised for such systems to ensure, at the linear level, stability of the implementa- tion. More information regarding the motivation for this approach can be found in [2, 27, 28]. By way of notation, we use letters from the beginning of the alphabet ( a , b , c ) for spacetime indices, while letters from the middle of the alphabet ( i , j , k ) range over spatial components. We adopt geometric units where c = G = 1. However, as discussed in Section III, when appropriate, we rescale the value of G to achieve improved accuracy in the con- servative to primitive variable conversion stage. The Einstein equations can be written as a system of ten nonlinear partial differential equations for the space- time metric g$_{ab}$ . The harmonic formulation of the Ein- stein equations exploits the fact that the coordinates x a can be chosen satisfying the generalized harmonic condi- tion [29, 30] where Γ a ≡ g $^{bc}$Γ$_{abc}$ are the contracted Christoffel sym- bols. The arbitrary source functions H $^{a}$( t, x $^{i}$) deter- mine the coordinate freedom of Einstein equations. The This GH formulation includes a number of constraints that must be satisfied for consistency. On one hand, Within these definitions we can write our evolution equa- tions in our GH formalism in the following way [31] We adopt a first order reduction of the second order differential equations represented in Eqs. (2). This re- duction can be achieved by introducing new independent variables related to the time and space derivatives of the fields Indices on spacetime quantities are raised and lowered with the 4-metric, g$_{ab}$ , and its inverse, while the 3-metric h$_{ij}$ and its inverse are used to raise and lower indices on spatial quantities. The spacetime can be foliated into hypersurfaces of constant coordinate time x 0 ≡ t = const. On these spacelike hypersurfaces, one defines a spatial 3-metric h$_{ij}$ = g$_{ij}$ . A vector normal to the hypersurfaces is given by n$_{a}$ ≡ −∇ $_{a}$t/ ||∇ $_{a}$t || , and coordinates defined on neigh- boring hypersurfaces can be related through the lapse, α , and shift vector, β $^{i}$. With these definitions, the space- time differential element can then be written as The matter is coupled to the geometry by means of the stress energy tensor T$_{ab}$ and its trace T ≡ g $^{ab}$T$_{ab}$ , which will be dictated by the particular model of magnetized fluid under consideration, detailed in the next subsection. original harmonic coordinates correspond to the case H $^{a}$( t, x $^{i}$) = 0, which is the choice here. The Einstein equations can be expressed in their generalized harmonic form [29], in particular 2 II. FORMULATION AND EQUATIONS OF MOTION A. The Einstein equations ∇ $^{c}$∇ $_{c}$x a = − Γ$_{a}$ = H $^{a}$( t, x $^{i}$) , (1) Q$_{ab}$ ≡ − n $^{c}$∂$_{c}$ g$_{ab}$ , D$_{iab}$ ≡ ∂$_{i}$g$_{ab}$ . (4) ∂$_{t}$g$_{ab}$ = β k D$_{kab}$ − α Q$_{ab}$, (5) k ∂$_{t}$Q$_{ab}$ = β ∂$_{k}$Q$_{ab}$ αh $^{ij}$∂$_{i}$ D$_{jab}$ − α ∂$_{a}$H$_{b}$ − α ∂$_{b}$H$_{a}$ + 2 α Γ$_{cab}$ H c + 2 α g cd ( h $^{ij}$D$_{ica}$ D$_{jdb}$ − Q$_{ca}$Q$_{db}$ − g $^{ef}$Γ$_{ace}$ Γ$_{bdf}$) − α 2 n $^{c}$n$^{d}$Q$_{cd}$ Q$_{ab}$ − α h $^{ij}$D$_{iab}$ Q$_{jc}$n c − 8 π α (2 T$_{ab}$ − g$_{ab}$T ) 2 σ$_{0}$ α [ n$_{a}$Z$_{b}$ + n$_{b}$Z$_{a}$ − g$_{ab}$n $^{c}$Z$_{c}$ ] + σ$_{1}$ β $^{i}$( D$_{iab}$ − ∂$_{i}$g$_{ab}$ ) , (6) ∂$_{t}$D$_{iab}$ = β $^{k}$∂$_{k}$ D$_{iab}$ α ∂$_{i}$Q$_{ab}$ + α 2 n $^{c}$n$^{d}$D$_{icd}$ Q$_{ab}$ + α h $^{jk}$n$^{c}$D$_{ijc}$ D$_{kab}$ − σ$_{1}$ α ( D$_{iab}$ − ∂$_{i}$g$_{ab}$ ) . (7) g $^{cd}$∂$_{cd}$ g$_{ab}$ + ∂$_{a}$H$_{b}$ + ∂$_{b}$H$_{a}$ = − 16 π ( T$_{ab}$ − T 2 g$_{ab}$ ) +2 Γ$_{cab}$ H c + 2 g $^{cd}$gef ( ∂$_{e}$g$_{ac}$ ∂$_{f}$ g$_{bd}$ − Γ$_{ace}$ Γ$_{bdf}$ ) . (2) d s 2 = g$_{ab}$ d x $^{a}$d x b = − α $^{2}$d t 2 + h$_{ij}$ ( d x i + β $^{i}$d t )( d x j + β j d t ) . (3)