© ESO 2020

1. Introduction

The recent observation of gravitational and electromagnetic radiation coming from the merger of a binary neutron star system (Abbott et al. 2017a) has given us a new opportunity to test the theory of general relativity (GR) in the strong-field regime (Will 2014), beyond the vacuum case of binary black hole (BH) mergers (Abbott et al. 2016), and to probe the physics of compact objects in unprecedented detail (Abbott et al. 2017b,c), fostering a renewed interest in neutron stars (NSs) as possible probes of new gravitational physics.

Indeed, it is well known that our understanding of the Universe lacks an explanation for what is called the “dark sector”. While a possible solution is to assume the existence of dark matter and dark energy (Trimble 1987; Peebles & Ratra 2003), a different approach is to consider the possibility that GR is not the definitive theory of gravity. Moreover, on a more theoretical basis, a consistent quantum theory of GR does not yet exist (Bars & Pope 1989; Deser 2000). This led to the development of theoretical frameworks [e.g. the hypothesis of Strings (Green et al. 1988)] which try to give an explanation of the fundamental interactions which is different from the mainstream one, leading to a modification, in the low energy limit, of GR itself.

Many attempts have been made to extend GR to account for such issues, giving rise to many alternative theories of gravity (Capozziello & de Laurentis 2011). Among the most studied are: f(R) theories (Buchdahl 1970; De Felice & Tsujikawa 2010), where deviations from GR are introduced by modifying the functional dependence of the gravitational Lagrangian on the Ricci scalar R; Gauss-Bonnet gravity (Lovelock 1971), which increases the dimensionality of the spacetime; scalar-tensor theories (STTs) – to some extent equivalent to f(R) theories (Sotiriou 2006) – which modify gravity with respect to GR replacing the gravitational constant G with a dynamical scalar field. STTs have been widely studied in the past (Brans & Dicke 1961; Nordtvedt 1970; Wagoner 1970; Matsuda & Nariai 1973; Damour & Esposito-Farèse 1993, 1996; Novak 1998b; Fujii & Maeda 2003; Faraoni 2004; Shibata et al. 2014; Langlois et al. 2018; Gong et al. 2018; Quiros 2019; Zhang et al. 2019) and are among the most promising alternatives to GR. This is due to a number of reasons: they are the most simple extensions of GR (Sect. 1.2 Papantonopoulos 2015); they are predicted to be the low-energy limit of some possible theories of Quantum Gravity (Damour et al. 2002); most of them respect the weak equivalence principle (WEP) – that is they are metric theories of gravity (Will 2014) – which has been extremely well tested (Touboul et al. 2017). They also seem to be free of some of the pathologies affecting other extensions of GR (De Felice et al. 2006; De Felice & Tanaka 2010; Bertolami & Páramos 2016). On the other hand, STTs violate the strong equivalence principle (SEP), which means that tests using self-gravitating bodies are ideal to constrain them (Barausse 2017).

The foundations of STTs were laid by Brans & Dicke (1961) in a seminal paper, in which the authors modified the Einstein-Hilbert action of GR attempting to bring it in conformity with Mach’s principle by replacing the gravitational constant G by a scalar field non-minimally coupled to the spacetime metric, giving birth to the Jordan-Fierz-Brans-Dicke theory (BD). Unfortunately, observational tests in the Solar System seem to have proved BD wrong, unless its only parameter is precisely fine tuned, in contrast with the principle of naturalness (Schärer et al. 2014). In this scenario, the study of NSs is especially important, because since the first work on massless mono-scalar STTs (Damour & Esposito-Farèse 1993), a non-perturbative strong field effect has been predicted, allowing the scalar field to exponentially grow in magnitude inside compact material objects. Even generalisations of STTs to massive scalar fields and other gravitational theories have been shown to be subject to a similar phenomenon (Salgado et al. 1998; Ramazanoğlu & Pretorius 2016; Ramazanoğlu 2017; Silva et al. 2018; Andreou et al. 2019). Scalarisation can happen in various contexts: binary systems of merging NSs can undergo a “dynamical scalarisation” process (Barausse et al. 2013), in which the initially non-scalarised NSs become scalarised once they get closer to each other; again, in a binary NS system, one scalarised star can prompt an “induced scalarisation” on its non-scalarised companion (Barausse et al. 2013); or even in an isolated NS system, where “spontaneous scalarisation” can develop (this was the first discovered non-perturbative strong field effect in STTs, Damour & Esposito-Farèse 1993). The importance of scalarisation is that STTs which include such effects predict strong deviations from GR only inside compact objects, while allowing the tight observational constraints in the weak-gravity regime to be fulfilled (Shao et al. 2017). As of today, the strongest limit on the strength of spontaneous scalarisation for massless STTs comes from observations of pulsars in binary systems, in particular in systems characterised by a large mass difference between the two stars, where STTs predict the emission of dipole scalar waves, potentially observable in the dynamics of the inspiral (Freire et al. 2012; Will 2014; Shao et al. 2017; Anderson et al. 2019). These, however, are binaries with large separations and the constraints do not apply in the case of screening (Yazadjiev et al. 2016; Doneva & Yazadjiev 2016).

Scalarisation modifies the relation between the mass and radius of the NS and its central density. In general, scalarised NSs have larger radii and higher maximum masses than the corresponding GR solutions computed with the same equation of state (EoS). Moreover, scalarisation is more effective at higher compactness. The presence of a strong scalar charge could, in principle, have important consequences on the phenomenology of NSs, even if many of these effects might be degenerate with the EoS. A different dependence of the mass and radius from the central density could lead to appreciable changes in the thermal evolution of NSs (Dohi et al. 2020), given the dependence of many cooling processes on the density itself (Yakovlev et al. 2005). Changes in radii could potentially be observable in the distribution function of millisecond pulsars (Papitto et al. 2014). The same holds for the distribution of NS masses, and the expected maximum mass (the recent measure of a 13 km radius for a 1.44 M_⊙ NS by NICER, Miller et al. 2019, suggests larger NS radii than previously thought, Özel & Freire 2016). Spontaneous scalarisation might impact the dynamics and evolution of the post merger remnant of binary NS coalescence (Raithel et al. 2018; Abbott et al. 2017). Indeed, there is some observational evidence suggesting the presence of long lived NSs powering the X-ray afterglow of Short-GRBs (Rowlinson et al. 2013), suggesting values of the maximum NS mass ≳2.2 M_⊙ (Gao et al. 2016; Margalit & Metzger 2017). Scalar fields can affect the deformability of NSs (Doneva et al. 2013, 2018), leaving an imprint in the pre-merger inspiral, and in the spin-down history of millisecond proto-magnetar as possible engines of GRBs (Dall’Osso et al. 2009). Scalarised NSs differ in the frequency of their normal modes (Sotani & Kokkotas 2005). On top of this STTs predicts also a new scalar wave emission, potentially detectable with future gravitational waves (GWs) observatories (Gerosa et al. 2016; Hagihara et al. 2020).

So far, only non-magnetised models of NSs have been studied in STTs in the full non-linear regime (see e.g. Suvorov 2018 for a perturbative approach to the magnetised scenario). Most of them focus on static (Damour & Esposito-Farèse 1993; Harada 1998; Novak 1998a; Taniguchi et al. 2015; Anderson & Yunes 2019; Doneva & Yazadjiev 2020) or slowly rotating (Damour & Esposito-Farèse 1996; Sotani 2012; Pani & Berti 2014; Silva et al. 2015) stars, while recently some work has been done for rapidly (Doneva et al. 2013; Doneva & Yazadjiev 2016; Pappas et al. 2019) and differentially (Doneva et al. 2018) rotating models. NSs have also been studied beyond the massless limit, and in the presence of a screening potential (Doneva & Yazadjiev 2016, 2020; Yazadjiev et al. 2016; Brax et al. 2017; Staykov et al. 2018, 2019). However, NSs are known to contain extremely powerful magnetic fields, inferred to be in the range 10⁸⁻¹²G for normal pulsars and up to 10¹⁶G at the surface of magnetars, while newly formed proto-NSs are hypothesised to store magnetic fields as high as 10¹⁷⁻¹⁸G in their core (Bonanno et al. 2003; Rheinhardt & Geppert 2005; Burrows et al. 2007; Spruit et al. 2009; Ferrario et al. 2015; Popov 2016; see also Price & Rosswog 2006; Kawamura et al. 2016; Ciolfi et al. 2019 for simulations showing remnants of binary NSs merger harboring such intense magnetic fields). These magnetic fields substantially affect the electromagnetic phenomenology of NSs, can act as a potentially detectable source of deformation, can modify the torsional oscillations of NSs and can also alter the cooling properties of the crust. This shows that an accurate modelling of the magnetised structure of NSs is fundamental for a correct understanding of their properties.

In GR, the first magnetised model of NS dates back to Chandrasekhar & Fermi (1953). Throughout the years, many magnetised models were proposed (Ferraro 1954; Roberts 1955; Prendergast 1956; Woltjer 1960; Monaghan 1965, 1966; Roxburgh 1966; Ostriker & Hartwick 1968; Miketinac 1975), up to more recent works (Tomimura & Eriguchi 2005; Yoshida et al. 2006; Fujisawa & Eriguchi 2015). Due to the non-linearity of the general-relativistic magnetohydrodynamic (GRMHD) equations, an accurate study of the structure of NSs must be done in a numerical way, and only recently numerical results in the full GR regime have appeared. Many of these models focus on either purely toroidal (Kiuchi & Yoshida 2008; Kiuchi et al. 2009; Frieben & Rezzolla 2012) or purely poloidal (Bocquet et al. 1995; Konno 2001; Yazadjiev 2012) magnetic field configurations [see also Pili et al. (2014, 2017)]. However, such models are shown to develop an instability which causes the magnetic field to rearrange in a mixed configuration, called Twisted Torus, which is roughly axisymmetric (Prendergast 1956; Tayler 1973; Wright 1973; Braithwaite & Nordlund 2006; Braithwaite & Spruit 2006; Braithwaite 2009; Lasky et al. 2011). Twisted Torus configurations have been studied only very recently (Ciolfi & Rezzolla 2013; Pili et al. 2014; Uryū et al. 2014; Bucciantini et al. 2015; Uryū 2019), because they require to solve a large set of coupled non-linear elliptic PDEs, which can be numerically unstable.

In this paper, we present the first numerical computations of a magnetised NS in a STT of gravity in the full non-linear regime. We wish to investigate how the mutual interplay of a strong magnetic field and a scalar field modifies both the magnetic properties of NSs, with respect to GR, and their scalarisation properties with respect to the un-magnetised case. For this reason we are going to provide a characterisation as complete as possible of our equilibrium configurations, including a parametrisation of their deformation, and to carry a comparison with GR, not just in terms of global quantities but also in the specific internal distribution of density and magnetic field. The purpose is to quantify for example how much the presence of a scalar field affects the magnetic deformability of NSs, which is a key parameter to evaluate the relative importance of GW vs electromagnetic dipole emission in the early spin-down of proto-NSs (Dall’Osso et al. 2009), and to assess the validity of the millisecond-magnetar model for Long GRBs (Metzger et al. 2011). On the other hand, we also want to evaluate if the presence of a magnetic field favours or disfavours the scalarisation of NSs, and how it changes the scalarisation range, or the maximum NS mass. For this reason we limit our analysis only to the two extreme cases of purely poloidal or purely toroidal magnetic fields, neglecting rotation. In this sense our work is both an extension of the existing literature on magnetised models of NSs in GR, and of un-magnetised models in STTs.

We also take the opportunity to introduce a computational strategy, which, for the sake of simplicity, we discuss here just in the case of non-rotating NSs, but that can easily be generalised to rotating and even dynamical regimes and that allows a straightforward extension of well established algorithms for GRMHD to handle MHD in STTs. Our algorithm is an extension of the well-tested XNS solver (Pili et al. 2014, 2017) to the case of a generic STT. It is based on the eXtended Conformally Flat Condition (XCFC) for the metric (Wilson et al. 1996; Wilson & Mathews 2003; Cordero-Carrión et al. 2009; Bucciantini & Del Zanna 2011), which, even if not formally exact, has proved to be highly accurate for rotating NSs (Camelio et al. 2019). We wish to point here that the accuracy of the solution with respect to full GR depends on which parameter, that is the central rotation rate or the surface ellipticity, is held fixed in the comparison (larger deviations have been found for differentially rotating models having the same surface ellipticity Iosif & Stergioulas 2014). The XCFC system has several advantages from a numerical point of view. These, as we are going to show, are retained also in STTs, and that can easily be adapted to the more complex case of time dependent dynamical evolution.

This paper is structured as follows. In Sect. 2 we introduce MHD within STTs, both from a Lagrangian point of view, and within the 3+1 formalism. In Sect. 3 we show how the formalism developed to model magnetised NSs in GR can be extended to STTs. In Sect. 4 we present the new version of XNS for STTs. In Sect. 5 we illustrate and discuss our results for various magnetic configurations and choice of STT and, finally, we conclude in Sect. 6.

2. Scalar-tensor theories and 3+1

In the following we assume a signature { − , + , + , + } for the spacetime metric and use Greek letters μ, ν, λ,…(running from 0 to 3) for 4D spacetime tensor components, while Latin letters i, j, k,…(running from 1 to 3) are employed for 3D spatial tensor components. Moreover, we use the dimensionless units where c = G = M_⊙ = 1, and we absorb the $\sqrt{4\pi }$ factors in the definition of the electromagnetic quantities. Variables denoted with a tilde, ·̃, are calculated in the Jordan frame, while quantities denoted with a bar, ·̄, are expressed in the Einstein frame.

2.1. STT frames and ideal MHD

