  Neutron stars harbour extremely powerful magnetic fields, leading to their
shape being deformed. Their magnetic deformation depends both on the geometry -
and strength - of their internal magnetic field and on their composition,
encoded by the equation of state. However, both the details of the internal
magnetic structure and the equation of state of the innermost part of neutron
stars are mostly unkown. We performed a study of numerical models of
magnetised, static, axisymmetric neutron stars in general relativity and in one
of its most promising extensions, scalar-tensor theories. We did so by using
several realistic equations of state currently allowed by observational and
nuclear physics constraints, considering also those for strange quark stars. We
show that it is possible to find simple relations among the magnetic
deformation of a neutron star, its Komar mass, and its circumferential radius
in the case of purely poloidal and purely toroidal magnetic configurations
satisfying the equilibrium criterion in the Bernoulli formalism. These
relations are quasi-universal, in the sense that they mostly do not depend on
the equation of state. Our results, being formulated in terms of potentially
observable quantities, could help to understand the magnetic properties of
neutron stars interiors and the detectability of continuous gravitational waves
by isolated neutron stars, independently of their equation of state. In the
case of scalar-tensor theories, these relations depend also on the scalar
charge of the neutron stars, thus potentially providing a new way to set
constraints on the theory of gravitation.
