Gravitational collapse involving electric charge in the decoupling limit of the dilatonic Gauss--Bonnet gravity

The paper discusses a collapse of a self-interacting electrically charged scalar field in the decoupling limit of the dilatonic Gauss-Bonnet gravity. The emerging spacetimes were either non-singular or singular and containing black holes of a Schwarzschild type. A size of the region, in which the gravitational dynamics was observed was controlled by an absolute value of the Gauss-Bonnet coupling. Dependencies of characteristics of the forming black holes on the dilatonic and Gauss-Bonnet parameters turned out to be similar in the case of black hole masses and radii as well as their time of formation in terms of retarded time. In the cases of masses and radii minima were observed, while in the remaining case a maximum existed. The electric charge of the emerging black holes possessed a maximum when measured versus the dilatonic coupling constant and was strictly decreasing with the Gauss-Bonnet coupling. All the characteristics changed monotonically with the field self-interaction strengths. The times of formation and charges of black holes decreased, while masses and radii increased with the self-interaction strengths of the dynamical fields. Values of the energy density, radial pressure, pressure anisotropy and the collapsing scalar fields were the biggest along the null hypersurface of propagation of the initial peaks of the scalar fields. For big values of the Gauss-Bonnet coupling constant, an increase in their values was also observed in the vicinity of the central singularity within the whole range of advanced time. Non-zero values of the dilaton field outside the black hole event horizon may indicate a formation of a hairy black hole. The local temperature calculated along the apparent horizon was increasing for late times of the evolution, that is in non-dynamical spacetime regions, and exhibited extrema in areas, where the dynamics of the gravity-matter system was observed.


Introduction
Dynamical gravitational collapse is a fully non-linearly and non-perturbatively described process in which a mutual evolution of the spacetime geometry and matter content of the evolving spacetime is tracked. It means that these two features affect each other in the course of the process and hence the described outcomes are closer to reality, but at the same time more difficult to obtain than in the case of matter evolutions on predefined backgrounds, even if backreaction is considered. When an electrically charged scalar field is included in the examined matter-geometry system, the whole system can be regarded as simulating real astrophysical evolutions. Thus obtained spacetime structures satisfactorily resemble the structures anticipated in the case of a non-charged, but rotating collapse, which takes place during real astrophysical events [1,2].
Up to now, structures of spacetimes emerging from the dynamical gravitational collapse were extensively studied first in the simplest cases of sole self-interacting neutral [3] and electrically charged [4][5][6][7][8] scalar fields. More sophisticated theories, complicated in either or both gravitational and matter sectors were studied in [9][10][11][12][13][14][15][16][17][18][19][20]. Apart from investigating structures of spacetimes forming in the process of interest, additional problems related to gravitational dynamics, such as pair creation during collapse and subsequent evaporation of black holes, issues of the cosmic censorship conjecture and the information loss problem, the impact of changing the number of dimensions and the role of various topologies, time measurements as well as observables, were discussed in [1, 2, 7-9, 11, 21-27]. The up-to-date summary of the sketched researches can be found in [28].
The Einstein-dilaton-Gauss-Bonnet (EdGB) theory is a beyond general relativity formulation, which is a scalar-tensor theory including higher order terms in curvature. It has been extensively studied in various cosmological large-scale contexts such as dark energy explanations, inflation or early universe bouncing models [29,30]. It has been also involved in small-scale astrophysical analyses like investigating black holes and their mergers or creation and stability of wormholes [31][32][33][34]. Experimental predictions of the EdGB theory are consistent with current general relativity tests, but the theory yields different results for mergers of the stellar mass black holes [35], what makes it interesting in the context of analyses of the recently extremely viable gravitational waves observations. So far, a spherical gravitational collapse in the EdGB theory was considered within its shift symmetric version in [36]. The research was conducted in coordinates which do not penetrate the horizon, hence the collapse results were described up to the forming horizon and its interior was not analysed. This successful attempt to describe the outcomes of a fully non-linear process within the EdGB theory was an introductory step towards analysis of the whole spacetime, including interiors of emerging singular objects. Our current studies present another approach to this problem, namely the engaged coordinates cover the whole spacetime, yet their usage also requires introducing a simplification of the full EdGB theory.
The decoupling limit of the EdGB theory is a version of the full theory, in which the Gauss-Bonnet term is not included in the derivation of the Einstein equations, only in the equations describing the matter sector of the theory. This means that the scalar sector of the theory does not backreact on geometry [33]. This simplified formulation of the full EdGB theory has been recently employed in research on black hole hair formation and scalar modes [37][38][39]. Investigating the collapse within the truncated version of the EdGB theory will be performed in double null coordinates, which enable to trace the gravitational evolution within the whole spacetime, that is from approximately null infinity up to the central singularity. This is the first step towards investigations of the course and results of gravitational evolutions in the full version of the EdGB theory in the whole forming spacetimes, not only regions exterior to the nascent singular objects.
The paper is organized as follows. In section 2 the considered model was presented. Section 3 contains details of numerical calculations and results analysis. The obtained results were described and discussed in sections 4-6. Section 7 summarizes the undertaken research. Technical details of the numerical code preparation and its tests are presented in appendix A.