The most general action S_J that describes the mutual interplay of an ideal magnetised fluid at thermodynamic equilibrium with a gravitational space-time containing one scalar field φ non-minimally coupled to the metric ${\stackrel{\sim }{g}}_{\mu \nu }$, is invariant under space-time diffeomorphisms, is at most quadratic in the derivatives of the fields, and which satisfies the WEP, can be written as the sum of two terms. The first term, encoding the information about the gravitational fields, ${\stackrel{\sim }{S}}_{\mathrm{g}}[{\stackrel{\sim }{g}}_{\mu \nu },\phi ]$, according to the “Bergmann-Wagoner formulation” (Bergmann 1968; Wagoner 1970; Berti et al. 2015) is

$$\begin{array}{c}\hfill {\stackrel{\sim }{S}}_{\mathrm{g}}=\frac{1}{16\pi }{\displaystyle \int d^4}x\sqrt{-\stackrel{\sim }{g}}[\phi \stackrel{\sim }{R}-\frac{\omega (\phi )}{\phi }{\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }\phi {\stackrel{\sim }{\mathrm{\nabla }}}^{\mu }\phi -U(\phi )],\end{array}$$(1)

where $\stackrel{\sim }{g}$ is the determinant of the spacetime metric ${\stackrel{\sim }{g}}_{\mu \nu }$, ${\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }$ its associated covariant derivative, $\stackrel{\sim }{R}$ its Ricci scalar, while ω(φ) and U(φ) are, respectively, the coupling function and the potential of the scalar field φ. The second term ${\stackrel{\sim }{S}}_{\mathrm{p}}[{\stackrel{\sim }{g}}_{\mu \nu },{\stackrel{\sim }{N}}^{\mu },{\stackrel{\sim }{A}}^{\mu },\stackrel{\sim }{\epsilon },\stackrel{\sim }{s}]$ contains information on the other physical fields and it is a function of the mass current density ${\stackrel{\sim }{N}}^{\mu }=\stackrel{\sim }{\rho }{\stackrel{\sim }{u}}^{\mu }$, expressed as a function of the rest mass density $\stackrel{\sim }{\rho }$ and four-velocity ${\stackrel{\sim }{u}}^{\mu }$, the specific entropy $\stackrel{\sim }{s}$, the internal energy density $\stackrel{\sim }{\epsilon }(\stackrel{\sim }{\rho },\stackrel{\sim }{s})$, and the electromagnetic four-potential ${\stackrel{\sim }{A}}^{\mu }$. For an ideal fluid neglecting polarisation, magnetisation (Chatterjee et al. 2015; Franzon et al. 2016), dynamo or resistivity (Bucciantini & Del Zanna 2013; Del Zanna et al. 2016; Del Zanna & Bucciantini 2018; Tomei et al. 2020), it is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\stackrel{\sim }{S}}_{\mathrm{p}}& ={\displaystyle \int d^4}x\sqrt{-\stackrel{\sim }{g}}[\stackrel{\sim }{\epsilon }({\stackrel{\sim }{N}}^{\mu }{\stackrel{\sim }{N}}_{\mu },\stackrel{\sim }{s})+\zeta {\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }{\stackrel{\sim }{N}}^{\mu }+\eta {\stackrel{\sim }{N}}^{\mu }{\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }\stackrel{\sim }{s}\hfill \\ \hfill & {\stackrel{\sim }{F}}_{\mu \nu }{\stackrel{\sim }{F}}^{\mu \nu }+{\tau }_{\nu }{\stackrel{\sim }{N}}_{\mu }{\stackrel{\sim }{F}}^{\mu \nu }],\hfill \end{array}\end{array}$$(2)

where ${\stackrel{\sim }{F}}_{\mu \nu }:={\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }{\stackrel{\sim }{A}}_{\nu }-{\stackrel{\sim }{\mathrm{\nabla }}}_{\nu }{\stackrel{\sim }{A}}_{\mu }$ is the Faraday tensor, and ζ, η, τ_ν, are Lagrangian multipliers that enforce mass conservation, entropy conservation, and the ideal MHD condition ${\stackrel{\sim }{u}}_{\mu }{\stackrel{\sim }{F}}^{\mu \nu }=0$ respectively (Hawking & Ellis 1973; Brown 1993; Bekenstein & Oron 2001).

The frame where the action reads $S_{\mathrm{J}}={\stackrel{\sim }{S}}_{\mathrm{g}}+{\stackrel{\sim }{S}}_{\mathrm{p}}$ is called the “Jordan frame” (J-frame). Variation of the action with respect to the various fields (and Lagrangian multipliers) leads to the Euler-Lagrange field equations (and to the constraints). In particular, variations with respect to the four potential ${\stackrel{\sim }{A}}_{\mu }$ lead to the Maxwell equation:

$$\begin{array}{c}\hfill \delta S_{\mathrm{J}}/\delta {\stackrel{\sim }{A}}_{\mu }=0\Rightarrow {\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }{\stackrel{\sim }{F}}^{\mu \nu }=-{\stackrel{\sim }{J}}^{\nu },\end{array}$$(3)

where ${\stackrel{\sim }{J}}^{\nu }$ is the electromagnetic four-current. Variations with respect to the matter four-current ${\stackrel{\sim }{N}}^{\mu }$ lead, ultimately, to the fluid Euler equation and to the momentum-energy conservation law:

$$\begin{array}{c}\hfill \delta S_{\mathrm{J}}/\delta {\stackrel{\sim }{N}}^{\mu }=0\Rightarrow {\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }{\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }=0,\end{array}$$(4)

where the energy momentum tensor is

$$\begin{array}{c}\hfill {\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }=[\stackrel{\sim }{\rho }+\stackrel{\sim }{\epsilon }+\stackrel{\sim }{p}]{\stackrel{\sim }{u}}^{\mu }{\stackrel{\sim }{u}}^{\nu }+\stackrel{\sim }{p}{\stackrel{\sim }{g}}^{\mu \nu }+{\stackrel{\sim }{F}}_{\lambda }^{\mu }{\stackrel{\sim }{F}}^{\nu \lambda }-\frac{1}{4}{\stackrel{\sim }{F}}^{\lambda \kappa }{\stackrel{\sim }{F}}_{\lambda \kappa }{\stackrel{\sim }{g}}^{\mu \nu }\end{array}$$(5)

and $\stackrel{\sim }{p}$ is the pressure. Given that the scalar field does not enter ${\stackrel{\sim }{S}}_{\mathrm{p}}$, the equations describing the behaviour of the physical quantities are unaffected by the presence of the scalar field. Introducing the Hodge dual of the Faraday tensor ${\stackrel{\sim }{F}}^{\star \mu \nu }=\frac{1}{2}{\stackrel{\sim }{ϵ}}^{\mu \nu \lambda \kappa }{\stackrel{\sim }{F}}_{\lambda \kappa }$, where ${\stackrel{\sim }{ϵ}}_{\mu \nu \lambda \kappa }=-{(-\stackrel{\sim }{g})}^{1/2}[\mu \nu \lambda \kappa ]$ is the Levi-Civita pseudo-tensor and [μνλκ] is the alternating Levi-Civita symbol, one can write the energy momentum tensor of ideal MHD in terms of the comoving magnetic field ${\stackrel{\sim }{b}}^{\mu }={\stackrel{\sim }{u}}_{\nu }{\stackrel{\sim }{F}}^{\star \mu \nu }$ as

$$\begin{array}{c}\hfill {\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }=(\stackrel{\sim }{\rho }\stackrel{\sim }{h}+{\stackrel{\sim }{b}}^2){\stackrel{\sim }{u}}^{\mu }{\stackrel{\sim }{u}}^{\nu }-{\stackrel{\sim }{b}}^{\mu }{\stackrel{\sim }{b}}^{\nu }+(\stackrel{\sim }{p}+\frac{1}{2}{\stackrel{\sim }{b}}^2){\stackrel{\sim }{g}}^{\mu \nu },\end{array}$$(6)

where ${\stackrel{\sim }{b}}^2={\stackrel{\sim }{b}}_{\mu }{\stackrel{\sim }{b}}^{\mu }$ and $\stackrel{\sim }{h}=1+(\stackrel{\sim }{\epsilon }+\stackrel{\sim }{p})/\stackrel{\sim }{\rho }$ is the specific enthalpy. On the other hand, variations of the action with respect to the metric lead to the generalisation of Einstein’s field equations:

$$\begin{array}{c}\hfill \delta S_{\mathrm{J}}/\delta {\stackrel{\sim }{g}}_{\mu \nu }=0\Rightarrow {\stackrel{\sim }{G}}^{\mu \nu }+{\stackrel{\sim }{G}}_{\mathrm{s}}^{\mu \nu }=8\pi {\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu },\end{array}$$(7)

where ${\stackrel{\sim }{G}}^{\mu \nu }$ is the standard Einstein tensor, while ${\stackrel{\sim }{G}}_{\mathrm{s}}^{\mu \nu }$ contains the contribution from the non-minimally coupled scalar field. Tensor ${\stackrel{\sim }{G}}_{\mathrm{s}}^{\mu \nu }$ contains higher-order derivatives of the scalar field, and its associated energy density is not positively defined (Santiago & Silbergleit 2000). As a consequence, in the J-frame the generalisation of Einstein’s field equations has a different mathematical structure than in GR, implying that standard solution techniques and algorithms developed for GR cannot be naively applied. However, it is possible to show (Santiago & Silbergleit 2000) that, by performing a conformal transformation of the metric,

$$\begin{array}{c}\hfill {\overline{g}}_{\mu \nu }\phi {\stackrel{\sim }{g}}_{\mu \nu },\end{array}$$(8)

and introducing a new scalar field χ related to φ according to

$$\begin{array}{c}\hfill \frac{\mathrm{d}\chi }{\mathrm{d}ln\phi }\sqrt{\frac{\omega (\phi )+3}{4}},\end{array}$$(9)

the gravitational part of the action becomes

$$\begin{array}{c}\hfill {\overline{S}}_{\mathrm{g}}=\frac{1}{16\pi }{\displaystyle \int d^4}x\sqrt{-\overline{g}}[\overline{R}-2{\overline{\mathrm{\nabla }}}_{\mu }\chi {\overline{\mathrm{\nabla }}}^{\mu }\chi -V(\chi )],\end{array}$$(10)

where $\overline{g}$ is the determinant of the spacetime metric ${\overline{g}}_{\mu \nu }$, ${\overline{\mathrm{\nabla }}}_{\mu }$ its associated covariant derivative, $\overline{R}$ its scalar curvature, and V(χ) = U(φ)/φ² the potential of the scalar field χ. The frame where the gravitational part of the action reads as in Eq. (10) is known as the “Einstein frame” (E-frame). In the E-frame the generalisation of Einstein’s field equations reads

$$\begin{array}{c}\hfill {\overline{G}}^{\mu \nu }=8\pi ({\overline{T}}_{\mathrm{s}}^{\mu \nu }+{\overline{T}}_{\mathrm{p}}^{\mu \nu }),\end{array}$$(11)

where ${\overline{T}}_{\mathrm{p}}^{\mu \nu }={\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }/{\phi }^3$ and the contribution from the scalar field now has the form

$$\begin{array}{c}\hfill {\overline{T}}_{\mathrm{s}}^{\mu \nu }=\frac{1}{4\pi }[{\overline{\mathrm{\nabla }}}^{\mu }\chi {\overline{\mathrm{\nabla }}}^{\nu }\chi -\frac{1}{2}{\overline{g}}^{\mu \nu }{\overline{\mathrm{\nabla }}}_{\lambda }\chi {\overline{\mathrm{\nabla }}}^{\lambda }\chi ].\end{array}$$(12)

It is evident that in the E-frame the metric field equations are equivalent to those of GR, and the scalar field acts only as an extra energy-momentum source term. However, this conformal transformation affects also the physical part of the action ${\overline{S}}_{\mathrm{p}}$, redefining both the metric (and its covariant derivative) but also the physical fields (e.g. the E-frame energy density is now a function also of the scalar field χ). As a result ${\overline{\mathrm{\nabla }}}_{\mu }{\overline{T}}^{\mu \nu }\ne 0$, and ${\overline{\mathrm{\nabla }}}_{\mu }(\overline{\rho }{\overline{u}}^{\mu })\ne 0$. Interestingly Maxwell equations retain their form, as expected from their pre-metric nature (Cartan 1986; van Dantzig & Dirac 1934; Delphenich 2005). As a consequence, in the E-frame standard methods, techniques, and algorithms developed in MHD, based on the conserved nature of the various physical quantities, and the locality of the EoS, cannot be naively applied.

In the E-frame one can also derive an equation for the scalar field by varying the action with respect to χ,

$$\begin{array}{c}\hfill {\overline{\mathrm{\nabla }}}_{\mu }{\overline{\mathrm{\nabla }}}^{\mu }\chi =-4\pi {\alpha }_{\mathrm{s}}{\overline{T}}_{\mathrm{p}},\end{array}$$(13)

where ${\overline{T}}_{\mathrm{p}}:={\overline{g}}_{\mu \nu }{\overline{T}}_{\mathrm{p}}^{\mu \nu }=3\overline{p}-\overline{\epsilon }-\overline{\rho }$ and

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{s}}(\chi )-\frac{\mathrm{d}ln\phi (\chi )}{2\mathrm{d}\chi }·\end{array}$$(14)

We note that the only direct sources of a massless scalar field are those physical fields with a non-vanishing trace of the energy-momentum tensor; as such, the EM field is not a direct source of the scalar field, and for the same reason purely metric black holes in STTs are undistinguishable from those in GR (Hawking 1972; Berti et al. 2015). Analogously, in the ultra-relativistic asymptotically free regime, ε + ρ = 3p and the same considerations apply.

This suggests that a simultaneous use of the E-frame, to compute the metric and scalar field, and of the J-frame, to compute the physical field, by performing the conformal transformations between the two whenever necessary, will enable us to easily extend the standard techniques of GRMHD to the case of STTs.

2.2. 3+1 decomposition

