Issue 
A&A
Volume 640, August 2020



Article Number  A44  
Number of page(s)  24  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202037918  
Published online  10 August 2020 
Axisymmetric equilibrium models for magnetised neutron stars in scalartensor theories
^{1}
Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, 50019 Sesto F. no (Firenze), Italy
^{2}
INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
email: soldateschi@arcetri.astro.it
^{3}
INFN – Sezione di Firenze, Via G. Sansone 1, 50019 Sesto F. no (Firenze), Italy
Received:
10
March
2020
Accepted:
25
May
2020
Among the possible extensions of general relativity that have been put forward to address some longstanding issues in our understanding of the Universe, scalartensor theories have received a lot of attention for their simplicity. Interestingly, some of these predict a potentially observable nonlinear phenomenon, known as spontaneous scalarisation, in the presence of highly compact matter distributions, as in the case of neutron stars. Neutron stars are ideal laboratories for investigating the properties of matter under extreme conditions and, in particular, they are known to harbour the strongest magnetic fields in the Universe. Here, for the first time, we present a detailed study of magnetised neutron stars in scalartensor theories. First, we showed that the formalism developed for the study of magnetised neutron stars in general relativity, based on the “extended conformally flat condition”, can easily be extended in the presence of a nonminimally coupled scalar field, retaining many of its numerical advantages. We then carried out a study of the parameter space considering the two extreme geometries of purely toroidal and purely poloidal magnetic fields, varying both the strength of the magnetic field and the intensity of scalarisation. We compared our results with magnetised generalrelativistic solutions and unmagnetised scalarised solutions, showing how the mutual interplay between magnetic and scalar fields affect the magnetic and the scalarisation properties of neutron stars. In particular, we focus our discussion on magnetic deformability, maximum mass, and range of scalarisation.
Key words: gravitation / stars: magnetic field / stars: neutron / magnetohydrodynamics (MHD) / methods: numerical / relativistic processes
© 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 strongfield 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; GaussBonnet gravity (Lovelock 1971), which increases the dimensionality of the spacetime; scalartensor 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 & EspositoFarè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 lowenergy 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 selfgravitating 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 EinsteinHilbert action of GR attempting to bring it in conformity with Mach’s principle by replacing the gravitational constant G by a scalar field nonminimally coupled to the spacetime metric, giving birth to the JordanFierzBransDicke 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 monoscalar STTs (Damour & EspositoFarèse 1993), a nonperturbative 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 nonscalarised 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 nonscalarised companion (Barausse et al. 2013); or even in an isolated NS system, where “spontaneous scalarisation” can develop (this was the first discovered nonperturbative strong field effect in STTs, Damour & EspositoFarè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 weakgravity 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 Xray afterglow of ShortGRBs (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 premerger inspiral, and in the spindown history of millisecond protomagnetar 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 nonmagnetised models of NSs have been studied in STTs in the full nonlinear regime (see e.g. Suvorov 2018 for a perturbative approach to the magnetised scenario). Most of them focus on static (Damour & EspositoFarèse 1993; Harada 1998; Novak 1998a; Taniguchi et al. 2015; Anderson & Yunes 2019; Doneva & Yazadjiev 2020) or slowly rotating (Damour & EspositoFarè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^{8−12}G for normal pulsars and up to 10^{16}G at the surface of magnetars, while newly formed protoNSs are hypothesised to store magnetic fields as high as 10^{17−18}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 nonlinearity of the generalrelativistic 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 nonlinear 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 nonlinear 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 unmagnetised 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 spindown of protoNSs (Dall’Osso et al. 2009), and to assess the validity of the millisecondmagnetar 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 unmagnetised 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 nonrotating 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 welltested 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; CorderoCarrió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. Scalartensor 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 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 spacetime containing one scalar field φ nonminimally coupled to the metric , is invariant under spacetime 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, , according to the “BergmannWagoner formulation” (Bergmann 1968; Wagoner 1970; Berti et al. 2015) is
where is the determinant of the spacetime metric , its associated covariant derivative, its Ricci scalar, while ω(φ) and U(φ) are, respectively, the coupling function and the potential of the scalar field φ. The second term contains information on the other physical fields and it is a function of the mass current density , expressed as a function of the rest mass density and fourvelocity , the specific entropy , the internal energy density , and the electromagnetic fourpotential . 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
where is the Faraday tensor, and ζ, η, τ_{ν}, are Lagrangian multipliers that enforce mass conservation, entropy conservation, and the ideal MHD condition respectively (Hawking & Ellis 1973; Brown 1993; Bekenstein & Oron 2001).
The frame where the action reads is called the “Jordan frame” (Jframe). Variation of the action with respect to the various fields (and Lagrangian multipliers) leads to the EulerLagrange field equations (and to the constraints). In particular, variations with respect to the four potential lead to the Maxwell equation:
where is the electromagnetic fourcurrent. Variations with respect to the matter fourcurrent lead, ultimately, to the fluid Euler equation and to the momentumenergy conservation law:
where the energy momentum tensor is
and is the pressure. Given that the scalar field does not enter , 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 , where is the LeviCivita pseudotensor and [μνλκ] is the alternating LeviCivita symbol, one can write the energy momentum tensor of ideal MHD in terms of the comoving magnetic field as
where and 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:
where is the standard Einstein tensor, while contains the contribution from the nonminimally coupled scalar field. Tensor contains higherorder derivatives of the scalar field, and its associated energy density is not positively defined (Santiago & Silbergleit 2000). As a consequence, in the Jframe 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,
and introducing a new scalar field χ related to φ according to
the gravitational part of the action becomes
where is the determinant of the spacetime metric , its associated covariant derivative, its scalar curvature, and V(χ) = U(φ)/φ^{2} 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” (Eframe). In the Eframe the generalisation of Einstein’s field equations reads
where and the contribution from the scalar field now has the form
It is evident that in the Eframe the metric field equations are equivalent to those of GR, and the scalar field acts only as an extra energymomentum source term. However, this conformal transformation affects also the physical part of the action , redefining both the metric (and its covariant derivative) but also the physical fields (e.g. the Eframe energy density is now a function also of the scalar field χ). As a result , and . Interestingly Maxwell equations retain their form, as expected from their premetric nature (Cartan 1986; van Dantzig & Dirac 1934; Delphenich 2005). As a consequence, in the Eframe 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 Eframe one can also derive an equation for the scalar field by varying the action with respect to χ,
where and
We note that the only direct sources of a massless scalar field are those physical fields with a nonvanishing trace of the energymomentum 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 ultrarelativistic asymptotically free regime, ε + ρ = 3p and the same considerations apply.
This suggests that a simultaneous use of the Eframe, to compute the metric and scalar field, and of the Jframe, 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 socalled “Eulerian observer”, n_{μ}n^{μ} = −1). The threemetric induced on Σ_{t} is γ_{μν} = g_{μν} + n_{μ}n_{ν} (and the induced rank3 LeviCivita pseudo tensor is ϵ^{ijk} = ϵ^{ijkμ}n_{μ}). Calling x^{μ} = [t, x^{i}] the coordinates adapted to the foliation, the generic line element takes the form
where α is the lapse function and β^{i} is the shift vector. If β^{i} = 0, the spacetimes is said to be static. n_{μ} and γ_{μν} allow one to project any tensor according to the foliation. The relation between the Eframe Eulerian observer and Jframe one is: , , where we have introduced the conformal function coupling the two frames.
The standard 3+1 decomposition of any vector is
where U_{∥} = −n_{μ}U^{μ} and , while any rank2 symmetric X^{μν} and antisymmetric A^{μν} tensor can be written as
where n_{μ}Z^{μ} = 0 = n_{μ}W^{μν} and n_{μ}C^{μ} = n_{μ}D^{μ} = 0. Recalling that the relations between the Jframe and Eframe physical energymomentum and Faraday tensors are and , one can easily recover the following relations among the various projections:
showing, for example, that the Lorentz factor Γ is the same in the two frames. The energy conservation law in Jframe, , together with the mass conservation and Maxwell equations, can be cast into a system for the evolution of the projected quantities , 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 Eframe, 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 Eframe, given that it is not needed in the Jframe, according to:
where Q^{μ} is purely spatial and Eq. (13) can also be cast into a set of evolutionary equations for P and Q^{i} (Salgado 2006; Salgado et al. 2008).
From now on, for the sake of clarity, and for ease of reading, we will drop the notation. All quantities referring either to the metric or the scalar field are assumed to be taken in the Eframe, while the MHD and fluid ones are to be considered in the Jframe. 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 sphericallike 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 twoplane Π = Vect(t^{μ}, ϕ^{μ}). Any vector V^{μ} ∈ Π is said to be toroidal, and takes the form V^{μ} = c_{t}t^{μ} + c_{ϕ}ϕ^{μ}; instead, it is said to be poloidal if it lies in the spacelike twoplane orthogonal to Π. Given the generalised Einstein’s equations for the metric, Eq. (11), if both the scalar and physical energymomentum tensors obey the relations
where the square brackets mean antisymmetrisation 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 sphericallike coordinates, the line element simplifies to
where is the quasiisotropic radius and ψ is the conformal factor. A metric in the form of Eq. (31) is said to be “quasiisotropic”. Stationarity and axisymmetry are enough to ensure that satisfies Eq. (30). However they are not enough to ensure the same for the physical part . Given that the energymomentum tensor of the E and Jframe 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 energymomentum tensor as in Eq. (6), on top of stationarity and axisymmetry, circularity requires the fourvelocity 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 massshedding limit, and magnetic fields as strong as 10^{19}G. Moreover it can be also shown that in GR the difference R_{qi} − ψ^{2}rsinθ is at most of order of 10^{−3} (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
where we have a common factor multiplying all flatspace metric terms in spherical coordinates.
From now on we shall restrict our analysis to static configurations alone, that is to the case of nonrotating stars, for which v^{i} = 0 and β^{i} = 0 (see Appendix A for a discussion on rotators). As a consequence, the idealMHD electric field and S^{i} = 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 Poissonlike elliptic equations for ψ and α:
where and 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
Under the same conditions, it can be shown that Eq. (13) reduces to
where ∂f∂g:=∂_{r}f∂_{r}g + (∂_{θ}f∂_{θ}g)/r^{2} and T_{p} = 3p − ε − ρ is the trace of the Jframe energy momentum tensor.
We note that the Poissonlike equations Eqs. (33) and (34) for ψ and αψ have the form Δu = su^{q}. In GR (𝒜 = 1, Q^{i} = 0) they satisfy the criterion for local uniqueness, sq ≥ 0. In STTs (Q^{i} ≠ 0), this is no longer true; in fact, the source term in Eq. (34), s = 2πψ^{−2}{𝒜^{4}[ε + 6p + 3B^{2}/2]−Q^{2}/8π} includes an additional factor −Q^{2}/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πψ^{4}T_{p} > 0. This implies that if α_{s}(χ)𝒜^{4} 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 GradShafranov 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 conformallyflat metric
and we recall that all metric terms are in the Eframe. Function A_{ϕ} is also called the magnetic flux function, and its isosurfaces A_{ϕ} = const, called magnetic surfaces, contain the magnetic poloidal field lines.
The Euler equation describing the static MHD equilibrium is
where J^{i} = 𝒜^{2}α^{−1}ϵ^{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”^{1}
where the magnetisation function ℳ(A_{ϕ}) defines the Lorentz force through
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^{ϕ}:
where . Given that, from Eq. (40), J^{ϕ} = ρh(dℳ/dA_{ϕ}), we can obtain the GradShafranov equation
where and . 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
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
where ℛ^{2} = α^{2}ψ^{4}r^{2}sin^{2}θ. This equation becomes integrable if we assume that the last term can be written as the gradient of a scalar function. Defining
this becomes possible if
It is customary to assume a barotropic expression for ℐ (Kiuchi & Yoshida 2008; Frieben & Rezzolla 2012):
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
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 nonbarotropic NSs (Camelio et al. 2019).
Given that the equations for the scalar quantities ψ, αψ, χ involve the Δ operator, and that the GradShafranov equation can be reduced to a nonlinear vector Poisson equation for , the solutions for u(r, θ) = ψ, αψ, χ are found as a sum of spherical harmonics Y_{l}(θ) with coefficients A_{l}(r) according to
and similarly for the vector potential,
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 nonlinear nature of the various elliptic equations, these are solved iteratively. If the source terms do not satisfy localuniqueness, 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 nonuniqueness of the solutions. In the iterative scheme used to solve Eq. (36) we opted to keep fixed the trace of the energymomentum tensor in the Jframe, and not in the Eframe. Fixing the trace in the Eframe leads to a source term of the form , which can be shown to violate local uniqueness for all values of χ. Fixing it in the Jframe instead leads to a source term of the form −4πψ^{4}α_{s}𝒜^{4}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 nonrotating and unmagnetised 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 “STOV” 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:

Given a distribution of the physical and scalar fields, Eqs. (33) and (34) for a new spacetime metric in the Eframe are solved in sequence;

Using the new metric in the Eframe and the old physical fields, scalar field Eq. (36) is solved, allowing one to define a new metric in the Jframe;

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 Jframe. 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 Jframe.

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 rdirection, 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^{−3}, 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 underrelaxed 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 & EspositoFarèse (1993), and used in many subsequent works (Novak 1998b; Mendes & Ortiz 2016), we adopt the choice of an exponential coupling function:
where α_{0} and β_{0} are parameters whose values are constrained by observations. It can be shown (Ramazanoğlu & Pretorius 2016) that, if , 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 β_{0}. Instead, α_{0} is constrained by weakfield observations (Will 2014), and has no role in this instability. It is customary (and we will follow this choice) to choose STTs with β_{0} < 0, because for most EoSs of NSs in the literature . We note that, in principle, spontaneous scalarisation can happen for positive values of β_{0} if the EoS predicts a strongly interacting behaviour of matter in the NS core, such that (Mendes & Ortiz 2016). The most stringent constraints to this day require that, for massless scalar fields, α_{0}≲3 × 10^{−3} and β_{0} ≳ −4.5. However, for a scalar field with mass, screening effects come into play and much more negative values of β_{0} are in principle allowed (Yazadjiev et al. 2016). We chose α_{0} = −2 × 10^{−4} and varied β_{0} 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 β_{0} = −6.
The global quantities used in the following are defined in Appendix C. It can be shown that in the Eframe the Komar and ADM masses have the same value, while in the Jframe 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 Eframe (). 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 Jframe. Moreover, since the metric field equations in the Eframe have the same mathematical structure as in GR, it is most natural to provide the quadrupole deformations in the Eframe, 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 oddfunction (e.g. if α_{0} = 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} ≠ 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 α_{0}Q_{s} > 0, while the other two scalarised branches split into two “strongly scalarised” solutions: one, , with α_{0}Q_{s} > 0, the other, , with α_{0}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} < ρ_{1} we have M[] < M[] < M[𝒮_{w}];

for ρ_{1} < ρ_{c} < ρ_{2} we have M[] < M[] < M[𝒮_{w}];

for ρ_{2} < ρ_{c}, < ρ_{3} we have M[] < M[𝒮_{w}] < M[];

for ρ_{3} < ρ_{c}, < ρ_{t} we have M[𝒮_{w}] < M[] < M[].
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 and . 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 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} = −0.05 and β_{0} = −6, for spherically symmetric unmagnetised and nonrotating NSs. Such a nonphysical high value of α_{0} was chosen in order to enhance the differences between the and branches. We found that, in terms of the net scalar charge, Q_{s}[𝒮_{w}] < Q_{s}[] < Q_{s}[], and similarly in terms of the NS circumferential radius R_{c}[𝒮_{w}] < R_{c}[] < R_{c}[]. In this sense the 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 Eframe, W:=M − M_{p}. We find that is the one with the smallest compactness and highest gravitational binding energy.
Values of various physical quantities describing the solutions 𝒮_{w}, and , for α_{0} = −0.05 and β_{0} = −6, and for selected values of the central density ρ_{c} (in the Jframe), corresponding from top to bottom to: ρ_{b} < ρ_{c} < ρ_{1}, ρ_{3} < ρ_{c} < ρ_{t}, ρ_{c} = ρ_{1}, ρ_{2}, ρ_{3}.
If we interpret spontaneous scalarisation as an effective phasetransition (Damour & EspositoFarèse 1996), then the difference in binding energy between the 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 branch, which is the one with the lowest W. Indeed we find that our code always selects the solution [we note that for α_{0} = 0, XNS always selects the GR solution, and that α_{0} ≠ 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 “descalarised” or “GRlike”.
5.2. Toroidal field models with β_{0} = −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 , 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^{14}g cm^{−3}, and the same maximum value of the magnetic field, B_{max} = 6.134 × 10^{17} G, as in Pili et al. (2014), for α_{0} = −2 × 10^{−4} and β_{0} = −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 , of the density ρ and of the scalar field χ for a model with a toroidal magnetic field of maximum strength B_{max} = 6.134 × 10^{17}G and central density ρ_{c} = 8.440 × 10^{14} 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_{0} = 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 ρ_{2}, on the 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 (see Appendix C).
It is meaningful to compare our reference model also to an unmagnetised model in STT with the same central density, which is characterised in Table 2 as T_{0}. 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 unmagnetised 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 unmagnetised model T_{0}, clearly showing the magnetic induced deformation on the density profile, that affects mostly the lowdensity 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 unmagnetised case, in the outer part of the star and outside it they are both shallower than in the unmagnetised 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 B_{max} (dashed lines), and with the density of the scalarised and unmagnetised model at the same ρ_{c}, T_{0} (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 unmagnetised model T_{0} (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_{0} = 1.68 M_{⊙}, we plot the deviations Δ of ρ_{c}, M, R_{c} and e with respect to the unmagnetised 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
except for e, in which case we just plot its value, since e(0, M_{0})=0. The results are compared with the GR sequence having the same baryonic mass.
Fig. 4. Variation, with respect to the unmagnetised model, of various quantities along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙} for purely toroidal magnetic field. From left to right, top to bottom: central density ρ_{c}, Komar mass M_{k}, circumferential radius R_{c} 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 descalarised 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 unmagnetised 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^{18} G, almost twice with respect to the one of the equivalant GR sequence, 6.13 × 10^{17} 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^{17} 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 descalarise (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^{17} 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 descalarises. 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 outwardpointing 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^{39}(B_{max}/10^{18}G)^{2}erg up to B_{max} ≃ 10^{18} G. As the magnetisation rises beyond the point where B_{max} = 1.03 × 10^{18} G, the magnetic field energy, in the scalarised part, reaches a maximum of ℋ = 1.65 × 10^{39} erg at B_{max} = 9.7 × 10^{17} G, finally relaxing to the GR profile when the sequence descalarises around B_{max} = 0.58 × 10^{17} G. The scalar charge, instead, drops with increasing magnetisation, being about 10% smaller at B_{max} = 1.05 × 10^{18} G. Beyond this point, the scalar charge drops substantially until the NS completely descalarises.
Fig. 5. Scalar charge Q_{s}, normalised to its value for the unmagnetised model (left panel), and magnetic field energy ℋ (right panel) as functions of B_{max} along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙} and purely toroidal magnetic field. The cyan dotted line highlights the descalarised 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_{0} (middle panel). The lower bound for scalarised models, ρ_{b}, moves to higher densities from ρ_{b} = 5 × 10^{14} g cm^{−3} for Φ = 0 to ρ_{b} = 7.5 × 10^{14} g cm^{−3} for Φ = 2.55 × 10^{30} G cm^{2}, 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 descalarised models. Analogously, the upper bound ρ_{t} for scalarised models increases from ρ_{t} = 3.5 × 10^{15} g cm^{−3} for Φ = 0 to ρ_{t} = 4 × 10^{15} g cm^{−3} for Φ = 1.46 × 10^{30} G cm^{2}, 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 unmagnetised one. However, as the density increases, the trend is reversed and we find magnetised models having a lower mass than the unmagnetised configuration at the same central density. This is reversed again once the density exceeds 2.72 × 10^{15} g cm^{−3} 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^{15} g cm^{−3} the Komar mass is independent of the magnetisation. Quantitatively, the density at which the maximum is reached always increases from ρ_{c} = 2.55 × 10^{15} g cm^{−3} for Φ = 0 to ρ_{c} = 2.95 × 10^{15} g cm^{−3} for Φ = 2.55 × 10^{30} G cm^{2}, while the value of the maximum mass drops initially from 2.08 M_{⊙} to 2.04 M_{⊙} for Φ = 1.46 × 10^{30} G cm^{2} and then rises again to 2.08 M_{⊙} for Φ = 2.55 × 10^{30} G cm^{2}. The full characterisation of the models at maximum mass is given in Table 3.
Fig. 6. Massdensity 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 unmagnetised case (red line). The dotted magenta lines represent the limit for spontaneous scalarisation. Dots mark the position of the maximum mass models UM_{0} (red), TM_{1} (light blue) and TM_{2} (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 unmagnetised one. 
Fig. 7. Sequences for the models with purely toroidal magnetic field and β_{0} = −6. Left panel: massdensity relation computed at fixed B_{max} (blue lines) compared with the unmagnetised sequence (red line). Middle panel: massdensity relation computed at fixed e (green lines) compared with the unmagnetised sequence (red line). Right panel: On top, scalar charge computed at fixed Φ (blue lines) compared with the unmagnetised sequence (red line); on bottom, trace quadrupole deformation e_{s}. 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^{30} G cm^{2}, while the density at which this maximum is reached increases from 2.09 × 10^{15} g cm^{−3} to 2.46 × 10^{15} g cm^{−3}. 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^{15} g cm^{−3} the scalar charge is independent of the magnetisation.
5.3. Poloidal field models with β_{0} = −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^{14} g cm^{−3} and the same maximum value of the magnetic field B_{max} = 6.256 × 10^{17} 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 ρ_{2}, on the sequence . In Fig. 8, we show the distribution of the magnetic field strength , 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 , of the density ρ and of the scalar field χ for a model with a poloidal magnetic field of maximum strength B_{max} = 6.256 × 10^{17} G and central density ρ_{c} = 5.15 × 10^{14} 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_{0} = 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 unmagnetised model with the same central density, characterised in Table 2 under the name P_{0}. 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 unmagnetised model P_{0}. We also plot the profiles of χ, normalised to its maximum value χ_{max}, for the models P and P_{0}. 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 B_{max} (dashed), and with the density of the scalarised and unmagnetised 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 unmagnetised 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^{17} G. In STT we found instead that up to values or order of 1 × 10^{18} 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^{18} 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 descalarise. 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 outwardpointing 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^{18} G, the total magnetic field energy ℋ scales with a good approximation as ℋ = 0.55 × 10^{39}(B_{max}/10^{18}G)^{2} erg, and the scalar charge increases by about 2% with respect to the unmagnetised case. We also found that the magnetic dipole scales as μ = 1.5 × 10^{35}(B_{max}/10^{18}G) erg G^{−1}, 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 unmagnetised model, of various quantities along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙} for purely poloidal magnetic field. From left to right, top to bottom: central density ρ_{c}, Komar mass M_{k}, circumferential radius R_{c} 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_{0} (top panel). The lower bound ρ_{b} for scalarised models now moves to lower densities  from ρ_{b} = 5 × 10^{14} g cm^{−3} for μ = 0 to ρ_{b} = 4.3 × 10^{14} g cm^{−3} for μ = 1.57 × 10^{35} erg G^{−1} – 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 descalarised models. Analogously, the upper bound ρ_{t} for scalarised models decreases  from ρ_{t} = 3.5 × 10^{15} g cm^{−3} for μ = 0 to ρ_{t} = 3.42 × 10^{15} g cm^{−3} for μ = 0.45 × 10^{35} erg G^{−1} – 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 unmagnetised 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^{15} g cm^{−3} 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^{35} erg G^{−1}, while the density at which the maximum is reached drops to 1.96 × 10^{15} g cm^{−3}. 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^{15} g cm^{−3} the scalar charge is independent on the magnetisation.
Fig. 11. Massdensity 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 unmagnetised case (red line). The dotted magenta lines represent the limit for spontaneous scalarisation. Dots mark the position of the maximum mass models UM_{0} (red), PM_{1} (light blue) and PM_{2} (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 unmagnetised one. 
Fig. 12. Sequences for models with purely poloidal magnetic field and β_{0} = −6. Left panel: massdensity relation computed at fixed B_{max} (blue lines) compared with the unmagnetised sequence (red line). Middle panel: massdensity relation computed at fixed e (green lines) compared with the unmagnetised sequence (red line). Right panel: on top, scalar chargedensity relation computed at fixed μ (blue lines) compared with the unmagnetised sequence (red line); on bottom, trace quadrupole deformation e_{s}. 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 β_{0} = −5
In order to understand how our results depend on the specific choice of the STT parameter β_{0}, we have computed equilibrium configurations also for β_{0} = −5 and β_{0} = −4.5, closer to the limit for spontaneous scalarisation, both in the case of pure toroidal and purely poloidal magnetic fields. For β_{0} = −5, the unmagnetised 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 unmagnetised values, at fixed baryonic mass M_{o} = 1.680 M_{⊙}, as was done for β_{0} = −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 β_{0} = −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^{17} G, before it descalarises, and then reaches a new maximum corresponding to the GR value of 6.13 × 10^{17} G. We can conclude that in STTs with β_{0} > −5 the upper limit to B_{max} is reached after the solution descalarises, while for β_{0} < −5 it is reached for scalarised configurations. On the other hand in models with a purely poloidal magnetic field, we observe no evidence for descalarisation with increasing k_{pol}. However, there seems to be an asymptote to a maximum value of B_{max} of ≃7.5 × 10^{17} G, slightly higher than in GR for β_{0} = −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 nonlinear way on the values of β_{0}.
Fig. 13. Value of the quadrupole deformation e along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙}, as a function of B_{max}, for β_{0} = −5 (blue lines) vs GR (red lines). The cyan dotted line highlights the unscalarised 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_{0}. The region of descalarisation ρ_{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^{14} g cm^{−3} for Φ = 0 to ρ_{c} = 1.06 × 10^{15} g cm^{−3} for Φ = 2 × 10^{30} G cm^{2}, 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 descalarised models. Analogously, the upper bound ρ_{t} for scalarised models increases, from ρ_{b} = 2.65 × 10^{15} g cm^{−3} for Φ = 0 to ρ_{c} = 3.05 × 10^{15} g cm^{−3} for Φ = 2 × 10^{30} G cm^{2}, 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^{30}G cm^{2}. In this case we also see that on sequences with Φ ≥ 1.64 × 10^{30} G cm^{2} the mass of equilibrium models is always larger than the relative unmagnetised 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 unmagnetised case (red line); on bottom, mass difference with respect to the unmagnetised case. Right panel: on top, scalar charge on sequences at fixed Φ; on bottom, trace quadrupole e_{s} 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 β_{0} = −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 descalarised 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 unmagnetised case (red line); on bottom, mass difference with respect to the unmagnetised case. Right panel: on top, scalar charge on sequences at fixed μ; on bottom, trace quadrupole e_{s} 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 β_{0} = −4.5
We consider here the case β_{0} = −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 unmagnetised models, ρ_{b} = 9.3 × 10^{14} g cm^{−3} and ρ_{t} = 2.0 × 10^{15} g cm^{−3}, 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 massdensity relation becomes progressively less evident: already at Φ = 0.9 × 10^{30} G cm^{2} 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 unmagnetised 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 unmagnetised 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 descalarisation.
This difference, in part already present at lower β_{0}, can be understood if one recalls that spontaneous scalarisation can be seen, from a dynamical point of view, as an instability (Damour & EspositoFarèse 1996), which can be excited only if the minimum wavelength of unstable modes (a function of β_{0}) 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 β_{0} ≈ −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 nonaxisymmetric 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 nonrotating and unmagnetised NS is
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 fluxfreezing 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
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 descalarised 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 descalarised part and one to the scalarised one. In the massdensity diagram there is a region where models are gravitationally unstable. Moreover, for any given value of Φ, there is a range of masses where both descalarised and scalarised solutions are possible. On the other hand, there is a lower limit to the values of the magnetic flux that can support descalarised 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 descalarised 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 descalarised configuration with M_{0} = 1.68 M_{⊙} can only exist for Φ > 2.06 × 10^{30} G cm^{2} and ρ_{c} < 7.12 × 10^{14} g cm^{−3}; 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^{15} g cm^{−3}, 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 β_{0} and found that this effect already disappears at β_{0} = −5. However, for β_{0} = −4.5 we found that two configurations, one scalarised and the other unscalarised, 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 unmagnetised solution. From bottom to top: other curves are computed at Φ = [2.55, 2.06, 1.46, 0.91]×10^{30} G cm^{2}. The various parts are: gravitationally stable descalarised 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 descalarised configuration with M_{0} = 1.68 M_{⊙} and Φ = 2.06 × 10^{30} G cm^{2}, while the arrow points to the blue dot where the configuration is expected to jump because of magneticallyinduced 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 descalarised 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 “magneticallyinduced 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 scalartensor 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 nonlinear 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 Einsteinframe (where the metric equation have the same mathematical structure as in GR, and the same numerical schemes can be applied) and the Jordanframe (where the magnetofluid equations retain their conservative, quasihyperbolic 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} ≠ 0. In particular, we verified that our numerical algorithm always selects the branch. We also showed that the solutions are not always the ones with the largest deviation from GR in the massdensity 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 ψ^{6} 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 Eframe, 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 unmagnetised 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 descalarise. This effect was evaluated for various values of β_{0} showing that there is a strong dependency.
We have then shown, using various parametrisations, how the massdensity 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 descalarised 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 counterbalance 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 β_{0}, 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 descalarisation, 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 β_{0} 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 ‘magneticallyinduced 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 β_{0} = −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 β_{0}, 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 nonminimally 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 nonlinear current terms in magnetised NSs (Pili et al. 2014), and spontaneous magneticvectorisation 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.
In analogy with the nonrelativistic case, the relativistic Bernoulli integral can be defined, in hydrodynamics, from the conservation law of hu_{t} along the trajectories of a stationary flow (see Friedman & Stergioulas 2013). This is a special case of the global first integral of Euler’s equation for isoentropic flows which, for stationary cases, reduces to Eq. (39). This is the reason why we refer to Eq. (39) as the generalised Bernoulli integral.
Acknowledgments
The authors acknowledge financial support from the “Accordo Attuativo ASIINAF n. 201714H.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
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 119, 161101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 848, L12 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJ, 848, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 851, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Alcubierre, M. 2008, Introduction to 3+1 Numerical Relativity, International Series of Monographs on Physics (Oxford: OUP) [CrossRef] [Google Scholar]
 Anderson, D., & Yunes, N. 2019, CQG, 36, 165003 [CrossRef] [Google Scholar]
 Anderson, D., Freire, P., & Yunes, N. 2019, CQG, 36, 225009 [CrossRef] [Google Scholar]
 Andreou, N., Franchini, N., Ventagli, G., & Sotiriou, T. P. 2019, Phys. Rev. D, 99, 124022 [CrossRef] [Google Scholar]
 Barausse, E. 2017, Proceedings of The 3rd International Symposium on Quest for the Origin of Particles and the Universe (SISSA Medialab), 294, 029 [CrossRef] [Google Scholar]
 Barausse, E., Palenzuela, C., Ponce, M., & Lehner, L. 2013, Phys. Rev. D, 87, 081506 [CrossRef] [Google Scholar]
 Bars, I., & Pope, C. N. 1989, Gen. Relativ. Gravity, 21, 545 [CrossRef] [Google Scholar]
 Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509 [NASA ADS] [CrossRef] [Google Scholar]
 Bekenstein, J. D., & Oron, A. 2001, Found. Phys., 31, 895 [CrossRef] [Google Scholar]
 Bergmann, P. G. 1968, Int. J. Theor. Phys., 1, 25 [CrossRef] [Google Scholar]
 Berti, E., Barausse, E., Cardoso, V., et al. 2015, CQG, 32, 243001 [CrossRef] [Google Scholar]
 Bertolami, O., & Páramos, J. 2016, Gen. Relativ. Gravity, 48, 34 [CrossRef] [Google Scholar]
 Bocquet, M., Bonazzola, S., Gourgoulhon, E., & Novak, J. 1995, A&A, 301, 757 [NASA ADS] [Google Scholar]
 Bonanno, A., Rezzolla, L., & Urpin, V. 2003, A&A, 410, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Braithwaite, J. 2009, MNRAS, 397, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Braithwaite, J., & Nordlund, A. 2006, A&A, 450, 1077 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Braithwaite, J., & Spruit, H. C. 2006, A&A, 450, 1097 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brans, C., & Dicke, R. H. 1961, Phys. Rev., 124, 925 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Brax, P., Davis, A.C., & Jha, R. 2017, Phys. Rev. D, 95, 083514 [CrossRef] [Google Scholar]
 Brown, J. D. 1993, CQG, 10, 1579 [CrossRef] [Google Scholar]
 Bucciantini, N., & Del Zanna, L. 2011, A&A, 528, A101 [CrossRef] [EDP Sciences] [Google Scholar]
 Bucciantini, N., & Del Zanna, L. 2013, MNRAS, 428, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Bucciantini, N., Pili, A. G., & Zanna, L. D. 2015, MNRAS, 447, 1 [CrossRef] [Google Scholar]
 Buchdahl, H. A. 1970, MNRAS, 150, 1 [NASA ADS] [Google Scholar]
 Burrows, A., Dessart, L., Livne, E., Ott, C. D., & Murphy, J. 2007, ApJ, 664, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Camelio, G., Dietrich, T., Marques, M., & Rosswog, S. 2019, Phys. Rev. D, 100, 123001 [CrossRef] [Google Scholar]
 Capozziello, S., & de Laurentis, M. 2011, Phys. Rep., 509, 167 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cartan, É. 1986, On Manifolds with an Affine Connection and the Theory of General Relativity, Monographs and Textbooks in Physical Science (Bibliopolis) [Google Scholar]
 Carter, B. 1969, J. Mat. Phys., 10, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, B. 1970, Commun. Mat. Phys., 17, 233 [CrossRef] [Google Scholar]
 Carter, B. 2009, Gen. Relativ. Gravity, 41, 2873 [CrossRef] [Google Scholar]
 Carter, B. 2010, Gen. Relativ. Gravity, 42, 653 [CrossRef] [Google Scholar]
 Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 116 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chatterjee, D., Elghozi, T., Novak, J., & Oertel, M. 2015, MNRAS, 447, 3785 [NASA ADS] [CrossRef] [Google Scholar]
 Ciolfi, R., & Rezzolla, L. 2013, MNRAS, 435, L43 [CrossRef] [Google Scholar]
 Ciolfi, R., Kastaun, W., Kalinani, J. V., & Giacomazzo, B. 2019, Phys. Rev. D, 100, 023005 [CrossRef] [Google Scholar]
 CorderoCarrión, I., CerdáDurán, P., Dimmelmeier, H., et al. 2009, Phys. Rev. D, 79, 024017 [NASA ADS] [CrossRef] [Google Scholar]
 Dall’Osso, S., Shore, S. N., & Stella, L. 2009, MNRAS, 398, 1869 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., & EspositoFarèse, G. 1993, Phys. Rev. Lett., 70, 2220 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., & EspositoFarèse, G. 1996, Phys. Rev. D, 54, 1474 [NASA ADS] [CrossRef] [Google Scholar]
 Damour, T., Piazza, F., & Veneziano, G. 2002, Phys. Rev. Lett., 89, 081601 [Google Scholar]
 Das, U., & Mukhopadhyay, B. 2015, J. Cosomology Astropart. Phys., 2015, 016 [Google Scholar]
 De Felice, A., & Tanaka, T. 2010, Prog. Theor. Phys., 124, 503 [CrossRef] [Google Scholar]
 De Felice, A., & Tsujikawa, S. 2010, Living Rev. Relativ., 13, 3 [CrossRef] [Google Scholar]
 De Felice, A., Hindmarsh, M., & Trodden, M. 2006, J. Cosomology Astropart. Phys., 2006, 005 [Google Scholar]
 Del Zanna, L., & Chiuderi, C. 1996, A&A, 310, 341 [NASA ADS] [Google Scholar]
 Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, A&A, 473, 11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Del Zanna, L., Papini, E., Landi, S., Bugli, M., & Bucciantini, N. 2016, MNRAS, 460, 3753 [NASA ADS] [CrossRef] [Google Scholar]
 Delphenich, D. 2005, 6th International Conference on Symmetry in Nonlinear Mathematical Physics (SNMP 05) Kiev, Ukraine, June 20–26, 2005 [Google Scholar]
 Del Zanna, L., & Bucciantini, N. 2018, MNRAS, 479, 657 [NASA ADS] [Google Scholar]
 Deser, S. 2000, Ann. Phys. (Berl.), 9, 299 [CrossRef] [Google Scholar]
 Dimmelmeier, H., Stergioulas, N., & Font, J. A. 2006, MNRAS, 368, 1609 [NASA ADS] [CrossRef] [Google Scholar]
 Dohi, A., Kase, R., Kimura, R., & Yamamoto, K. 2020, ArXiv eprints [arXiv: 2003.12571] [Google Scholar]
 Doneva, D. D., & Yazadjiev, S. S. 2016, J. Cosomology Astropart. Phys., 2016, 019 [Google Scholar]
 Doneva, D. D., & Yazadjiev, S. S. 2020, Phys. Rev. D, 101, 064072 [CrossRef] [Google Scholar]
 Doneva, D. D., Yazadjiev, S. S., Stergioulas, N., & Kokkotas, K. D. 2013, Phys. Rev. D, 88, 084060 [CrossRef] [Google Scholar]
 Doneva, D. D., Yazadjiev, S. S., Stergioulas, N., & Kokkotas, K. D. 2018, Phys. Rev. D, 98, 104039 [CrossRef] [Google Scholar]
 Faraoni, V. 2004, Cosmology in ScalarTensor Gravity, Fundamental Theories of Physics (Netherlands: Springer) [CrossRef] [Google Scholar]
 Ferrario, L., Melatos, A., & Zrake, J. 2015, Space Sci. Rev., 191, 77 [CrossRef] [Google Scholar]
 Ferraro, V. C. A. 1954, ApJ, 119, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Franzon, B., Dexheimer, V., & Schramm, S. 2016, Phys. Rev. D, 94, 044018 [NASA ADS] [CrossRef] [Google Scholar]
 Freire, P. C. C., Wex, N., EspositoFarèse, G., et al. 2012, MNRAS, 423, 3328 [NASA ADS] [CrossRef] [Google Scholar]
 Fricke, K. 1969, A&A, 1, 388 [NASA ADS] [Google Scholar]
 Frieben, J., & Rezzolla, L. 2012, MNRAS, 427, 3406 [NASA ADS] [CrossRef] [Google Scholar]
 Friedman, J. L., & Stergioulas, N. 2013, Rotating Relativistic Stars, Cambridge Monographs on Mathematical Physics (Cambridge University Press) [CrossRef] [Google Scholar]
 Fujii, Y., & Maeda, K.I. 2003, The ScalarTensor Theory of Gravitation [Google Scholar]
 Fujisawa, K., & Eriguchi, Y. 2015, PASJ, 67, 53 [CrossRef] [Google Scholar]
 Gao, H., Zhang, B., & Lü, H.J. 2016, Phys. Rev. D, 93, 044065 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., Sperhake, U., & Ott, C. D. 2016, CQG, 33, 135002 [CrossRef] [Google Scholar]
 Gong, Y., Papantonopoulos, E., & Yi, Z. 2018, Eur. Phys. J C, 78, 738 [CrossRef] [EDP Sciences] [Google Scholar]
 Gourgoulhon, É. 2012, 3+1 Formalism in General Relativity: Bases of Numerical Relativity: Lecture Notes in Physics (Berlin Heidelberg: Springer) [CrossRef] [Google Scholar]
 Green, M. B., Schwartz, J. H., & Witten, E. 1988, Astron. Nachr., 309, 297 [CrossRef] [Google Scholar]
 Hagihara, Y., Era, N., Iikawa, D., Takeda, N., & Asada, H. 2020, Phys. Rev. D, 101, 041501 [CrossRef] [Google Scholar]
 Harada, T. 1998, Phys. Rev. D, 57, 4802 [CrossRef] [Google Scholar]
 Hawking, S. W. 1972, Commun. Mat. Phys., 25, 167 [CrossRef] [MathSciNet] [Google Scholar]
 Hawking, S. W., & Ellis, G. F. R. 1973, The Large Scale Structure of SpaceTime, Cambridge Monographs on Mathematical Physics (Cambridge University Press) [CrossRef] [Google Scholar]
 Heisenberg, L. 2014, J. Cosomology Astropart. Phys., 2014, 015 [CrossRef] [Google Scholar]
 Hellings, R. W., & Nordtvedt, K. 1973, Phys. Rev. D, 7, 3593 [CrossRef] [Google Scholar]
 Iosif, P., & Stergioulas, N. 2014, Gen. Relativ. Gravity, 46, 1800 [Google Scholar]
 Isenberg, J. A. 2008, Int. J. Mod. Phys. D, 17, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Just, K. 1959, Z. Naturforsch. A, 14, 751 [CrossRef] [Google Scholar]
 Kase, R., Minamitsuji, M., & Tsujikawa, S. 2018, Phys. Rev. D, 97, 084009 [CrossRef] [Google Scholar]
 Kase, R., Minamitsuji, M., & Tsujikawa, S. 2020, Phys. Rev. D, 102, 024067 [CrossRef] [Google Scholar]
 Kawamura, T., Giacomazzo, B., Kastaun, W., et al. 2016, Phys. Rev. D, 94, 064012 [NASA ADS] [CrossRef] [Google Scholar]
 Kiuchi, K., & Yoshida, S. 2008, Phys. Rev. D, 78, 044045 [NASA ADS] [CrossRef] [Google Scholar]
 Kiuchi, K., Kotake, K., & Yoshida, S. 2009, ApJ, 698, 541 [CrossRef] [Google Scholar]
 Konno, K. 2001, A&A, 372, 594 [CrossRef] [EDP Sciences] [Google Scholar]
 Kundt, W., & Trümper, M. 1966, Z. Phys. A, 192, 419 [CrossRef] [Google Scholar]
 Langlois, D., Saito, R., Yamauchi, D., & Noui, K. 2018, Phys. Rev. D, 97, 061501 [CrossRef] [Google Scholar]
 Lasky, P. D., Zink, B., Kokkotas, K. D., & Glampedakis, K. 2011, ApJ, 735, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Lovelock, D. 1971, J. Mat. Phys., 12, 498 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Margalit, B., & Metzger, B. D. 2017, ApJ, 850, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuda, T., & Nariai, H. 1973, Prog. Theor. Phys., 49, 1195 [CrossRef] [Google Scholar]
 Mendes, R. F., & Ortiz, N. 2016, Phys. Rev. D, 93, 124035 [CrossRef] [Google Scholar]
 Metzger, B. D., Giannios, D., Thompson, T. A., Bucciantini, N., & Quataert, E. 2011, MNRAS, 413, 2031 [NASA ADS] [CrossRef] [Google Scholar]
 Miketinac, M. J. 1975, Ap&SS, 35, 349 [CrossRef] [Google Scholar]
 Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24 [Google Scholar]
 Monaghan, J. J. 1965, MNRAS, 131, 105 [CrossRef] [Google Scholar]
 Monaghan, J. J. 1966, MNRAS, 134, 275 [CrossRef] [Google Scholar]
 Nordtvedt, K. 1970, ApJ, 161, 1059 [CrossRef] [MathSciNet] [Google Scholar]
 Novak, J. 1998a, Phys. Rev. D, 58, 064019 [NASA ADS] [CrossRef] [Google Scholar]
 Novak, J. 1998b, Phys. Rev. D, 57, 4789 [CrossRef] [Google Scholar]
 Oppenheimer, J. R., & Volkoff, G. M. 1939, Phys. Rev., 55, 374 [NASA ADS] [CrossRef] [Google Scholar]
 Oron, A. 2002, Phys. Rev. D, 66, 023006 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, J. P., & Hartwick, F. D. A. 1968, ApJ, 153, 797 [NASA ADS] [CrossRef] [Google Scholar]
 Ott, C. D., Dimmelmeier, H., Marek, A., et al. 2007, CQG, 24, S139 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., & Freire, P. 2016, ARA&A, 54, 401 [Google Scholar]
 Pani, P., & Berti, E. 2014, Phys. Rev. D, 90, 024025 [CrossRef] [Google Scholar]
 Papantonopoulos, E. 2015, Modifications of Einstein’s Theory of Gravity at Large Distances, Lecture Notes in Physics (Springer International Publishing) [Google Scholar]
 Papitto, A., Torres, D. F., Rea, N., & Tauris, T. M. 2014, A&A, 566, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pappas, G., Doneva, D. D., Sotiriou, T. P., Yazadjiev, S. S., & Kokkotas, K. D. 2019, Phys. Rev. D, 99, 104014 [CrossRef] [Google Scholar]
 Peebles, P. J. E., & Ratra, B. 2003, Rev. Mod. Phys., 75, 559 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Pili, A. G., Bucciantini, N., & Del Zanna, L. 2014, MNRAS, 439, 3541 [CrossRef] [Google Scholar]
 Pili, A. G., Bucciantini, N., & Del Zanna, L. 2015, MNRAS, 447, 2821 [CrossRef] [Google Scholar]
 Pili, A. G., Bucciantini, N., & Del Zanna, L. 2017, MNRAS, 470, 2469 [CrossRef] [Google Scholar]
 Pili, A. G., Bucciantini, N., Drago, A., Pagliara, G., & Del Zanna, L. 2016, MNRAS, 462, L26 [CrossRef] [Google Scholar]
 Popov, S. B. 2016, A&AT, 29, 183 [Google Scholar]
 Prendergast, K. H. 1956, ApJ, 123, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., & Rosswog, S. 2006, Science, 312, 719 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Quiros, I. 2019, Int. J. Mod. Phys. D, 28, 1930012 [NASA ADS] [CrossRef] [Google Scholar]
 Raithel, C. A., Özel, F., & Psaltis, D. 2018, ApJ, 857, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Ramazanoğlu, F. M. 2017, Phys. Rev. D, 96, 064009 [CrossRef] [Google Scholar]
 Ramazanoğlu, F. M., & Pretorius, F. 2016, Phys. Rev. D, 93, 064005 [CrossRef] [Google Scholar]
 Rheinhardt, M., & Geppert, U. 2005, A&A, 435, 201 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roberts, P. H. 1955, ApJ, 122, 508 [CrossRef] [Google Scholar]
 Rowlinson, A., O’Brien, P. T., Metzger, B. D., Tanvir, N. R., & Levan, A. J. 2013, MNRAS, 430, 1061 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W. 1966, MNRAS, 132, 347 [NASA ADS] [Google Scholar]
 Salgado, M. 2006, CQG, 23, 4719 [CrossRef] [Google Scholar]
 Salgado, M., Martínez del Río, D., Alcubierre, M., & Núñez, D. 2008, Phys. Rev. D, 77, 104010 [CrossRef] [Google Scholar]
 Salgado, M., Sudarsky, D., & Nucamendi, U. 1998, Phys. Rev. D, 58, 124003 [CrossRef] [Google Scholar]
 Santiago, D. I., & Silbergleit, A. S. 2000, Gen. Relativ. Gravitation, 32, 565 [Google Scholar]
 Schubert, G. 1968, ApJ, 151, 1099 [CrossRef] [Google Scholar]
 Schärer, A., Angélil, R., Bondarescu, R., Jetzer, P., & Lundgren, A. 2014, Phys. Rev. D, 90, 123005 [NASA ADS] [CrossRef] [Google Scholar]
 Shao, L., Sennett, N., Buonanno, A., Kramer, M., & Wex, N. 2017, Phys. Rev. X, 7, 041025 [Google Scholar]
 Shibata, M., & Sekiguchi, Y.I. 2005, Phys. Rev. D, 72, 044014 [NASA ADS] [CrossRef] [Google Scholar]
 Shibata, M., Taniguchi, K., Okawa, H., & Buonanno, A. 2014, Phys. Rev. D, 89, 084005 [CrossRef] [Google Scholar]
 Silva, H. O., Macedo, C. F. B., Berti, E., & Crispino, L. C. B. 2015, CQG, 32, 145008 [CrossRef] [Google Scholar]
 Silva, H. O., Sakstein, J., Gualtieri, L., Sotiriou, T. P., & Berti, E. 2018, Phys. Rev. Lett., 120, 131104 [CrossRef] [Google Scholar]
 Sotani, H. 2012, Phys. Rev. D, 86, 124036 [CrossRef] [Google Scholar]
 Sotani, H., & Kokkotas, K. D. 2005, Phys. Rev. D, 71, 124038 [CrossRef] [Google Scholar]
 Sotiriou, T. P. 2006, CQG, 23, 5117 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Spruit, H. C. 2009, in Cosmic Magnetic Fields: From Planets, to Stars and Galaxies, eds. K. G. Strassmeier, A. G. Kosovichev, & J. E. Beckman, IAU Symp., 259, 61 [Google Scholar]
 Staykov, K. V., Popchev, D., Doneva, D. D., & Yazadjiev, S. S. 2018, Eur. Phys. J C, 78, 586 [CrossRef] [EDP Sciences] [Google Scholar]
 Staykov, K. V., Doneva, D. D., Popchev, D., & Yazadjiev, S. S. 2019, Am. Inst. Phys. Conf. Ser., 2075, 040006 [Google Scholar]
 Suvorov, A. G. 2018, Phys. Rev. D, 98, 084026 [CrossRef] [Google Scholar]
 Taniguchi, K., Shibata, M., & Buonanno, A. 2015, Phys. Rev. D, 91, 024033 [CrossRef] [Google Scholar]
 Tayler, R. J. 1973, MNRAS, 161, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Tolman, R. C. 1939, Phys. Rev., 55, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Tomei, N., Del Zanna, L., Bugli, M., & Bucciantini, N. 2020, MNRAS, 491, 2346 [Google Scholar]
 Tomimura, Y., & Eriguchi, Y. 2005, MNRAS, 359, 1117 [NASA ADS] [CrossRef] [Google Scholar]
 Touboul, P., Métris, G., Rodrigues, M., et al. 2017, Phys. Rev. Lett., 119, 231101 [NASA ADS] [CrossRef] [Google Scholar]
 Trimble, V. 1987, ARA&A, 25, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Uryū, K., Gourgoulhon, E., Markakis, C. M., et al. 2014, Phys. Rev. D, 90, 101501 [NASA ADS] [CrossRef] [Google Scholar]
 Uryū, K., b. o., Yoshida, S., Gourgoulhon, E., et al. 2019, Phys. Rev. D, 100, 123019 [CrossRef] [Google Scholar]
 van Dantzig, D., & Dirac, P. A. M. 1934, Math. Proc. Camb. Philos. Soc., 30, 421 [CrossRef] [Google Scholar]
 Wagoner, R. V. 1970, Phys. Rev. D, 1, 3209 [NASA ADS] [CrossRef] [Google Scholar]
 Will, C. M. 2014, Living Rev., Relativ, 17 [Google Scholar]
 Wilson, J. R., & Mathews, G. J. 2003, Relativistic Numerical Hydrodynamics, Cambridge Monographs on Mathematical Physics (Cambridge University Press) [CrossRef] [Google Scholar]
 Wilson, J. R., Mathews, G. J., & Marronetti, P. 1996, Phys. Rev. D, 54, 1317 [NASA ADS] [CrossRef] [Google Scholar]
 Woltjer, L. 1960, ApJ, 131, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, G. E. 1973, MNRAS, 162, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Yakovlev, D. G., Gnedin, O. Y., Gusakov, M. E., et al. 2005, Nucl. Phys. A, 752, 590 [NASA ADS] [CrossRef] [Google Scholar]
 Yazadjiev, S. S. 2012, Phys. Rev. D, 85, 044030 [CrossRef] [Google Scholar]
 Yazadjiev, S. S., Doneva, D. D., & Popchev, D. 2016, Phys. Rev. D, 93, 084038 [CrossRef] [Google Scholar]
 Yoshida, S., Yoshida, S., & Eriguchi, Y. 2006, ApJ, 651, 462 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, X., Niu, R., & Zhao, W. 2019, Phys. Rev. D, 100, 024038 [CrossRef] [Google Scholar]
Appendix A: XCFC for a rotating NS
We show here how the standard techniques of XNS, based on the XCFC approach to the solution of the metric functions, can be adapted to take into account the presence of a scalar field. For simplicity we are going to consider here only unmagnetised rotators. The generalisation to magnetised ones is trivial and strictly follows what was done in Pili et al. (2017). In the Eframe the standard set of XCFC equations is:
where f_{ij} is the flat 3metric, is the flat covariant derivative (), and is the usual Laplacian operator in flat 3space. Δ_{L} is defined as:
and
The source terms come from the 3+1 decomposition of the energymomentum tensor in the Eframe:
For stationary (∂_{t} = 0) and axisymmetric (∂_{ϕ} = 0) configurations, for the metric given by Eq. (32) (where the only non vanishing component of the shift vector is β^{ϕ}), assuming that the only non vanishing component of the velocity is v^{ϕ}, it can be shown that
where e, p, and are all in the Jframe and is instead in the Eframe.
If on a time slice the values of the physical quantities are provided as well as the scalar field, then the XCFC set of equations can be solved for the metric component in the Eframe. It is evident that the XCFC scheme retains its main interesting property of decoupling the various equations, allowing to solve them separately, one after the other. This holds also in the more general time dependent case. In fact 3+1 schemes for GRHD and MHD evolve the conserved quantities in the Jframe and . Combined with a scheme that evolves also the 3+1 components of the scalar field P and Q^{i}, the XCFC equations can then be used to solve for the metric.
Appendix B: The STOV system
The STOV system of equations can be derived setting B^{i} = 0 in Eqs. (33), (34), (36), and (38):
These must be supplemented by a barotropic EoS p = p(ρ), ε = ε(ρ). This system can be solved, given the value at r = 0 of the density ρ_{c}, the conformal factor ψ_{c} and the scalar field χ_{c} (recalling that all radial derivatives of scalar quantities vanish in r = 0). The value of the lapse function at the center, α_{c}, is irrelevant to the solution per se, since only its derivative appears in Eqs. (B.1)–(B.6). This means that the lapse function is derived minus an arbitrary constant, which is then chosen in order to satisfy the correct asymptotic behaviour at r → ∞.
The correct STT solution satisfies the following requirements:

the ratio C = αQ_{r}/2∂_{r}α must be constant outside the NS, because it can be shown that it is equal to the ratio Q_{s}/2M between the net scalar charge and twice the Komar mass in the Eframe;

in vacuum α and ψ must behave like the Just metric (Just 1959) in isotropic coordinates.
Given that the Just metric in isotropic coordinates has no analytical form, we provide here an approximation that proves to be accurate with a precision ∼10^{−4}, already at a couple of NS radii. If one writes the metric terms ψ and α, outside of the NS surface, as
one finds that the first values of m_{i} for i > 0 are:
and n_{i} = (− 1)^{i}m_{i}. When Q_{s} = 0 one finds m_{0} = n_{0} and m_{i} = n_{i} = 0 for i > 0, recovering the GR solution.
Appendix C: Global quantities
In this appendix we list the main global quantities used in the paper. We give their general form, valid also in the case of a nonstatic, but stationary, spacetime.
The Komar mass in the Eframe is
where Σ_{t} is a spacelike hypersurface of constant coordinate time t and ξ^{ν} is the timelike Killing vector associated to the stationarity of the spacetime. In our static case it reduces to
The baryonic mass, which is the same in the Eframe and in the Jframe, is
which in our case is
The proper mass in the Eframe is
The scalar charge of the star in the Eframe, , is defined as the monopole component of the scalar field at asymptotically large radii:
By integrating Eq. (13) over a spherical volume of asymptotically large radius, using Stokes’ theorem and using the fact that outside the star’s surface, we obtain
where T_{p} = 3p − ε − ρ. The circumferential radius in the Jframe is
The magnetic energy in the Jframe is
which reduces to
The binding energy of the star in the Eframe is defined as
The flux of the toroidal magnetic field, which is the same in the Eframe and in the Jframe, is
The magnetic dipole moment in the Jframe is
The value converges already for r ≃ 5−10R_{c}.
The quadrupole deformation of the star in the Eframe is defined as
where the physical and scalar field moments of inertia around the polar axis z and the x axis are, respectively,
We note that we defined the moments of inertia of the scalar field in the same way as the usual physical ones: as integrals of the energy density .
The quadrupolar deformation of the trace is related to the quadrupolar and monopolar distributions of the scalar field at asimptotically large radii, and is defined as
We note that the denominator is the Newtonian equivalent of .
All Tables
Values of various physical quantities describing the solutions 𝒮_{w}, and , for α_{0} = −0.05 and β_{0} = −6, and for selected values of the central density ρ_{c} (in the Jframe), corresponding from top to bottom to: ρ_{b} < ρ_{c} < ρ_{1}, ρ_{3} < ρ_{c} < ρ_{t}, ρ_{c} = ρ_{1}, ρ_{2}, ρ_{3}.
All Figures
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 and . 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. 

In the text 
Fig. 2. From left to right: meridional distribution of the magnetic field strength , of the density ρ and of the scalar field χ for a model with a toroidal magnetic field of maximum strength B_{max} = 6.134 × 10^{17}G and central density ρ_{c} = 8.440 × 10^{14} 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 the text 
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 B_{max} (dashed lines), and with the density of the scalarised and unmagnetised model at the same ρ_{c}, T_{0} (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 unmagnetised model T_{0} (dotted purple). 

In the text 
Fig. 4. Variation, with respect to the unmagnetised model, of various quantities along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙} for purely toroidal magnetic field. From left to right, top to bottom: central density ρ_{c}, Komar mass M_{k}, circumferential radius R_{c} 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 descalarised 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 unmagnetised model in GR. The arrows show the direction of increasing magnetisation. 

In the text 
Fig. 5. Scalar charge Q_{s}, normalised to its value for the unmagnetised model (left panel), and magnetic field energy ℋ (right panel) as functions of B_{max} along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙} and purely toroidal magnetic field. The cyan dotted line highlights the descalarised configurations. The arrows show the direction of increasing magnetisation. 

In the text 
Fig. 6. Massdensity 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 unmagnetised case (red line). The dotted magenta lines represent the limit for spontaneous scalarisation. Dots mark the position of the maximum mass models UM_{0} (red), TM_{1} (light blue) and TM_{2} (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 unmagnetised one. 

In the text 
Fig. 7. Sequences for the models with purely toroidal magnetic field and β_{0} = −6. Left panel: massdensity relation computed at fixed B_{max} (blue lines) compared with the unmagnetised sequence (red line). Middle panel: massdensity relation computed at fixed e (green lines) compared with the unmagnetised sequence (red line). Right panel: On top, scalar charge computed at fixed Φ (blue lines) compared with the unmagnetised sequence (red line); on bottom, trace quadrupole deformation e_{s}. In all panels, the dotted magenta lines represent the limit for spontaneous scalarisation and the yellow square represents model T of Fig. 2. 

In the text 
Fig. 8. From left to right: meridional distribution of the magnetic field strength , of the density ρ and of the scalar field χ for a model with a poloidal magnetic field of maximum strength B_{max} = 6.256 × 10^{17} G and central density ρ_{c} = 5.15 × 10^{14} 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 the text 
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 B_{max} (dashed), and with the density of the scalarised and unmagnetised 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 unmagnetised model P0 (dotted purple). 

In the text 
Fig. 10. Variation, with respect to the unmagnetised model, of various quantities along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙} for purely poloidal magnetic field. From left to right, top to bottom: central density ρ_{c}, Komar mass M_{k}, circumferential radius R_{c} 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 the text 
Fig. 11. Massdensity 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 unmagnetised case (red line). The dotted magenta lines represent the limit for spontaneous scalarisation. Dots mark the position of the maximum mass models UM_{0} (red), PM_{1} (light blue) and PM_{2} (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 unmagnetised one. 

In the text 
Fig. 12. Sequences for models with purely poloidal magnetic field and β_{0} = −6. Left panel: massdensity relation computed at fixed B_{max} (blue lines) compared with the unmagnetised sequence (red line). Middle panel: massdensity relation computed at fixed e (green lines) compared with the unmagnetised sequence (red line). Right panel: on top, scalar chargedensity relation computed at fixed μ (blue lines) compared with the unmagnetised sequence (red line); on bottom, trace quadrupole deformation e_{s}. In all panels, the dotted magenta lines represent the limit for spontaneous scalarisation and the yellow square represents model P of Fig. 8. 

In the text 
Fig. 13. Value of the quadrupole deformation e along the equilibrium sequence with constant M_{0} = 1.68 M_{⊙}, as a function of B_{max}, for β_{0} = −5 (blue lines) vs GR (red lines). The cyan dotted line highlights the unscalarised configurations. Left panel: purely toroidal magnetic field; right panel: purely poloidal magnetic field. The arrows show the direction of increasing magnetisation. 

In the text 
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 unmagnetised case (red line); on bottom, mass difference with respect to the unmagnetised case. Right panel: on top, scalar charge on sequences at fixed Φ; on bottom, trace quadrupole e_{s} on the same sequences. The dotted magenta lines represent the limit for spontaneous scalarisation. 

In the text 
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 unmagnetised case (red line); on bottom, mass difference with respect to the unmagnetised case. Right panel: on top, scalar charge on sequences at fixed μ; on bottom, trace quadrupole e_{s} on the same sequences. The dotted magenta lines represent the limit for spontaneous scalarisation. 

In the text 
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 unmagnetised 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 unmagnetised case (red line). Bottom panel: value of the scalar charge on the same sequences at fixed μ. 

In the text 
Fig. 17. Sequences at fixed magnetic flux Φ, computed in the case β_{0} = −6. The red curve is the unmagnetised solution. From bottom to top: other curves are computed at Φ = [2.55, 2.06, 1.46, 0.91]×10^{30} G cm^{2}. The various parts are: gravitationally stable descalarised 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 descalarised configuration with M_{0} = 1.68 M_{⊙} and Φ = 2.06 × 10^{30} G cm^{2}, while the arrow points to the blue dot where the configuration is expected to jump because of magneticallyinduced spontaneous scalarisation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.