Numerical simulations of stellar collapse in scalar-tensor theories of gravity

We present numerical-relativity simulations of spherically symmetric core collapse and compact-object formation in scalar-tensor theories of gravity. The additional scalar degree of freedom introduces a propagating monopole gravitational-wave mode. Detection of monopole scalar waves with current and future gravitational-wave experiments may constitute smoking gun evidence for strong-field modifications of general relativity. We collapse both polytropic and more realistic pre-supernova profiles using a high-resolution shock-capturing scheme and an approximate prescription for the nuclear equation of state. The most promising sources of scalar radiation are protoneutron stars collapsing to black holes. In case of a galactic core collapse event forming a black hole, Advanced LIGO may be able to place independent constraints on the parameters of the theory at a level comparable to current solar-system and binary-pulsar measurements. In the region of the parameter space admitting spontaneously scalarised stars, transition to configurations with prominent scalar hair before black-hole formation further enhances the emitted signal. Although a more realistic treatment of the microphysics is necessary to fully investigate the occurrence of spontaneous scalarisation of neutron star remnants, we speculate that formation of such objects could constrain the parameters of the theory beyond the current bounds obtained with solar-system and binary-pulsar experiments.


Introduction
General Relativity (GR) is currently assumed to be the standard theory of gravity, and has so far passed all experimental tests with flying colours [1][2][3][4].Theoretical and observational evidence, however, suggests that some modifications of GR may be inevitable.Cosmological and astrophysical observations require most of the energy content of the Universe to be present in the form of dark energy and dark matter [5][6][7].On theoretical grounds, GR is expected to represent the low-energy limit of a more fundamental (quantum) theory [8].Presently considered candidates for such theories predict modifications of GR at higher energies which also provide means to circumvent the formation of mathematical singularities inevitable in GR [9].
Attempts to generalise GR in these directions often involve additional fields that mediate the gravitational interaction together with the spacetime metric.The simplest class of such models is that of scalar-tensor (ST) theories, where one scalar field is included in the gravitational sector of the action.Ever since the pioneering work of Jordan, Fierz, Brans, and Dicke [10][11][12], ST theories have received a great deal of attention, both from a theoretical and a phenomenological point of view (see e.g.[13][14][15][16][17] and references therein).This class of theories is simple enough to allow for detailed predictions to be worked out, but also complicated enough to introduce a richer phenomenology leading to potentially observable deviations from GR. ST theories make predictions in the largely untested strong-field regime, while remaining compatible with the weak-field constraints imposed on GR by Solar System experiments (cf.Sec.3.2 below).
Black hole (BH) spacetimes might at first glance appear to represent an ideal testing ground for strong-field effects.The classical no-hair theorems, first proven for Brans-Dicke theory [18][19][20] and later extended to a wider range of ST theories (see [21,22] for reviews), however, strongly constrain the potential for deviations of BH spacetimes in ST theory from their GR counterparts.At leading post-Newtonian (PN) order, for example, the dynamics of a BH binary in Brans-Dicke theory is indistinguishable from the GR case [23].Indeed, considering the ST field equations given below as Eqs.(2.5)-(2.8),one immediately sees that vacuum solutions of GR are also solutions to the ST equations with a constant scalar field.Non-trivial BH dynamics can still be obtained by relaxing some of the fundamental ingredients of the no-hair theorems as for example a non-vanishing potential term [24] or non-asymptotic flatness [25].Due to the additional coupling introduced by the energy momentum tensor in the ST equations, however, compact matter sources such as neutron stars (NSs) and collapsing protoneutron stars forming BHs appear to be more promising objects for exploring observational signatures of ST theories.
Guided by this expectation, we shall focus in this paper on the formation of compact objects through gravitational collapse.Gravitational collapse is the expected evolutionary endpoint of stars of zero-age main sequence (ZAMS) mass of 10M M 130M [26][27][28].After exhausting their available fuel, the star's central core (mostly made of iron group nuclei) collapses under the strength of gravity as it exceeds its effective Chandrasekhar mass [29].Collapse proceeds until mass densities become comparable to those of nuclear matter.Thereafter, the increasingly repulsive character of the nuclear interactions leads to core bounce, which results in an outgoing hydrodynamic shock.The outgoing shock soon stalls because of dissociation of nuclei and neutrino emission in the post-shock region, and must be revived to successfully drive a supernova (SN) explosion [29].The physical mechanism responsible for the shock revival is still a topic of active research (see e.g.[30] and references therein).Multi-dimensional fluid instabilities and neutrino interactions are generally believed to play a crucial role in driving most SN explosions with the possible exception of hyper-energetic ones [31,32].One single core-collapse SN provides photon luminosities comparable to those of an entire galaxy and outshines all stars in the Universe in neutrinos.If the explosion is successful, a NS is left behind.If the explosion fails or is very weak, continued accretion will push the central NS over its maximum mass of 2 − 3 M and lead to the formation of a BH.The details of BH formation depend on the structure of the progenitor star and on the nuclear equation of state (EOS) [26].
ST theories may play a crucial role in this picture of NS and BH formation.A peculiar non-linear effect called "spontaneous scalarisation" [33,34] -somewhat similar to spontaneous magnetisation in ferromagnets -represents a particularly strong form of non-trivial scalar-field dynamics leading to additional branches of stationary NS families (see also [35][36][37] for dynamical scalarisation in binary NS systems).Moreover, ST theories provide a new channel for emission of gravitational waves (GWs) in stellar collapse.Whereas in GR conservation of mass and momentum exclude monopole and dipole radiation, monopole waves are permitted in ST theories in the form of scalar radiation, the so-called breathing mode.Detection of this breathing mode generated by a galactic SN would constitute smoking-gun evidence for a deviation from GR in the strong-field regime.Such tests of GR represent a major scientific goal [2] of the new era of GW astronomy initiated with the recent breakthrough detection of GW150914 [38], and thus add to the enormous scientific potential of exploring the physics of stellar collapse with GWs (see Ref. [39] for a comprehensive review on the topic).
The impact of ST theories on the equilibrium structure of NSs has been extensively studied in the literature (see, e.g., [33,[40][41][42][43][44]). Surprisingly few studies, however, have been devoted to their formation processes.Following pioneering numerical relativity simulations in Brans-Dicke theory [45], early studies have been devoted to dust-fluid collapse [46][47][48][49].The collapse of NSs into BHs [50] and the transition between different static NS branches [51] was first addressed by Novak using pseudo-spectral methods.To the best of our knowledge, the only published simulations of NS formation in ST theories have been presented by Novak and Ibáñez in Ref. [52], who combined pseudo-spectral techniques and high-resolution shock-capturing to study core collapse.The only other study we are aware of is Ref. [53], which numerically models the collapse of spherically symmetric fluids with a Γ-law EOS in Brans-Dicke theory and finds the monopole radiation to dominate at frequencies near the GW detectors' maximum sensitivity regime f ∼ 100 Hz, independently of the Brans-Dicke coupling parameter.The systematic exploration of GW emission from core collapse in ST theories thus represents a largely uncharted area in SN research.The dawning age of observational GW physics makes the filling of this gap a timely task, the first step of which is the main goal of this paper.
For this purpose, we have extended the open-source code gr1d of O'Connor and Ott [54] to ST theory and performed numerical simulations of NS and BH formation following core collapse to address the detectability of the monopole GWs with Advanced LIGO [55,56] and the proposed Einstein Telescope [57].We tackle the following questions.
• Are non-trivial scalar-field profiles and correspondingly large amplitudes in the scalar radiation naturally triggered in compact remnants following stellar collapse?• Can future GW observations of core collapse provide smoking gun evidence for deviations from GR in the framework of ST theories?
This paper is organised as follows.The action and the evolution equations of the theory are presented in Sec. 2. Additional physical ingredients entering our simulations are given in Sec. 3. Our numerical procedure is described in Sec. 4. We present our results on core collapse dynamics and monopole GW emission in Sec. 5. We summarise our findings in Sec. 6. Supporting material is provided online at Ref. [58].Throughout the paper, we generally use geometrical units c = G = 1, but occasionally restore factors of G for clarity of presentation.

Evolution equations
In this Section, we first review different ways to formulate ST theories and then arrive at the equations for the metric, scalar field, and matter sector in general covariant form (Sec. 2.1).Next, we derive the hydrodynamic equations for the matter sources, the metric and scalar field for the specific case of radial-gauge, polar-slicing coordinates (Sec.2.2).

A tale of two formulations
In ST theories, gravity is mediated by the spacetime metric g µν and an additional scalar field φ.The most general action which (i) involves a single scalar field coupled non-minimally to the metric, (ii) is invariant under space-time diffeomorphisms, (iii) is at most quadratic in derivatives of the field, and (iv) satisfies the weak equivalence principle can be written in the form [1,13,59] Here, d 4 x is the standard coordinate volume element, R is the Ricci scalar built from g µν , g = det g µν and the symbol ψ m collectively denotes all non-gravitational fields.
The theory has only two free functions of the scalar field: the potential V = V (φ) and the coupling function F = F (φ) * .If the potential V is a slowly varying function of φ -as expected on cosmological grounds, see [60] -it causes negligible effects on the propagation of φ on stellar scales.For simplicity, we thus set V = 0 throughout this paper; GR is then recovered for F = 1.Details on the choice of the coupling function F are postponed to Sec. 3.2.The weak equivalence principle -which has been verified experimentally to very high precision [3] -is guaranteed to hold as long as the matter part of the action S m does not couple to the scalar field, and its motion is therefore governed by the geodesics of the metric g µν .In this formulation, the scalar field does not interact with ordinary matter directly, but influences the motion of particles exclusively through its coupling with the spacetime metric.
The theory described by the action (2.1) is said to be formulated in the Jordan frame [10].Probably the most famous case of a ST theory, though by now severely constrained by solar-system tests [61], is Brans-Dicke theory [12]: the specific theory obtained by setting F (φ) = 2πφ 2 /ω BD where ω BD is constant [59].
Alternatively to the above Jordan-frame description, ST theories can also be formulated in the so-called Einstein frame.Here, one considers the conformal transformation and the action of Eq. (2.1) becomes The Ricci scalar R is now built from the Einstein metric ḡµν and ϕ is a redefinition of the scalar field φ through [59,62], The key advantage resulting from this conformal transformation is a minimal coupling between the conformal metric and the scalar field, evident at the level of the action.
The fact that such a redefinition of the theory exists has an important consequence for attempts to constrain the theory's parameters through observations of compact objects: BHs are less suitable to obtain such constraints because the action (2.3) in vacuum (S m = 0) reduces to the Einstein-Hilbert action of GR with a minimally coupled scalar field.In the action as well as the field equations given further below, it is evident that matter sources represent an additional and more straightforward channel to couple the metric and scalar sectors.
The equations of motion in the Jordan frame can be obtained by varying the action (2.1) with respect to the spacetime metric g µν and the scalar field φ: Combining the Bianchi identities with the field equations can be shown to imply that the matter part of the energy momentum tensor, is conserved on its own, i.e. ∇ µ T µν = 0 . (2.10) This feature makes the Jordan frame particularly suitable for studying stellar collapse: the matter equations, which are expected to develop shocks, do not need to be modified from their GR counterparts (cf.Sec.2.2.3).The drawback of this choice is that the scalar field is not minimally coupled to the metric, i.e. the Hilbert term in the action (2.1) acquires a φ-dependent factor.This factor F (φ) leads to the term T F µν on the right-hand side of Eq. (2.5) additionally to the minimally coupling term T φ µν and the standard matter sources T µν .

Equation of motions
We now restrict the equations of motion to spherical symmetry in radial-gauge, polarslicing coordinates [63].The line element in the Jordan frame is where the metric functions α = α(t, r) and X = X(t, r) can be more conveniently rewritten in terms of the metric potential, and the enclosed mass, (2.13) Note that in Eq. (2.11) we multiplied the angular part of the metric dΩ 2 by a factor 1/F , thus effectively imposing the radial gauge in the Einstein frame.In this formulation, the (Jordan-frame) areal radius is given by r/ √ F .This choice allows for comparisons with Refs.[50][51][52], where the analysis is entirely carried out in the Einstein frame.Likewise, Φ and m are Einstein-frame variables and their definition in terms of the Jordan metric components in Eqs.(2.12), (2.13) acquires factors of F .
Following [54], we assume ideal hydrodynamics as described by an energymomentum tensor of the form and the matter current density Here ρ is the baryonic density, P is the fluid pressure, h is the specific enthalpy (which is related to the specific internal energy and the pressure P by h = 1 + + P/ρ), and u µ is the 4-velocity of the fluid.Spherical symmetry implies where v = v(t, r).
The equations of motion can be reformulated in flux conservative form using conserved variables and thus become amenable to a numerical treatment using highresolution shock-capturing schemes [64,65].These conserved variables D, S r and τ are related to the to the primitive variables ρ, , v and P by The definitions above generalise Eq. ( 8) in Ref. [54] to ST theory.We take advantage of the Einstein-frame scalar-field redefinition φ → ϕ of Eq. (2.4) because it simplifies the wave equation (2.8).Moreover, the space of ST theories and the weak-field experimental constraints are traditionally described in terms of ϕ (cf.Sec.3.2).Following Refs.[50][51][52], we introduce auxiliary variables for the derivatives of the scalar field defined by (2.20) 2.2.1.Metric equations.The evolution equations (2.5-2.7) for the metric potential Φ and the mass function m expressed in terms of the conserved variables read These equations are not independent; the last equation for ∂ t m directly follows from the other two combined with the conservation of the energy momentum tensor (2.10).
For convenience, we follow standard practice and compute the metric functions using the constraints (2.21), (2.22) and discard the time evolution equation for m.From Eq. (2.21), we further notice that the metric potential Φ is determined only up to an additive constant.In GR, this freedom is commonly used to match the outer edge of the computational domain to an external Schwarzschild metric.This cannot be done in ST theories, as such theories do not obey a direct analogue of the Birkhoff theorem [66,67].We therefore specify a boundary condition for Φ using the method put forward by Novak [50]: Φ is constrained on the outer boundary of the computational domain by requiring that is approximately constant in the weak-field regime, far away from the star.K is then evaluated for the initial profile (cf.Sec.3.2), fixed to be constant during the evolution and determines Φ on the outer edge of the grid r = R out by inverting (2.24) Note that the Birkhoff theorem in GR corresponds to the case K = 1.The error incurred from this procedure can be estimated by comparing results obtained for different extents of the computational domain.We obtain variations of order |∆ϕ/ϕ| ∼ 10 −3 at the radius of extraction when the grid extent is decreased by a factor 2 (cf.Sec.4.1 for more details on our numerical setup).Similar errors are detected in the collapse of a ST polytrope if K is set to 1, rather than evaluated from the initial profile.
2.2.2.Scalar-field equations.The wave equation for the scalar field (2.8) can be written as a first-order system using the definitions (2.20) and the identity In order to prescribe the behaviour of ϕ at the outer boundary, we consider the asymptotic behaviour of the scalar field at spatial infinity [13] where ϕ 0 = const and and ω denotes the scalar charge of the star.Physically, we require that no radiation enters the spacetime from infinity and therefore impose an outgoing boundary condition [68] at spatial infinity where f is a free function of retarded time.This condition can be translated into the following differential expressions for η and ψ, and the scalar field ϕ is directly obtained from Eq. (2.26).As shown in more detail below [see Eq. (3.8) and the following discussion], the value of ϕ 0 is degenerate with one of the parameters used to describe the coupling function and is set to zero in our study without loss of generality.In practice, our computational domain extends to large but finite radii and we approximate the physical boundary conditions by imposing Eqs.(2.31), (2.32) at the outer edge rather than at infinity.As already mentioned, we have tested the influence of the outer boundary location on our results and observe only tiny variations of the order of |∆ϕ/ϕ| ∼ 10 −3 in the extraction region when comparing with simulations performed with R out twice as large.

2.2.3.
Matter equations in flux-conservative form.The evolution equations (2.5-2.8)can be conveniently written in flux-conservative form [65,69], where U is the vector of the conserved variables U = [D, S r , τ ] defined in Eqs.(2.17 ) ) Note that Ref. [52] misses a factor 1/a (in their notation) inside the argument of the radial derivative in their Eq.(11).Inclusion of this factor and pulling the term proportional to ηv in Eq. (2.37) out of the radial derivative enables us to cast the evolution equation for D in the same form (2.33) as the other matter equations.For the integration of the evolution equation for D, we therefore do not need the additional considerations described in Sec.2.1 of [52].
The hyperbolic structure of the system of equations (2.33) is dictated by the Jacobian matrix of the fluxes [70], (2.40) The characteristic speeds associated with the propagation of the matter fields are the eigenvalues λ of J U , Here c s = (dP/dρ) S /h (S is the entropy) is the local speed of sound given for our choice of EOS of the form P = P (ρ, ) by The characteristic speeds are therefore exactly the same as in GR, since they do not depend on the conformal factor F .The high-resolution shock-capturing scheme implemented in gr1d for GR [54] can therefore be used in ST theories as well, provided the conserved variables U and their fluxes f(U) are generalised using the expressions presented above.

