Gravitational dynamics in the Higgs-dark matter sector toy model -- the field theoretic perspective

The objective of the paper was to examine the course and results of gravitational collapse in the toy model of the Higgs--dark matter sector. The real part of the Higgs doublet written in the unitary gauge was modelled by a neutral scalar field. Two dark matter candidates were introduced to the model. One of them is the dark photon, which can be associated with one of the included $U(1)$ gauge fields. The other one is a complex scalar field, charged under the second $U(1)$ gauge field that represents the Maxwell field. Additionally, non-minimal couplings of both scalars to gravity were taken into account. There were two coupling channels between the ordinary and dark matter sectors, that is a kinetic mixing between the $U(1)$ gauge fields and the Higgs portal coupling among the scalars. Numerical, fully non-linear simulations of the investigated gravitational process were performed within the model of interest and its truncated version, with scalars minimally coupled to gravity. The structures of emerging singular spacetimes were analyzed via the behavior of dynamical horizons forming in them. The features of dynamical black holes appearing in the spacetimes were described as functions of the parameters of the model. A set of quantities associated with an observer moving with the collapsing matter was proposed and calculated for the dynamical spacetimes.


Introduction
Considering dynamically formed spacetimes containing complex matter systems, their geometric structures and features of objects existing within them, such as black holes, white holes or wormholes, are in most cases impossible to describe analytically. Moreover, even for relatively simple gravity-matter systems, such as for example a gravitationally selfinteracting electrically charged scalar field, structures of dynamically formed spacetimes differ significantly from their non-dynamical counterparts. The dynamical Reissner-Nordström spacetime is considerably different from the static one [1]. There exists a central spacelike singularity and a null Cauchy horizon, which is also singular due to the the mass inflation phenomenon. In the static case, the central singularity is timelike and the Cauchy horizon is non-singular. For the above reasons, numerical investigations are usually performed to obtain a description of spacetimes, which were formed or evolving during dynamical processes, such as gravitational collapse or matter accretion.
A set of various matter-geometry systems were hitherto considered for establishing outcomes of dynamical processes within them. The simplest ones are self-interacting scalar fields, neutral [2] and electrically charged [3][4][5][6][7]. The gravitational evolution of the inflaton field was tracked in [8,9] and within the f (R) gravity in [10,11]. The role of dilatonic and phantom couplings during the collapse within the Einstein-Maxwell theory was elaborated in [12][13][14][15]. The gravitational evolution of scalar fields in the Brans-Dicke theory was inspected in [16][17][18]. The influence of the dark sector, i.e., the presence of dark matter and dark energy in spacetime, on the collapse of an electrically charged scalar field was examined in [19].
Investigating gravitational dynamics within various physical models is often accompanied by inspecting additional issues, related to the dynamics. Pair creation during collapse and subsequent evaporation of black holes, i.e., quantum effects accompanying the collapse, were described in [6,20,21]. The issues of the cosmic censorship conjecture and the information loss problem were raised in [6][7][8][9]. The impact of changing the number of dimensions on the course and outcomes of the collapse was described in [22], while the role of various topologies in [23,24]. The problem of measuring time using a collapsing scalar field was discussed in [25][26][27].
The existence and features of dark matter and the Higgs field have been a viable topic of both experimental and theoretical research recently. How these two matter components influence dynamical gravitational processes is an intriguing issue, which we decided to raise in the current paper. To describe the matter sector we opted to use the model containing two scalar fields and two U (1) gauge fields. One of the gauge fields describes electromagnetism and the second one can be associated with the dark photon [28][29][30] (see also [31] for a modern review). Among the scalars, one is real and not charged under either of the U (1) gauge fields. This field may represent a real part of the Higgs doublet written in the unitary gauge. The second one, which is complex scalar charged under a U (1) field may also represent a stable dark matter candidate [30,32,33]. Additionally, we also considered a possibility that both scalars possess non-minimal couplings to gravity. Moreover, there are two coupling channels among the ordinary matter sector, which consists of the real scalar and electromagnetic field, and the dark sector, composed of the complex scalar and an additional U (1) field. These channels are given by a kinetic mixing between the U (1) gauge fields and the Higgs portal coupling among the scalars.
The double null formalism in investigating dynamical gravitational processes is beneficial as it allows to track their course from an asymptotically flat region situated close to past null infinity up to the forming singularity. It means that both the external and internal geometrical structures of the forming objects can be resolved. Formation of additional spacetime special regions, such as event and Cauchy horizons, dynamical horizons, wormhole throats, vacuum bubbles, etc., can be followed within this formalism. A comprehensive review on employing the double null formalism in investigations of gravitational collapse of matter can be found in [34].
The aim of the current paper is investigating the course and results of gravitational dynamics within the toy model of the Higgs-dark matter sector using the double null formalism. The structures of emerging spacetimes will be analyzed via locations of dynamically formed horizons and singularities. The dependence of the characteristics of forming objects, precisely dynamical black holes, on the parameters of the model of interest, will be presented. Additionally, a set of observables will be proposed to describe the outcomes of the processes.
The paper is organized as follows. In section 2 the theoretical model of the mattergravity system of interest is presented. Section 3 contains basic information on solving the derived evolution equations and particulars of the results interpretation. In sections 4 and 5 the course and outcomes of the investigated process are discussed. The summary of the obtained results is placed in section 6. A presentation of the numerical setup for solving the equations of motion within the examined model and tests of the code constitute appendix A.