According to the 3+1 formalism (Alcubierre 2008; Gourgoulhon 2012), any globally hyperbolic spacetime admits a foliation with a family of spacelike hypersurfaces Σ_t with normal timelike vector n^μ (which is, by definition, the velocity of the so-called “Eulerian observer”, n_μn^μ = −1). The three-metric induced on Σ_t is γ_μν = g_μν + n_μn_ν (and the induced rank-3 Levi-Civita pseudo tensor is ϵ^ijk = ϵ^ijkμn_μ). Calling x^μ = [t, xⁱ] the coordinates adapted to the foliation, the generic line element takes the form

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-{\alpha }^2\mathrm{d}{t}^2+{\gamma }_{\mathit{ij}}(\mathrm{d}{x}^i+{\beta }^i\mathrm{d}t)(\mathrm{d}{x}^j+{\beta }^j\mathrm{d}t),\end{array}$$(15)

where α is the lapse function and βⁱ is the shift vector. If βⁱ = 0, the spacetimes is said to be static. n_μ and γ_μν allow one to project any tensor according to the foliation. The relation between the E-frame Eulerian observer and J-frame one is: ${\stackrel{\sim }{n}}_{\mu }=\mathcal{A}{\overline{n}}_{\mu }$, ${\stackrel{\sim }{\gamma }}_{\mu \nu }={\mathcal{A}}^2{\overline{\gamma }}_{\mu \nu }$, where we have introduced the conformal function $\mathcal{A}=1/\sqrt{\phi (\chi )}$ coupling the two frames.

The standard 3+1 decomposition of any vector is

$$\begin{array}{c}\hfill U^{\mu }=U_{\Vert }{n}^{\mu }+U_{\perp }^{\mu },\end{array}$$(16)

where U_∥ = −n_μU^μ and $n_{\mu }{U}_{\perp }^{\mu }=0$, while any rank-2 symmetric X^μν and antisymmetric A^μν tensor can be written as

$$\begin{array}{cc}\hfill X^{\mu \nu }& =Y{n}^{\nu }{n}^{\nu }+Z^{\mu }{n}^{\nu }+Z^{\nu }{n}^{\mu }+W^{\mu \nu },\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill A^{\mu \nu }& =C^{\mu }{n}^{\nu }+C^{\nu }{n}^{\mu }+{ϵ}^{\mu \nu \lambda \kappa }{D}_{\lambda }{n}_{\kappa },\hfill \end{array}$$(18)

where n_μZ^μ = 0 = n_μW^μν and n_μC^μ = n_μD^μ = 0. Recalling that the relations between the J-frame and E-frame physical energy-momentum and Faraday tensors are ${\mathcal{A}}^6{\stackrel{\sim }{T}}^{\mu \nu }={\overline{T}}^{\mu \nu }$ and ${\mathcal{A}}^4{\stackrel{\sim }{F}}^{\mu \nu }={\overline{F}}^{\mu \nu }$, one can easily recover the following relations among the various projections:

$$\begin{array}{cc}& \stackrel{\sim }{\mathrm{\Gamma }}=-{\stackrel{\sim }{n}}_{\mu }{\stackrel{\sim }{u}}^{\mu }=-{\overline{n}}_{\mu }{\overline{u}}^{\mu }=\overline{\mathrm{\Gamma }},\hfill \end{array}$$(19)

$$\begin{array}{cc}& \mathcal{A}{\stackrel{\sim }{v}}^j={\stackrel{\sim }{\gamma }}_{\mu }^j\mathcal{A}{\stackrel{\sim }{u}}^{\mu }={\overline{\gamma }}_{\mu }^j{\overline{u}}^{\mu }={\overline{v}}^j,\hfill \end{array}$$(20)

$$\begin{array}{cc}& {\mathcal{A}}^4{\stackrel{\sim }{E}}_{\mathrm{p}}={\mathcal{A}}^{-2}{\stackrel{\sim }{n}}_{\mu }{\stackrel{\sim }{n}}_{\nu }{\mathcal{A}}^6{\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }={\overline{n}}_{\mu }{\overline{n}}_{\nu }{\overline{T}}_{\mathrm{p}}^{\mu \nu }={\overline{E}}_{\mathrm{p}},\hfill \end{array}$$(21)

$$\begin{array}{cc}& {\mathcal{A}}^5{{\stackrel{\sim }{S}}^j}_{\mathrm{p}}=-{\mathcal{A}}^{-1}{\stackrel{\sim }{n}}_{\mu }{\stackrel{\sim }{\gamma }}_{\nu }^j{\mathcal{A}}^6{\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }=-{\overline{n}}_{\mu }{\overline{\gamma }}_{\nu }^j{\overline{T}}_{\mathrm{p}}^{\mu \nu }={{\overline{S}}^j}_{\mathrm{p}},\hfill \end{array}$$(22)

$$\begin{array}{cc}& {\mathcal{A}}^6{{\stackrel{\sim }{W}}^{\mathit{ij}}}_{\mathrm{p}}={\stackrel{\sim }{\gamma }}_{\mu }^i{\stackrel{\sim }{\gamma }}_{\nu }^j{\mathcal{A}}^6{\stackrel{\sim }{T}}_{\mathrm{p}}^{\mu \nu }={\overline{\gamma }}_{\mu }^i{\overline{\gamma }}_{\nu }^j{\overline{T}}_{\mathrm{p}}^{\mu \nu }={{\overline{W}}^{\mathit{ij}}}_{\mathrm{p}},\hfill \end{array}$$(23)

$$\begin{array}{cc}& {\mathcal{A}}^3{\stackrel{\sim }{B}}^{\mu }={\mathcal{A}}^4{\stackrel{\sim }{F}}^{\star \mu \nu }{\stackrel{\sim }{n}}_{\nu }{\mathcal{A}}^{-1}={\overline{F}}^{\star \mu \nu }{\overline{n}}_{\nu }={\overline{B}}^{\mu },\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\mathcal{A}}^3{\stackrel{\sim }{E}}^{\mu }={\mathcal{A}}^4{\stackrel{\sim }{F}}^{\mu \nu }{\stackrel{\sim }{n}}_{\nu }{\mathcal{A}}^{-1}={\overline{F}}^{\mu \nu }{\overline{n}}_{\nu }={\overline{E}}^{\mu },\hfill \end{array}$$(25)

showing, for example, that the Lorentz factor Γ is the same in the two frames. The energy conservation law in J-frame, ${\stackrel{\sim }{\mathrm{\nabla }}}_{\mu }{\stackrel{\sim }{T}}^{\mu \nu }=0$, together with the mass conservation $\stackrel{\sim }{\mathrm{\nabla }}(\stackrel{\sim }{\rho }{\stackrel{\sim }{u}}^{\mu })=0$ and Maxwell equations, can be cast into a system for the evolution of the projected quantities ${\stackrel{\sim }{E}}_{\mathrm{p}},{{\stackrel{\sim }{S}}^j}_{\mathrm{p}},{\stackrel{\sim }{B}}^{\mu },{\stackrel{\sim }{E}}^{\mu }$, once an EoS and a closure for the electromagnetic currents (e.g. the Ideal MHD conditions) are provided, according for example to Del Zanna et al. (2007) and Bucciantini & Del Zanna (2011). Then, the above equations allow to rescale those quantities to the E-frame, where they are used to solve the 3+1 evolutionary equations for the metric and the scalar field. For this purpose one needs also the 3+1 projection of the latter. This is only done in the E-frame, given that it is not needed in the J-frame, according to:

$$\begin{array}{cc}& {\overline{\mathrm{\nabla }}}^{\mu }\chi =P{\overline{n}}^{\mu }+Q^{\mu },\hfill \end{array}$$(26)

$$\begin{array}{cc}& {\overline{E}}_{\mathrm{s}}={\overline{n}}_{\mu }{\overline{n}}_{\nu }{\overline{T}}_{\mathrm{s}}^{\mu \nu }=Q^2+P^2,\hfill \end{array}$$(27)

$$\begin{array}{cc}& {{\overline{S}}^j}_{\mathrm{s}}=-{\overline{n}}_{\mu }{\overline{\gamma }}_{\nu }^j{\overline{T}}_{\mathrm{s}}^{\mu \nu }=P{Q}^j,\hfill \end{array}$$(28)

$$\begin{array}{cc}& {{\overline{W}}^{\mathit{ij}}}_{\mathrm{s}}={\overline{\gamma }}_{\mu }^i{\overline{\gamma }}_{\nu }^j{\overline{T}}_{\mathrm{s}}^{\mu \nu }=Q^i{Q}^j+(Q^2+P^2){\overline{\gamma }}^{\mathit{ij}},\hfill \end{array}$$(29)

where Q^μ is purely spatial and Eq. (13) can also be cast into a set of evolutionary equations for P and Qⁱ (Salgado 2006; Salgado et al. 2008).

From now on, for the sake of clarity, and for ease of reading, we will drop the $\overline{·},\stackrel{\sim }{·}$ notation. All quantities referring either to the metric or the scalar field are assumed to be taken in the E-frame, while the MHD and fluid ones are to be considered in the J-frame. Whenever necessary, in case of possible ambiguity, the bar and tilde notation will be restored to specify the frame of reference for the given quantity.

3. Static magnetised configurations

For the problem we are interested in, we chose spherical-like coordinates x^μ = [t, r, θ, ϕ] and considered only configurations that are stationary and axisymmetric. This means that there exist two commuting Killing vectors, the timelike t^μ = (∂_t)^μ and the spacelike ϕ^μ = (∂_ϕ)^μ (Carter 1970, 2009, 2010). These two vectors span a timelike two-plane Π = Vect(t^μ, ϕ^μ). Any vector V^μ ∈ Π is said to be toroidal, and takes the form V^μ = c_tt^μ + c_ϕϕ^μ; instead, it is said to be poloidal if it lies in the spacelike two-plane orthogonal to Π. Given the generalised Einstein’s equations for the metric, Eq. (11), if both the scalar and physical energy-momentum tensors obey the relations

$$\begin{array}{c}\hfill \begin{array}{cc}& t_{\mu }{\overline{T}}^{\mu [\nu }{t}^{\kappa }{\varphi }^{\lambda ]}=0,\hfill \\ \hfill & {\varphi }_{\mu }{\overline{T}}^{\mu [\nu }{t}^{\kappa }{\varphi }^{\lambda ]}=0,\hfill \end{array}\end{array}$$(30)

where the square brackets mean anti-symmetrisation with respect to the enclosed indices, then the spacetime has the additional property of being “circular” (Kundt & Trümper 1966; Carter 1969). In this case, β^r = β^θ = 0, γ_rϕ = γ_θϕ = γ_rθ = 0 and all the remaining metric components depend solely on r and θ.

In case of circular spacetimes and spherical-like coordinates, the line element simplifies to

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-{\alpha }^2\mathrm{d}{t}^2+{\psi }^4(\mathrm{d}{r}^2+r^2\mathrm{d}{\theta }^2)+R_{\mathrm{qi}}^2{(\mathrm{d}\varphi +{\beta }^{\varphi }\mathrm{d}t)}^2,\end{array}$$(31)

where $R_{\mathrm{qi}}:=\sqrt{{\gamma }_{\varphi \varphi }}$ is the quasi-isotropic radius and ψ is the conformal factor. A metric in the form of Eq. (31) is said to be “quasi-isotropic”. Stationarity and axisymmetry are enough to ensure that ${\overline{T}}_{\mathrm{s}}^{\mu \nu }$ satisfies Eq. (30). However they are not enough to ensure the same for the physical part ${\overline{T}}_{\mathrm{p}}^{\mu \nu }$. Given that the energy-momentum tensor of the E and J-frame are related by a simple conformal transformation, and the same holds for the Killing vectors and the metric, the conditions that ensure circularity in one of them will also ensure it in the other. For an ideal plasma, having an energy-momentum tensor as in Eq. (6), on top of stationarity and axisymmetry, circularity requires the four-velocity to be toroidal, u^r = u^θ = 0, and the magnetic field b^μ to be either purely toroidal or purely poloidal (in this latter case, rotation must also be uniform). On the contrary, even if the configuration is static and axisymmetric, for a magnetic field with a mixed configuration, Eq. (30) does not hold, and in principle the metric of Eq. (31) is no longer correct. However, even in this case it has been shown in GR (Oron 2002; Shibata & Sekiguchi 2005; Dimmelmeier et al. 2006; Ott et al. 2007; Bucciantini & Del Zanna 2011; Pili et al. 2014, 2017) that Eq. (31) provides a good approximation of the correct metric, and leads to small errors in the structure of rotating stars, mostly in the outer layers close to the surface, even in the extreme cases of a rotation at the mass-shedding limit, and magnetic fields as strong as 10¹⁹G. Moreover it can be also shown that in GR the difference R_qi − ψ²rsinθ is at most of order of 10⁻³ (Pili et al. 2017). Thus, to a good level of accuracy, the metric can be further simplified to the conformally flat (CFC) approximation (Wilson & Mathews 2003; Isenberg 2008), for which

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-{\alpha }^2\mathrm{d}{t}^2+{\psi }^4[\mathrm{d}{r}^2+r^2\mathrm{d}{\theta }^2+r^2{sin}^2\theta {(\mathrm{d}\varphi +{\beta }^{\varphi }\mathrm{d}t)}^2],\end{array}$$(32)

where we have a common factor multiplying all flat-space metric terms in spherical coordinates.

From now on we shall restrict our analysis to static configurations alone, that is to the case of non-rotating stars, for which vⁱ = 0 and βⁱ = 0 (see Appendix A for a discussion on rotators). As a consequence, the ideal-MHD electric field $E_i=-{\stackrel{\sim }{ϵ}}_{\mathit{ijk}}{\mathit{v}}^j{B}^k=-{\mathcal{A}}^{-3}{\overline{ϵ}}_{\mathit{ijk}}{\mathit{v}}^j{B}^k=0$ and Sⁱ = 0. Then, it can be shown that the extrinsic curvature K_ij = 0, which means that maximal slicing, K = 0, holds (see Gourgoulhon 2012 for a discussion of the interesting properties of this kind of slicing). Under these assumptions, Einstein’s equations reduce to a system of two Poisson-like elliptic equations for ψ and α:

$$\begin{array}{cc}& \mathrm{\Delta }\psi =[-2\pi \widehat{E}]{\psi }^{-1},\hfill \end{array}$$(33)

$$\begin{array}{cc}& \mathrm{\Delta }\left(\alpha \psi \right)=\left[2\pi (\widehat{E}+2\widehat{S}){\psi }^{-2}\right]\left(\alpha \psi \right),\hfill \end{array}$$(34)

where $\mathrm{\Delta }=f^{\mathit{ij}}{\widehat{\mathrm{\nabla }}}_i{\widehat{\mathrm{\nabla }}}_j$ and ${\widehat{\mathrm{\nabla }}}_i$ are, respectively, the 3D Laplacian and nabla operator of the flat space metric f_ij. We note that the two equations are decoupled, such that Eq. (33) can be solved before Eq. (34). The source terms take the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{E}& ={\psi }^6\{{\mathcal{A}}^4[e+\frac{1}{2}{B}^2]+\frac{1}{8\pi }{Q}^2\},\hfill \\ \hfill \widehat{S}& ={\psi }^6\{{\mathcal{A}}^4[3p+\frac{1}{2}{B}^2]-\frac{1}{8\pi }{Q}^2\}·\hfill \end{array}\end{array}$$(35)

Under the same conditions, it can be shown that Eq. (13) reduces to

$$\begin{array}{c}\hfill \mathrm{\Delta }\chi =-4\pi {\psi }^4{\alpha }_{\mathrm{s}}(\chi ){\mathcal{A}}^4{T}_{\mathrm{p}}-\partial ln\left(\alpha {\psi }^2\right)\partial \chi ,\end{array}$$(36)

where ∂f∂g:​​=∂_rf∂_rg + (∂_θf∂_θg)/r² and T_p = 3p − ε − ρ is the trace of the J-frame energy momentum tensor.

We note that the Poisson-like equations Eqs. (33) and (34) for ψ and αψ have the form Δu = su^q. In GR (𝒜 = 1, Qⁱ = 0) they satisfy the criterion for local uniqueness, sq ≥ 0. In STTs (Qⁱ ≠ 0), this is no longer true; in fact, the source term in Eq. (34), s = 2πψ⁻²{𝒜⁴[ε + 6p + 3B²/2]−Q²/8π} includes an additional factor −Q²/8π such that it cannot be excluded that in particular conditions, when the scalar field is extremely strong one has s <  0. However we verified that this does not happen in any of the many configurations we computed, not even the most compact ones. Still, it remains to be verified that this holds also in the case of collapse to BH. Concerning instead Eq. (36) at first order in χ, neglecting the higher order second term on the right, it has the form Δu = sf(u). It can be shown that the condition for local uniqueness is s(df/du) ≥ 0. Now s = −4πψ⁴T_p >  0. This implies that if α_s(χ)𝒜⁴ is a decreasing function of χ, as it happens to be for STTs with spontaneous scalarisation, Eq. (36) will not satisfy local uniqueness, and multiple solutions are expected. This will be further investigated and discussed in Sect. 5.1

3.1. Purely poloidal configuration

We begin by showing how the Grad-Shafranov formalism used in GR (Del Zanna & Chiuderi 1996; Pili et al. 2017), for the case of equilibrium configurations with a purely poloidal magnetic field, can be extended to the case of STTs. The solenoidal condition of the magnetic field allows us to write it as a function of the ϕ-component of the vector potential, A_ϕ. In conformally-flat metric

$$\begin{array}{c}\hfill B^r=\frac{{\partial }_{\theta }{A}_{\varphi }}{{\mathcal{A}}^3{\psi }^6{r}^{2}sin\theta },B^{\theta }=-\frac{{\partial }_r{A}_{\varphi }}{{\mathcal{A}}^3{\psi }^6{r}^{2}sin\theta },\end{array}$$(37)

and we recall that all metric terms are in the E-frame. Function A_ϕ is also called the magnetic flux function, and its iso-surfaces A_ϕ = const, called magnetic surfaces, contain the magnetic poloidal field lines.

The Euler equation describing the static MHD equilibrium is

$$\begin{array}{c}\hfill {\partial }_{i}p+(\epsilon +p){\partial }_{i}ln(\mathcal{A}\alpha )={ϵ}_{\mathit{ijk}}{J}^i{B}^k/{\mathcal{A}}^3=L_i,\end{array}$$(38)

where Jⁱ = 𝒜²α⁻¹ϵ^ijk∂_j(𝒜αB_k) and L_i is the Lorentz force.

NSs are often assumed to be well described by a barotropic EoS, that is ε = ε(ρ) and p = p(ρ). Then, also h = h(ρ) and Eq. (38) becomes (Pili et al. 2014) the “generalised Bernoulli integral”¹

$$\begin{array}{c}\hfill ln\left(\frac{h}{h_{\mathrm{c}}}\right)+ln\left(\frac{\mathcal{A}\alpha }{{\mathcal{A}}_{\mathrm{c}}{\alpha }_{\mathrm{c}}}\right)-\mathcal{M}=0,\end{array}$$(39)

where the magnetisation function ℳ(A_ϕ) defines the Lorentz force through

$$\begin{array}{c}\hfill L_i=\rho h\frac{\mathrm{d}\mathcal{M}}{\mathrm{d}{A}_{\varphi }}{\partial }_i{A}_{\varphi },\end{array}$$(40)

and h_c, α_c, and 𝒜_c are the values of h, α and 𝒜 at the center of the star, respectively (we have assumed ℳ_c = 0). By working out the derivatives of the poloidal components of the magnetic field, one can find an equation for J^ϕ:

$$\begin{array}{c}\hfill J^{\varphi }=-\frac{1}{{\mathcal{A}}^4{\psi }^8{r}^2{sin}^2\theta }[{\mathrm{\Delta }}_{\ast }{A}_{\varphi }+\partial A_{\varphi }\partial ln\left(\alpha {\psi }^{-2}\right)],\end{array}$$(41)