Physical setup
In this Section, we discuss in more detail the physical ingredients entering our simulations.We discuss the EOS for the fluid used in our work (Sec.3.1), the various choices for the coupling function that relate the physical metric to its conformally rescaled counterpart (Sec.3.2) and the initial stellar profiles used in our study (Sec.3.3).We also provide information on the quantities used to compare GW signals and detector sensitivities in the context of monopole waves (Sec.3.4).

Equation of state
An EOS is required to close the hydrodynamical system of equations.Specifically, it provides a prescription for the pressure P and other thermodynamic quantities as a function of the mass density, internal energy (or temperature), and possibly the chemical composition.In this paper we study stellar collapse using the so-called hybrid EOS.This EOS was introduced in Ref. [71] and qualitatively captures in closed analytic form the expected stiffening of the nuclear matter EOS at nuclear density and includes nonisentropic (thermal) effects to model the response of shocked material.
The hybrid EOS was widely used in early multi-dimensional core-collapse simulations (e.g.[72,73]), and the results of simulations using a hybrid EOS have been compared in detail with those obtained with modern finite-temperature EOS; see e.g.Ref. [74,75].
The hybrid EOS consists of a cold and a thermal part: The cold component P c is modelled in piecewise polytropic form with adiabatic indices Γ 1 and Γ 2 , This expression models both the pressure contribution from relativistic electrons, which dominates at ρ ≤ ρ nuc , and the stiffening at nuclear density due to the repulsive character of the nuclear force.The two components are matched at "nuclear density" which we set to ρ nuc = 2 × 10 14 g/cm 3 following [73].We set K 1 = 4.9345 × 10 14 [cgs], as predicted for a relativistic degenerate gas of electrons with electron fraction Y e = 0.5 [76], while is then obtained from requiring continuity in P at ρ = ρ nuc .The specific internal energy follows from the first law of thermodynamics applied to the case of adiabatic processes where the integration constant E 3 is determined by continuity at ρ = ρ nuc .The thermal contribution P th is described by a Γ-law with adiabatic index Γ th , where th = − c is the thermal contribution to the internal energy, computed from the primitive variable .The flow is adiabatic before bounce, implying that c and the total pressure is described by considering only its cold contribution.At core bounce, however, the hydrodynamic shock results in nonadiabatic flow and thus triggers the onset of a non-negligible thermal contribution to the EOS.
We consider a hybrid EOS characterised by three parameters: Γ 1 , Γ 2 and Γ th .The physical range of these adiabatic indices has been explicitly studied in Refs.[74,75], where 2+1 GR simulations of core collapse were used to compute the effective adiabatic index of the finite-temperature EOS of Lattimer and Swesty [77,78] and Shen et al. [79,80].In the collapse phase, electron capture decreases the effective adiabatic index below the value Γ 1 = 4/3 predicted for a relativistic gas of electrons.More precisely, comparisons with more detailed simulations yields a range from Γ 1 1.32 to Γ 1 1.28 [74,75,81].In particular, lower values of Γ 1 are found when deleptonisation is taken into account because electron capture onto nuclei before neutrino trapping decreases Y e for given ρ, thus softening the EOS.Collapse is stopped by the stiffening of the EOS at nuclear density which raises the effective adiabatic index Γ 2 above 4/3.Reference [75] finds Γ 2 3.0 for the Shen et al.EOS and Γ 2 2.5 for the Lattimer-Swesty EOS.Finally, the thermal adiabatic index Γ th models a mixture of relativistic and non-relativistic gas, and is therefore physically bounded to 4/3 < Γ th < 5/3.We select fiducial values Γ 1 = 1.3, Γ 2 = 2.5, Γ th = 1.35 for our code tests presented in Sec. 4, and explore a more extended parameter grid around this model in Sec. 5.