Theoretical model of the evolution
The action functional for our theory consists of three additive parts and is given by where L G = 1 16πG R is the Einstein-Hilbert action and the matter part is In the above formulas B µν ≡ 2∇ [µ B ν] is the strength tensor for the dark photon field B µ , F µν ≡ 2∇ [µ A ν] is the Maxwell tensor for an ordinary electromagnetic field A µ . In the employed representation the gauge fields are coupled through the kinetic mixing term α 2 F µν B µν . The strength of the mixing is controlled by the coupling constant α, whose upper value is experimentally constrained to be less than 10 −3 [35].
The remaining matter is represented by two scalar fields. The real neutral scalar field h may be regarded as the neutral part of the Higgs field (for an appropriate choice of the λ h and m 2 h parameters, roughly speaking λ h ∼ 0.13 and m 2 h 0). The complex scalar field X is charged under a dark U (1) field (the covariant derivative is d µ = ∇ µ + ieB µ , where e is its charge) and may be a model of a dark matter candidate, see for example [36][37][38][39] and references therein. The fields h and X interact through the Higgs portal term λ hX h 2 |X| 2 . An influence of non-minimal couplings of the scalars to gravity is described by the L GM part of the action.
Before we proceed further it is convenient to diagonalize the gauge sector. For this purpose we define C µ ≡ B µ + αA µ and write In this representation the X field is also charged under the ordinary electromagnetic U (1) group, but its charge is suppressed by the small α parameter (in the dark matter literature this type of field is called millicharged dark matter, provided that m 2 X > 0). On the other hand, if m 2 X < 0 the appropriate dark matter candidate is the dark photon which also becomes massive due to the existence of the non-zero vacuum expectation value for the X field.
To write the equations of motion in a fully explicit form we need to specify the metric ansatz. We will be using the double null (u, v, Θ, Φ) spherically symmetric line element [40] where u and v null coordinates are called retarded and advanced time, respectively, and dΩ 2 = dΘ 2 + sin 2 ΘdΦ 2 is the line element of the unit sphere, where Θ and Φ are angular coordinates. Such a coordinate choice determines the double null spacetime foliation [41]. From now on, partial derivatives with respect to the null coordinates will be denoted by ,u and ,v . Vector and tensor elements indexes will be marked by adequate sub-and superscripts. Before we present the equations of motion it is convenient to introduce a set of dimensionless quantities. For this purpose we make the following rescaling: where M 2 P = 1 8πG is the reduced Planck mass squared. In what follows we will be using the rescaled dimensionless fields and coordinates, and moreover for the purpose of making the notation more readable we will drop the tilde above them. To recapitulate, the Lagrangian densities of the action (2.1) written in terms of the rescaled variables are The equations of motion for the evolving fields were obtained via the variation of the action (2.1). The equations for the Higgs field h, the complex scalar field X and gauge fields A µ and C µ are 14) The Einstein equations derived by varying the action (2.1) with respect to gravitational field complement the above set of equations, which describes the dynamics of the inspected dynamical system.
Regarding spherical symmetry, in the chosen coordinate system, the only non-vanishing components of the field tensors are F uv , F vu , C uv and C vu . Due to the gauge freedom A u → A u + ∇ u θ and C u → C u + ∇ u θ , where θ = A v dv and θ = C v dv, the only nonzero four-vector components are A u and C u . They are functions of retarded and advanced time.
The evolution equations of the fields from within the examined theoretical setup expressed in double null coordinates, i.e., using the assumed line element (2.5), are the following. The dynamics of the Higgs field described covariantly by (2.10) is governed by the equation The equations of motion of the scalar field X (2.11) and its complex conjugate X * (2.12) are given by The equations of motion of the only non-zero components of the respective four-vectors of the U (1) gauge fields define the gauge fields related charges T and Q, whose evolution equations are The above pairs (2.18), (2.20) and (2.19), (2.21) stem from a separation of the v-components of the gauge fields equations (2.13) and (2.14), respectively, into two first-order differential equations. Both T and Q depend on retarded and advanced time. They correspond to charges associated with the C µ and A µ fields contained within a sphere of a radius r(u, v), on a spacelike hypersurface containing the point (u, v). In the latter case, the charge is simply the electric charge. The stress-energy tensor for the considered theory is (2.24) Its non-vanishing components expressed in double null coordinates are composed of with |X| = √ XX * . Finally, the gravitational field equations are obtained using the adequate components of the Einstein tensor resulting from the line element (2.5) and the stressenergy tensor (2.22) components The set (2.15)-(2.21), (2.33)-(2.36) forms a complete set of equations of motion which describe the dynamics of the system of interest. In order to solve the obtained system of equations of motion, a set of auxiliary variables and quantities defined as was introduced. The substitutions enabled us to rewrite the second-order differential equations from the set (2.15)-(2.21), (2.33)-(2.36) as first-order ones. Additionally, the real fields X 1 and X 2 were introduced instead of conjugate fields X and X * simply according to X = X 1 + iX 2 and X * = X 1 − iX 2 . These relations result in (2.39) The final system of equations of motion, which governs the investigated evolution yields