where ${\mathrm{\Delta }}_{*}={\partial }_r^2+r^{-2}{\partial }_{\theta }^2-r^{-2}{(}^{\tan }{\partial }_{\theta }$. Given that, from Eq. (40), J^ϕ = ρh(dℳ/dA_ϕ), we can obtain the Grad-Shafranov equation

$$\begin{array}{c}\hfill {\check{\mathrm{\Delta }}}_3{\check{A}}_{\varphi }+\frac{\partial A_{\varphi }\partial ln\left(\alpha {\psi }^{-2}\right)}{rsin\theta }+{\mathcal{A}}^4{\psi }^{8}rsin\theta \left(\rho h\frac{\mathrm{d}\mathcal{M}}{\mathrm{d}{A}_{\varphi }}\right)=0,\end{array}$$(42)

where ${\check{A}}_{\varphi }:=A_{\varphi }/(rsin\theta )$ and ${\check{\mathrm{\Delta }}}_3{\check{A}}_{\varphi }={\mathrm{\Delta }}_{\ast }{A}_{\varphi }/(rsin\theta )$. Equation 42 allows one to find the magnetic field and current components once the metric (α and ψ) is known and the free function ℳ has been chosen. The simplest choice, found for example in Pili et al. (2014), is

$$\begin{array}{c}\hfill \mathcal{M}=k_{\mathrm{pol}}{A}_{\varphi },\end{array}$$(43)

where k_pol is the poloidal magnetisation constant. This leads to dipolar magnetic field configurations and guarantees that the currents are confined within the star.

3.2. Purely toroidal configuration

For a purely toroidal magnetic field, ℳ in Eq. (39) is no longer a function of A_ϕ and L_i = ρh∂_iℳ. Deriving the generalised Bernoulli integral and writing the Lorentz force in terms of the magnetic field components, we obtain

$$\begin{array}{c}\hfill {\partial }_{i}lnh+{\partial }_{i}ln(\mathcal{A}\alpha )+\frac{\mathcal{A}\alpha B_{\varphi }{\partial }_i\left(\mathcal{A}\alpha B_{\varphi }\right)}{\rho h{\mathcal{A}}^4{\mathcal{R}}^2}=0,\end{array}$$(44)

where ℛ² = α²ψ⁴r²sin²θ. This equation becomes integrable if we assume that the last term can be written as the gradient of a scalar function. Defining

$$\begin{array}{c}\hfill \mathcal{G}=\rho h{\mathcal{A}}^4{\mathcal{R}}^2,\end{array}$$(45)

this becomes possible if

$$\begin{array}{c}\hfill B_{\varphi }=\frac{\mathcal{I}(\mathcal{G})}{\mathcal{A}\alpha },\mathrm{and}\mathcal{M}(\mathcal{G})=-{\displaystyle \int \frac{\mathcal{I}}{\mathcal{G}}}\frac{\mathrm{d}\mathcal{I}}{\mathrm{d}\mathcal{G}}\mathrm{d}\mathcal{G}·\end{array}$$(46)

It is customary to assume a barotropic expression for ℐ (Kiuchi & Yoshida 2008; Frieben & Rezzolla 2012):

$$\begin{array}{c}\hfill \mathcal{I}=k_{\mathrm{tor}}{\mathcal{G}}^{\mathrm{m}}\mathrm{and}\mathcal{M}=-\frac{m{k}_{\mathrm{tor}}^2}{2m-1}{\mathcal{G}}^{2m-1},\end{array}$$(47)

where k_tor is the toroidal magnetisation constant and m ≥ 1 is the toroidal magnetisation index. This form of ℐ ensures that the magnetic field is confined within the star and that its configuration is symmetric with respect to the equatorial plane. The generalised Bernoulli integral then becomes

$$\begin{array}{c}\hfill ln\left(\frac{h}{h_{\mathrm{c}}}\right)+ln\left(\frac{\mathcal{A}\alpha }{{\mathcal{A}}_{\mathrm{c}}{\alpha }_{\mathrm{c}}}\right)+\frac{m{k}_{\mathrm{tor}}^2}{2m-1}{\left(\rho h{\mathcal{A}}^4{\mathcal{R}}^2\right)}^{2m-1}=0·\end{array}$$(48)

4. The XNS code

The XNS code (Bucciantini & Del Zanna 2011; Pili et al. 2014, 2015, 2017) solves the coupled equations for the metric, scalar field, and MHD structure of a NS under the assumptions of stationarity and axisymmetry, adopting conformal flatness and maximal slicing. It is based on an iterative scheme, which computes the various quantities separately. It has been applied also to the case of white dwarves (Das & Mukhopadhyay 2015) and to non-barotropic NSs (Camelio et al. 2019).

Given that the equations for the scalar quantities ψ, αψ, χ involve the Δ operator, and that the Grad-Shafranov equation can be reduced to a non-linear vector Poisson equation for ${\check{A}}_{\varphi }$, the solutions for u(r, θ) = ψ, αψ, χ are found as a sum of spherical harmonics Y_l(θ) with coefficients A_l(r) according to

$$\begin{array}{c}\hfill u(r,\theta )={\displaystyle \sum _{l=0}^{\infty }}\left[A_l(r)Y_l(\theta )\right],\end{array}$$(49)

and similarly for the vector potential,

$$\begin{array}{c}\hfill {\check{A}}_{\varphi }(r,\theta )={\displaystyle \sum _{l=0}^{\infty }}\left[C_l(r){\partial }_{\theta }{Y}_l(\theta )\right].\end{array}$$(50)

This choice leads to a series of radial, second order boundary values ODEs for the coefficients A_l(r) and C_l(r), which are solved using a tridiagonal matrix inversion. The decomposition in terms of spherical harmonics ensures the correct behaviour of the solutions on the symmetry axis, and allows us to enforce the proper boundary conditions at r = 0, where A_l(r) and C_l(r) go to zero with parity (−1)^l, and at the outer radial boundary, where we assume that A_l(r) and C_l(r) go to zero as r^{−(l + 1)}.

Given the non-linear nature of the various elliptic equations, these are solved iteratively. If the source terms do not satisfy local-uniqueness, iterative schemes might fail to converge. This issue is particularly relevant for Eq. (36) for the scalar field. As we discussed, the very nature of spontaneous scalarisation is tied to the non-uniqueness of the solutions. In the iterative scheme used to solve Eq. (36) we opted to keep fixed the trace of the energy-momentum tensor in the J-frame, and not in the E-frame. Fixing the trace in the E-frame leads to a source term of the form $-4\pi {\psi }^4{\alpha }_{\mathrm{s}}{\overline{T}}_{\mathrm{p}}$, which can be shown to violate local uniqueness for all values of χ. Fixing it in the J-frame instead leads to a source term of the form −4πψ⁴α_s𝒜⁴T_p, and it can be shown that local uniqueness is violated only in a finite range of values for χ. This ensures at least the boundedness of the solution.

In the following we briefly describe the flow structure of XNS. The code computes at the beginning the solution for a spherically symmetric non-rotating and un-magnetised NS in isotropic coordinates, at the desired central density ρ_c, solving the generalisation of the Tolman–Oppenheimer–Volkoff (TOV) equations (Tolman 1939; Oppenheimer & Volkoff 1939) to STTs (the “S-TOV” system, see Appendix B). This is achieved with a nested shooting technique requiring that in the final solution the ratio Q_r/∂_rα is constant outside the NS, and that the conformal factor ψ corresponds to the Just metric (Just 1959) in isotropic coordinates. Then, starting with an initial guess, the XNS code performs iteratively the following steps until a converged solution is found:

1. Given a distribution of the physical and scalar fields, Eqs. (33) and (34) for a new space-time metric in the E-frame are solved in sequence;

2. Using the new metric in the E-frame and the old physical fields, scalar field Eq. (36) is solved, allowing one to define a new metric in the J-frame;

3. If the magnetic field is purely toroidal, Eq. (48) is solved, and new values of the physical fields, including the magnetic field components through Eq. (46), are found in the J-frame. If the magnetic field is purely poloidal, first the equation for the vector potential Eq. (42) and then Eq. (39) are solved, determining the new physical fields in the J-frame.

4. Convergence is checked and, if not reached, the new physical metric and scalar fields are used to define a new starting model.

5. Results

In this section, we present various equilibrium configurations, analysing how the global quantities that parametrise the resulting models depend on the strength and geometry of the magnetic field. All our models, unless otherwise specified, have been computed on a 2D grid in spherical coordinates extending over the range r = [0, 100] in dimensionless units, corresponding to a range of ∼150 km, and θ = [0, π]. The grid has 400 points in the r-direction, with the first 200 points equally spaced, and covering the range r = [0, 20], and the remaining 200 points logarithmically spaced (Δr_i/Δr_i − 1 = const), and 200 equally spaced point in the angular direction. For the reference models shown in Sects. 5.2 and 5.3, the radial resolution was doubled. We have verified that at these resolutions our results have an accuracy of the order of 10⁻³, that the radius of the outer edge is far enough not to affect the solution, and the same holds for the choice of a stretched grid. In all cases the elliptic solvers use 20 spherical harmonics. We found that in order to avoid strongly oscillatory behaviours in the relaxation scheme of XNS, iterations over the various quantities Q had to be under-relaxed according to: Q_new = [Q_new + Q_old]/2.

For the ease of comparison, and in line with previous literature in GR (Bocquet et al. 1995; Kiuchi & Yoshida 2008; Frieben & Rezzolla 2012; Pili et al. 2014) we adopted a simple polytropic EoS p = K_aρ^γ_a, with an adiabatic index γ_a = 2 and a polytropic constant K_a = 110 (in dimensionless units). Concerning the magnetic field structure, for purely toroidal magnetic fields we chose a magnetic barotropic law, Eq. (47), with toroidal magnetisation index m = 1, while for purely poloidal magnetic fields we opted for the simplest choice Eq. (43) (for more complex choices see Pili et al. 2014).

The coupling function 𝒜(χ) is the only free function of a STT with zero potential. As introduced in Damour & Esposito-Farèse (1993), and used in many subsequent works (Novak 1998b; Mendes & Ortiz 2016), we adopt the choice of an exponential coupling function:

$$\begin{array}{c}\hfill \mathcal{A}\left(\chi \right)exp[{\alpha }_0\chi +\frac{{\beta }_0}{2}{\chi }^2],\end{array}$$(51)

where α₀ and β₀ are parameters whose values are constrained by observations. It can be shown (Ramazanoğlu & Pretorius 2016) that, if ${\beta }_0{\overline{T}}_{\mathrm{p}}>0$, a tachyonic instability is triggered, and modes with wavelength smaller than the NS radius grow exponentially, leading to spontaneous scalarisation, which only depends on the value of the parameter β₀. Instead, α₀ is constrained by weak-field observations (Will 2014), and has no role in this instability. It is customary (and we will follow this choice) to choose STTs with β₀ <  0, because for most EoSs of NSs in the literature ${\overline{T}}_{\mathrm{p}}<0$. We note that, in principle, spontaneous scalarisation can happen for positive values of β₀ if the EoS predicts a strongly interacting behaviour of matter in the NS core, such that ${\overline{T}}_{\mathrm{p}}>0$ (Mendes & Ortiz 2016). The most stringent constraints to this day require that, for massless scalar fields, |α₀|≲3 × 10⁻³ and β₀ ≳ −4.5. However, for a scalar field with mass, screening effects come into play and much more negative values of β₀ are in principle allowed (Yazadjiev et al. 2016). We chose α₀ = −2 × 10⁻⁴ and varied β₀ in the range [−6, −4.5]. In order to enhance and highlight the effect of spontaneous scalarisation, a particular focus will be devoted to the case β₀ = −6.

The global quantities used in the following are defined in Appendix C. It can be shown that in the E-frame the Komar and ADM masses have the same value, while in the J-frame they differ by an amount proportional to the scalar charge. For this reason, in the following, when referring generically to the mass of the NS, we always mean the Komar mass in the E-frame ($M:={\overline{M}}_{\mathrm{k}}$). On the other hand, given that the circumferential radius is a potentially measurable quantity, when referring to it we always mean its value in the J-frame. Moreover, since the metric field equations in the E-frame have the same mathematical structure as in GR, it is most natural to provide the quadrupole deformations in the E-frame, as this is where GWs should be studied.

5.1. Uniqueness of scalarised NSs

It can be shown that, given a central density ρ_c, NSs in STTs admit multiple solutions. If 𝒜 is an even function of χ then α_s is an odd-function (e.g. if α₀ = 0 in Eq. (51)) and Eq. (36) is invariant under the transformation χ → −χ (the same holds for Eqs. (33), (34) and Eqs. (39), (48)). This implies that there are three possible NS solutions: one corresponding to χ = 0, identical to GR, and two with χ ≠ 0, that only differ by the sign of χ. If α_s is an arbitrary function of χ, this symmetry breaks. If α₀ ≠ 0 in Eq. (51), then these three solutions, split into three branches: the GR solution becomes a “weakly scalarised” solution 𝒮_w, where the total scalar charge Q_s is such that α₀Q_s >  0, while the other two scalarised branches split into two “strongly scalarised” solutions: one, ${\mathcal{S}}_{\text{s}}^{+}$, with α₀Q_s >  0, the other, ${\mathcal{S}}_{\text{s}}^{-}$, with α₀Q_s <  0.

In Fig. 1, we illustrate qualitatively how these three branches behave in terms of their mass M as a function of the central density ρ_c. The range of spontaneous scalarisation, ρ_b <  ρ_c <  ρ_t, can be divided into 4 subregions depending on the relative values of the masses of the branches:

• for ρ_b <  ρ_c <  ρ₁ we have M[${\mathcal{S}}_{\text{s}}^{-}$] < M[${\mathcal{S}}_{\text{s}}^{+}$] < M[𝒮_w];

• for ρ₁ <  ρ_c <  ρ₂ we have M[${\mathcal{S}}_{\text{s}}^{+}$] < M[${\mathcal{S}}_{\text{s}}^{-}$] < M[𝒮_w];

• for ρ₂ <  ρ_c, < ρ₃ we have M[${\mathcal{S}}_{\text{s}}^{+}$] < M[𝒮_w] < M[${\mathcal{S}}_{\text{s}}^{-}$];

• for ρ₃ <  ρ_c, < ρ_t we have M[𝒮_w] < M[${\mathcal{S}}_{\text{s}}^{+}$] < M[${\mathcal{S}}_{\text{s}}^{-}$].

Fig. 1. Qualitative behaviour of multiple solutions for NSs in STTs, in terms of the relation of their mass to the central density ρc . The black, orange and red sequences represent, respectively, the weakly scalarised solutions 𝒮w and the strongly scalarised solutions ${\mathcal{S}}_{\text{s}}^{+}$ and ${\mathcal{S}}_{\text{s}}^{-}$. Green diamonds mark the position with central densities ρc = ρ1, ρ2, ρ3 where two branches have the same mass; triangles select intermediate densities (see e.g. the values in Table 1); ρb and ρt (magenta circles) represent the lower and upper limits of the central density for which spontaneous scalarisation happens.

The densities ρ_1, 2, 3 correspond to the points two branches have the same mass. Almost always, the ${\mathcal{S}}_{\text{s}}^{-}$ branch is the one where the mass shows the largest deviation from the GR (or from 𝒮_w) and is also the one with the maximum mass. In Table 1, we report the values of global quantities characterising solutions of the three branches, for few selected values of the central density, assuming α₀ = −0.05 and β₀ = −6, for spherically symmetric un-magnetised and non-rotating NSs. Such a non-physical high value of α₀ was chosen in order to enhance the differences between the ${\mathcal{S}}_{\text{s}}^{-}$ and ${\mathcal{S}}_{\text{s}}^{+}$ branches. We found that, in terms of the net scalar charge, Q_s[𝒮_w] < Q_s[${\mathcal{S}}_{\text{s}}^{+}$] < Q_s[${\mathcal{S}}_{\text{s}}^{-}$], and similarly in terms of the NS circumferential radius R_c[𝒮_w] < R_c[${\mathcal{S}}_{\text{s}}^{+}$] < R_c[${\mathcal{S}}_{\text{s}}^{-}$]. In this sense the ${\mathcal{S}}_{\text{s}}^{-}$ solution is the one with the largest deviation from GR. One can compare the three branches also in terms of their compactness 𝒞:​​=M/R_c, or in terms of their gravitational binding energy, defined as the difference between the Komar and proper masses in the E-frame, W:​​=M − M_p. We find that ${\mathcal{S}}_{\text{s}}^{-}$ is the one with the smallest compactness and highest gravitational binding energy.

If we interpret spontaneous scalarisation as an effective phase-transition (Damour & Esposito-Farèse 1996), then the difference in binding energy between the ${\mathcal{S}}_{\mathrm{s}}^{\pm }$ and 𝒮_w branches can be though of as an effective latent heat that the appearance of a scalar field releases into the system, inflating the star and reducing |W|. Within this interpretation, it is reasonable to expect that NSs undergoing spontaneous scalarisation should settle in the ${\mathcal{S}}_{\text{s}}^{-}$ branch, which is the one with the lowest |W|. Indeed we find that our code always selects the ${\mathcal{S}}_{\text{s}}^{-}$ solution [we note that for α₀ = 0, XNS always selects the GR solution, and that α₀ ≠ 0 is required to get a scalarised one; see Bucciantini et al. (2015) for a discussion of this issue with relaxation schemes for elliptic equations]. It remains to be understood, in a dynamical evolving system, which branch is selected and under what physical conditions.

In the following, we will refer to strongly scalarised solutions, in the regime where spontaneous scalarisation leads to sizeable scalar charges, simply as “scalarised”, while weakly scalarised solutions or in general solutions showing a negligible scalar charge, will be referred to as “de-scalarised” or “GR-like”.

5.2. Toroidal field models with β₀ = −6

To illustrate how a purely toroidal magnetic field affects the properties of scalarised NSs, and to allow a comparison with GR, in Fig. 2 we show the distribution of the magnetic field strength $B=\sqrt{B^{\varphi }{B}_{\varphi }}$, of the density ρ, and of the scalar field χ, for a reference model chosen in order to have the same central density, ρ_c = 8.440 × 10¹⁴g cm⁻³, and the same maximum value of the magnetic field, B_max = 6.134 × 10¹⁷ G, as in Pili et al. (2014), for α₀ = −2 × 10⁻⁴ and β₀ = −6. Comparing Fig. 2 to the GR solution (Pili et al. 2014, Fig. 1), we see that the overall distribution of the magnetic field and of the density are very similar, both in their shape and in their values: as expected for a toroidal field, the magnetic field vanishes on the symmetry axis and reaches a maximum deep inside the star, close to its center. Again, as expected, the star displays a prolate shape in density, caused by the magnetic field stress, and the outer layers are inflated to large radii by the magnetic pressure. We note that this deformation is much more pronounced in the inner parts of the star compared to its outer layers, where the density isosurfaces show only a mild deviation from a spherical shape. On the other hand, we see that the effect of the magnetic stress on the shape of the scalar field is far less evident than on the density, and the scalar field isosurfaces show the same level of prolateness throughout the star.

Fig. 2. From left to right: meridional distribution of the magnetic field strength $B=\sqrt{B^{\varphi }{B}_{\varphi }}$, of the density ρ and of the scalar field χ for a model with a toroidal magnetic field of maximum strength Bmax = 6.134 × 1017G and central density ρc = 8.440 × 1014 g cm−3. The white curve represents the surface of the star. More quantitative details on this configuration can be found in Table 2, where it is named “model T”.

In Table 2, we give the values of various global quantities characterising this model (T). Its mass M = 1.460 M_⊙ is lower than that of its GR counterpart, 1.596 M_⊙, by roughly 10%. The same holds for the baryonic mass which now is M₀ = 1.520 M_⊙, lower than in the GR case where its value is 1.680 M_⊙. With reference to the regimes shown in Fig. 1, our reference model sits between ρ_b and ρ₂, on the ${\mathcal{S}}_{\mathrm{s}}^{-}$ sequence. Interestingly, the circumferential radius R_c = 20.59 km is just 2% higher than in GR. The “radius ratio” between the surface radial coordinate at the pole, r_p, and at the equator, r_e, is r_p/r_e = 1.15, not much higher than 1, and only marginally higher than the corresponding GR value. The same holds for the quadrupole deformation e (see Appendix C for its definition). This might seem counterintuitive, because the scalar field is known to make NSs more spherical (Doneva et al. 2013), in part because the contribution of the scalar field to the quadrupole deformation has the opposite sign with respect to the matter, in part because the scalar field pressure tends to counteract matter deformations. We also provide an estimate of the quadrupolar deformation of the scalar field through the quantity e_s, that corresponds to the quadrupolar deformation of the trace of ${\overline{T}}_{\mathrm{p}}^{\mu \nu }$ (see Appendix C).

It is meaningful to compare our reference model also to an un-magnetised model in STT with the same central density, which is characterised in Table 2 as T₀. The main differences to note are the lower values of both the Komar and baryonic mass, and of the circumferential radius with respect to the magnetised case. This gives a quantitative estimate of how strong the effects of the magnetic field are and, as in GR, it shows that the magnetic field can provide extra pressure support to sustain a larger total mass. On the other hand, the compactness is higher: 𝒞 = 0.09 without a magnetic field versus 𝒞 = 0.07 in the magnetised model. This reflects in the fact the scalar charge Q_s is higher in the un-magnetised model, by about one third.

To provide a more accurate comparison of model T with the corresponding GR one, in Fig. 3 we plot for both of them the profiles of B and ρ, normalised to their maximum value. In particular, we clearly see that the STT profiles are virtually coincident with the GR ones: only the polar radius gets slightly larger. This agrees with the fact that apart from integrated quantities, that differ at most ∼10%, all other quantities characterising those models are very close, suggesting that it is not the dynamical action of the scalar field that gives rise to the differences in mass, but more likely changes in the volume element, associated to small changes in the metric. In the same figure we also compare model T to the un-magnetised model T₀, clearly showing the magnetic induced deformation on the density profile, that affects mostly the low-density outer part of the NS, nearly doubling the star’s polar radius. We also compare the profiles of χ, normalised to its maximum value χ_max. While in the central part of the star, r ≲ 7 km, the equatorial and polar profiles are respectively steeper and shallower than in the un-magnetised case, in the outer part of the star and outside it they are both shallower than in the un-magnetised case, as expected for a lower total scalar charge.

Fig. 3. Upper panel: profile of the polar (solid blue lines) and equatorial (solid orange lines) density, and of the magnetic field strength at the equator (solid green lines), normalised to their maximum values, for the equilibrium model T (with purely toroidal magnetic field) of Table 2. These are to be compared to the corresponding GR model at the same ρc and Bmax (dashed lines), and with the density of the scalarised and un-magnetised model at the same ρc, T0 (dotted purple line). Lower panel: profile of the equatorial (orange line) and polar (blue line) scalar field, normalised to the maximum value, for the equilibrium model T (solid), compared to the un-magnetised model T0 (dotted purple).

In line with Pili et al. (2014), in order to characterize the interplay of the scalar and magnetic field, in Fig. 4, for equilibrium models having all the same baryonic mass M₀ = 1.68 M_⊙, we plot the deviations Δ of ρ_c, M, R_c and e with respect to the un-magnetised case, as functions of the maximum value of the magnetic field strength inside the star B_max. The deviation of a quantity f is defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }f≔\frac{f(B_{\max },M_0)-f(0,M_0)}{f(0,M_0)},\end{array}$$(52)