Coupling function
As introduced in Sec.2.1, ST theories with a single scalar field and vanishing potential are described by a single free function F (ϕ).The phenomenology of ST theories is simplified, however, by the fact that all modifications of gravity at first PN order depend only on two parameters.These are the asymptotic values of the first and second derivatives of ln F [13, 14, 34] * , The effective gravitational constant determining the attraction between two bodies as measured in a Cavendish experiment is where G is the bare gravitational constant entering the action.Furthermore, the Eddington Parameterised post-Newtonian parameters [82,83] can be expressed exclusively in terms of α 0 and β 0 through For an interpretation of these equations in terms of fundamental interactions, see Ref. [84].In consequence, weak-field deviations from GR are completely determined by the Taylor expansion of ln F to quadratic order about lim r→∞ ϕ = ϕ 0 .For these reasons, most of the literature on ST theories has focused on coupling functions of quadratic form [33,34] and we follow this approach by employing a coupling function The asymptotic value ϕ 0 does not represent an additional degree of freedom in the theory because it can be reabsorbed by a field redefinition ϕ → ϕ + ϕ 0 [60] and we therefore set ϕ 0 = 0 without loss of generality in the rest of this paper.* We can furthermore assume α 0 ≥ 0 because the sign of α 0 is degenerate with the field redefinition ϕ → −ϕ.Despite its apparent simplicity, this two-parameter family of ST theories is representative of all ST theories with the same phenomenology up to first PN order.Brans-Dicke theory [12] is a special case of Eq. (3.8) with the Brans-Dicke parameter [defined above Eq.(2. 2)] given by ω BD = (1 − 6α 2 0 )/2α 2 0 and β 0 = 0.It is worth mentioning here that theories with the coupling function (3.8) and strictly vanishing potential have been shown to exhibit non-viable cosmological evolutions [85,86]; however, this can be cured by introducing a suitable (sufficiently flat) potential which leaves the phenomenology on stellar scales unchanged [87,88].
It is well known that all deviations in the structure of spherically symmetric bodies in ST theory from their general relativistic counterparts are given in terms of a series of PN terms proportional to α 2 0 [13,89]; cf. also Eq. (3.7) above.Any ST theory with α 0 = 0 is therefore perturbatively equivalent to GR and current observations (see below) constrain α 0 to be very small.In 1993, however, Damour and Esposito-Farèse [33,34] discovered a remarkable non-perturbative effect called spontaneous scalarisation, which introduces macroscopic modifications to the structure of NSs even when α 0 is very small or vanishes [90].For certain values β 0 < 0, there exists a threshold in the compactness (M/R, where M is the total mass of the object and R is its radius) of stellar structure at which spherically symmetric equilibrium solutions develop significant scalar hair.One can find three distinct solutions in this regime: besides a weakly scalarised solution where the ratio between the scalar charge and the star's mass ω/M is of the order of α 0 , two strong-field solutions appear where this ratio is of order unity [51,91].If α 0 = 0, the weak-field solution is a GR star and the two strong field solutions coincide.Notably, scalarised solutions are present for compactness values of order M/R 0.2 [89], as realised in NSs.When present, scalarised neutron stars can be energetically favoured over their weak-field counterparts [33,34,92], allowing for the possibility of dynamical transitions between the two branches of solutions [51].Spontaneously scalarised stars have been found for β 0 −4.35 [51,91], but the exact value of this threshold depends on the EOS.
The (α 0 , β 0 ) parameter space of ST theories has been severely constrained by observations.Solar System probes include measurements of Mercury's perihelion shift [93], Lunar Laser Ranging [94], light deflection measured with Very-Long-Baseline Interferometry [95], and the impressive bound α 0 < 3.4 × 10 −3 obtained with the Cassini space mission [61].Timing of binary pulsars currently provides the tightest constraints in the β 0 direction of the parameter space [96].In particular, observations from pulsars PSR J1738+0333 [97] and PSR J0348+0432 [98] (both orbiting a white dwarf companion) rule out a wide range of theories exhibiting prominent spontaneous scalarisation.Current observational constraints are summarised in Fig. 3.1 where the shaded area is now excluded.Note, however, that the binary-pulsar constraints apply to the case of a single massless scalar field.Scalar-tensor theories with multiple scalar fields [43] or with one massive field [99] may still lead to spontaneously scalarised neutron stars over a wide range of the theories' parameters without coming into conflict with the binary pulsar observations.