Details of computer simulations and results analysis
The final system of evolution equations (2.40)-(2.57) was solved numerically. Appendix A presents the details of the numerical code and its tests. The equations were solved in the region of the (vu)-plane, which is shown on a Carter-Penrose diagram in figure 1. The background spacetime is the one containing a dynamical Reissner-Nordström black hole with a spacelike central singularity and a Cauchy horizon, both surrounded by an event horizon [3]. In all current simulations the computational domain spread from v = 0 to v = 7.5 and from u = 0 to u = 22.5. The only arbitrary input data were initial profiles of the evolving fields, posed on the null u = 0 hypersurface. The initial profile of the h field was Gaussian and the complex field was represented by the trigonometric profile The profiles were one-parameter families with free family parameters being the amplitudesp h andp X , which regulate the strength of the particular field gravitational selfinteraction [42]. The values of the remaining constants v 0 = 1.3, D = 0.21, v i = 0 and v f = 2.5 were chosen arbitrarily and were the same in all simulations. The parameter δ ph determining the amount of initial charge was equal to π 2 . The employed profiles (3.1) and (3.2) ensure that the spacetime slice at the initial null hypersurface is regular and hence their choice is representative for the examined evolutions as the outcomes of computations are independent of the profiles types [25,42].
The non-zero value of the electric coupling constant does not influence the results of the gravitational collapse [12]. Hence, it was invariable in all evolutions and equal to 1. Due to the argument relating the field h with the neutral part of the Higgs field raised in section 2, the parameter λ h was set as equal to 0.13 during computations. The value of the parameter α was set as equal to 10 −5 in all simulations, which is in agreement with its experimental constraint [31,35,43].
From equations (2.20) and (2.21) it can be seen that the difference between the absolute values of the U (1) charges T and Q is of the order of α 2 . Taking into account the value of α adopted in the computations, the difference is of the order of 10 −10 . It was confirmed in the course of numerical simulations and for this reason the text of subsequent sections of the paper will refer only to one of the charges.
There is a significant number of quantities to be set to completely define a particular collapse case (the parameters of the model λ X , λ hX , m 2 , ξ X , ξ h , α m and field amplitudes p h ,p X ). For this reason, a vast collection of simulations has been run, which involved possible configurations of the above parameters from within their admissible ranges. It should be emphasized that the spacetime structures as well as observables and fields spacetime distributions presented in sections 4 and 5 are representative for the particular case. Moreover, to make the presentation of the results as clear as possible, the field amplitudesp h andp X were set as equal to each other, with no loss of generality.
Spacetime structures stemming from the studied dynamical processes will be shown on Penrose diagrams. The diagrams present r = const. contour lines in the (vu)-plane. The outermost thick line refers to r = 0. The spacetime regions depicted on Penrose diagrams will be those essential for the performed analyses. The structures of spacetimes will be described qualitatively at first and then analyzed via the behavior of selected black hole features listed in the next section when the parameters of the model change.

Horizons
Since the spacetimes obtained as outcomes of the gravitational evolutions of interest are dynamical, global characteristics are of limited use in their structures interpretation. Instead, local notions ought to be employed in order to interpret the obtained structures. One of the most straightforward tools to describe causal relations within dynamical spacetimes are local horizons, which are lines of expansion The line r ,v = 0 will be marked on the spacetime diagrams as a green solid line. The horizon corresponding to the r ,u = 0 line is not present in neither of the obtained spacetimes. As null coordinates tend to infinity, the dynamics of the spacetime-matter system diminishes and hence in these regions local horizons coincide with the global ones. In particular, considering the placement of the numerical domain within spacetime, the location of the event horizon of the black hole forming during the gravitational collapse can be determined by inspecting the behaviour of a local horizon as v tends to infinity. For large values of advanced time, the local horizon settles along a hypersurface u = const. indicating the location of the null event horizon in spacetime.
For purposes of interpretation of the results and numerical checks of the code, a set of quantities characterizing the emerging black hole will be inspected. These are the charges Q and T and mass contained within the event horizon. The mass will be calculated as the quasi-local Hawking mass for two coupled gauge fields A µ and C µ [19] It describes mass contained in a sphere of a radius r (u, v) on a spacelike hypersurface containing the point (u, v). The values of Q, T and m H will be obtained at a local horizon in the area where it becomes null, i.e., on the boundary of the numerical domain, where v = v f . Additionally, u-locations of the black hole event horizons and their radii will be also read out at the same point and will be used in the results interpretation.
The above quantities characterizing emerging black holes will be investigated within the model parameters ranges which, apart from the scalar-gravity couplings, are limited by the model determinants. In the case of the constants ξ X and ξ h , their ranges were established such that the apparent horizons of the nascent black holes are separated from the initial data hypersurface. This is equivalent to ensuring that the calculations begin in the region outer with respect to the black hole horizon which can be regarded as nearly flat.

