Introduction
Neutron stars are not only the densest objects in the Universe, but in some cases they possess a magnetic field billions of times larger than the strongest magnet on Earth. These so called magnetars have magnetic fields that surpass the quantum electrodynamic threshold of $4\times10^{13} G$ and are responsible for many exotic phenomena, such as vacuum birefringence, photon splitting, and the distortion of atoms. They are invoked in order to explain the large bursts of gamma-rays and X-rays in soft-gamma repeaters or the related anomalous X-ray pulsars.
Despite the large amount of research in analytical and perturbative magnetohydrodynamics, self-consistent general relativistic solutions of the Einstein-Maxwell-Euler system are still in their infancy. Self-consistent equilibria have been obtained with only poloidal or only toroidal magnetic fields. Solutions with mixed poloidal and toroidal magnetic fields were also obtained, but with the price of greatly reducing the number of Einstein equations solved. On the other hand equilibrium solutions are not necessarily stable, and indeed, the first general relativistic MHD simulations with either purely toroidal magnetic fields or purely poloidal magnetic fields confirmed the unstable nature of these solutions predicted decades ago. Such simulations either adopt the Cowling approximation or are performed in axisymmetry.
As a first step we go beyond the previous works above by constructing rotating, magnetized equilibria with mixed ultrastrong poloidal and toroidal components and evolve them in full GRMHD in order to assess their evolutionary fate. Our initial data are constructed with the magnetized, rotating NS libraries of the COCAL code, where the $\textit{whole}$ set of Einstein-Maxwell-Euler system is solved to construct models under the assumptions of perfect conductivity, stationarity, and axisymmetry. These models are then evolved using the Illinois GRMHD code in full general relativity.
The behavior of our magnetized neutron stars can be broadly described by the following characteristics:
(i) Large radial density oscillations, especially for the rapidly rotating magnetars.
(ii) Uniform rotation is destroyed in the core of the stars, which at instances becomes counterrotating.
(iii) The NSs remain axisymmetric to a large degree.
(iv) The toroidal magnetic field is the first to become unstable.
(v) The timescale of the instability is longer for smaller toroidal magnetic field strengths, although the degree of instability is comparable in all cases. Models with the strongest toroidal B-field exhibit a density dip inside the star, are the most unstable and develop a "gear"-like shape at the NS surface.
(vi) All our models develop the "varicose" and "kink" instabilities.