Initial profiles
We perform simulations of stellar collapse starting from two types of initial data: (i) polytropic models generated in the static limit of the ST theory equations and (ii) "realistic" SN progenitors obtained from stellar evolutionary computations performed by Woosley & Heger [100].
(i) In the static limit, the evolution equations presented in Sec.2.2.1-2.2.3 reduce to (cf.[33,43]) ∂ r ϕ = Xη , (3.12) which generalise the Tolman-Oppenheimer-Volkoff [101,102] equations to ST theory.As in GR, the equation for the metric potential Φ decouples from the remainder and we need an EOS P = P EOS (ρ) to close the system.
In practice, we integrate the system (3.9)-(3.13)outwards starting at the origin where boundary and regularity conditions require Here, P c (or, alternatively, ρ c ) is a free parameter determining the overall mass and size of the star and the central value of the scalar field ϕ c is related through the integration to the value of ϕ at infinity.In our case, the boundary condition for the scalar field is ϕ(r → ∞) = ϕ 0 = 0 and the task is to identify the "correct" central value ϕ c that satisfies the outer boundary condition.From a numerical point of view, this task represents a two point boundary value problem [103] and we use a shooting algorithm to solve it.For this purpose, we note that the integration terminates at the stellar surface r s defined as the innermost radius where P = 0. From this radius r s , we could in principle continue the integration to infinity by setting the matter sources to zero and switching to a compactified radial coordinate such as y ≡ 1/r.We have found such a scheme to work successfully [43], but here we implement an equivalent, but conceptually simpler algorithm.The numerical solution computed for r ≤ r s can be matched to a vacuum solution at r > r s to relate the scalar field at the stellar surface ϕ s to its asymptotic value ϕ 0 at r = ∞ [33]: where the subscript s denotes quantities evaluated at r s .The shooting algorithm starts the integration of Eqs.(3.9-3.13) with some initial guess ϕ(0), obtains the corresponding ϕ s and then iteratively improves the choice of ϕ(0) until it leads to a ϕ s that satisfies Eq. (3.15) within some numerical tolerance (10 −10 for the absolute difference in our case).
The central density or pressure can be freely chosen and parameterises the family of static solutions for a given ST theory (α 0 , β 0 ) in the same way as it does in GR.
The members of this one-parameter family of solutions are often characterised by their total gravitational mass which is given by [33] (3.16)All polytropic initial profiles used in this work are generated using a polytropic EOS P = Kρ Γ with K = 4.9345 × 10 14 [cgs], Γ = 4/3 and central mass density ρ c = 10 10 g cm −3 ; these parameters are considered qualitatively reasonable approximations to model iron cores supported by the degeneracy pressure of relativistic electrons [76].In particular, the choice ρ c = 10 10 g cm −3 results in stars with baryonic mass ∼ 1.44M , slightly below the Chandrasekhar limit [29].(ii) We also perform core collapse simulations using more realistic pre-SN models.
Woosley and Heger [100] evolved non-rotating single stars up to the point of iron core collapse [104,105].Here, we consider two specific models of their catalogue obtained from the evolution of stars with zero-age-main-sequence (ZAMS) mass M ZAMS = 12M and 40M .We refer to these models as WH12 and WH40 respectively.Model WH12 has a steep density gradient outside its iron core, and M ZAMS = 40M (WH40) pre-SN models of Woosley and Heger [100] while the solid lines show three Γ = 4/3 polytropes generated in ST theories with (α 0 , β 0 ) = (0.001, −3), (0.003, 0), (0.01, −5).The mass-density distributions of all three ST polytropes are indistinguishable from their GR counterparts.The more realistic models WH12 and WH40 mostly differ from the polytropic ones through the presence of outer low-density layers.Note that we cut the WH models at rs = 2 × 10 4 km and pad them with an artificial atmosphere of ρatm = 1 g cm −3 .The scalar field is more pronounced in models with higher α 0 , but the low compactness of these models prevents spontaneous scalarisation.The scalar field is initialised to zero when the WH models are evolved.
which results in a low accretion rate after bounce.Even if no explosion occurs, the delay to BH formation would be multiple seconds and no BH forms over the time we simulate.Model WH40, on the other hand, has a very shallow density gradient, resulting in a high accretion rate after bounce.This pushes the protoneutron star over its maximum mass and leads to BH formation within a few hundred milliseconds of bounce (cf.[26]).Hence, we use model WH12 to explore ST theory for a scenario in which core collapse results in a stable NS and model WH40 for a scenario in which the protoneutron star collapses to a BH.Since WH12 and WH40 are Newtonian models, we initialise the scalar-field variables ϕ, ψ and η to 0. An unfortunate consequence of this artificial (but unavoidable) approximation is that no scalar-field dynamics occur at all if α 0 = 0: all source terms on the right-hand side of Eqs.(2.26-2.28)vanish at all times, and the evolution proceeds exactly as in GR.We overcome this problem, by using small but non-zero values for α 0 , which triggers a brief initial transient in the scalar field that afterwards settles down into a smooth but non-trivial configuration eventually leading to significant scalar field dynamics as the collapse progresses through increasingly compact stages of the core.
For both classes of initial data there remains one degree of freedom that we need to specify: the metric function Φ is determined by Eq. (3.9)only up to an additive constant.While our integration in case (i) starts with Φ(0) = 0, we can trivially shift the profile of Φ(r) by a constant (leaving all other variables unchanged) and still have a solution of the system of Eqs.(3.9)-(3.13).We use this freedom to enforce that the physical metric component g tt = 1 as r → ∞, so that coordinate time is identical to the proper time measured by an observer at infinity.In practice, this is achieved by using very large grids and fitting Φ = Φ 0 + Φ 1 /r on the outer parts.Φ 0 is then the constant we subtract from the entire profile Φ(r).The realistic initial models of case (ii) above are calculated without a scalar field and in that case our procedure is equivalent to the standard matching in GR based on the Birkhoff theorem.
For illustration, we show in Fig. 3.2 some of the initial profiles used in this study.Because of the low compactness of iron cores, the polytropic profiles for all values of α 0 ≤ 0.01 present very similar mass-density distributions which also very closely resemble their GR counterpart.The magnitude of the scalar field inside the star increases as larger values are chosen for α 0 (cf.right panel of Fig. 3.2) while outside the star ϕ rapidly approaches the 1/r behaviour of Eq. (2.29).In the left panel of Fig. 3.2, we also see that the realistic pre-SN models WH12 and WH40 are well approximated by a Γ = 4/3 polytrope in their central regions r 10 3 km; outer less degenerate layers of lighter elements, however, substantially broaden the mass-density distribution outside the iron core.
In order to overcome instabilities arising in our numerical scheme due to zero densities ρ [54], we add an artificial atmosphere outside the stellar surface r s .More specifically, we pad the polytropic profiles with a layer of constant mass density ρ atm = 1 g cm −3 .The WH models are cut at r s = 2 × 10 4 km (cf.Fig. 3.2) and padded with an artificial atmosphere of ρ atm = 1 g cm −3 .By comparing evolutions using different values for the atmospheric density, we find the atmosphere to be completely irrelevant to the dynamics of the star, but we observe that significantly larger values than ρ atm = 1 g cm −3 unphysically affect the propagation of the scalar wave signal such that the wave signal does not converge in the limit of large extraction radii.We estimate the resulting error for our choice by comparison with otherwise identical simulations using instead ρ atm = 10 g cm −3 ; the observed differences are |∆h(t)|/h(t) ∼ 0.3% in the extracted waveform [cf.Eq. (3.23) below].

GW extraction and detector sensitivity curves
The output of a GW detector s(t) = n(t) + h(t) is the sum of noise n(t) and signal h(t).For quadrupole GWs, as present in GR, h(t) is related to the metric perturbation h +,× in the transverse traceless gauge through the beam pattern functions A +,× : [106].Monopole GWs are present in ST theory and are related to the dynamics of the scalar field ϕ.In this case, the detector response h(t) = A • h • (t) is given by the metric perturbation h • (t) weighted by the correspondent beam pattern A • [107,108] * .If we denote by h(f ) and ñ(f ) the Fourier transform of h(t) and n(t), respectively, the (one-sided) noise power spectral density S n (f ) is defined as where • denotes a time average for stationary stochastic noise.The signal-to-noise ratio is defined as (see [109] where the numerical factor is derived; see also [110]) The characteristic strains for noise and signal are defined as such that ρ 2 can be written as the squared ratio between signal and noise: The most common convention used to display detector sensitivity curves involves plotting the square root of the power-spectral density and the analogous quantity [109] which characterises the GW signal.* In the following we will use sensitivity curves S n (f ) for: (i) the Advanced LIGO detectors [55,56] in their zero-detuned high-power configuration, as anticipated in [111]; (ii) the proposed Einstein Telescope [57], using the analytic fit reported in [110].
Scalar waves are also promising sources for future GW experiments targeting the deci-Hertz regime, such as the proposed space mission DECIGO [112].
In contrast to GR, ST theories admit gravitational radiation in spherical symmetry, specifically in the form of a radiative monopole of the scalar field or, equivalently, a so-called breathing mode when considering the Jordan frame.The metric perturbation of a monopole scalar wave in ST theory is [13] h where D is the distance between the detector and the source and the scalar field ϕ is evaluated at radius r.The factor α 0 is due to the coupling between the scalar field and the detector and limits the potential of GW observations to constrain ST theories [113].Throughout this paper we consider optimally oriented sources, such that In analysing our simulations, we proceed as follows.At a given radius r ext , we extract ϕ(t) and compute h(t).In order to eliminate the brief unphysical transient (cf.Sec.4.1), we truncate this early part from the time domain waveform h(t).We then obtain h(f ) numerically using a Fast Fourier Transform (FFT) algorithm.To reduce spectral leakage, the FFT algorithm is applied to data h(t) mirrored about the latest timestep available and the result h(f ) is normalised accordingly [114].This confines spectral leakage to frequencies f 200 Hz (cf. the tails in Fig. 5.4 and 5.5) where the signal is very weak.Finally, we compute S h (f ) from Eq.(3.22) and compare it with the detectors' sensitivity curves S n (f ).

Numerical implementation
In this Section, we provide details of our numerical scheme, stressing the modifications needed in ST theories with respect to the GR version of the code (Sec 4.1).We present the convergence properties in Sec.4.2.