Observables
Another way of interpreting the results of gravitational collapse within the model of interest is inspecting a set of local spacetime quantities, whose derivation is the following.
The tetrad connected to the observer moving with the medium, associated with the double null metric (2.5) is [44,45] e 3 = r sin θdφ, (3.8) or, alternatively, may be written as The energy densityρ, energy fluxesf , radialp r and tangential pressuresp t as seen by the observer can be defined asρ Moreover, defining two null vectors the surface gravity for a dynamic black hole is [44,45] The quantities used for the results interpretation were the energy densityρ, radial pressurep r , pressure anisotropyp a ≡p t −p r and local temperature T l = κ l 2π with the surface gravity as defined in [44,45]. It is worth emphasizing that there is no unique definition of surface gravity for dynamical black holes, as dynamical spacetimes within which they emerge do not admit a Killing vector field [46]. The first three of the above quantities are related with the stress-energy tensor components viaρ and, taking into account (2.22), for the inspected gravity-matter model considered in double null coordinates are given by the following relations: . (3.29) 4 Collapse dynamics within a model involving scalars minimally coupled to gravity To begin with, a set of results of gravitational collapse within a simplified version of the investigated model will be presented. The non-minimal couplings of scalars to gravity will be omitted, that is ξ X = ξ h = 0.

Spacetime structures
The structures of singular spacetimes resulting from the dynamical evolution of fields are presented in figure 2. All the spacetimes contain a spacelike central singularity along r = 0 surrounded by an apparent horizon that is either spacelike or null, when situated along a u = const. hypersurface. For positive values of m 2 the resulting spacetime is a typical dynamical Schwarzschild spacetime. The apparent horizon which is spacelike for small values of advanced time, settles along a null hypersurface of constant retarded time as v tends to infinity thus indicating the location of the event horizon in spacetime. When m 2 is negative, the apparent horizon has two spacelike sections, for small and large values of advanced time, separated by a null section. Such a behavior has been already observed during collapses within models containing more than one scalar [13,14].

Black hole characteristics
The characteristics of emerging black holes were investigated for the case with the following parameters: λ X = λ hX = 0.1, m 2 = 0.25, α m = 0 andp X =p h = 0.05. The spacetime structure for this collection is shown in figure 2(a). While the dependence on the particular parameter of the model is presented, the remaining ones are as above.
The dependence of the u-locations of the event horizons, radii and masses of black holes formed during the gravitational collapse within the model of interest as functions of m 2 , α m , λ X and λ hX is presented in figure 3. The u-locations of the black hole event horizons increase as m 2 and λ hX increase, decreases with α m and is not affected by a value of λ X . The changes of the black hole radii and masses are the same qualitatively, that is they decrease with an increase of m 2 and α m and increase with increasing quartic couplings.  The changes of u eh , r eh and m eh H against m 2 become much less dynamical as the mass parameter becomes bigger than −0.15 and the dependence on λ X is very weak within its whole range. Their changes with λ X , m 2 and α m are monotonic, while in the dependency on λ hX there is a discontinuity for λ hX ≈ 0.0515.
The black hole charges related to the U (1) gauge fields, Q eh and T eh , versus m 2 , α m and λ X , λ hX are depicted in figure 4. An increase of Q eh with the mass parameter is small when the parameter exceeds −0.05. When its value is smaller, the dependence is not monotonic, the charge slightly decreases for small values of m 2 and after reaching a minimum around m 2 = −0.4, starts to increase. Q eh weakly decreases with α m and both quartic couplings with a characteristic discontinuity around λ hX = 0.0515, which was also observed for the black hole characteristics discussed above, as shown in figure 3(d).   tionally, lines of constant r were plotted on the graphs in order to visualize their behavior on a Penrose diagram when the spacetime is flat. The plots will serve as a comparison for further analyses of singular spacetimes. In a spacetime that does not contain a singularity, maximal absolute values of the energy density, radial pressure, pressure anisotropy, as well as the neutral scalar field and the moduli of the complex scalar field are distributed along a null direction of constant advanced time. This is a direction in which a propagation of peaks of the field functions, whose profiles (3.1) and (3.2) were imposed on the initial hypersurface u = 0, advances as the dynamical collapse proceeds.

Observables and fields
Observables and fields distributions in a spacetime containing a black hole resulting from a collapse with the same model parameters as the non-singular one described above and higher values of field amplitudes are presented in figures 7 and 8. The structure of the spacetime is presented in figure 2(b). As in the case of a flat spacetime, peaks in absolute values of the energy density, radial pressure, pressure anisotropy and both neutral and complex scalar fields are observed along the null direction which indicates the field propagation in spacetime. Moreover, an increase of absolute values of the energy density, radial pressure and pressure anisotropy is visible in a close vicinity of the central singularity. A similar behavior is observed for the absolute values of the field h, while the complex scalar field with the modulus |X| does not display such a behavior as its values increase less considerably as the singularity is approached. All the quantities of interest remain finite within the whole spacetime, up to a close vicinity of the central singularity, as far as it is reachable during numerical computations. The energy density is positive within the whole spacetime. The sign of the radial pressure and pressure anisotropy functions change as the singularity is approached, in comparison to the region of their high values located along a constant value of advanced time. The local temperature calculated along the black hole apparent horizon using the definition of surface gravity for dynamical spacetimes adopted in the current paper (3.23) is positive and increases monotonically with advanced time. Such a behavior as v → ∞ is intuitive, as the black hole radius and mass also increase with advanced time.