except for e, in which case we just plot its value, since e(0, M₀)=0. The results are compared with the GR sequence having the same baryonic mass.

Fig. 4. Variation, with respect to the un-magnetised model, of various quantities along the equilibrium sequence with constant M0 = 1.68 M⊙ for purely toroidal magnetic field. From left to right, top to bottom: central density ρc, Komar mass Mk, circumferential radius Rc and quadrupole deformation e. The blue lines represent our STT results, to be compared to the red lines, describing the GR models in Pili et al. (2014, Fig. 2). The cyan dotted lines highlight the de-scalarised configurations; it is connected by the black dashed segments to the magenta dotted lines, which represent the same STT deviations when calculated with respect to the un-magnetised model in GR. The arrows show the direction of increasing magnetisation.

It is immediately evident that the qualitative trends are unchanged. The sequence shows that at a fixed baryonic mass there is a limit to the strength of the magnetic field that a NS can host. We find that in out STT models this value is 1.05 × 10¹⁸ G, almost twice with respect to the one of the equivalant GR sequence, 6.13 × 10¹⁷ G. As the magnetisation parameter k_m increases, so does at the beginning also B_max, until it reaches its limiting value. A further increase of k_m leads to a reduction of the magnetic field. The central density first rises with k_m, reaching a value about 10% larger at B_max ≃ 9 × 10¹⁷ G and then beginning to decrease. For weak magnetisations, we find that, for the same B_max, the deviation is about one fourth than in GR. However, once the magnetisation parameter k_m increases beyond the point where the limiting magnetic field is reached, the deviation of our STT models becomes about a factor two higher than GR. We also find that, as the magnetisation increases even farther, solutions de-scalarise (cyan dotted line), becoming equivalent to GR. When looking at ΔM or ΔR_c, one recovers similar trends, with deviations that are smaller than in GR for weak magnetic fields. Interestingly, along the scalarised part of our sequence, there seems to be a maximum value of ΔM = 0.05 at B_max = 8 × 10¹⁷ G, a behaviour not present in GR. Similarly, the quadrupolar deformation e is about one fourth than that of GR for weak magnetisations and, again, GR is recovered at high magnetisations, when the NS de-scalarises. Just focusing on the weakly magnetised part of the sequence, before the limiting magnetic field is reached, we found that the same deviations are usually achieved at twice the value of B_max with respect to GR. This indicates that NSs in STTs are far less deformable than their GR counterparts of the same baryonic mass. The origin of this behaviour is to be looked for in the effective pressure support provided by the scalar field. A purely toroidal magnetic field exerts a stress on the star that leads to a prolate matter distribution. This, as a consequence, acting as a source for the scalar field, leads to a prolate distribution of the scalar field itself. Given that the effective pressure of χ depends on its gradient, a prolate distributions leads, with respect to a spherically symmetric one, to an increased outward-pointing force along the equator and a decreased one along the polar axis (see e.g. the scalar field profiles on a prolate system shown in Fig. 3). This might seem to contradict what was found before, where we showed only marginal differences between STT and GR. But while previously the comparison was done at the same central density, here is instead done at the same baryonic mass.

In Fig. 5, we show how the magnetic energy ℋ and the scalar charge Q_s change with B_max. As the magnetisation parameter k_m rises, the magnetic energy scales with good approximation as ℋ = 1.25 × 10³⁹(B_max/10¹⁸G)²erg up to B_max ≃ 10¹⁸ G. As the magnetisation rises beyond the point where B_max = 1.03 × 10¹⁸ G, the magnetic field energy, in the scalarised part, reaches a maximum of ℋ = 1.65 × 10³⁹ erg at B_max = 9.7 × 10¹⁷ G, finally relaxing to the GR profile when the sequence de-scalarises around B_max = 0.58 × 10¹⁷ G. The scalar charge, instead, drops with increasing magnetisation, being about 10% smaller at B_max = 1.05 × 10¹⁸ G. Beyond this point, the scalar charge drops substantially until the NS completely de-scalarises.

Fig. 5. Scalar charge Qs, normalised to its value for the un-magnetised model (left panel), and magnetic field energy ℋ (right panel) as functions of Bmax along the equilibrium sequence with constant M0 = 1.68 M⊙ and purely toroidal magnetic field. The cyan dotted line highlights the de-scalarised configurations. The arrows show the direction of increasing magnetisation.

In Fig. 6, we show how the Komar mass changes with central density holding fixed the magnetic flux Φ (top panel) or the baryonic mass M₀ (middle panel). The lower bound for scalarised models, ρ_b, moves to higher densities from ρ_b = 5 × 10¹⁴ g cm⁻³ for Φ = 0 to ρ_b = 7.5 × 10¹⁴ g cm⁻³ for Φ = 2.55 × 10³⁰ G cm², while the corresponding Komar (baryonic) mass changes from 1.25 M_⊙ (1.33 M_⊙) to 1.75 M_⊙ (1.81 M_⊙). We find no evidence suggesting the existence of an upper bound to the mass of the possible de-scalarised models. Analogously, the upper bound ρ_t for scalarised models increases from ρ_t = 3.5 × 10¹⁵ g cm⁻³ for Φ = 0 to ρ_t = 4 × 10¹⁵ g cm⁻³ for Φ = 1.46 × 10³⁰ G cm², while the corresponding Komar (baryonic) mass changes from 1.60 M_⊙ (1.73 M_⊙) to 1.62 M_⊙ (1.71 M_⊙). Contrary to GR, where it is found that the maximum mass of sequences at fixed Φ increases with the magnetic flux while the central density of the related models first rises and then drops (Pili et al. 2014, Fig. 4), in our STT sequences we found that the behaviour is more complex. At densities just above ρ_b, the mass of magnetised models is found to be always larger than the un-magnetised one. However, as the density increases, the trend is reversed and we find magnetised models having a lower mass than the un-magnetised configuration at the same central density. This is reversed again once the density exceeds 2.72 × 10¹⁵ g cm⁻³ as a consequence of the shift of the position of the maximum mass. This trend is also evident by looking at configurations at fixed baryonic mass and when sequences are parametrised at fixed values of B_max or at fixed e, in Fig. 7. It is interesting to notice that close to ρ_c ≃ 2.72 × 10¹⁵ g cm⁻³ the Komar mass is independent of the magnetisation. Quantitatively, the density at which the maximum is reached always increases from ρ_c = 2.55 × 10¹⁵ g cm⁻³ for Φ = 0 to ρ_c = 2.95 × 10¹⁵ g cm⁻³ for Φ = 2.55 × 10³⁰ G cm², while the value of the maximum mass drops initially from 2.08 M_⊙ to 2.04 M_⊙ for Φ = 1.46 × 10³⁰ G cm² and then rises again to 2.08 M_⊙ for Φ = 2.55 × 10³⁰ G cm². The full characterisation of the models at maximum mass is given in Table 3.

Fig. 6. Mass-density sequences for models with purely toroidal magnetic field and β0 = −6. Upper panel: sequences computed at fixed values of the magnetic flux Φ (blue lines), compared with the un-magnetised case (red line). The dotted magenta lines represent the limit for spontaneous scalarisation. Dots mark the position of the maximum mass models UM0 (red), TM1 (light blue) and TM2 (dark blue) of Table 3. The yellow square represents model T of Fig. 2. Middle panel: sequences computed at fixed baryonic mass (green lines). Lower panel: mass difference of sequences at fixed Φ with respect to the un-magnetised one.

Fig. 7. Sequences for the models with purely toroidal magnetic field and β0 = −6. Left panel: mass-density relation computed at fixed Bmax (blue lines) compared with the un-magnetised sequence (red line). Middle panel: mass-density relation computed at fixed e (green lines) compared with the un-magnetised sequence (red line). Right panel: On top, scalar charge computed at fixed Φ (blue lines) compared with the un-magnetised sequence (red line); on bottom, trace quadrupole deformation es. In all panels, the dotted magenta lines represent the limit for spontaneous scalarisation and the yellow square represents model T of Fig. 2.

In a similar way, in Fig. 7, we have also analysed how the scalar charge Q_s changes with magnetisation. The maximum of the scalar charge goes from Q_s = 1.16 M_⊙ at Φ = 0, to Q_s = 1.14 M_⊙ when Φ = 2.55 × 10³⁰ G cm², while the density at which this maximum is reached increases from 2.09 × 10¹⁵ g cm⁻³ to 2.46 × 10¹⁵ g cm⁻³. Globally, this appears as a shift to higher density of the sequences. The maximum of the scalar charge is always reached before the maximum of the mass. Analogously to the mass, we find that close to ρ_c ≃ 2.33 × 10¹⁵ g cm⁻³ the scalar charge is independent of the magnetisation.

5.3. Poloidal field models with β₀ = −6

As it was done in the toroidal case, also for purely poloidal magnetic fields, our reference model was chosen in order to have the same central density ρ_c = 5.15 × 10¹⁴ g cm⁻³ and the same maximum value of the magnetic field B_max = 6.256 × 10¹⁷ G, as in Pili et al. (2014). Analogously to the previous toroidal case, this model sits in the part of Fig. 1 between ρ_b and ρ₂, on the sequence ${\mathcal{S}}_{\mathrm{s}}^{-}$. In Fig. 8, we show the distribution of the magnetic field strength $B=\sqrt{B^r{B}_r+B^{\theta }{B}_{\theta }}$, of the density ρ and of the scalar field χ for this model. Comparing them to the GR ones in Pili et al. (2014, Fig. 5), we see that, even for a purely poloidal magnetic field, the overall distributions of the various quantities are very similar to GR, both in their shape and in their values. As expected for a poloidal field, the magnetic field reaches a maximum at the center of the star, and vanishes in an equatorial ring located at r ≃ 12 km. The star displays an oblate shape in density, caused by the magnetic field stress, with an equatorial density profile which is almost flat close to the center. As in GR, increasing farther the magnetic field strength produces configurations where the density maximum is no longer at the center (analogously to Pili et al. 2014, Fig. 6). Again, we see that the effect of the magnetic stress on the shape of the scalar field is far less pronouced than on the density.

Fig. 8. From left to right: meridional distribution of the magnetic field strength $B=\sqrt{B^r{B}_r+B^{\theta }{B}_{\theta }}$, of the density ρ and of the scalar field χ for a model with a poloidal magnetic field of maximum strength Bmax = 6.256 × 1017 G and central density ρc = 5.15 × 1014 g cm−3. The white curve represents the surface of the star. The light white lines on the left panel represent magnetic surfaces. More quantitative details on this configuration can be found in Table 2, where it is named “model P”.

In Table 2, we give the values of various global quantities characterizing this model (P). The Komar mass M = 1.360 M_⊙, is lower than the GR mass, 1.597 M_⊙ by roughly 15%, and the same holds for the baryonic mass which is M₀ = 1.42 M_⊙, compared to the value of the GR counterpart, 1.680 M_⊙. The radius ratio r_p/r_e = 0.67 is instead marginally smaller than the GR value of 0.69. On the other hand, its circumferential radius R_c = 16.71 km is less than 1% smaller than the GR one. The quadrupole deformation e is the same as in GR. As before, it seems that the presence of a scalar field, at the same central density and for the same maximum magnetic field, does not affect the distribution of fluid quantities. Moreover, we provide an estimate of the quadrupolar deformation of the scalar field through the quantity e_s, which is comparable in strength to the quadrupole deformation e.

We can also make a comparison to the un-magnetised model with the same central density, characterised in Table 2 under the name P₀. The main differences are the values of the masses and of the circumferential radius, that are smaller for B = 0. Also the compactness is slightly lower: 𝒞 = 0.0795 without a magnetic field versus 𝒞 = 0.0814 in the magnetised model. Differently than in the toroidal case, the scalar charge Q_s is much higher in the magnetised model.

In Fig. 9, we show the profiles of the magnetic field B and densityρ, normalised to their maximum value, for the model P (solid lines) and for the corresponding GR model (dashed lines) with the same B_max and ρ_c together with the un-magnetised model P₀. We also plot the profiles of χ, normalised to its maximum value χ_max, for the models P and P₀. Again, the STT profiles are almost coincident with the GR ones: only the equatorial radius gets marginally increased. This is slightly different than the effect of the magnetic field, which changes the density profile and decreases the star’s polar radius and increases the equatorial one. The profile of the scalar field reflects the oblateness of the matter distribution, showing deviations that are somewhat smaller than the toroidal case. The same conclusions drawn in the toroidal case apply here too.