Second-order finite differences and high-resolution shock capturing
Our numerical code is built on top of gr1d, an open-source spherically-symmetric Fortran 90+ code developed by O'Connor and Ott [54].gr1d has been applied to a range of problems in stellar collapse and BH formation (e.g., [26,115]).Its most recent GR version is available at [116] and includes energy-dependent neutrino radiation transport [117].
As in the GR case, the constraint equations (2.21) and (2.22) for the metric functions Φ and m are integrated using standard second-order quadrature.In the scalar-field equations (2.26-2.28), the source terms are discretised using centred second-order stencils.Due to the potential formation of shocks in the matter variables, their evolution is handled with a high-resolution shock capturing scheme as described in detail in Sec.2.1 of [54].For our evolutions in ST theory, we extended the flux and source terms of gr1d in accordance with our Eqs.(2.34)-(2.39).Integration in time of the evolution equations for the matter and scalar fields is performed using the method of lines with a second-order Runge-Kutta algorithm.One significant difference from the GR case arises from the presence of the scalar field as a dynamical degree of freedom with the characteristic speed of light, whereas in spherical symmetry in GR we only have to consider the characteristic speed of sound for the matter degrees of freedom (cf.Sec.2.2.3).In order to satisfy the Courant-Friedrichs-Lewy stability condition we therefore determine the timestep using the speed of light instead of the speed of sound, which results in smaller values for the allowed timestep as compared with the corresponding evolutions in GR.
As discussed in Sec.2.2.3, a key ingredient in the implementation of shock capturing methods is the use of a system of evolution equations in flux conservative form which is available in terms of the conserved variables D, S r , τ but not in the primitive variables ρ, v and .The primitive variables appear in the constraint equations for the metric, the flux terms of the shock-capturing scheme, in the EOS, and also form convenient diagnostic output.Conversion between the two sets of variables is thus required at each timestep.This process is straightforward for the direction primitive → conserved; cf.Eqs.(2.17-2.19).The reverse transformation, however, is non-trivial because of the presence of the pressure P which is an intrinsic function of ρ and given by the EOS.This conversion is performed iteratively using a Newton-Raphson algorithm: given an initial guess P for the pressure from the previous timestep, we first calculate in this order Then we compute an updated estimate for the pressure from the EOS P = P (ρ, ), and iterate this procedure until convergence.The evolution of the scalar field turns out to be susceptible to numerical noise near the origin r = 0, which represents a coordinate singularity.In order to obtain longterm stable evolutions, we add artificial dissipation terms of Berger-Oliger type [118] to the scalar evolution equations.Specifically, we add a dissipation term of the form D × ∆r 4 × ∂ 4 u/∂r 4 to the right-hand-side of Eqs.(2.26)-(2.28),where u stands for either of the scalar-field variables, ∆r is the width of the grid cell, and D is a dissipation coefficient.In practice, we obtain good results using D = 2.
In all our simulations, the grid functions exhibit much stronger spatial variation in the central region of the star than in the wave zone.In order to accommodate these space dependent requirements on the resolution of our computational domain, we use a numerical domain composed of an inner grid with constant and an outer grid with logarithmic spacing.This setup enables us to capture the dynamics of the inner core with high accuracy while maintaining a large grid for GW extraction at tolerable computational cost.Unless specified otherwise, we use the following grid setup.The outer edge of the grid is placed at R out = 1.8 × 10 5 km and the two grid components are matched at R match = 40 km.The cell width of the inner grid is ∆r = 0.25 km.The total number of zones is set to N zones = 5000, such that 160 (4840) zones are present in the inner (outer) grid.4 ghost cells are added at both r = 0 and r = R out for implementing symmetry and boundary conditions.GW signals are extracted at r ext = 3 × 10 4 km, which is well outside the surface of the star but sufficiently far from the outer edge of the grid R out to avoid contamination from numerical noise from the outer boundary.We simulate the evolution for 0.7 s to allow for the entire GW signal to cross the extraction region.
Radial gauge, polar slicing coordinates are not well adapted to BH spacetimes: as the star approaches BH formation, the lapse function α tends to zero in the inner region [119] and inevitably introduces significant numerical noise.The stellar evolution, however, is effectively frozen as α → 0. Following Novak [50], we handle BH formation by explicitly stopping the evolution of the matter variables while we let the scalar field propagate outwards.In practice, we freeze the matter evolution whenever the central value of α becomes smaller than α T = 5 × 10 −3 .By varying the threshold α T over two orders of magnitude, we verified this procedure introduces a negligible error |∆ϕ|/ϕ 1% on the extracted wave signal in case of BH formation.
A final note on the numerical methods concerns the time window used for the wave extraction.As mentioned in Sec.3.2, our initial data for the realistic progenitor models require us to trigger scalar dynamics by using a small but non-zero value for α 0 that induces a brief transient in the wave signal.This transient is not part of the physical signal we are interested in and is removed by calculating waveforms in an interval starting not at zero retarded time, but shortly afterwards: we use for this purpose the time window [t i , t f ] with t i = r ext /c + 0.006 s to t f = r ext /c + 0.6 s from the beginning of the simulation.This provides us with waveforms of total length ∆t ∼ 0.6 s corresponding to a lower bound f ∼ 1.7 Hz in the frequency domain.Note that our waveforms are significantly longer than those obtained in previous studies of collapse in ST theories [52].A polytropic core is collapsed using the fiducial hybrid EOS for three different resolutions (see text for details).The top panels show the evolution of the mass density ρ (left) and the scalar-field ϕ (right) for the highest-resolution run a t = 7 (grey), 28 (yellow), 37 (green), 38 (blue), 57 (red) ms after starting the simulations.Bounce happens at t ∼ 38 ms and the shock reaches the surface of the star at t ∼ 131 ms.The bottom panels show the self-convergence properties of the gravitational mass m (left) and the scalar field ϕ (right) at the same times.As detailed in the text, solid and dotted (dashed) lines are expected to coincide for second-(first-) order convergence.We initially observe second order convergence which decreases to first order due to (i) the shock capturing scheme when a discontinuity forms at bounce and (ii) numerical noise in the scalar field propagating in from the outer boundary.

Self-convergence test
Here we present the convergence properties of our dynamical code.Given three simulations of increasing resolutions with grid spacings ∆r 1 > ∆r 2 > ∆r 3 , the selfconvergence factor Q of a quantity q is defined by where q i indicates the quantity q obtained at resolution ∆r i and n is the convergence order of the implemented numerical scheme.We collapse a Γ = 4/3 polytropic core in ST theory with α 0 = 10 −4 and β 0 = −4.35using the hybrid EOS with Γ 1 = 1.3, Γ 2 = 2.5 and Γ th = 1.35.This model is evolved for three uniformly spaced * grids of size R out = 2 × 10 3 km with N = 6000, 12000 and 24000 grid cells.For these grids, we expect Q = 2 (Q = 4) for first-(second-) order convergence.The bottom panels of Fig. 4.1 show the convergence properties of the gravitational mass m and the scalar field ϕ at various timesteps.Solid lines show the difference between the coarse and the medium resolution runs q 1 − q 2 ; dashed (dotted) lines show the difference between the medium and the fine resolution runs q 2 − q 3 multiplied by the expected first-(second-) order self-convergence factor Q = 2 (Q = 4).Second-order convergence is achieved if the solid and dotted lines coincide, while the code is only first-order convergent if the solid and dashed lines coincide.The evolution of ρ and ϕ is displayed in the top panels for orientation.
The enclosed gravitational mass m shows good second-order convergence properties before bounce t 38 ms, while convergence deteriorates to first order as the shock propagates outwards at t 38 ms.This is a characteristic feature of high-resolution shock-capturing schemes; they are second-order (or higher) schemes for smooth fields, but become first-order accurate in the presence of discontinuities [54].Note that the behaviour of the total gravitational and baryonic mass is more complex than in the GR limit where both are conserved because of the absence of gravitational radiation in spherical symmetry and the vanishing of the source term s D in Eq. (2.37).The convergence properties of the scalar field are more complicated.While evolved with a second-order accurate scheme, we observe that the scalar field's convergence may deteriorate for the following two reasons: (i) the drop to first-order convergence of the matter fields which source the scalar dynamics; (ii) numerical noise generated at the outer boundary, especially during the early transient (note that in this convergence analysis the outer boundary is located much closer to the core than in our production runs because of the limit imposed by a strictly uniform grid).The observed convergence in the scalar field bears out these effects.Initially convergent at second order, we note a drop to roughly first order after one light crossing time R out /c ∼ 7 ms.As the noise is gradually dissipated away, the convergence increases back towards second order, but drops once more to first order at the time the shock forms in the matter profile around 38 ms.
We also tested the convergence of the scalar waveform ϕ(t) extracted at finite radius in these simulations and observe first-order convergence which we attribute to the relatively small total computational domain such that the outer boundary effects discussed above causally affect the extraction radius early in the simulation.The resulting uncertainty in the waveform is obtained by comparing the finite resolution result with the Richardson extrapolated (see, e.g., [120]) waveform.We find a relative error of 10% which we regard as a conservative estimate since the production runs have much larger computational domains.