Spacetime and matter dynamics within the Higgs-dark matter toy model
The outcomes of simulations of gravitational dynamics within the non-truncated version of the model of interest, that includes the non-minimal scalars-gravity couplings, will be presented twofold. The first set of results will refer to both non-zero coupling constatnts ξ X and ξ h . The second set will comprise the evolutions with one non-vanishing scalar-gravity coupling, either ξ X or ξ h .

Spacetime structures
Spacetime structures stemming from gravitational evolutions with ξ X = 0 and ξ h = 0 which lead to a formation of a black hole are presented in figure 9. For the cases with positive values of m 2 , presented on the plots 9(a) and 9(b), the spacetimes are dynamical Schwarzschild spacetimes with a spacelike singular r = 0 line surrounded by an apparent horizon that settles along a u = const. hypersurface as v → ∞. The non-vanishing coupling constants between scalar fields and gravity influence the behavior of the apparent horizon in the region where it changes its character from spacelike to null (between v equal to 1.2 and 2.3). There are several timelike, spacelike and null portions of the horizon, in contrast to a smooth behavior which was observed in the minimally coupled case, in figure 2. When the parameter m 2 is negative, the collapse results in the dynamical Reissner-Nordström-type spacetime. It contains a spacelike central singularity along the vanishing r line. The singularity is surrounded by two apparent horizons. The outer one is spacelike for small values of advanced time and becomes null as v tends to infinity, indicating the location of the event horizon. Similarly to the case discussed above, there are several timelike, spacelike and null portions of the apparent horizon in the region where it transforms from spacelike to null. The inner apparent horizon is spacelike within the whole domain of computations. An appearance of additional apparent horizons has been already observed during a collapse involving exotic matter, precisely phantom scalar fields [14]. There also exists a Cauchy horizon in the emerging spacetime. It is located at an infinite-v hypersurface. Its existence is manifested by a collection of r = const. lines that settle along constant-u and are located between the apparent horizon and the singularity [3].
Penrose diagrams of spacetimes stemming from the gravitational collapse within the investigated model with one non-zero scalar-gravity coupling constant are shown in figure 10. In all cases the spacetimes are dynamical Schwarzschild spacetimes. Central singularities along r = 0 are spacelike within them and surrounded by apparent horizons, which are spacelike for small values of v and become null as v increases. The behavior of the apparent horizons in the regions where they transform from spacelike to null again involves several changes of their character between timelike, spacelike and null. These changes are more significant when ξ X = 0, in comparison with evolutions with a non-vanishing ξ h .

Black hole characteristics
In the case of both non-zero scalar-gravity coupling constants, the characteristics of forming black holes were examined for the evolution with the following parameters: λ X = 0, λ hX = While the dependence on the particular parameter of the model is presented, the remaining ones are as above.
The dependence of the u-locations of the event horizons, radii and masses of black holes formed during the gravitational collapse within the model of interest as functions of m 2 , α m , λ X , λ hX and ξ X , ξ h is presented in figures 11 and 12 for the cases of two and one non-zero scalar-gravity couplings, respectively. The dependencies of all these quantities on the mass parameter are qualitatively the same in both cases. After a significant increase of u eh and decrease of r eh and m eh H up to m 2 equal approximately −0.15 their values vary much less. This behavior is similar to the one observed in the minimally coupled case discussed in section 4.2. All the discussed characteristics decrease with α m in both examined cases, except u eh when ξ X = 0, whose value, after a decrease, begins to increase for large values of α m . The decrease was also observed in the simplified version of the model with ξ X = ξ h = 0. In the case of two non-vanishing scalar-gravity coupling constants, the relations between u eh , r eh , m eh H and the quartic couplings are qualitatively the same, i.e., the first of the quantities increases and the remaining ones decrease as λ X and λ hX increase. In the dependencies on λ X in the case when ξ X = 0, which are qualitatively the same as the ones discussed above, two discontinuities around λ X equal to 0.013 and 0.0815 appear. The u-locations of the event horizons do not change with λ hX , while the radii and masses of nascent black holes increase and decrease with the parameter, respectively. The changes of u eh , r eh and m eh H with the quartic self-interaction coupling constant of the Higgs field are qualitatively the same in both cases of two and one non-zero scalar-gravity couplings. There is a maximal value of u eh and minima of r eh and m eh H , present at ξ h = 0. The relations between the black hole features and the quartic self-interaction coupling constant of the complex scalar field in the case of ξ X = 0 and ξ h = 0 also display a maximum in the case of the u-location of the event horizon and maxima in the cases of black hole radii and masses, but these shallow extrema are located at ξ X = −0.5.
The black hole charges related to the U (1) gauge fields, Q eh and T eh , as functions of m 2 , α m , λ X , λ hX and ξ X , ξ h are depicted in figures 13 and 14. When both non-minimal scalargravity couplings do not vanish, the relation between the charges and m 2 is similar to the one observed for the minimally coupled case, described in section 4.2. It possesses a minimum around m 2 = −0.35 and its variations become small when the parameter becomes bigger than −0.05. For the case with ξ X = 0, Q eh decreases within the whole m 2 range and the changes become less significant when the mass parameter exceeds 0.2. The dependency of the charge on α m displays a minimum at 0.06 when both non-minimal couplings are   non-zero and is monotonically increasing when only ξ h = 0. The charge increases as λ X and λ hX increase when both scalar-gravity couplings do not vanish. It decreases with λ hX and increases with two discontinuities at 0.013 and 0.0815 as λ X increases in the case with ξ X = 0. The dependencies against ξ X and ξ h possess local minima, situated at ξ X = −0.25 and ξ h equal to −0.1 and −0.45 for the cases of two and one non-minimal scalar-gravity couplings, respectively.