Fig. 9. Top panel: profile of the polar (solid blue lines) and equatorial (solid orange lines) density, and of the magnetic field strength (solid green lines) at the equator, normalised to their maximum values, for the equilibrium model P (with purely poloidal magnetic field) of Table 2. These are to be compared to the corresponding GR model at the same ρc and Bmax (dashed), and with the density of the scalarised and un-magnetised model at the same ρc, P0 (dotted purple line). Bottom panel: profile of the equatorial (orange line) and polar (blue line) scalar field, normalised to their maximum value, for the equilibrium model P (solid), compared to the un-magnetised model P0 (dotted purple).

In Fig. 10, we show the deviations Δ as it was done in Fig. 4. The qualitative trends are the same as in GR, and do not show the complexity of the toroidal case. In GR there was some evidence indicating that the maximum magnetic field for a NS of 1.68 M_⊙ could not exceed ≈6.2 × 10¹⁷ G. In STT we found instead that up to values or order of 1 × 10¹⁸ G there is no evidence of a saturation or limit of the maximum value of the magnetic field, which does not rule out the possibility that it might exist above 10¹⁸ G. The behaviour of all quantities appears to be monotonic in B_max: the central density decreases, while the mass, the circumferential radius and the quadrupole deformation rise. As in the toroidal case, for a given value of B_max the deviation appears to be about one fourth than in GR, while the same deviation is reached for values of B_max about twice higher than in GR. There is no evidence that the sequence would de-scalarise. As in the poloidal case, this trend can again be understood based on the effective pressure support provided by the scalar field. A purely poloidal magnetic field exerts a stress on the star that leads to an oblate matter distribution. This leads to an oblate distribution of the scalar field itself which, in turn, increases the outward-pointing force along the pole and decreases the one along the equator with respect to a spherically symmetric model. We found that, up to B_max ≈ 10¹⁸ G, the total magnetic field energy ℋ scales with a good approximation as ℋ = 0.55 × 10³⁹(B_max/10¹⁸G)² erg, and the scalar charge increases by about 2% with respect to the un-magnetised case. We also found that the magnetic dipole scales as μ = 1.5 × 10³⁵(B_max/10¹⁸G) erg G⁻¹, about 30% less than in GR. Given that the dipole moment is ultimately a measure of the net toroidal current, this can be considered a kind of global measure of a quantity integrated throughout the NS; as such, even in this case strongly affected by variations in the value of the volume element, related to the metric itself.

Fig. 10. Variation, with respect to the un-magnetised model, of various quantities along the equilibrium sequence with constant M0 = 1.68 M⊙ for purely poloidal magnetic field. From left to right, top to bottom: central density ρc, Komar mass Mk, circumferential radius Rc and quadrupole deformation e. The blue line represents our STT results, to be compared to the red line, describing the GR models of Pili et al. (2014, Fig. 7). The arrows show the direction of increasing magnetisation.

In Fig. 11, we show how the Komar mass changes with central density holding fixed the magnetic dipole moment μ or the baryonic mass M₀ (top panel). The lower bound ρ_b for scalarised models now moves to lower densities - from ρ_b = 5 × 10¹⁴ g cm⁻³ for μ = 0 to ρ_b = 4.3 × 10¹⁴ g cm⁻³ for μ = 1.57 × 10³⁵ erg G⁻¹ – while the corresponding Komar (baryonic) mass rises, going to 1.31 M_⊙ (1.38 M_⊙). Contrary to the toroidal case, we see from Fig. 12 (left panel) that, for purely poloidal magnetic fields, above a Komar mass of 1.34 M_⊙ there are no de-scalarised models. Analogously, the upper bound ρ_t for scalarised models decreases - from ρ_t = 3.5 × 10¹⁵ g cm⁻³ for μ = 0 to ρ_t = 3.42 × 10¹⁵ g cm⁻³ for μ = 0.45 × 10³⁵ erg G⁻¹ – while the Komar mass remains almost unchanged. As in GR, it is found that the maximum mass of sequences at fixed μ increases with the magnetic dipole moment, and the central density at which the maximum is reached drops. The characterisation of the models at maximum mass is given in Table 3. Similarly to GR we found that, at a given central density, the mass of equilibrium configurations is always above the un-magnetised case (in the stable part of the sequence). This same trend is also evident when sequences are parametrised at fixed values of B_max or at fixed e, in Fig. 12. Again, close to ρ_c ≃ 2.72 × 10¹⁵ g cm⁻³ the Komar mass is independent on the magnetisation. We have also analysed in Fig. 12 how the scalar charge changes with magnetisation. The maximum of the scalar charge changes from Q_s = 1.16 M_⊙ to Q_s = 1.21 M_⊙ when μ = 0.54 × 10³⁵ erg G⁻¹, while the density at which the maximum is reached drops to 1.96 × 10¹⁵ g cm⁻³. Globally, this appears as a shift to lower density of the sequences. Analogously to the mass, we find that close to ρ_c ≃ 2.33 × 10¹⁵ g cm⁻³ the scalar charge is independent on the magnetisation.

Fig. 11. Mass-density sequences for models with purely poloidal magnetic field and β0 = −6. Upper panel: sequences computed at fixed values of the magnetic dipole moment μ (blue lines) and at fixed baryonic mass (green lines), compared with the un-magnetised case (red line). The dotted magenta lines represent the limit for spontaneous scalarisation. Dots mark the position of the maximum mass models UM0 (red), PM1 (light blue) and PM2 (dark blue) of Table 3. The yellow square represents the model of Fig. 8. Lower panel: mass difference of sequences at fixed μ with respect to the un-magnetised one.

Fig. 12. Sequences for models with purely poloidal magnetic field and β0 = −6. Left panel: mass-density relation computed at fixed Bmax (blue lines) compared with the un-magnetised sequence (red line). Middle panel: mass-density relation computed at fixed e (green lines) compared with the un-magnetised sequence (red line). Right panel: on top, scalar charge-density relation computed at fixed μ (blue lines) compared with the un-magnetised sequence (red line); on bottom, trace quadrupole deformation es. In all panels, the dotted magenta lines represent the limit for spontaneous scalarisation and the yellow square represents model P of Fig. 8.

5.4. Magnetised models with β₀ = −5

In order to understand how our results depend on the specific choice of the STT parameter β₀, we have computed equilibrium configurations also for β₀ = −5 and β₀ = −4.5, closer to the limit for spontaneous scalarisation, both in the case of pure toroidal and purely poloidal magnetic fields. For β₀ = −5, the un-magnetised model with baryonic mass M_o = 1.680 M_⊙ is scalarised. It is then possible to compute deviations of various quantities with respect to their un-magnetised values, at fixed baryonic mass M_o = 1.680 M_⊙, as was done for β₀ = −6. In Fig. 13, we show how the quadrupole deformation e changes with the maximum strength of the magnetic field B_max. Again, we find that the scalarised part of the sequence shows a lower quadrupole deformation than in GR, but now this difference is not as strong as for β₀ = −6. In general e is about 2/3 of the value of the corresponding GR counterpart at the same B_max, both in the toroidal and poloidal magnetic field case. For purely toroidal magnetic fields, there is some indication that the scalarised part reaches a maximum value B_max ≃ 5.8 × 10¹⁷ G, before it de-scalarises, and then reaches a new maximum corresponding to the GR value of 6.13 × 10¹⁷ G. We can conclude that in STTs with β₀ >  −5 the upper limit to B_max is reached after the solution de-scalarises, while for β₀ <  −5 it is reached for scalarised configurations. On the other hand in models with a purely poloidal magnetic field, we observe no evidence for de-scalarisation with increasing k_pol. However, there seems to be an asymptote to a maximum value of B_max of ≃7.5 × 10¹⁷ G, slightly higher than in GR for β₀ = −5. The same conclusions can be found looking at the deviations of other variables. What we see is that changes with respect to GR depend in a strongly non-linear way on the values of β₀.

Fig. 13. Value of the quadrupole deformation e along the equilibrium sequence with constant M0 = 1.68 M⊙, as a function of Bmax, for β0 = −5 (blue lines) vs GR (red lines). The cyan dotted line highlights the un-scalarised configurations. Left panel: purely toroidal magnetic field; right panel: purely poloidal magnetic field. The arrows show the direction of increasing magnetisation.

In Fig. 14, we repeat the same analysis of Fig. 6, for purely toroidal fields. We show how the Komar mass and scalar charge change with central density holding fixed the magnetic flux Φ, and the Komar mass for fixed values of the baryonic mass M₀. The region of de-scalarisation ρ_c = [ρ_b, ρ_t] is smaller, but the behaviour of the lower and upper bounds with magnetisation is the same. The lower bound ρ_b moves to higher densities, from ρ_b = 7.07 × 10¹⁴ g cm⁻³ for Φ = 0 to ρ_c = 1.06 × 10¹⁵ g cm⁻³ for Φ = 2 × 10³⁰ G cm², and the corresponding Komar (baryonic) mass from 1.461 M_⊙ (1.57 M_⊙) to 1.75 M_⊙ (1.84 M_⊙). Again we find no evidence suggesting the existence of an upper bound to the mass of the possible de-scalarised models. Analogously, the upper bound ρ_t for scalarised models increases, from ρ_b = 2.65 × 10¹⁵ g cm⁻³ for Φ = 0 to ρ_c = 3.05 × 10¹⁵ g cm⁻³ for Φ = 2 × 10³⁰ G cm², and the corresponding Komar (baryonic) mass from 1.67 M_⊙ (1.83 M_⊙) to 1.77 M_⊙ (1.85 M_⊙). Again, we find that for toroidal magnetic fields the density at which the maximum is reached increases, and the value of the maximum mass first remains almost constant at 1.81 M_⊙, and then rises to 1.86 M_⊙ for Φ = 2 × 10³⁰G cm². In this case we also see that on sequences with Φ ≥ 1.64 × 10³⁰ G cm² the mass of equilibrium models is always larger than the relative un-magnetised counterpart at the same central density.

Fig. 14. Models with purely toroidal magnetic field and β0 = −5. Left panel: on top, sequences computed at fixed values of the magnetic flux Φ (blue lines) and at fixed baryonic mass (green lines), compared with the un-magnetised case (red line); on bottom, mass difference with respect to the un-magnetised case. Right panel: on top, scalar charge on sequences at fixed Φ; on bottom, trace quadrupole es on the same sequences. The dotted magenta lines represent the limit for spontaneous scalarisation.

For poloidal magnetic fields, we observe in Fig. 15 a more regular trend, similar to the case with β₀ = −6, where the maximum mass initially seems to remain unchanged to then rises at higher magnetisation. We find that, for poloidal fields, above a Komar mass of 1.7 M_⊙ there are no de-scalarised models.

Fig. 15. Models with purely poloidal magnetic field and β0 = −5. Left panel: on top, sequences computed at fixed values of the magnetic dipole moment μ (blue lines) and at fixed baryonic mass (green lines), compared with the un-magnetised case (red line); on bottom, mass difference with respect to the un-magnetised case. Right panel: on top, scalar charge on sequences at fixed μ; on bottom, trace quadrupole es on the same sequences. The dotted magenta lines represent the limit for spontaneous scalarisation.

It is evident that now the magnetic field plays a more dominant role that the scalar field, and the general trends of the various sequences tend to approach what was found in GR. However, in the region where the scalar charge reaches its maximum, the trends are still in line with more scalarised configurations.

5.5. Magnetised models with β₀ = −4.5

We consider here the case β₀ = −4.5, which is close to the upper limit on massless STTs set by binary pulsar constraints (Freire et al. 2012; Shao et al. 2017; Anderson et al. 2019). In Fig. 16, we show how the Komar mass changes holding fixed the magnetic flux Φ for configurations with a purely toroidal magnetic field. The scalarised range is now strongly reduced. For the un-magnetised models, ρ_b = 9.3 × 10¹⁴ g cm⁻³ and ρ_t = 2.0 × 10¹⁵ g cm⁻³, with a Komar mass that changes from 1.58 M_⊙ to 1.71 M_⊙. As the magnetic flux increases, the typical scalarised trend in the mass-density relation becomes progressively less evident: already at Φ = 0.9 × 10³⁰ G cm² the seuqence is almost indistinguishable from GR. This is made even more evident looking at the scalar charges in Fig. 16, where we observe simultaneously both a reduction of Q_s and of the scalarisation range.

Fig. 16. Left figure: models with purely toroidal magnetic field and β0 = −4.5. Upper panel: sequences computed at fixed values of the magnetic flux Φ (blue lines) and at fixed baryonic mass (green lines), compared with the un-magnetised case (red line). Bottom panel: value of the scalar charge on the same sequences at fixed Φ. Right figure: models with purely poloidal magnetic field and β0 = −4.5. Upper panel: sequences computed at fixed values of the magnetic dipole moment μ (blue lines) and at fixed baryonic mass (green lines), compared with the un-magnetised case (red line). Bottom panel: value of the scalar charge on the same sequences at fixed μ.

In case of a purely poloidal magnetic field, the trend is instead quite different, as can be seen in Fig. 16. Increasing the magnetic flux Φ, both the scalar charge and the scalarisation range increase, with ρ_b moving to lower values. The maximum mass rises, and there is no evidence for the de-scalarisation.

This difference, in part already present at lower β₀, can be understood if one recalls that spontaneous scalarisation can be seen, from a dynamical point of view, as an instability (Damour & Esposito-Farèse 1996), which can be excited only if the minimum wavelength of unstable modes (a function of β₀) is smaller or of the order of the typical highscale of the matter distribution (roughly the size of the compact star). Detailed calculations set this limit for NSs around β₀ ≈ −4.2, −4.0. It is obvious, that close to this threshold limit, any process that modifies the distribution of matter in compact stars can have deep consequances on their spontaneous scalarisability. A strong toroidal magnetic field leads to a prolate distribution of density, that on average corresponds to a reduction of the typical highscale of the matter distribution, potentially pushing the NS below the threshold for spontaneous scalarisation. On the other way a strong poloidal magnetic field leads to an oblate distribution of density, corresponding to an increase of the typical highscale of the matter distribution, potentially pushing the NS above the threshold for spontaneous scalarisation.