Results and discussion
In this Section, we present the results of our simulations.After illustrating the main features of stellar collapse in ST theories (Sec.5.1), we present our predictions for monopole gravitational radiation (Sec.5.2).All waveforms presented in this section are publicly available at [58].
Since Γ 1 < 4/3, the initial iron cores are not equilibrium solutions of the evolution equations and collapse is triggered dynamically.While the polytropic profile collapses smoothly from the very beginning of the simulation, a brief transient in the scalarfield evolution is present in the early stages of the collapse of both the WH12 and WH40 models.As already mentioned in Sec.3.3, this is due to the fact that these initial models are Newtonian and their initially vanishing scalar profiles are not fully consistent with the ST theory used in the evolution.This transient generates a pulse in the scalar field propagating outwards at the speed of light.The scalar field quickly settles down in the stellar interior while the spurious pulse reaches the outer edge of the grid at R out /c ∼ 0.6 s where it is absorbed by the outgoing boundary condition.
As the collapse proceeds in either of the three models, the central mass density increases from its initial value up to beyond nuclear densities ρ nuc 2 × 10 14 g cm −3 .The EOS suddenly stiffens to an effective adiabatic index Γ 2 4/3 and the inner core bounces after t b ∼ 38, 39, and 84 ms from the beginning of the simulations, for the ST polytrope, WH12 and WH40 profile, respectively * .Core bounce launches a hydrodynamic shock into the outer core.Due to the steep density profile of the polytrope, the shock explodes the polytrope promptly, reaching its surface at ∼ 130 ms from the start of the simulation.Since we set Γ th = 1.35 to qualitatively account for reduced pressure due to nuclear dissociation and neutrino losses, the pressure behind the shock is not sufficient to lead to a prompt explosion of the more realistic WH12 and WH40 progenitors.The shock stalls and only secularly moves to larger radii as the accretion rate decreases.Core bounce is paralleled by a small reversal in the scalar field amplitude.For example, in the collapse of the polytropic model shown in Fig. 5.1, the central value of ϕ reaches a minimum ∼ −2.6 × 10 −5 at bounce before settling down to ∼ −2.3 × 10 −5 .A more detailed description of the scalar field dynamics is postponed to Sec. 5.2.* The WH40 profile takes longer to reach ρnuc because of its lower initial central density (cf.Fig. 3.2).The inner regions of the promptly exploding polytropic model settle down to a stable compact remnant with compactness m/r ∼ 0.053 (measured from the metric potential at r = 10 km).While simulations with model WH12 show that the shock stalls and then only slowly shifts to larger radii, the low accretion rate in this model does not increase its compactness above the values that we find for the polytropic model.In both models, the scalar charge ω evolves from ∼ −10 −4 M to ∼ −2 × 10 −4 M during the entire evolution and thus remains of the order of ω/M ∼ α 0 , as predicted for weakly scalarised NS solutions (cf.Sec.3.2).In both simulations, the NSs do not evolve to strongly spontaneously scalarised solutions because the compactness of the core remains lower than the threshold at which multiple solutions appear (m/r ∼ 0.2 [89]).
On the other hand, the WH40 model forms a protoneutron star that subsequently collapses to a BH within t BH ∼ 0.46 s from the beginning of the simulation (∼ 0.38 s from bounce).The high accretion rate in this model quickly increases the central compactness.As BH formation is approached, our gauge choice causes the lapse function α to collapse to zero near the origin (Fig. 5.2) and the dynamics of the inner region effectively freezes.In this regime, spontaneously scalarised NS solutions are not only present but energetically favourable [33,34,92].While collapsing towards a BH, the core first transits through a spontaneously scalarised NS.BH formation generates strong scalar-field excitation, enhanced in this case (β 0 = −4.35)by spontaneous scalarisation (cf.Sec.5.2.2).The central value ϕ c , which through collapse and bounce remains close to values of the order of ∼ −10 −5 , increases in magnitude to ∼ −2×10 −3 .This signal propagates outwards at the speed of light, rapidly leaves the region of the stalled shock, and reaches the extraction radius after about ∼ 0.56 s from the beginning of the simulation.
Our gauge choice does not allow us to follow the evolution of the inner region of the star beyond BH formation.Following [50], we terminate the evolution of the matter variables at the onset of BH formation in order to ensure numerical stability.At this point, the inner core has reached a compactness of ∼ 0.466, close to the BH value of 0.5.We are still able to gain insight into the late-time behaviour of the scalar field, however, by solving the wave equation (2.8) on the now frozen background (cf.Fig. 5.1).We observe in these evolutions that, as the NS (now spontaneously scalarised) collapses to a BH, the scalar field slowly relaxes to a flat profile as predicted by the no-hair theorems [18][19][20][21]121].
Overall, the entire dynamics strongly resembles GR.The scalar field is mostly driven by the matter evolution, which in turn is largely independent of the scalar field propagation.This point is illustrated in Fig. 5

Monopole gravitational-wave emission
Unlike GR, ST theories of gravity admit propagating monopole GWs.This breathing mode of the scalar field is potentially detectable with current and future GW interferometers which have therefore the potential of constraining the parameters of the theory.We now analyse the scalar GW signal extracted from our numerical simulations, separately discussing the effects of the EOS and the ST parameters.

Effect of the equation of state
As detailed in Sec.3.1, the hybrid EOS is a simplified EOS model that qualitatively approximates more sophisticated microphysical EOS in the core collapse context (e.g., [77]).The hybrid EOS is characterised by three adiabatic indices for the pre-bounce dynamics (Γ 1 ), the repulsion at nuclear densities (Γ 2 ), and the response of the shocked material (Γ th ).
The effect of the EOS on the emitted GW waveforms is explored in Fig. (ii) As the core collapses, the scalar field signal significantly grows in amplitude.
Although the first ∼ 0.02 s of the waveform appear to be rather insensitive to the EOS, the adiabatic indices strongly affect the total amount of time the star spends in the collapse phase before bounce.In the hybrid EOS, this is controlled by Γ 1 .Collapse is triggered by Γ 1 ≤ 4/3 and the smaller Γ 1 , the more rapid the collapse and the smaller the mass of the inner core that collapses homologously (e.g., [72,74]).We note that in reality (and in simulations using more realistic microphysics), Γ 1 is not a parameter.Instead, the effective adiabatic index is a complex function of the thermodynamics and electron capture during collapse [74,75].Figure 5.3 shows that core bounce occurs in our simulations at retarded time t − r ext /c ∼ 0.03, 0.04, 0.08 s (0.06, 0.07, 0.14) for model WH12 (WH40) and Γ 1 = 1.28, 1.3, 1.32 respectively.The bounce itself is a rapid process with a duration of ∆t ∼ 1 − 2 ms.
(iii) After bounce, the scalar field in the inner core settles down to a non-trivial profile, as illustrated in Fig. 5.1.The post-bounce value of ϕ at r ext , hence the value of h(t), encodes information about all three adiabatic indices.In particular, larger values of Γ 1 and smaller values of Γ th both produce stronger wave signals h(t), which can be intuitively understood as follows.Larger values of Γ 1 result in a more massive inner core in the pre-bounce stage because the speed of sound is larger and, hence, more matter remains in sonic contact in the central region.At bounce, this implies a more compact core and a correspondingly larger amplitude in the scalar wave.Smaller values of Γ th imply lower pressure in the shocked material, and, therefore, material that accretes through the shock settles faster onto the protoneutron star.In terms of microphysical processes, this effect is driven by neutrino cooling [74,75], which is not included in our simulations.In contrast, we find that Γ 2 has a relatively minor effect on this phase of the wave signal: the scalar field profile is only slightly more pronounced for lower values of Γ 2 , which result in a deeper bounce and a more compact postbounce configuration.Note that in the waveforms shown in Fig. 5.3, both Γ th and Γ 2 only affect the wave signal at and after bounce.This is expected, since Γ th only plays a role in the presence of shocked material and Γ 2 affects only the high density regime of the EOS not encountered in the collapse evolution prior to bounce.
(iv) Two of the simulations shown in Fig. 5.3 (namely Γ 1 = 1.28, 1.3; Γ 2 = 2.5; Γ th = 1.35) for the WH40 profile collapse to BHs.BH formation is triggered when the protoneutron star exceeds its maximum mass and is therefore facilitated and accelerated by smaller values of the adiabatic indices.BH formation generates a very large pulse in the scalar field which dominates the entire GW signal.Spontaneous scalarisation (marginally allowed for the value β 0 = −4.35chosen here) before BH formation further enhances the signal.The amplitude of the scalar field signal from this phase is more than an order of magnitude larger than the bounce signal in the absence of BH formation.We expect BH-forming collapse events to be the most promising source of monopole GWs in the context of ST theory.