Observables and fields
The (vu)-distributions of observables defined in section 3.2 and the evolving scalar fields along with the relation between local temerature along the apparent horizon and v will be depicted and discussed for selected spacetimes whose structures were shown in section 5.1. Two representative cases with spacetimes resulting from gravitational evolutions with two non-zero scalar-gravity coupling constants will be presented. The behavior of the quantities of interest within spacetimes, which stem from collapses when one of the non-minimal couplings vanishes are qualitatively the same.  Their signs in the spacetime region further away from the singularity vary, which is a result of non-minimal couplings of scalars to gravity. It can be deduced when compared with the behavior of observables in spacetimes in the minimally coupled case, which was discussed in section 4.3. The changes of the black hole local temperature, which is positive, along the apparent horizon as v increases are not monotonic. It increases for small values of advanced time, reaches a maximum in the v-range, within which the apparent horizon changes its character from spacelike to null and then, after a slight decrease. the temperature increases along the null segment of the horizon.
The spacetime distributions of observables and field functions stemming from a dynamical evolution, which results in a Reissner-Nordström-type spacetime are shown in figures 17 and 18, respectively. The structure of the spacetime was depicted in figure 9(d). As in the cases discussed above, a region of high absolute values of the energy density, radial pressure, pressure anisotropy and the neutral scalar field is situated along the direction of initial peaks propagation along constant v. This increase in values does not appear for the moduli of the complex scalar field. The values of all discussed quantities increase considerably beyond the inner apparent horizon, where the r = const. lines settle down along null hypersurfaces of constant retarded time indicating the location of the Cauchy horizon at infinite advanced time, and diverge as the central singularity is approached. The local temperature in the studied case is positive and behaves in the same manner as in both cases discussed above. It increases for small values of advanced time, where the apparent horizon is spacelike. Then, it reaches an extremum in the region, where the character of the apparent horizon transforms from spacelike to null. Afterwards, after a slight decrease, it begins to increase along the null part of the horizon.