5.6. On the stability of magnetised equilibrium models

It is well known that NSs endowed with either a purely toroidal or a purely poloidal magnetic field are unstable against non-axisymmetric perturbations (Braithwaite & Nordlund 2006; Braithwaite & Spruit 2006; Braithwaite 2009). This is due to a magnetofluid instability that, on a typical Alfvénic timescale, leads to a reconfiguration of the magnetic field geometry toward a more tangled structure. Magnetic stability requires mixed configurations, with comparable amount of energy in the poloidal and toroidal components of the magnetic field.

With respect to axisymmetric perturbations, on the other hand, it is found that, purely poloidal magnetic fields are stable, while the stability of purely toroidal magnetic fields, against interchange modes, depends on their stratification. Toroidal configurations with m = 1 are found to be stably stratified (Schubert 1968; Fricke 1969).

Independently of their magnetofluid stability, we are going to show that in STTs, NSs with purely toroidal magnetic fields, are also gravitationally unstable against spontaneous scalarisation. The criterion for gravitational instability for non-rotating and un-magnetised NS is

$$\begin{array}{c}\hfill \frac{\partial M_0}{\partial {\rho }_{\mathrm{c}}}\le 0,\end{array}$$(53)

where the equality defines the maximum mass. We note that in GR and STTs it is the baryonic mass that formally enters the criterion, and not the Komar mass, given that the former is the dynamically conserved quantity. However in GR and STTs the Komar mass is always a monotonically increasing function of the baryonic mass and one can safely use it to evaluate stability. This criterion can be generalised to magnetic configurations. Recalling that the flux-freezing condition of ideal MHD, ensures that the magnetic flux Φ is conserved in axisymmetry, one has that NSs with a purely toroidal magnetic field are unstable when

$$\begin{array}{c}\hfill \frac{\partial M_0}{\partial {\rho }_{\mathrm{c}}}{|}_{\mathrm{\Phi }}\le 0.\end{array}$$(54)

In Fig. 17, we plot how the baryonic mass of various equilibrium configurations change with density at fixed values of the magnetic flux Φ. It is immediately evident that each sequence shows four parts:

• a gravitationally stable de-scalarised GR part;

• a gravitationally unstable scalarised part;

• a gravitationally stable scalarised part (up to the density of the model of maximum mass for the entire sequence);

• a gravitationally unstable scalarised part (beyond the density of the model of maximum mass for the entire sequence).

This is in sharp contrast to GR, where only two parts are found (stable and unstable), separated by the model with maximum mass. In principle now we can have two maxima for the mass of NSs with purely toroidal magnetic fields: one corresponding to the de-scalarised part and one to the scalarised one. In the mass-density diagram there is a region where models are gravitationally unstable. Moreover, for any given value of Φ, there is a range of masses where both de-scalarised and scalarised solutions are possible. On the other hand, there is a lower limit to the values of the magnetic flux that can support de-scalarised configurations of a given baryonic mass. Lowering the magnetic flux beyond this limit could lead to a gravitational instability where the star jumps from the de-scalarised branch to the scalarised one. This is a gravitational instability, unrelated to rearrangements of the magnetic field geometry, that will take place on a typical scalarisation timescale, of the order of the light crossing time of the NS. For example, with reference to Fig. 17, a de-scalarised configuration with M₀ = 1.68 M_⊙ can only exist for Φ >  2.06 × 10³⁰ G cm² and ρ_c <  7.12 × 10¹⁴ g cm⁻³; below this limiting value of the magnetic flux, the NS will jump at the same baryonic mass but with a central density ρ_c >  1.32 × 10¹⁵ g cm⁻³, and a scalar charge Q_s = 0.8 M_⊙. Interestingly, these two limiting configurations have not just the same baryonic mass, and magnetic flux, but also the same Komar mass M_k = 1.62 M_⊙. We have repeated this analysis also for higher values of β₀ and found that this effect already disappears at β₀ = −5. However, for β₀ = −4.5 we found that two configurations, one scalarised and the other un-scalarised, with the same baryonic mass still exist, but in this case they have the same central density.

Fig. 17. Sequences at fixed magnetic flux Φ, computed in the case β0 = −6. The red curve is the un-magnetised solution. From bottom to top: other curves are computed at Φ = [2.55, 2.06, 1.46, 0.91]×1030 G cm2. The various parts are: gravitationally stable de-scalarised branch (solid magenta); gravitationally unstable scalarised branch (black dashed); gravitationally stable scalarised branch (solid blue). The yellow region corresponds to gravitationally unstable models. The black dot represents the de-scalarised configuration with M0 = 1.68 M⊙ and Φ = 2.06 × 1030 G cm2, while the arrow points to the blue dot where the configuration is expected to jump because of magnetically-induced spontaneous scalarisation.

Independently of the specific choice of magnetic field distribution, that in our case is dictated by the request of an integrable form for the generalised Bernoulli equation, our results have shown that a strong toroidal magnetic field can support de-scalarised configurations, and that, in principle, if such magnetic field drops below a limiting value (for example because of non ideal processes or magnetic instabilities) such configuration can undergo a rapid “magnetically-induced spontaneous scalarisation”. In the case of purely poloidal configurations, the quantity that is dynamically conserved for axisymmetric perturbations is the net flux of the toroidal current J^ϕ. This can be equivalently parametrised by the magnetic dipole moment. If we repeat the same analysis done in the toroidal case, considering sequences at fixed magnetic dipole moment, we see no evidence for the presence of an unstable part.

6. Conclusions

On the one hand, a proper understanding of the role of magnetic fields is fundamental in the physics and phenomenology of NSs. Magnetic fields affect virtually all of their observational properties and can modify their structure to the point of affecting also their gravitational behaviour. On the other hand, several existing issues in our understanding of the physical Universe have led many theorists to postulate extensions of GR, some of which make interesting predictions on the structure of NSs. In particular, a class of theories known as scalar-tensor theories allow a phenomenon called spontaneous scalarisation, that in principle can lead to sizeable deviations in the structure of NSs from GR. Here, for the first time, we modelled and investigated the properties of magnetised NSs in STTs subject to spontaneous scalarisation, in the full non-linear regime, assuming either purely toroidal or purely poloidal magnetic fields. This is an extension and improvement of our previous work on magnetised NSs in GR.

We have shown how to develop a strategy, within the framework of the 3+1 formalism, to extend standard techniques developed for GRMHD to the case of STTs, by making simultaneous use of the Einstein-frame (where the metric equation have the same mathematical structure as in GR, and the same numerical schemes can be applied) and the Jordan-frame (where the magnetofluid equations retain their conservative, quasi-hyperbolic form, and thus are amenable to be treated with standard finite volume or finite difference conservative schemes for fluid dynamics). In particular, for simplicity and ease of discussion, we have focused on the case of static configurations, illustrating how the equations that describe the density and magnetic field distribution change in the presence of a scalar field, and how the effects of a scalar field can be fully encapsulated in the conformal scaling factor 𝒜.

Our formalism is based on the so called eXtended Conformally Flat Approximation (XCFC), which has proved to be very accurate in GR – even for strongly deformed NSs – and in the fully dynamical regime, also for systems undergoing collapse to black hole, as long as one is not interested in the GW emission. The XCFC approach has several advantages in GR: the source terms of the metric equations are the same conserved variables evolved by the conservative algorithm for the fluid dynamics; the equations are decoupled and can be solved sequentially; local uniqueness is satisfied. We have verified that in STTs, the XCFC approach retains these properties. Even if, in principle, because of the sign of the scalar field term in the equation for αψ, local uniqueness could be violated, we have checked that, practically, this is never the case, even for the most scalarised of our configurations.

We have shown that spontaneous scalarisation leads to multiple solutions for NSs, either weakly or strongly scalarised, and we have shown and characterised how the symmetry of the strongly scalarised branch is broken if one chooses a value α₀ ≠ 0. In particular, we verified that our numerical algorithm always selects the ${\mathcal{S}}_{\mathrm{s}}^{-}$ branch. We also showed that the ${\mathcal{S}}_{\mathrm{s}}^{-}$ solutions are not always the ones with the largest deviation from GR in the mass-density diagram, but are always the ones with the largest scalar charge and smallest compactness.

In this paper, we carry out a detailed study of the properties of magnetised NSs in STT with spontaneous scalarisation, trying to characterise them as completely as possible, not just in term of their masses or radii, but also considering how the interplay of the magnetic and scalar fields affect their internal structure and deformation. We also tried to characterise the deformation of the scalar field, and introduced the parameter e_s related to the emission of quadrupolar scalar waves.

In general, we find that the action of different configurations of the magnetic field on the overall structure of a NS leads to qualitatively similar results: a toroidal magnetic field produced prolate configurations, while a poloidal field leads to oblate one. However, significative changes are found when we proceed to a quantitative comparison.

When comparing STT to GR models, computed at the same central density ρ_c and maximum value of the magnetic field B_max, we found that the distribution of density and magnetic field vary less than few percent. This suggests that GR models can be used as good proxy for the internal structure of magnetised NSs in STT. On the other hand, when for the same models we compare global integrated quantities like the mass, or the quadrupole deformation, we found deviations from GR up to 10–20%. This difference can be easily understood recalling that while the distributions of density and magnetic field depend on the ratio α/α_c (i.e. on relative changes of the metric terms), the value of integrated quantities depends on the conformal factor ψ⁶ through the volume element (i.e. on the absolute values of the metric terms). On top of this, the quadrupole deformation e, used to estimate the possible emission of GWs from deformed system, is properly computed in the E-frame, where the metric equations have the same mathematical structure of GR.

We have also investigated sequences at fixed baryonic mass, which is the conserved quantity from a dynamical and evolutionary perspective, and compared typical trends with those of GR for the same baryonic mass. We found that, in general, the presence of a scalar field reduces the deformability of NSs and tends to reduce the typical deviations from the spherically symmetric un-magnetised configuration. This also implies that with respect to GR, NSs at the same baryonic mass can host stronger magnetic fields. For configurations with purely toroidal magnetic fields we also showed that as the magnetisation rises the models de-scalarise. This effect was evaluated for various values of β₀ showing that there is a strong dependency.

We have then shown, using various parametrisations, how the mass-density relation changes with the magnetisation of the system, revealing both how this affects the region of spontaneous scalarisation and the location of the configuration with maximum mass, together with its value. In particular, we have shown that while for toroidal magnetic fields there is a de-scalarised region, for purely poloidal magnetic fields there is a limiting mass above which only scalarised solutions are possible. We have also shown that contrary to GR, where the maximum mass is always an increasing function of magnetisation, in STTs, for purely toroidal magnetic fields, the maximum mass decreases with increasing magnetisation for systems with B_max lower than a threshold magnetic field, and then rises. We verified that the quadrupolar term arising from magnetic deformations in the source of the scalar field equation is of the same order of the one in Einstein’s equations, suggesting comparable levels of gravitational losses in tensor and scalar waves.

In general, we found that for weakly magnetised models the presence of a scalar field dominates the properties of NSs, and its effect is to counter-balance the magnetic stresses, either by reducing the deformation, or leading to saturation of the values of the maximum mass. We verified, by changing the value of β₀, that when scalarisation effects become smaller the typical trends of GR tend to be recovered, with the significative difference that while for purely toroidal fields a rise in magnetisation leads to de-scalarisation, for purely poloidal magnetic fields, on the contrary, it increases the total scalar charge. Depending on its geometry, the magnetic field can either favour or suppress spontaneous scalarisation when β₀ is close to the threshold limit on the range of this effect.

Finally, we have also shown that the mutual interplay of a scalar and toroidal magnetic field, in the presence of strong scalarisation effects, leads to unstable configurations and potentially to events of spontaneous scalarisation due to the loss of magnetic support - a ‘magnetically-induced spontaneous scalarisation’.

This paper is mostly devoted to a global study of the properties of magnetised NSs in STT, with a particular focus on the comparison with their respective GR counterparts. For this reason, we adopted a simple polytropic EoS and considered only the two extreme cases of purely toroidal and purely poloidal magnetic fields, focusing the discussion on the case β₀ = −6 to enhance and highlight the main differences. We plan to investigate in more detail, in a future work, how the deformability of NSs in STT depends on the choice of β₀, and on the EoS (Pili et al. 2016), and how it scales with the mass, radius, and compactness of NSs to see if it is possible to derive scaling laws that can parametrise the magnetic deformability, in a similar way to what has been previously done in GR (Pili et al. 2017).

We conclude by recalling that STTs are just a subset of a more extended class of alternative theories of gravity, TeVeS (Bekenstein 2004), which predict also the possible existence of non-minimally coupled vector fields. As STTs, even theories with vector fields can present phenomena of “spontaneous vectorisation” (Hellings & Nordtvedt 1973; Heisenberg 2014; Kase et al. 2018, 2020). Interestingly the mathematics behind spontaneous vectorisation is not dissimilar to the one used to model non-linear current terms in magnetised NSs (Pili et al. 2014), and spontaneous magnetic-vectorisation has already been treated and discussed within the framework of the standard techniques that we have illustrated here (Bucciantini et al. 2015). This shows that the algorithms and approaches we have introduced, even if developed in the context of the specific case of magnetic fields, have a far larger applicability to vector fields in general.

Acknowledgments

The authors acknowledge financial support from the “Accordo Attuativo ASI-INAF n. 2017-14-H.0 Progetto: on the escape of cosmic rays and their impact on the background plasma” and from the INFN Teongrav collaboration. We also thanks the referee for pointing us a few mistakes in the original text.

References