Effect of the ST parameters.
As introduced in Sec.3.2, PN deviations from GR in ST theories only depend on two parameters, α 0 and β 0 .While α 0 mainly controls the perturbative deviation from GR, β 0 is responsible for non-linear effects such as spontaneous scalarisation.Our primary interest in this section is to explore the effect of these parameters on the detectability of signals with current and future GW detectors and, in particular, comparison with their sensitivity curves.Figures 5.4 and 5.5 show frequency domain waveforms S h (f ) compared with the expected (design) sensitivity curves S n (f ) of Advanced LIGO and the Einstein Telescope.We use the WH12 and WH40 initial profiles, together with the hybrid EOS with fiducial values Γ 1 = 1.3, Γ 2 = 2.5, and Γ th = 1.35 (cf.Sec.3.1).To better disentangle the effect of the two ST parameters, Fig. 5.4 (5.5) presents a series of simulations where only α 0 (β 0 ) varies while the other parameter is kept fixed at β 0 = 0 (α 0 = 3×10 −3 ).These two parameter sets overlap at α 0 = 3×10 −3 and β 0 = 0 and this specific simulation is shown in both figures.The location of our runs in the (α 0 , β 0 ) parameter space is shown in Fig. 3.1.Throughout our analysis, we consider optimally oriented sources placed at a fiducial distance of D = 10 kpc, i.e. within the Milky Way.
As mentioned above, the most pronounced feature in the emitted waveform arises from the collapse of the protoneutron star to a BH.As a consequence, the GW strains emitted during collapse of the M ZAMS = 40M profile WH40 are over an order of magnitude larger than the corresponding signals obtained from the collapse of the WH12 profile.BH formation (possibly enhanced by spontaneous scalarisation -see below) following the protoneutron star phase is the most promising signature of monopole GWs in the context of ST theory.
Simulations presented in Fig. 5.4 are performed in ST theory with β 0 = 0 and various values of α 0 , equivalent to Brans-Dicke theory with ω BD = (1 − 6α 2 0 )/2α 2 0 .Since spontaneously scalarised stars are not permitted in this regime, this set of simulations illustrates the effect of perturbative deviations from GR.In ST theory with α 0 ∼ 3 × 10 −3 , just compatible with the Cassini bound, GW signals generated by BH formation in our Galaxy, are marginally detectable by second-generation ground-based detectors and fall well into the sensitivity range of future experiments like the Einstein Telescope.Observation of a BH forming core collapse event with Advanced LIGO therefore has the potential of independently constraining ST theory at a level comparable with the most stringent present tests.Future third-generation observatories, on the other hand, will be able to push the constraint to new levels: α 0 ×10 −4 corresponding to |γ PPN − 1| 2 × 10 −8 ; cf.Eq. (3.7).On the other hand, our present results suggest that core collapse forming NSs (such as in our WH12 model) will at best allow for an independent confirmation of existing bounds, even when observed with third-generation observatories.
By analysing the curves in Fig. 5.4 quantitatively, we observe that the amplitude of the GW signal scales approximately as α 2 0 .One factor of α 0 is evidently due to the local coupling between the scalar field and the detector [see Eq. (3.23)].In our simulations, however, we find that the amplitude of the emitted scalar field ϕ also depends (roughly linearly) on α 0 .This second factor of α 0 is entirely due to the source dynamics and therefore separate from that arising in the coupling between the wave and the detector at the moment of observation.Even though the dynamics in the matter variables only mildly deviates from the GR case (cf.Fig. 5.2), such perturbative deviations from GR of the order α 0 can leave a significant imprint on the generation of monopole GWs.
The strongest effect of β 0 on the structure of NSs is that of allowing for spontaneously scalarised stars in the range β 0 −4.35.In fact, it is precisely the strength of this effect that enables binary pulsar observations to severely constrain β 0 as displayed in Fig. 3.1.For our simulations using values of β 0 significantly above the spontaneous scalarisation threshold of about −4.35, we only identify a relatively minor variation of the scalar wave with β 0 .This is well illustrated by the curves in Fig. 5.5 corresponding to β 0 = 0 and −2.Deviations of this kind can become particularly pronounced for β 0 −4.35 if the strongly non-linear effects of spontaneous scalarisation are triggered [33,34,51,91]; cf.Fig. 5.5.Whether this effect is triggered in our simulations and, in consequence, the shape and magnitude of the resulting waveform, critically depends on the stellar progenitor.
• If the core collapse leads to a protoneutron star that subsequently collapses to a BH, spontaneous scalarisation can be triggered by the high compactness reached shortly before BH formation, leading to a large enhancement of the GW signal.
In the bottom panel of Fig. 5.5, we compare the frequency-domain GW signals for the BH-forming WH40 progenitor.Spontaneous scalarisation occurs in a very strong way for the model with β 0 = −5 (already ruled out by current constraints) and leads to an enhancement of ∼two orders of magnitude in the amplitude compared to models that do not exhibit this strong non-linear behaviour (cases with β 0 = 0 and β 0 = −2).The waveform of the model with β = −4 (still allowed) is also somewhat enhanced by non-linear scalar field dynamics.Given the quantitative differences between these waveforms, present and future detectors have the potential of either observing scalar waves from BH-forming core collapse events or use their absence in the data stream to constrain the parameter β 0 beyond current limits.This will, however, require that other uncertain parameters such as the distance to the source etc. can be determined with high precision.• None of our simulations of the progenitor model WH12 leads to BH formation in the time simulated.This is so because this moderate-mass progenitor has a steep density gradient outside its core and thus a lower postbounce accretion rate.If no explosion is launched, a BH would still result, though on a timescale of O(10) s [26].Furthermore, we do not observe any signature of spontaneous scalarisation in the waveform or in the protoneutron star of the WH12 model, even for the extreme case β 0 = −5 (cf.Fig. 5.5, top panel).An analogous conclusion holds for collapse of ST polytropes, cf.Sec.3.3.The reason for this absence of spontaneous scalarisation in these models lies in the insufficient compactness of their protoneutron stars.At the end of our simulations, the protoneutron star in model WH12 has a compactness of m/r ∼ 0.05 (at r = 10 km), significantly lower than the threshold of ∼ 0.2 at which multiple families of stationary solutions appear [89].
The final compactness reached by NS remnants is naturally model dependent and the microphysics implemented in our analysis is greatly simplified by the use of the hybrid EOS.The possibility of triggering spontaneous scalarisation in stellar core collapse forming NSs (as opposed to BHs) clearly requires further exploration with more realistic finite-temperature EOS, which is left to future work.Dissociation of accreting heavy nuclei at the shock and neutrino cooling act indeed in the direction of lowering the effective adiabatic index in the postshock region, thus facilitating a more rapid increase in the protoneutron star's mass and compactness [122].We probe this expectation within our current framework by evolving the WH12 model with an adiabatic index Γ th artificially lowered to 1.25.With such a low value of Γ th , the shock stalls at a small radius and material accreting through the shock quickly settles onto the protoneutron star, driving up its mass and compactness.At the end of our simulation, at 0.7 s, the protoneutron star in this model has reached a compactness of ∼ 0.18 and is spontaneously scalarised.Configurations with non-trivial scalar-field profiles are energetically favoured over their weak-field counterparts and the dynamical evolutions naturally settle there.The GW strain S h (f ) increases by roughly two orders of magnitude when compared to runs performed using the more realistic value Γ th = 1.35.Galactic signals from spontaneously scalarised NSs, if formed in core collapse, will likely be detectable by Advanced LIGO even beyond the Cassini bound α 0 = 3 × 10 −3 .Given its observational potential, this topic definitely merits further investigation with more realistic microphysics.