Conclusions
Gravitational dynamics within a toy model which involves a Higgs field and two dark matter candidates was investigated. The theoretical setup consisted of two scalar fields and two U (1) gauge fields. The gauge fields described electromagnetism and the dark photon, which was one of the dark matter candidates. One of the scalars was not charged under any of the U (1) gauge fields and represented a real part of the Higgs doublet written in the unitary gauge. The second scalar was complex and charged under a U (1) gauge field. It represented a second stable dark matter candidate. A possibility that both scalars possess non-minimal couplings to gravity was also included. The coupling channels among the ordinary matter sector, which consisted of the real scalar and electromagnetic field, and the dark sector, composed of the complex scalar and an additional U (1) field, were given by a kinetic mixing between the U (1) gauge fields and the Higgs portal coupling among the scalars.
The course and results of gravitational collapse within the setup described above was simulated numerically for two cases with different characters of couplings between the scalars and gravity. One of them was a minimal coupling with vanishing scalar-gravity coupling constants ξ X = ξ h = 0. The other one was the full version of the model which involves the non-minimal couplings ξ X = 0 and ξ h = 0.
The evolutions led to either non-singular spacetimes or spacetimes containing black holes. In the case of a model with scalars minimally coupled to gravity the outcome of a singular collapse is a dynamical Schwarzschild spacetime. There is a central spacelike singularity within it, surrounded by an apparent horizon, which is spacelike in the dynamical spacetime region, that is for small values of advanced time. For positive values of the parameter m 2 the horizon settles down along a null u = const. hypersurface as v tends to infinity and indicates the location of the event horizon in the spacetime. When the mass parameter is negative, there is a null part of the apparent horizon, which is followed by its second spacelike part visible for large values of advanced time.
Inclusion of the non-minimal couplings between scalars and gravity results in a formation of more complex spacetime structures. In all cases, the spacetimes contain central spacelike singularities. When both of the scalar-gravity coupling constants are non-zero and the mass parameter is positive and when only one of the couplings does not vanish, irrespectively of the sign of m 2 , the emerging spacetimes are dynamical Schwarzschild spacetimes. The apparent horizon surrounding the singularity is spacelike and null for small and large values of advanced time, respectively. Between these two parts there exist also timelike portions of the horizon. It seems characteristic for the non-minimality of the scalar-gravity couplings, as was not observed when only minimally coupled scalars were involved in the process, neither in the current paper nor in previous works on the subject (see, e.g., [34] and references therein).
For both non-zero couplings ξ X and ξ h and negative values of the parameter m 2 dynamical Reissner-Nordström spacetimes form. The behavior of the outer apparent horizon is qualitatively similar to the one observed in the cases of non-minimally coupled scalars decribed in the previous paragraph. The second, inner apparent horizon is spacelike within the whole spacetime. There also exists a Cauchy horizon at infinite advanced time. The type of the spacetime which forms during the gravitational process involving non-minimal couplings of scalars to gravity depends on the mass parameter of the complex scalar field. A similar dependency was observed during a collapse with both dark energy and dark matter present in spacetime [19].
The dependencies of the u-locations of the event horizons, radii and masses of black holes emerging from the investigated process on m 2 in both minimally and non-minimally coupled cases are qualitatively the same. As the mass parameter increases, the black holes form later in terms of retarded time and both their radii and masses become smaller. Qualitatively similar relations were observed within the outcomes of gravitational collapse of an electrically charged scalar field accompanied by dark matter and, possibly, additionally by dark energy [19]. The changes are significant up to about m 2 = −0.15 and become much less visible for its larger values. The dependence of the absolute values of black hole U (1) charges on the mass parameter posesses a minimum for its negative values for the cases of both couplings ξ X , ξ h either vanishing or not equal to zero. When m 2 becomes bigger than approximately −0.05, the charges increase become less significant. When only one scalargravity coupling constant is non-zero, Q eh and |T eh | decrease with m 2 and the changes become less dynamical as the mass parameter increases.
The black holes form earlier in terms of retarded time and their radii, masses and charges become smaller as α m increases. Exceptions from this rule are observed in the case of only one non-vanishing scalar-gravity coupling, for which an increase of the u-locations of the event horizons and absolute values of black hole charges is observed for large values of the parameter α m .
Black holes form later in terms of u, their radii and masses decrease while their charges increase with both λ X and λ hX model parameters in the case when ξ X = 0 and ξ h = 0. When either both or only one scalar-gravity coupling vanishes, the quartic couplings do not influence the time of the event horizon formation. With their increase, it turns out that the black holes possess bigger radii and masses and smaller charges. The dependencies of the u-locations of the event horizons, the radii, masses and charges of forming black holes on λ hX in the minimally coupled case and on λ X when only one scalar-gravity coupling does not vanish are complicated as they exhibit discontinuities.
The dependencies of the black hole features on the scalar-gravity coupling ξ h possess extrema. There exists a value of ξ h for which the nascent black hole forms the latest in terms of retarded time and simultaneously possesses the smallest radius and mass. What is interesting, there is also a minimum in the relation between the black hole charges and ξ h , but it does not overlap with the one for u eh , r eh and m eh H . The dependencies on the scalar-gravity coupling ξ X are qualitatively the same with the difference that the observed extrema are shallower.
In all the investigated cases, both non-singular and singular, an increase of the energy density, radial pressure, pressure anisotropy and values of the evolving scalar fields was observed along a null direction of the propagation of the maxima of initially imposed field profiles in spacetime. When dynamical Schwarzschild black holes formed, another increase in values of the quantities measured by an observer moving with the collapsing matter was visible in a close vicinity of the emerging singularity. For dynamical Reissner-Nordström spacetimes the increase was also considerable in the region, in which the r = const. lines settled at null hypersurfaces of constant retarded time, indicating the existence of the Cauchy horizon at infinite advanced time. It seems that calculating the proposed observables within the dynamical spacetimes may allow to distinguish between their types. The increase of these quantities is much bigger as the spacelike singularity is approached in dynamical Reissner-Nordström spacetimes containing the Cauchy horizon, when compared to dynamical spacetimes of the Schwarzschild type.
The local temperature calculated along the apparent horizon of the nascent black holes using the definition for dynamical black holes adopted in the current paper increases for large values of advanced time for all investigated cases. This is a region, in which the apparent horizons are situated along null hypersurfaces. In dynamical regions of spacetimes, where the horizons are either spacelike or timelike, the changes of the values of local temperature are monotonic and increasing in the minimally coupled case and non-monotonic when at least one non-minimal coupling is involved.

A.1 Algorithm setup
The dynamics of the physical system of interest described by equations (2.40)-(2.57) was resolved numerically. The set of equations of motion involves the following functions: d, z 1 , z 2 , y, X 1 , X 2 , h, a, w 1 , w 2 , x, r, f , g, Q, β, T , γ, each of which is a function of two null coordinates, i.e., advanced and retarded times. The dynamics of functions d, z 1 , z 2 and y was followed along the u-coordinate according to equations E4, X (Re) , X (Im) and H, respectively. The remaining functions, X 1 , X 2 , h, a, w 1 , w 2 , x, r, f , g, Q, β, T and γ, evolved along the v-coordinate in line with the respective equations P 6, P 7, P 2, X (Re) , X (Im) , H, P 4, E3, E2, D2, D1, G2 and G1.
The system was solved within a bounded region of the (vu)-plane, which is presented in figure 1 in section 3. A null hypersurface u = const. was taken as an initial data surface. The boundary conditions were posed on a hypersurface of constant v. For purposes of numerics, the two surfaces were marked as u = 0 and v = 0, respectively. Initial conditions involve arbitrary profiles of the evolving fields functions, X 1 , X 2 and h, which were assigned according to (3.2) and (3.1). The initial values of z 1 , z 2 and y were calculated analytically using the relations P 6 and P 7. Since in the employed setup the distribution of matter is shell-shaped, the boundary is not affected by it and the field functions X 1 , X 2 and h vanish there. The boundary values of z 1 , z 2 and y were obtained via integration of equations X (Re) , X (Im) and H, respectively.
Within the investigated setup there remains gauge freedom to choose the initial and boundary profiles of the r function. We chose r (0, 0) to be equal to 7.5 for numerical purposes. The initial and boundary values of g and f , respectively, which determine the distances between the null lines were chosen to be constant, namely g (0, v) = 1 2 and f (v, 0) = − 1 2 . Their precise values are justified by the fact that mass (3.4) should vanish at the central point (0, 0). The r values along the initial null segment were obtained with the use of the equation P 4 and along the boundary using the equation P 3. The initial and boundary profiles of f and g, respectively, were calculated via the integration of the equation E3.