Theoretical model of the evolution
The studied theoretical model of the dynamical collapse involves an electrically charged scalar field collapsing within the decoupling limit of the dilaton-Gauss-Bonnet gravity. The general form of the action written, due to the fact that the whole construction involves the string theory concepts, in the string frame is the following: where φ stands for the dilaton field, α is the dilatonic coupling constant and γ determines the coupling between the Gauss-Bonnet contribution and the dilaton. The geometrized units system, in which 8πG = c = 1, was used in the computations. The Lagrangian of the electrically charged scalar field ψ is given by the expression where F βσ is the Maxwell field strength tensor. The covariant derivative isD β =∇ β +ieA β , where e is the electric coupling constant, A β is the four-potential and i denotes the imaginary unit. The Gauss-Bonnet Lagrangian is defined in a standard waŷ 3) The variation of the constructed action (2.1) with respect to matter fields, namely the dilaton field φ, the Maxwell field A µ and the complex scalar field ψ results in the following set of evolution equations: The derivation of the above equations of motion involved a conversion of the string frame into the Einstein frame, which are related via the conformal transformation where g µν andĝ µν denote metric tensors in the Einstein and string frames, respectively [40]. The transformation between the two frames preserves causality. Since the obtained results will be mainly elaborated on employing notions related to causal structures of emerging spacetimes and features of the intrinsic dynamical objects, we do not expect the transformation to influence the ultimate conclusions and hence they can be regarded as physically relevant. From now on, all variables and quantities are written in the Einstein frame. The gravitational field equations derived by varying the action (2.1) with respect to the metric tensor, taking the truncation of the theory into account, will complement the above set of equations describing the examined dynamical system. For the current studies, the double null spherically symmetric line element [41] is chosen where u and v are retarded and advanced time null coordinates, respectively, and dΩ 2 = dΘ 2 + sin 2 ΘdΦ 2 is the line element of the unit sphere, where Θ and Φ are angular coordinates. This choice of coordinates determines the spacetime foliation for conducting computations, which becomes 2+2 [42]. The double null formalism has been widely and successfully employed in dynamical gravitational collapse researches, e.g., [1,2,4,[6][7][8].
Its advantage is that it enables to follow the dynamical evolution from approximately past null infinity, through the formation of possible horizons up to the final central singularity of singular spacetimes.
In spherical symmetry the only non-vanishing components of the electromagnetic field tensor are F uv and F vu . Due to the gauge freedom A u → A u + ∇ u θ , where θ = A v dv, the only non-zero four-vector component is A u . It is a function of retarded and advanced time.
The dilaton field equation of motion (2.4) in the chosen coordinate system is given by where we set Q is a function of retarded and advanced time and determines the electric charge within a sphere of a radius r(u, v) on a spacelike hypersurface containing the point (u, v). Partial derivatives with respect to the null coordinates are marked as ,u and ,v . Concerning the assumed line element (2.9) and the definition of electric charge (2.11), the v-component of Maxwell equations (2.5) can be separated into two first-order differential equations. The first one governs the evolution of the only non-zero component of the four-vector of the Maxwell field A u,v − Qa 2 2r 2 = 0, (2.12) and the second one describes the behavior of the quantity Q during the evolution The equations for the complex scalar field (2.6)-(2.7) become The stress-energy tensor for the considered electrically charged scalar field evolving within the decoupling limit of the dilaton-Gauss-Bonnet theory is the following: Its non-vanishing components written in double null coordinates are Combining the adequate components of the Einstein tensor resulting from the metric (2.9) and the obtained stress-energy tensor components, the Einstein equations yield They complement the matter equations presented above in order to obtain the complete set of equations describing the dynamics of the examined system, which are (2.10)-(2.15), without the relation (2.11), which is a definition of the physical quantity Q.
The following set of auxiliary variables is introduced to prepare the obtained set of dynamical equations to be solved numerically: It is supplemented by the quantities The above variables make it possible to rewrite the second-order differential equations (2.10), (2.14)-(2.15) and (2.21)-(2.24) as first-order ones. Additionally, real fields ψ 1 and ψ 2 are introduced instead of conjugate fields ψ and ψ * according to relations After introducing the above substitutions, the final system of equations of motion describing the gravitational collapse of interest can be written as P 1 : a ,u = ac, (2.28) (2.44)

Details of numerical simulations and results analysis
The obtained complex coupled system of differential equations (2.28)-(2.44) does not possess an analytic solution and needs to be solved numerically. The prepared code and tests run to confirm its correctness are presented in appendix A. The equations describing the dynamics of the considered systems were solved in the bounded region of the (vu)-plane. It is shown in figure 1 on the background of a dynamical Schwarzschild spacetime, whose Carter-Penrose diagram does not differ from the static case. The borders of the computational domain were marked for numerical purposes as 0 and 7.5 in the v-direction, 0 and 15 in the u-direction in all conducted simulations. The only arbitrary data of the numerical computations were initial profiles of the evolving fields, posed on the null hypersurface denoted as u = 0. In order to describe the behavior of the real and complex scalar fields properly, their initial profiles were of the following Gaussian and trigonometric types, respectively, [3,4,6,43]: The profiles were one-parameter families with amplitudesp h andp s being the free family parameters. The amplitudes determine the strength of the gravitational self-interaction of the particular field [44]. The remaining arbitrary constants were invariable during computations, precisely c 1 = 1.3, c 2 = 0.21 and the parameter of the amount of initial charge δ = π 2 . The final value of advanced time was v f = 2.5. The initial conditions are representative for the conducted evolutions, as their outcomes do not depend of the profiles types provided that they are regular, what means that they result in a regular spacetime slice at the initial u = const hypersurface. This condition is fulfilled by the above profiles (3.1) and (3.2).
When the value of the electric coupling constant is not equal to zero, it does not affect the results of the collapse [13]. It was confirmed for the investigated cases. Thus e was set as equal to 1 in all evolutions. The case of a vanishing electric coupling constant is beyond the scope of the current research, as it comes down to an analysis of two neutral fields instead of one electrically charged in the spacetime.
The values of the remaining model parameters, i.e., α and γ were assigned as follows. The dilatonic coupling constant was equal to − √ 3, −1 and 0, while γ was taken from within the range −1, 1 . The values of α corresponded to the dimensionally reduced Kaluza-Klein theory, dilaton gravity and the Einstein-Maxwell theory. The range of γ was consistent with observational constraints (for a thorough analysis see [33] and references therein).
The spacetime structures resulting from the dynamical evolutions of interest will be presented on Penrose diagrams, which contain contours of r = const lines plotted in the (vu)-plane. The outermost line refers to r = 0, which in all presented cases will be a part of a central singularity. The lines indicating the vanishing expansion will be also presented on the diagrams and they will indicate the location of an apparent horizon in the spacetime. On the spacetime diagrams they will be marked as red solid lines. On the remaining plots they will be marked in black. One of physical quantities which was employed in the interpretation of the obtained results is the quasi-local Hawking mass [45]. Its value calculated for the spherically symmetric spacetime with a gauge field A µ is the following: It describes the amount of mass contained within a sphere of a radius r(u, v) on a spacelike hypersurface containing the particular point (u, v). The mass of a singular object can be expressed by a value of the above expression calculated at the event horizon in the region of the emerging spacetime which is not dynamical, that is for large values of advanced time.
The outcomes of the gravitational collapse within the discussed model were also inspected with the use of a set of local spacetime quantities that may be interpreted as observables related with an observer moving with the evolving matter. The analysed quantities were energy densityρ, radial pressurep r and pressure anisotropyp a ≡p t −p r . Their precise covariant derivation for any general case can be found in [27]. Moreover, local temperature T l = κ l 2π with the surface gravity defined as κ l = a ,u a −2 [46,47] was also calculated. For the studied gravity-matter model considered in double null coordinates the first three of the above observables are given by relationŝ

Dynamical spacetime structures
The spacetimes which emerge during the investigated collapse are either non-singular, for small values of the field self-interaction strengths, or singular, containing in all cases a dynamical black hole of a Schwarzschild type, for sufficiently big self-interactions represented by values of field amplitudes. The singular spacetimes resulting from dynamical evolutions   Each resulting dynamical Schwarzschild spacetime contains a central spacelike singularity located along the r = 0 line. It is indicated by the outermost r = const line on the Penrose diagram. The singularity is surrounded by a single apparent horizon located along the r ,v = 0 line, whose behavior divides the whole spacetime into two regions, i.e., a dynamical one, within which the actual collapse proceeds, and a non-dynamical one, remaining after the dynamical evolution. The former corresponds to the range of small values of advanced time, in which the apparent horizon is spacelike and changes its location in the u-coordinate. The latter region refers to the v-range of big values, where v → ∞ and the apparent horizon settles along a null hypersurface of constant retarded time, thus indicating the location of the event horizon in spacetime.
The analysis of the behavior of the apparent horizons in the presented cases allows to

Black hole characteristics
The examined characteristics of forming black holes were the u-locations of the event horizons, radii, masses and electric charges of black holes formed during the gravitational evolu-  tions within the model of interest. They were examined as functions of the model couplings α − √ 3, √ 3 and γ −1, 1 , as well as field amplitudesp s =p h 0.01, 0.09 for a selected evolution characterized by the following parameters: α = −1, γ = 1,p s =p h = 0.04. When the dependence on the specific parameter of the model is presented, the remaining ones are as listed above. In the case of the dependency on the field self-interaction strengths, several combinations of the model parameters were considered.
The dependencies of the abovementioned nascent black holes characteristics on the model parameters are presented in figures 6 and 7. The dependencies of the u-locations of the event horizons, radii and masses on both α and γ are in both cases qualitatively the same, that is there exist a maximum for u eh at α = 0.25 and γ = −0.2 and minima at α = 0.1 and γ = −0.1 for the remaining characteristics. The black hole electric charge increases for small values of the dilatonic coupling constant up to an extremum α = 0.05 and decreases for its larger values. Q eh decreases monotonically with an increasing γ. Figure 8 presents the discussed black hole characteristics as functions of the collapsing fields amplitudes for several combinations of the couplings present within the studied model. In all cases, for increasing field self-interaction strengths the changes of the black hole characteristics are monotonic. The u-locations of the event horizons and black hole electric charges decrease, while radii and masses of the forming singular objects increase with increasingp s =p h . The changes of u eh , r eh and m eh H become smaller with increasing field amplitudes. An opposite tendency is observed for Q eh .

Observables and fields
The (vu)-distributions of observables presented in section 3 and the evolving scalar fields, along with the relation between local temperature along the apparent horizon and the vcoordinate will be shown and discussed for selected spacetimes, the structures of which were shown in section 4. Figures 9 and 10 present the considered spacetime distributions for the spacetime whose structure was shown in figure 2(b), with the model parameters α = − √ 3 and γ = −0.01. The highest absolute values of the discussed observables as well as both the neutral and the moduli of the complex scalar field functions are located along a constant-v null direction, which is a direction of peaks of initially imposed field functions propagation. A considerable increase in their values appears as the central singularity is approached. This increase is observed only within the range of advanced time, in which the highest field values were imposed initially. The energy density, radial pressure and the moduli of the  The distributions of observables and field functions resulting from the evolution leading to the spacetime with the structure shown in figure 3(b), obtained with α = −1 and γ = −1, are presented in figures 11 and 12, respectively. As in the case discussed above, the energy Figure 11.
(color online) The (vu)distribution of (a) energy density,ρ, (b) radial pressure,p r , and (c) pressure anisotropy, p a , and (d) local temperature along the black hole apparent horizon, T l , as a function of advanced time for a dynamical evolution characterized by parameters α = −1, γ = −1 and p s =p h = 0.04 (the same as in figure 3(b)).
density and the moduli of the complex scalar field are positive, while the pressure anisotropy and the neutral scalar field function are negative within the whole domain of integration. What distinguishes this case from the other discussed is the fact that non-zero energy density and radial pressure persist also for large values of advanced time. Moreover, pressure anisotropy is also non-zero for large v. This implies that there is a persisting non-trivial matter distribution around the central singularity. This observation finds its confirmation in figure 12(a). The dilaton field does not collapse completely but instead attains a nonzero value even outside the black hole event horizon. This may indicate that in this case a formation of a hairy black hole is observed. The local temperature in the examined case is positive and increases monotonically with advanced time along the apparent horizon with a slight inclination in the region, where an inclination of the apparent horizon is also visible. Unlike the case depicted in figure 9(d), the late-lime increase is not linear.
The distributions of the discussed quantities, that is observables and field functions, for the evolution with the emerging spacetime structure presented in figure 4(b), resulting from the evolution with α = −1 and γ = 1, are shown in figures 13 and 14. The signs of the particular quantities are the same as in the cases presented above. An increase of absolute values of the observables is visible in the region corresponding to advanced time where the field values was the highest initially. The energy density, radial pressure and pressure anisotropy also display an increase in values nearby the central singularity, as in the latter of the above cases. The local temperature along the black hole apparent horizon behaves like in the first case above, namely it increases for small values of v, then, after reaching a maximum it slightly decreases and the decrease turns into the late-time nearly linear increase as v → ∞.  Figure 12. (color online) The (vu)-distribution of (a) the neutral scalar field, h, and (b) the moduli of the complex scalar field, |s|, for the same parameters and field amplitudes as in figure 11. Figure 13.
(color online) The (vu)distribution of (a) energy density,ρ, (b) radial pressure,p r , and (c) pressure anisotropy,p a , and (d) local temperature along the black hole apparent horizon, T l , as a function of advanced time for a dynamical evolution characterized by parameters α = −1, γ = 1 andp s =p h = 0.04 (the same as in figure 4(b)).
The spacetime distribution of the analysed observables and fields for the evolution with α = 0 and γ = −0.1, with the structure of spacetime shown in figure 5(b), presented in figures 15 and 16, is analogous to the one presented above for the case of α = − √ 3 and γ = −0.01. The energy density, radial pressure and the moduli of the complex scalar field are positive, while the pressure anisotropy and the neutral scalar field function are negative within the dynamical spacetime region covered by numerical calculations. A sole increase in absolute values in all the distributions is visible along the null v = const direction of propagation of the highest initial values of the imposed field functions. The changes in values of the local temperature calculated along the black hole apparent horizon tend to a late-time increase close to linear, after reaching a shallow extrema, precisely a maximum  Figure 14. (color online) The (vu)-distribution of (a) the neutral scalar field, h, and (b) the moduli of the complex scalar field, |s|, for the same parameters and field amplitudes as in figure 13. Figure 15.
(color online) The (vu)distribution of (a) energy density,ρ, (b) radial pressure,p r , and (c) pressure anisotropy, p a , and (d) local temperature along the black hole apparent horizon, T l , as a function of advanced time for a dynamical evolution characterized by parameters α = 0, γ = −0.1 and p s =p h = 0.04 (the same as in figure 5(b)). and a minimum for small values of advanced time in the range where an inclination of the apparent horizon appears.

Conclusions
The course and outcomes of gravitational evolutions involving electric charge in a decoupling limit of the dilatonic Gauss-Bonnet gravity was investigated. The system consisted of two scalar fields, namely a neutral scalar dilaton field and a complex scalar field coupled with a U (1) gauge field, i.e., electrically charged. The dilaton was coupled both with the matter and Gauss-Bonnet sectors of the theory. The course of the dynamical collapse was traced numerically. The formation of dynamical spacetimes and emerging black holes was observed and analysed via description of spacetime structures, inspection of the black hole characteristics and values of gravitational observables and evolving fields within the forming spacetimes.
The collapse resulted in either non-singular spacetimes formed for small self-interaction strengths of the fields, that is small initial values of their amplitudes, or singular spacetimes containing black holes. The particular spacetime contained a central spacelike singularity along r = 0 surrounded by a single apparent horizon, which was spacelike in the dynamical region, i.e., for small values of advanced time, and became null as v → ∞. Disregarding the values of the model parameters α and γ, the black holes were of a Schwarzschild type, despite the fact that the evolving scalar field was electrically charged. Such a situation was observed in the case of a dynamical gravitational collapse of an electrically charged scalar field in dilaton gravity, but only for non-zero values of the dilatonic coupling constant [13]. Taking this into account, a conclusion that the Gauss-Bonnet term in the gravitational sector suppresses the tendency to form a Reissner-Nordström spacetime with an inner Cauchy horizon, whose appearance is characteristic for evolutions involving electric charge, can be drawn. The absolute value of the parameter γ influenced the v-range of the dynamical spacetime region, precisely the range of advanced time in which gravitational dynamics was observed was narrower for smaller values of |γ|.
The dependencies of the u-locations of the event horizons, radii and masses of black holes emerging from the investigated process on α and γ are qualitatively the same. As values of the model parameters increase, the black holes form later in terms of retarded time and both their radii and masses become smaller, up to an extremum, in which the tendencies reverse. Similar dependencies were observed in the case of non-minimal scalar-gravity couplings during dynamics investigations within the theory involving Higgs and dark matter sectors [27]. The only difference between the dependencies on the model parameters is that the changes are more significant for small values of α and large values of γ. The dependence of the values of black hole U (1) charge on α and γ do not overlap the ones describe above. In the case of the dilatonic coupling constant, Q eh increases significantly for its small values and decreases also significantly after reaching a maximum. In the case of the Gauss-Bonnet coupling, the black hole electric charge decreases monotonically within the whole parameter range.
The changes of the u-locations of the event horizons, radii, masses and electric charges of the nascent black holes with increasing self-interaction strengths of the evolving scalar fields are monotonic and qualitatively the same for various combinations of values of the model parameters α and γ. The features u eh and Q eh increase, while r eh and m eh H decrease with the parameters. The changes are more significant for the electric charge in comparison to the remaining black hole characteristics.
In all the investigated cases an increase of the energy density, radial pressure, pressure anisotropy and values of the collapsing scalar fields was observed along a null direction of propagation of the maxima of initially imposed field profiles in spacetime. For large absolute values of the γ parameter 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 also for large values of advanced time. This implies a persisting nontrivial matter distribution around the central singularity. Additionally, an observation that the dilaton field possesses non-zero values outside the black hole event horizon may indicate a formation of a hairy black hole in this case. The local temperature calculated along the apparent horizon of the emerging black holes shows a late-time monotonic increase for all the investigated cases. Within dynamical spacetime regions, where inclinations of the horizons appear, the changes of the values of local temperature are not monotonic and extrema are observed.
The ultimate aim of the undertaken research is to investigate gravitational evolution of scalar fields in the full version of the EdGB theory in double null coordinates. These coordinates cover the whole spacetime, that is both the exterior and interior of objects that may arise in dynamically formed spacetimes. The presented analysis within the truncated version of the underlying theory is a first step towards calculations in the full EdGB. Its equations require a more sophisticated treatment, as they are coupled in a complicated way. The issue has been only partially resolved by introducing a shift symmetry in the case of coordinates, which do not penetrate emerging horizons, i.e., when only the exterior of the forming objects is possible to be thus dealt with [36].

Algorithm setup
The evolution of the studied physical system is described by equations (2.28)-(2.44). It was resolved numerically. The set of equations of motion involves quantities d, q 1 , q 2 , y, s 1 , s 2 , h, a, p 1 , p 2 , x, r, f , g, Q, β. Each function depends on two null coordinates, namely advanced and retarded times. The dynamics of d, q 1 , q 2 and y was followed along u in line with the equations E4, S (Re) , S (Im) and D, respectively. The remaining quantities, s 1 , s 2 , h, a, p 1 , p 2 , x, r, f , g, Q and β, evolved along v according to the respective equations P 6, P 7, P 2, S (Re) , S (Im) , D, P 4, E3, E2, M 2 and M 1.
The system of evolution equations was solved in a bounded region of the (vu)-plane presented in figure 1 in section 3. An arbitrary null hypersurface of constant retarded time was taken as an initial data surface. The boundary conditions were posed on a hypersurface of constant advanced time. The two surfaces were marked as u = 0 and v = 0, respectively, for computational purposes.
Initial conditions are arbitrary profiles of the fields functions, s 1 , s 2 and h, which were posed according to (3.2) and (3.1). The initial values along v of q 1 , q 2 and y were calculated analytically using the relations P 6 and P 7. In the employed setup the distribution of matter is shell-shaped, hence the boundary is unaffected by it and the field functions s 1 , s 2 and h vanish there. The boundary values of q 1 , q 2 and y were obtained through integration of equations S (Re) , S (Im) and D, respectively.
A gauge freedom to choose initial and boundary profiles of the r function remains within the investigated setup. r (0, 0) was chosen to be equal to 7.5 for computational purposes. Initial and boundary values of g and f , respectively, determine the distances between the null lines and were chosen to be constant, that is g (0, v) = 1 2 and f (v, 0) = − 1 2 . These values are justified by the fact that mass (3.4) should vanish at the central point (0, 0). The r values on the initial null segment were obtained using the relation P 4 and along the boundary with the equation P 3. Initial and boundary profiles of f and g, respectively, were obtained via an integration of E3.
Initial values of the quantity d were calculated with the use of the equation E2, and its boundary values were obtained using E4. The spherical shell shape of the matter distribution justifies imposing the following boundary values: a (u, 0) = 1, Q (u, 0) = β (u, 0) = 0 and p 1 (u, 0) = p 2 (u, 0) = x (u, 0) = 0. Initial profiles of these functions were obtained using the equations P 2, M 2, M 1, S (Re) , S (Im) and D, respectively.

Employed schemes
The numerical code was written in Fortran from scratch. Integration along the u-coordinate involved the 2 nd order accurate Runge-Kutta method. Integration of the partial differential equations along advanced time was performed with the 2 nd order accurate Adams-Bashforth-Moulton method, except the first point, where the trapezoidal rule was used.
The double null coordinates selected for the analysis ensure regular behaviour of all the evolving quantities within the computational domain. Numerical difficulties arise as the event horizon is approached, as the function f diverges there. For this reason, a relatively dense numerical grid is indispensable to determine the horizon location and to examine the behaviour of fields beyond it, especially for large v-coordinate values. The efficiency of calculations was ensured by the use of an adaptive grid and performing integration with a smaller step in problematic regions. For the gravitational collapse investigations, a sufficient refinement algorithm is the one making the grid denser solely in the direction of retarded time [13]. The determination of the area of the computational grid, where it should be denser was made through a local error indicator. The quantity should be bounded with the evolving quantities and change its value significantly in adequate regions. The quantity ∆r r along the u-coordinate meets the requirements and indicates the numerically problematic surrounding of the event horizon in spacetime [6].

Tests of the code
Analytical solutions do not exist for the investigated process and hence the accuracy of the code has to be checked via numerical possibilities in this regard. The convergence tests were performed for two evolutions characterized by the following parameters. Evolution 1 was initiated with α = − √ 3, γ = −0.01, while Evolution 2 with α = −1, γ = −1, and in both casesp s =p h = 0.04. The corresponding spacetime structures are presented in figures 2(b) and 3(b), respectively.
To monitor the numerical outcomes convergence, the computations for Evolutions 1 and 2 were carried out on four grids with integration steps equal to multiples of δ = 10 −4 . An integration step of a particular grid was twice the size of a denser one. The convergence was examined on a hypersurface of constant retarded time selected arbitrarily with u = 1. The chosen hypersurface was located close to the forming event horizon in the region where the adaptive mesh on neither of the grids was active, what enabled a proper comparison of the results.
The evolving field functions along the chosen hypersurface of u = const from within the range of advanced time in which the functions are initially non-vanishing for all the examined integration steps are shown in figure 17. The maximal observed discrepancy between the functions calculated on the grids with the smallest and biggest steps was equal to 2.55 · 10 −4 %. Figure 18 confirms 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 4 · 10 −1 %. The errors decreased as the grid density increased. The overall analysis revealed that the expected convergence was achieved and both the algorithm and the numerical code were appropriate for solving the obtained system of equations (2.28)-(2.44) describing the dynamics of interest.  Figure 18. (color online) The convergence of the code. The differences between the scalar field functions, ∆h, and the moduli of complex scalar field, ∆|s|, calculated on grids with different integration steps (multiples of δ = 10 −4 ) and their multiples were obtained along the same hypersurfaces of constant u as in figure 17 for (a) Evolution 1 and (b) Evolution 2.