Conclusions
This paper presents an extension of the open-source code gr1d [54] to ST theories of gravity.The required additions to gr1d can be summarised as follows: (i) generalisation of the flux and source terms in the high-resolution shock capturing scheme according to (2.34)-(2.39)as well as the constraint equations (2.21), (2.22) for the metric components; (ii) implementation of the evolution of the scalar field according to Eqs. (2.26)-(2.28)using standard finite differencing; (iii) outgoing radiation boundary condition for the scalar field (2.31), (2.32).
The scalar field furthermore introduces a new radiative degree of freedom propagating at the speed of light, which requires a smaller numerical timestep.All presented time evolutions start from one of two types of initial data, (i) polytropic models obtained by solving the time independent limit of the evolution equations and (ii) more realistic pre-SN models from zero-age main-sequence stars of masses 12 and 40M [100].
In this framework, we have simulated a large number of collapse scenarios which are characterised by five parameters: the linear and quadratic coefficients α 0 and β 0 determining the coupling function of the ST theory and the adiabatic indices Γ 1 , Γ 2 and Γ th , characterising the phenomenological hybrid EOS used in the time evolution.We summarise our main observations as follows.
• The most prominent GW signals are detected from the collapse of progenitor stars that form BHs after a protoneutron star phase (such as the M ZAMS = 40M model of [100]), as opposed to collapse events forming long-term stable NSs.The collapse of protoneutron stars to BHs is the most promising dynamical feature for monopole gravitational radiation in the context of ST theories.• The dynamical features in the matter fields (density, mass function, pressure) resemble closely those obtained in the general relativistic limit α 0 = β 0 = 0.In other words, the effect of the scalar field on the matter dynamics is weak.• The opposite is not true.The scalar radiation or GW breathing mode is sensitive to the specifics of the collapse dynamics as well as the choice of ST parameters α 0 and β 0 .The observed dependencies are of the kind one would intuitively expect.EOS and progenitors giving rise to more compact post-collapse configurations result in stronger radiation and the amplitude of the scalar wave sensitively increases with α 0 with approximately a quadratic dependence.• The ST parameter β 0 is known to generate strongly non-linear effects in the scalar field for β 0 −4.35, the so-called spontaneous scalarisation [33,34].For progenitors collapsing to BHs after a protoneutron star phase, transition of the central core to a spontaneously scalarised configuration before BH formation further enhances the outgoing GW signal.For progenitors forming NSs but not BHs, we do not find spontaneously scalarised configurations for physically plausible values of the adiabatic indices in our hybrid EOS.We attribute this to the stellar compactness achieved in those collapse scenarios being insufficient to trigger spontaneous scalarisation.This observation may be an artefact of our simplistic treatment of microphysics in our simulations • We have extracted waveforms from a large set of simulations and compared their amplitude for the case of a fiducial distance D = 10 kpc with the sensitivity curves of Advanced LIGO and the Einstein Telescope.Given the present constraints from the Cassini probe, α 0 3 × 10 −3 , scalar radiation may be marginally detectable from galactic sources.This offers the possibility of providing constraints on ST theory with GW observations in case of a favourable event occurring in the Milky Way.Considerable power is emitted at low frequencies f 10 Hz, thus making core collapse in ST theory ideal sources for future experiments such as DECIGO [112].
The impact of more realistic microphysic, as for example nuclear dissociation at the shock and neutrino cooling, on the compactness of the core and, thus, its degree of scalarisation, represents one key extension left for future work.Our analysis has shown that the massive increase in wave amplitudes due to spontaneous scalarisation and BH formation has the potential to drastically increase the range for detection.The O(10 3 ) waveforms generated for this work were completed in less than a week using O(10 2 ) CPU cores simultaneously.A moderate increase in the computational resources will make simulations using tabulated finite-temperature EOS feasible.Since the matter fields' dynamics is very similar to the GR case, one may perhaps take advantage of existing GR simulations and simulate the scalar field evolution using such GR results as backgrounds.A further numerical improvement may consist in computing (perhaps iteratively) approximate initial conditions for the scalar field from existing pre-collapse stellar models (such as WH12 and WH40 used in this paper), in order to reduce the brief unphysical transient in the GW signal.
Aside from the treatment of the microphysics, our study offers further scope for extension.The effects of multiple scalar fields in ST theories on gravitational collapse remains largely unknown in spite of some early studies [13] (see [43] for an analysis of static NS solutions in this framework), but represents a relatively minor addition to our code.The same holds for ST theories with non-vanishing potential, as for example massive fields [99,123].
As GW physics and astronomy are ushering in a new era, the community will be offered a wealth of unprecedented opportunities to observationally test generalisations of GR.Stellar collapse clearly offers a vast potential for such fundamental tests.

Figure 3 . 1 .
Figure 3.1.Experimental constraints on the ST-theory parameters (α 0 , β 0 ) entering the coupling function F .The shaded area is currently ruled out by observations; GR lies at α 0 = β 0 = 0.The most stringent constraints on α 0 are provided by the Cassini space mission while the binary-pulsar experiments impose strong bounds on β 0 .Circles mark our choices for (α 0 , β 0 ) used in Sec. 5.This figure is produced using the data published in Fig.6.3 of Ref.[1].

Figure 3 . 2 .
Figure 3.2.Mass-density (left panel) and scalar-field (right panel) profiles for the initial data considered in this study.In particular, dashed and dotted lines show the M ZAMS = 12M (WH12)and M ZAMS = 40M (WH40) pre-SN models of Woosley and Heger[100] while the solid lines show three Γ = 4/3 polytropes generated in ST theories with (α 0 , β 0 ) = (0.001, −3), (0.003, 0), (0.01, −5).The mass-density distributions of all three ST polytropes are indistinguishable from their GR counterparts.The more realistic models WH12 and WH40 mostly differ from the polytropic ones through the presence of outer low-density layers.Note that we cut the WH models at rs = 2 × 10 4 km and pad them with an artificial atmosphere of ρatm = 1 g cm −3 .The scalar field is more pronounced in models with higher α 0 , but the low compactness of these models prevents spontaneous scalarisation.The scalar field is initialised to zero when the WH models are evolved.

Figure 4 . 1 .
Figure 4.1.Convergence test of stellar collapse in ST theory with α 0 = 10 −4 and β 0 = −4.35.A polytropic core is collapsed using the fiducial hybrid EOS for three different resolutions (see text for details).The top panels show the evolution of the mass density ρ (left) and the scalar-field ϕ (right) for the highest-resolution run a t = 7 (grey), 28 (yellow), 37 (green), 38 (blue), 57 (red) ms after starting the simulations.Bounce happens at t ∼ 38 ms and the shock reaches the surface of the star at t ∼ 131 ms.The bottom panels show the self-convergence properties of the gravitational mass m (left) and the scalar field ϕ (right) at the same times.As detailed in the text, solid and dotted (dashed) lines are expected to coincide for second-(first-) order convergence.We initially observe second order convergence which decreases to first order due to (i) the shock capturing scheme when a discontinuity forms at bounce and (ii) numerical noise in the scalar field propagating in from the outer boundary.

Figure 5 . 1 .
Figure 5.1.Collapse of a Γ = 4/3 polytrope (top), the 12M (centre) 40M (bottom) pre-SN profiles of [100] in ST theory with α 0 = 10 −4 and β 0 = −4.35,assuming Γ 1 = 1.3, Γ 2 = 2.5 and Γ th = 1.35.The evolution of the mass density ρ (left) and the scalar field ϕ (right) is shown as a function of the radius r at various timesteps t − t B = −0.01,−0.001, 0.001, 0.1, 0.3, 0.375, 0.38 s, measured from the bounce time t B .Timesteps are coloured from darker (early times) to lighter (late times) solid lines as labelled; initial profiles are shown with black dashed lines.The inset in the bottom right panel shows the wide variation of the scalar field when a BH is formed.An animated version of this figure is available online at Ref. [58].

Figure 5 . 2 .
Figure 5.2.Evolution of the central values of the mass density ρc (left panels) and lapse function αc (right panel) through collapse, bounce and late time evolution in ST theory with α 0 = 10 −4 and β 0 = −4.35.We use the hybrid EOS with fiducial parameters (Γ 1 = 1.3, Γ 2 = 2.5 and Γ th = 1.35) and three different initial profiles: ST polytrope (top), WH12 (centre) and WH40 (bottom).Gray dashed lines mark the bounce time t b ; the WH40 profile first collapses to a protoneutron star and then to a BH at t BH ∼ 0.46 s marked by grey dotted lines.Relative differences with analogous simulations performed in GR are shown in the lower subpanels (red lines).Deviations in the dynamics are very small: of the order of |∆ρc|/ρc ∼ 10 −5 and |∆αc|/αc ∼ 10 −6 .

. 2 ,
where the central values of the mass density ρ c and lapse function α c obtained in ST theory and in GR are compared.The relative differences between these two scenarios are about |∆ρ c |ρ c ∼ 10 −5 and |∆α c |/α c ∼ 10 −6 throughout collapse, bounce, and (eventually) BH formation.

Figure 5 . 3 .
Figure 5.3.Effect of the hybrid EOS adiabatic indices on the emitted monopole gravitational waveform h(t) ∝ rϕ [cf.Eq.(3.23)].The signal is plotted against the retarded time t − r/c at the extraction radius.Simulations are performed using the preSN initial profiles WH12 (top) and WH40 (bottom) in ST theory with α 0 = 10 −4 and β 0 = −4.35.The curves encode the value of Γ 1 in their brightness (colour) and the value of Γ th in their line style: Γ 1 = 1.28 (red), 1.3 (blue), 1.35 (green);Γ th = 1.35 (solid), 1.5 (dashed).For each of these combinations, two curves are present: circles mark simulations with Γ 2 = 2.5, while no symbols are shown for Γ 2 = 3.For some cases, these two curves overlap to such high precision that they become indistinguishable in the plot.The lower-case Roman labels refer to the key phases of the GW signal described in Sec.5.2.1: (i) initial pulse of the spurious radiation; (ii) collapse and bounce; (iii) NS configuration; (iv) BH formation.The bounce time is marked with vertical dotted lines following the same colour codes of the other curves.Note that Γ 1 is the only adiabatic index that has an effect on the bounce time.Waveforms presented in this figure are available at[58].
5.3, where we show time-domain monopole waveforms h(t) ∝ rϕ [cf.Eq. (3.23) with ϕ 0 = 0] for various choices of Γ 1 , Γ 2 , and Γ th .All simulations shown in Fig. 5.3 are performed in ST theory with α 0 = 10 −4 and β 0 = −4.35, the lower limit of β 0 compatible with present observations, using the WH12 and WH40 initial profiles.We plot the GW signals as a function of the retarded time t − r/c, such that the origin corresponds to a single light-crossing time at the extraction radius r ext = 3 × 10 4 km.The structure of the emitted signals displayed in Fig. 5.3 consists of the following four main stages.(i) The initial pulse of spurious radiation arises from the initialisation of the scalar field, as already pointed out in Secs.3.3 and 5.1.This pulse propagates outwards and leaves the extraction region after a retarded time of about 0.006 s.