A.2 Employed schemes
The numerical code was written from scratch in Fortran. The integration along the retarded time coordinate was conducted with the use of the 2 nd order accurate Runge-Kutta method. The integration of the partial differential equations along the v-coordinate was performed with the 2 nd order accurate Adams-Bashforth-Moulton method, apart from the first point, where the trapezoidal rule was applied. Adequate v-derivatives of functions z 1 , z 2 and y whose calculation was indispensable to perform the integration of (2.48) were obtained with 2 nd order accurate rules, symmetrical everywhere apart from the points at boundaries of the computational region. The function c is not explicitly involved in the integration of the set of evolution equations (2.40)-(2.57) described above, however, it is needed to compute the observables (3.27)- (3.29). For this purpose, the values of c were calculated according to its definition c = a,u a . The double null coordinates ensure regular behaviour of all the evolving quantities within the domain of integration. However, considerable numerical difficulties arise as the event horizon, where function f diverges, is approached. A relatively dense numerical grid is necessary to determine its location and to examine the behaviour of fields beyond it, especially for large values of advanced time. The efficiency of calculations requires using an adaptive grid and performing integration with a smaller step in particular regions. For the gravitational collapse investigations, the refinement algorithm making the grid denser solely in the u-direction, is sufficient [12]. In order to determine the area of the integration grid, where it should be denser, a local error indicator is needed. This quantity should be bounded with the evolving quantities and change its value significantly in the adequate region. The function ∆r r along the u-coordinate thus indicates the surrounding of the event horizon in spacetime and hence meets our requirements [3].  Figure 19. (color online) The convergence of field functions. The scalar field, h, and the moduli of complex scalar field, |X|, were plotted versus v for evolutions conducted with integration steps, which were multiples of δ = 10 −4 , along hypersurfaces of constant u equal to 1 for (a) Evolution 1 and (b) Evolution 2.

A.3 Tests of the code
The accuracy of the numerical code was checked indirectly, as no analytical solutions exist for the investigated process. The tests were performed for two evolutions initiated with the following parameters. Evolution 1 was conducted with λ X = λ hX = 0.1, m 2 = 0.25, ξ X = ξ h = 0, α m = 0,p X =p h = 0.05 and Evolution 2 with λ X = λ hX = 0.1, m 2 = 0.25, ξ X = ξ h = −0.5, α m = 0,p X =p h = 0.025. The spacetime structures are presented in figures 2(a) and 9(a), respectively.
The first test was based on checking the convergence of the code. In order to monitor the convergence, the computations for Evolutions 1 and 2 were conducted on four grids with integration steps being multiples of the basic value δ = 10 −4 . A step of a particular grid was twice the size of a denser one. The convergence was examined on a u = const. hypersurface chosen arbitrarily with the value of u equal to 1. The selected hypersurface was situated close to the emerging event horizon, but in the region where the adaptive mesh on neither of the grids was active, which was necessary for performing a proper comparison of the results.
The field functions along the selected hypersurface of constant u from within the v-range in which the functions are initially non-zero for the examined integration steps are shown in figure 19. The maximal observed discrepancy between the functions calculated on the finest and coarsest grids was equal to 3 · 10 −4 %. Figure 20 presents the 2 nd order convergence of the numerical code. The maximal divergence between the field profiles obtained on two grids with a quotient of integration steps equal to 2 and their respective quadruples was 10 −1 %. As expected, the errors decreased with an increase of the grid density. The overall error analysis revealed that the expected convergence was achieved and both the algorithm and the numerical code were appropriate for solving the system of equations (2.40)-(2.57).
The second test of the numerical code was based on checking the mass and charge conservation in spacetime. The Hawking mass (3.4) and charges related with the U (1) gauge fields defined by relations (2.19) and (2.18) as functions of retarded time along the line v = 7.5, which was a maximal value of advanced time achieved numerically, are presented in figure 21 for the investigated Evolutions. Since during the course of gravitational collapse the matter was scattered by the gravitational potential barrier when the collapsing shell approached its gravitational radius, the plotted physical quantities were not conserved during the whole evolution. The effect of the outgoing flux was negligible everywhere except for the vicinity of the event horizon. The deviation from the constancy increased with advanced time, as the horizon was approached. The maximal percentage deviations from the particular quantity conservation up to the value of u corresponding to the location of the event horizon were equal to 0.49%, 0.07% and 0.07% for the mass and charges Q and T , respectively. The analysis of mass and charge conservation in spacetime led to the conclusion that the behavior of matter investigated numerically was correct within the computational domain.